added pendulum

experiment/newton/pendulum.c

@@ -35,15 +35,25 @@ point point_sub(point lhs, point rhs)
   return res;
 }
 
+point point_neg(point lhs)
+{
+  return point_new(-lhs.x, -lhs.y);
+}
+
 point point_mul(double lhs, point rhs)
 {
   point res = {lhs * rhs.x, lhs * rhs.y};
   return res;
 }
 
+double point_scalar_product(point lhs, point rhs)
+{
+  return lhs.x * rhs.x + lhs.y * rhs.y;
+}
+
 double point_norm_sqr(point lhs)
 {
-  return lhs.x * lhs.x + lhs.y * lhs.y;
+  return point_scalar_product(lhs, lhs);
 }
 
 double point_norm(point lhs)
@@ -51,6 +61,11 @@ double point_norm(point lhs)
   return sqrt(point_norm_sqr(lhs));
 }
 
+point point_normalize(point lhs)
+{
+  return point_mul(1.0 / point_norm(lhs), lhs);
+}
+
 point compute_deriv(double dt, point p, point p_1)
 {
   return point_mul(1.0 / dt, point_sub(p, p_1));
@@ -58,10 +73,12 @@ point compute_deriv(double dt, point p, point p_1)
 
 point compute_force(point pos, double g, double mass)
 {
-  double angle = atan2(pos.y, pos.x);
+  // double angle = atan2(pos.y, pos.x);
   double amplitude = g * mass;
   point grav = point_new(0.0, -amplitude);
-  point string = point_from_polar(amplitude, angle);
+  point neg = point_neg(point_normalize(pos));
+  double prod = point_scalar_product(neg, grav);
+  point string = point_mul(-prod, neg);
   return point_add(grav, string);
 }
 
@@ -83,8 +100,8 @@ static void render(struct gfx_context_t *context, point pos, double length)
   int cx = context->width / 2;
   int cy = context->height / 2;
   render_point(context, cx, cy, COLOR_RED);
-  int x = cx + pos.x / length * context->width / 2;
-  int y = cy + pos.y / length * context->height / 2;
+  int x = cx + pos.x / length * context->width / 4;
+  int y = cy - pos.y / length * context->height / 4;
   render_point(context, x, y, COLOR_WHITE);
 
   gfx_present(context);
@@ -109,7 +126,7 @@ int main(int argc, char **argv)
   {
     mass = atof(argv[1]);
     length = atof(argv[2]);
-    debug = (0 == strcmp("-d", argv[2]));
+    debug = (0 == strcmp("-d", argv[3]));
   }
   else
   {
@@ -128,11 +145,12 @@ int main(int argc, char **argv)
     return EXIT_FAILURE;
   }
 
-  double angle = (double)rand() / (double)RAND_MAX * 2.0 * M_PI;
+  // double angle = (double)rand() / (double)RAND_MAX * 2.0 * M_PI;
+  double angle = (double)-M_PI/2.0+0.1;
   point p = point_from_polar(length, angle);
   point p_1 = p; // TODO: change that....
   point p_2 = p; // TODO: change that....
-  double dt = 0.01;
+  double dt = 0.001;
 
   struct timespec time = {0, dt * 1.0e9};
 
@@ -154,7 +172,7 @@ int main(int argc, char **argv)
     point vel = compute_deriv(dt, p, p_1);
     point vel_1 = compute_deriv(dt, p_1, p_2);
     point acc = compute_deriv(dt, vel, vel_1);
-    if (i > 0)
+    if (i >= 0)
     {
       if (debug)
       {
@@ -171,7 +189,7 @@ int main(int argc, char **argv)
     nanosleep(&time, NULL);
     point p_tmp = p;
     point new_acc = point_mul(1.0 / mass, compute_force(p, GRAV, mass));
-    p = verlet(p, p_1, new_acc, dt); // add some noise of the order of
+    p = verlet(p, p_1, new_acc, dt);
     p_2 = p_1;
     p_1 = p_tmp;
   }