added code

experiment/newton/.gitignore

 *.o
-newton
\ No newline at end of file
+newton
+pendulum
\ No newline at end of file

experiment/newton/Makefile

@@ -5,10 +5,13 @@ CFLAGS=-Wall -Wextra -std=gnu11 -O3
 LIBS=-lSDL2 -lm
  
 
-all: newton
+all: newton pendulum
 
 newton: newton.o gfx.o
 	$(C_GCC) $(CFLAGS) $^ -o $@ $(LIBS)
 
+pendulum: pendulum.o gfx.o
+	$(C_GCC) $(CFLAGS) $^ -o $@ $(LIBS)
+
 clean:
-	rm -f *.o gfx_example newton
+	rm -f *.o gfx_example newton pendulum

experiment/newton/pendulum.c

+#include "gfx.h"
+#include <math.h>
+#include <stdbool.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+
+#define GRAV 9.8
+
+typedef struct _point
+{
+  double x, y;
+} point;
+
+point point_new(double x, double y)
+{
+  point p = {x, y};
+  return p;
+}
+
+point point_from_polar(double r, double theta)
+{
+  return point_new(r * cos(theta), r * sin(theta));
+}
+
+point point_add(point lhs, point rhs)
+{
+  point res = {lhs.x + rhs.x, lhs.y + rhs.y};
+  return res;
+}
+
+point point_sub(point lhs, point rhs)
+{
+  point res = {lhs.x - rhs.x, lhs.y - rhs.y};
+  return res;
+}
+
+point point_mul(double lhs, point rhs)
+{
+  point res = {lhs * rhs.x, lhs * rhs.y};
+  return res;
+}
+
+double point_norm_sqr(point lhs)
+{
+  return lhs.x * lhs.x + lhs.y * lhs.y;
+}
+
+double point_norm(point lhs)
+{
+  return sqrt(point_norm_sqr(lhs));
+}
+
+point compute_deriv(double dt, point p, point p_1)
+{
+  return point_mul(1.0 / dt, point_sub(p, p_1));
+}
+
+point compute_force(point pos, double g, double mass)
+{
+  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);
+  return point_add(grav, string);
+}
+
+static void render_point(struct gfx_context_t *context, int x, int y, uint32_t color)
+{
+  for (int i = -2; i <= 2; ++i)
+  {
+    for (int j = -2; j <= 2; ++j)
+    {
+      gfx_putpixel(context, x + i, y + j, color);
+    }
+  }
+}
+
+/// @param context graphical context to use.
+static void render(struct gfx_context_t *context, point pos, double length)
+{
+  gfx_clear(context, COLOR_BLACK);
+  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;
+  render_point(context, x, y, COLOR_WHITE);
+
+  gfx_present(context);
+}
+
+point verlet(point p, point p_1, point a, double dt)
+{
+  return point_add(point_sub(point_mul(2.0, p), p_1), point_mul(dt * dt, a));
+}
+
+int main(int argc, char **argv)
+{
+  double mass = 0.0;
+  double length = 0.0;
+  bool debug = false;
+  if (argc == 3)
+  {
+    mass = atof(argv[1]);
+    length = atof(argv[2]);
+  }
+  else if (argc == 4)
+  {
+    mass = atof(argv[1]);
+    length = atof(argv[2]);
+    debug = (0 == strcmp("-d", argv[2]));
+  }
+  else
+  {
+    printf("Error: missing mass or length parameters.\n\
+                Normal use is: ./newton <mass> <length> -d\n\
+                where -d is for debug mode and is optional.\n\
+                The units are in IS.\n\
+                Example: ./newton 1.0 1.0\n");
+    return EXIT_FAILURE;
+  }
+  int width = 1000, height = 1000;
+  struct gfx_context_t *ctxt = gfx_create("The Pendulum", width, height);
+  if (!ctxt)
+  {
+    fprintf(stderr, "Graphics initialization failed!\n");
+    return EXIT_FAILURE;
+  }
+
+  double angle = (double)rand() / (double)RAND_MAX * 2.0 * M_PI;
+  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;
+
+  struct timespec time = {0, dt * 1.0e9};
+
+  int num = 1000;
+  if (debug)
+  {
+    printf("time position velocity acceleration\n");
+  }
+  else
+  {
+    printf("time position\n");
+  }
+
+  gfx_clear(ctxt, COLOR_BLACK);
+
+  point avg_acc = point_new(0.0, 0.0);
+  for (int i = 0; i < num; ++i)
+  {
+    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 (debug)
+      {
+        printf("t = %f, pos = (%f, %f), vel = (%f, %f), acc = (%f, %f)\n",
+               i * dt, p.x, p.y, vel.x, vel.y, acc.x, acc.y);
+      }
+      else
+      {
+        printf("%f, %f, %f\n", i * dt, p.x, p.y);
+      }
+      avg_acc = point_add(avg_acc, acc);
+    }
+    render(ctxt, p, length);
+    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_2 = p_1;
+    p_1 = p_tmp;
+  }
+  if (debug)
+  {
+    point a = point_mul(1.0 / (num - 1), avg_acc);
+    printf("Simulation ended: avg acceleration = (%f, %f)\n", a.x, a.y);
+  }
+
+  while (gfx_keypressed() != SDLK_ESCAPE)
+  {
+  }
+
+  gfx_destroy(ctxt);
+
+  return EXIT_SUCCESS;
+}
\ No newline at end of file