Fixed orbit

magnetism.c

@@ -6,19 +6,27 @@
 #define EPSILON 8.85418782E-12
 
 void magnetism(universe_t* uni){
-  mass_body_t body_1;
-  mass_body_t body_2;
-  Vector2 diff, force;
-  for (unsigned int i = 0; i < uni->num_bodies-1; i++) {
-    body_1 = uni->bodies[i];
-      for (unsigned int y = i + 1; y < uni->num_bodies; y++) {
-        body_2 = uni->bodies[y];
-        diff=vec2_diff(body_2.position,body_1.position);
-        force = vec2_scale(diff,body_1.charge*body_2.charge/(4*PI*EPSILON*vec2_len(diff)));
-        body_1.speed=vec2_sum(body_1.speed,vec2_scale(force,1/body_1.mass));
-        body_2.speed=vec2_diff(body_2.speed,vec2_scale(force,1/body_2.mass));
-        uni->bodies[y] = body_2;
+  for (unsigned int i = 0; i < uni->num_bodies; i++) {
+    mass_body_t *body_1 = &uni->bodies[i];
+    for (unsigned int y = i + 1; y < uni->num_bodies; y++) {
+      mass_body_t *body_2 = &uni->bodies[y];
+      if (body_1->charge == 0 || body_2->charge == 0) {
+        continue;
+      }
+      Vector2 diff = vec2_diff(body_2->position, body_1->position);
+      float distance_squared = vec2_dot(diff, diff);
+      if (distance_squared == 0) {
+        continue;
+      }
+      // Charges are in coulombs; convert acceleration to Earth radii/s^2.
+      float factor = PHYSICS_TIMESTEP * body_1->charge * body_2->charge /
+                     (4 * PI * EPSILON * EARTH_MASS * DISTANCE_UNIT_METERS *
+                      DISTANCE_UNIT_METERS * DISTANCE_UNIT_METERS * distance_squared);
+      Vector2 impulse = vec2_scale(vec2_normalize(diff), factor);
+      body_1->speed = vec2_diff(body_1->speed,
+                               vec2_scale(impulse, 1 / body_1->mass));
+      body_2->speed = vec2_sum(body_2->speed,
+                              vec2_scale(impulse, 1 / body_2->mass));
     }
-    uni->bodies[i] = body_1;
   }
 }