Add better physical simulation, solved 1. tas

ovm/ovm.c

@@ -3,6 +3,7 @@
 #include <glib/gstdio.h>
 
 #undef CMPEPS
+// #define CMPEPS
 
 void ovm_free(ovm_t *ovm) {
 	if (!ovm) return;
@@ -82,7 +83,7 @@ void ovm_print_c(ovm_t *ovm, const gchar *filename) {
 				break;
 			case SOP_CMPZ:
 				switch (instr_sop_cmp(oi[i])) {
-#if CMPEPS /* cmp with eps */
+#ifdef CMPEPS /* cmp with eps */
 				case CMP_LTZ:
 					g_string_printf(str, "\tovm_status = (v%u.d < -eps);\n", USED(instr_sop_r1(oi[i])));
 					break;
@@ -142,7 +143,7 @@ void ovm_print_c(ovm_t *ovm, const gchar *filename) {
 			g_string_printf(str, "\tv%u.d = v%u.d * v%u.d;\n", USED(i), USED(instr_dop_r1(oi[i])), USED(instr_dop_r2(oi[i])));
 			break;
 		case OP_DIV:
-#if CMPEPS
+#ifdef CMPEPS
 			g_string_printf(str, "\tv%u.d = (fabs(v%u.d) < eps) ? 0.0 : v%u.d / v%u.d;\n", USED(i), USED(instr_dop_r2(oi[i])), USED(instr_dop_r1(oi[i])), instr_dop_r2(oi[i]));
 #else
 			g_string_printf(str, "\tv%u.d = (v%u.d == 0) ? 0.0 : v%u.d / v%u.d;\n", USED(i), USED(instr_dop_r2(oi[i])), USED(instr_dop_r1(oi[i])), instr_dop_r2(oi[i]));
@@ -173,7 +174,7 @@ void ovm_print_c(ovm_t *ovm, const gchar *filename) {
 		"static gboolean ovm_status;\n"
 		"\n"
 	));
-#if CMPEPS
+#ifdef CMPEPS
 	write(f, CONST_STR_LEN(
 		"static const gdouble eps = 1e-300;\n"
 		"\n"

ovm/physics.h

@@ -9,5 +9,57 @@
 #define Me 6.0e24
 #define Re 6.357e6
 
+typedef struct {
+	gdouble x, y;
+} vector_t;
+
+typedef struct {
+	vector_t g, v;
+	gdouble v_abs;
+} move_data_t;
+
+typedef struct {
+	vector_t pos_old, pos;
+	gdouble rad2_old, rad_old, rad2, rad;
+	vector_t dv_old, dv;
+
+	gboolean b_old, b; /* whether the following was calculated */
+	move_data_t move_old, move;
+} satellite_t;
+
+static inline void satellite_update_pos(satellite_t *s, gdouble x, gdouble y) {
+	s->pos_old = s->pos;
+	s->rad2_old = s->rad2;
+	s->rad_old = s->rad;
+	s->dv_old = s->dv;
+	s->b_old = s->b;
+	s->move_old = s->move;
+
+	s->pos.x = x; s->pos.y = y;
+	s->rad2 = s->pos.x*s->pos.x + s->pos.y*s->pos.y;
+	s->rad = sqrt(s->rad2);
+	s->dv.x = 0.0; s->dv.y = 0.0;
+
+	s->b = FALSE;
+}
+
+static inline void satellite_update_move(satellite_t *s) {
+	if (!s->b_old) {
+		s->move_old.g.x = -(G*Me/(s->rad2_old*s->rad_old)) * s->pos_old.x;
+		s->move_old.g.y = -(G*Me/(s->rad2_old*s->rad_old)) * s->pos_old.y;
+		s->move_old.v.x = s->pos.x - s->pos_old.x - 0.5*(s->move_old.g.x + s->dv_old.x);
+		s->move_old.v.y = s->pos.y - s->pos_old.y - 0.5*(s->move_old.g.y + s->dv_old.y);
+		s->move_old.v_abs = sqrt(s->move_old.v.x*s->move_old.v.x+s->move_old.v.y*s->move_old.v.y);
+		s->b_old = TRUE;
+	}
+	if (!s->b) {
+		s->move.g.x = -(G*Me/(s->rad2*s->rad)) * s->pos.x;
+		s->move.g.y = -(G*Me/(s->rad2*s->rad)) * s->pos.y;
+		s->move.v.x = s->move_old.v.x + s->dv_old.x + 0.5*(s->move_old.g.x + s->move.g.x);
+		s->move.v.y = s->move_old.v.y + s->dv_old.y + 0.5*(s->move_old.g.y + s->move.g.y);
+		s->move.v_abs = sqrt(s->move.v.x*s->move.v.x+s->move.v.y*s->move.v.y);
+		s->b = TRUE;
+	}
+}
 
 #endif

ovm/task.h

@@ -9,6 +9,7 @@ typedef struct {
 	guint in_size, out_size;
 	gdouble *in, *out, *in_old;
 	int tracefd;
+	gboolean finished;
 } task_t;
 
 task_t* task_new(guint scenario);

ovm/task1.c

@@ -14,8 +14,8 @@ Address  Sensor
 
 #define SCORE (task->out[0])
 #define FUEL (task->out[1])
-#define S_X (task->out[2])
+#define S_X (-task->out[2])
-#define S_Y (task->out[3])
+#define S_Y (-task->out[3])
 #define RADIUS (task->out[4])
 
 #define DV_X (task->in[2])
@@ -30,44 +30,84 @@ static void debug(task_t * task) {
 		task->timestamp, SCORE, FUEL, S_X, S_Y, RADIUS, rad);
 }
 
-gdouble init_rad, last_x, last_y;
-gboolean running = FALSE, started = FALSE;
-gdouble dv, dv_tic;
-
 static void run(task_t *task, gpointer userdata) {
-	gdouble cur_rad;
+	static satellite_t sat;
+	static gdouble init_rad, target_v;
+	static gdouble dv, dv_tic;
+	static gboolean reached = FALSE, outoffuel = TRUE;
+	static guint32 reached_ts = 0;
 	UNUSED(userdata);
 
-	if (0 == task->timestamp) return;
+	if (task->finished) debug(task);
+	if (0 == task->timestamp || task->finished) return;
+
+	satellite_update_pos(&sat, S_X, S_Y);
 
-	cur_rad = sqrt(S_X*S_X+S_Y*S_Y);
 	DV_X = DV_Y = 0.0;
 
-	if (!running) {
+	if (!reached) {
-		running = TRUE;
+		if (1 == task->timestamp) {
-		started = TRUE;
+			debug(task);
-		init_rad = cur_rad;
+			init_rad = sat.rad;
+
+			/* v^2/r == G*M_e/r^2 => v^2 = mu/r */
+			target_v = sqrt(mu/RADIUS);
+
 			dv = sqrt(mu/init_rad)*(sqrt(2*RADIUS/(init_rad + RADIUS))-1);
 			dv_tic = sqrt(mu/RADIUS)*(1-sqrt(2*init_rad/(init_rad + RADIUS)));
-		printf("Orbit: %f -> %f; dv = %f, dv' = %f\n", init_rad, RADIUS, dv, dv_tic);
+			printf("-- Orbit: %f -> %f; dv = %f, dv' = %f\n", init_rad, RADIUS, dv, dv_tic);
-	} else if (started) {
+		} else if (2 == task->timestamp) {
-		gdouble vx = S_X - last_x, vy = S_Y - last_y, v = sqrt(vx*vx+vy*vy);
+			debug(task);
-		started = FALSE;
+			satellite_update_move(&sat);
-		DV_X = -vx * (dv / v);
+			if (sat.move.v.x * S_Y > 0.0 && sat.move.v.y * S_X < 0.0) target_v = -target_v;
-		DV_Y = -vy * (dv / v);
+			DV_X = sat.move.v.x * (dv / sat.move.v_abs);
-// 		printf("v: %f / %f\n", vx, vy);
+			DV_Y = sat.move.v.y * (dv / sat.move.v_abs);
+			printf("-- Leaving initial orbit at %u\n", (guint) task->timestamp);
+			debug(task);
+			printf("v: %f / %f\n", sat.move.v.x, sat.move.v.y);
 			printf("Dv: %f / %f\n", DV_X, DV_Y);
-	} else if (fabs(cur_rad - RADIUS) < 1000.0) {
+		} else if (fabs(sat.rad - RADIUS) < 100.0) {
-		gdouble vx = S_X - last_x, vy = S_Y - last_y, v = sqrt(vx*vx+vy*vy);
+			satellite_update_move(&sat);
-		running = FALSE;
+			debug(task);
-		printf("Reached target radius\n");
+			reached = TRUE;
-		DV_X = -vx * (dv_tic / v);
+			DV_X = sat.move.v.x * (dv_tic / sat.move.v_abs);
-		DV_Y = -vy * (dv_tic / v);
+			DV_Y = sat.move.v.y * (dv_tic / sat.move.v_abs);
+			printf("-- Reached target orbit at %u\n", (guint) task->timestamp);
+			printf("v: %f / %f\n", sat.move.v.x, sat.move.v.y);
 			printf("Dv: %f / %f\n", DV_X, DV_Y);
 		}
+	} else if (!outoffuel) {
+		/* target_v is clockwise (TODO) */
+		gdouble tar_dv_x = -target_v * S_Y / sat.rad;
+		gdouble tar_dv_y = target_v * S_X / sat.rad;
+		satellite_update_move(&sat);
+		DV_X = tar_dv_x - sat.move.v.x;
+		if (fabs(DV_X) < 0.5) DV_X = 0.0;
+		DV_Y = tar_dv_y - sat.move.v.y;
+		if (fabs(DV_Y) < 0.5) DV_Y = 0.0;
+		if (DV_X*DV_X + DV_Y*DV_Y > FUEL*FUEL) {
+			printf("-- Out of fuel\n");
+			printf("v: %f / %f\n", sat.move.v.x, sat.move.v.y);
+			printf("Dv: %f / %f\n", DV_X, DV_Y);
+			DV_X = 0.0; DV_Y = 0.0;
+			outoffuel = TRUE;
+		}
+	}
 
-	last_x = S_X;
+	if (reached_ts > 0 && fabs(sat.rad - RADIUS) >= 1000.0) {
-	last_y = S_Y;
+// 		debug(task);
+		printf("-- Left target radius at %u (duration: %u)\n", (guint) task->timestamp, (guint) task->timestamp - reached_ts);
+		reached_ts = 0;
+	} else if (reached_ts == 0 && fabs(sat.rad - RADIUS) < 1000.0) {
+// 		debug(task);
+		reached_ts = (guint) task->timestamp;
+		printf("-- Reached target radius at %u\n", (guint) task->timestamp);
+	}
+
+	if (fabs(DV_X) < 0.1) DV_X = 0.0;
+	if (fabs(DV_Y) < 0.1) DV_Y = 0.0;
+
+	sat.dv.x = DV_X; sat.dv.y = DV_Y;
 }
 
 static void run_debug(task_t *task, gpointer userdata) {

ovm/task_main.h

@@ -20,6 +20,7 @@ static inline void trace_simulate(task_t *task, gpointer userdata) {
 		guint32 i;
 		trace->pos++;
 		for (i = 0; i < te->count; i++) {
+			if (te->data[i].addr != 0x3e80)
 				task->in[te->data[i].addr] = te->data[i].value;
 		}
 	}
@@ -36,6 +37,8 @@ static inline gdouble task_run(task_t *task, task_app_t app, gpointer userdata,
 		if (dotrace) task_trace_step(task);
 		ovm_step(task->scenario, task->in, task->out);
 	}
+	task->finished = TRUE;
+	app(task, userdata);
 	if (dotrace) task_trace_end(task);
 	return task->out[0];
 }

ovm/trace.h

@@ -3,13 +3,15 @@
 
 #include "base.h"
 
+typedef struct {
+	guint32 addr;
+	gdouble value  __attribute__ ((__packed__));
+} trace_entry_point_t;
+
 typedef struct {
 	guint32 timestamp;
 	guint32 count;
-	struct {
+	trace_entry_point_t data[];
-		guint32 addr;
-		gdouble value;
-	} data[];
 } trace_entry_t;
 
 typedef struct {