update fixmes about root body simulation

also manual TODO

docs/TODO

@@ -13,8 +13,8 @@ If you see modifications to this file in git status, IGNORE them and do not comm
         frame. we need to test, and make a compromise on performance vs sim
         accuracy, and have this set as a config constant, MAX_STEP_SIZE, with
         a good discription nearby in the code.
-    - see FIXMEs in simulation.cpp, and physics.cpp about not simulating root
-        bodies anymore
+        - could also try an alternative step method with error function
+            eg) https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method
     - remove simulation dependancy from physics.h
     - fix moon tests with SOI coordinate transforms for planets
     - remove bodies on non closed orbits after they are a certain distance from

src/main.cpp

@@ -65,26 +65,26 @@ void run_gui_simulation(SimulationState* sim) {
             printf("Speed multiplier: %.1fx\n", speed_multiplier);
         }
 
-#if 1
+        if (!paused) {
+            int steps = (int)(physics_steps_per_frame * speed_multiplier);
+            for (int i = 0; i < steps; i++) {
+                update_simulation(sim);
+            }
+        }
+
+#if 0
         /****
         * Hacks ahead...
         */
-		int focus_index = 4;  // earth: 2, jupiter, 4
+        int focus_index = 4;  // earth: 2, jupiter, 4
         if (sim->body_count >= focus_index) {
-			Vector3 offset = {0, 100, 100};
+            Vector3 offset = {0, 2, 10};
             focus_camera(&render_state, &sim->bodies[focus_index], offset);
         }
 #else
         update_camera(&render_state);
 #endif
 
-        if (!paused) {
-            int steps = (int)(physics_steps_per_frame * speed_multiplier);
-            for (int i = 0; i < steps; i++) {
-                update_simulation(sim);
-            }
-        }
-
         render_simulation(sim, &render_state);
     }
 

src/physics.cpp

@@ -77,29 +77,17 @@ Vec3 evaluate_acceleration(Vec3 pos, Vec3 vel, AccelerationContext* ctx) {
 
     Vec3 total_force = {0.0, 0.0, 0.0};
 
-    // FIXME: we shouldn't be simulating root bodies anymore
-    if (temp_body.parent_index == -1) {
-        for (int j = 0; j < ctx->sim->body_count; j++) {
-            if (j == ctx->body_index) continue;
-            CelestialBody* other = &ctx->sim->bodies[j];
-            if (other->parent_index == -1) {
-                Vec3 force = calculate_gravity_force(&temp_body, other);
-                total_force = vec3_add(total_force, force);
-            }
-        }
-    } else {
-        if (temp_body.parent_index >= 0 && temp_body.parent_index < ctx->sim->body_count) {
-            CelestialBody* parent = &ctx->sim->bodies[temp_body.parent_index];
-
-            double distance = vec3_magnitude(pos);
-            if (distance < 1.0) {
-                distance = 1.0;
-            }
-
-            double force_magnitude = G * temp_body.mass * parent->mass / (distance * distance);
-            Vec3 direction = vec3_normalize(vec3_scale(pos, -1.0));
-            total_force = vec3_scale(direction, force_magnitude);
+    if (temp_body.parent_index >= 0 && temp_body.parent_index < ctx->sim->body_count) {
+        CelestialBody* parent = &ctx->sim->bodies[temp_body.parent_index];
+
+        double distance = vec3_magnitude(pos);
+        if (distance < 1.0) {
+            distance = 1.0;
         }
+
+        double force_magnitude = G * temp_body.mass * parent->mass / (distance * distance);
+        Vec3 direction = vec3_normalize(vec3_scale(pos, -1.0));
+        total_force = vec3_scale(direction, force_magnitude);
     }
 
     return calculate_acceleration(total_force, temp_body.mass);

src/simulation.cpp

@@ -76,22 +76,6 @@ void update_soi(CelestialBody* body, CelestialBody* parent, double semi_major_ax
 }
 
 void update_simulation(SimulationState* sim) {
-    for (int i = 0; i < sim->body_count; i++) {
-        CelestialBody* body = &sim->bodies[i];
-
-        // FIXME: we shouldn't be simulating root bodies anymore
-        if (body->parent_index == -1) {
-            AccelerationContext ctx;
-            ctx.sim = sim;
-            ctx.current_body = body;
-            ctx.body_index = i;
-
-            rk4_step(body, &ctx, sim->dt);
-        }
-    }
-
-    compute_global_coordinates(sim);
-
     for (int i = 0; i < sim->body_count; i++) {
         CelestialBody* body = &sim->bodies[i];
 
@@ -101,6 +85,8 @@ void update_simulation(SimulationState* sim) {
 
         int new_parent = find_dominant_body(sim, i);
         if (new_parent != body->parent_index && new_parent != -1) {
+            // FIXME: there's now a bug here when changing parent bodes
+            //  I think it's that the local position is never updated
             body->parent_index = new_parent;
         }
 
@@ -123,31 +109,18 @@ static void compute_orbital_velocity_from_vis_viva(CelestialBody* body,
                                                      CelestialBody* parent) {
     Vec3 r = vec3_sub(body->position, parent->position);
     double distance = vec3_magnitude(r);
-
-    if (distance < 1.0) {
-        body->velocity = {0.0, 0.0, 0.0};
-        return;
-    }
-
     double e = body->eccentricity;
     double a = body->semi_major_axis;
-
-    // FIXME: there's a bug here
-    //  could be unit related, mass is in units of kg, and distance is meters
-    //  insteam of km, velocity is also m/s
-    //double speed = (double) sqrt(G * parent->mass * (2.0 / distance - 1.0 / a));
     double v_squared = G * parent->mass * (2.0 / distance - 1.0 / a);
+    // FIXME: need error handling here
+    //  if v_squared is negative, the inital body parameters are bad
     assert(v_squared >= 0);
     double speed = (double) sqrt(v_squared);
-
-    if (e > 0.01) {
-        printf("  %s: eccentric orbit (e=%.2f, a=%.3e m), speed at r=%.3e m: %.3f km/s\n",
-               body->name, e, a, distance, speed / 1000.0);
-    }
-
     Vec3 z_axis = {0.0, 0.0, 1.0};
     Vec3 vel_dir = vec3_cross(r, z_axis);
 
+    // NOTE: I suppose this is for the case of a high inclination orbit?
+    // FIXME: make a test to see if this is necessary
     if (vec3_magnitude(vel_dir) < 0.01) {
         Vec3 x_axis = {1.0, 0.0, 0.0};
         vel_dir = vec3_cross(r, x_axis);
@@ -159,6 +132,7 @@ static void compute_orbital_velocity_from_vis_viva(CelestialBody* body,
     body->velocity = vec3_add(body->velocity, parent->velocity);
 }
 
+// FIXME: remove this function since we're already looping in calculate_initial_velocities
 static void set_child_bodies_velocity(SimulationState* sim) {
     for (int i = 0; i < sim->body_count; i++) {
         CelestialBody* body = &sim->bodies[i];
@@ -198,6 +172,8 @@ void calculate_soi_radii(SimulationState* sim) {
     }
 }
 
+// FIXME: actually all the above functions should be called inside this loop
+//  or the loop should be in config_loader
 void initialize_local_coordinates(SimulationState* sim) {
     for (int i = 0; i < sim->body_count; i++) {
         CelestialBody* body = &sim->bodies[i];