diff --git a/experiment/newton/.gitignore b/experiment/newton/.gitignore
index ba31bf0a27bc5a896a8bd1351e44ba31e8b69ac7..6ff2e4cd4de8ffb1d3a6427db87b796f71b0aa43 100644
--- a/experiment/newton/.gitignore
+++ b/experiment/newton/.gitignore
@@ -1,2 +1,3 @@
*.o
-newton
\ No newline at end of file
+newton
+pendulum
\ No newline at end of file
diff --git a/experiment/newton/Makefile b/experiment/newton/Makefile
index adec1dbfd2678a7fce0a8715aee9aab3d2503979..ca339506d540b1b9296677fdcc3e5b93b5e29af3 100644
--- a/experiment/newton/Makefile
+++ b/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
diff --git a/experiment/newton/pendulum.c b/experiment/newton/pendulum.c
new file mode 100644
index 0000000000000000000000000000000000000000..066f4b21ca0905ece935c628a2b847f86db1bc5c
--- /dev/null
+++ b/experiment/newton/pendulum.c
@@ -0,0 +1,191 @@
+#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