Add a particle simulation

include/quad.h

@@ -61,5 +61,6 @@ void node_remove_point(node_t* node, point_t* point);
 void node_merge(node_t* node);
 void qtree_update_point(quadtree_t* qtree, point_t* old_point, point_t* new_point);
 void qtree_delete(node_t* node);
+void qtree_graph(node_t* node);
 
 #endif // QUAD_H

main.c

@@ -3,8 +3,7 @@
 #include <limits.h>
 #include <stdlib.h>
 #include <time.h>
-// TODO: this should be less than what the visualiser can handle
-#define NPARTICLES 100
+#define NPARTICLES 200 // make sure it's less than the visualiser's capacity
 
 
 static size_t id_particle = 0;
@@ -15,47 +14,55 @@ typedef struct particle_t {
     point_t point;
 } particle_t;
 
-void reflect(particle_t* particle, rect_t* boundary, float dt) {
-    const int w = boundary->x1;
-    const int h = boundary->y1;
-    if (particle->point.x + particle->velx * dt > w) {
-        particle->point.x -= w - particle->velx * dt;
-        particle->velx *= -1;
-    } else if (particle->point.x + particle->velx * dt < boundary->x0) {
-        particle->point.x = -particle->velx * dt - particle->point.x;
-        particle->velx *= -1;
-    }
-    if (particle->point.y + particle->vely * dt > h) {
-        particle->point.y -= h - particle->vely * dt;
-        particle->vely *= -1;
-    } else if (particle->point.y + particle->vely * dt < boundary->y0) {
-        particle->point.y = -particle->vely * dt - particle->point.y;
-        particle->vely *= -1;
-    }
-}
-
 static particle_t particles[NPARTICLES];
-
 static size_t ids = 0;
+static float random_norm() {
+    return (float) random() / RAND_MAX;
+}
+static float accel = -5;
+static int niters = 300;
 
 int main() {
     srand(time(NULL));
     rect_t boundary = {.x0=0, .y0=0, .x1=400, .y1=400};
     point_t points[NPARTICLES];
     particle_t particles[NPARTICLES];
-    for (int i = 0; i < NPARTICLES; ++i) {
-        points[i].x = random() % boundary.x1;
-        points[i].y = random() % boundary.y1;
-        points[i].id = ids++;
-    }
     quadtree_t qtree;
     viz_init(boundary.x1, boundary.y1);
-    qtree.root = node_new(&boundary);
     for (int i = 0; i < NPARTICLES; ++i) {
-        node_insert(qtree.root, &points[i]);
+        point_t p = {random() % boundary.x1, random() % boundary.y1, ids++};
+        particles[i].point = p;
+        // TODO: clamp it
+        particles[i].velx = (random_norm() > 0.66) ?  5 + random_norm()*10 : -5 - random_norm()*10;
+        particles[i].vely = (random_norm() > 0.66) ? -5 - random_norm()*10 :  5 + random_norm()*10;
+        particles[i].accelerationy = accel;
+        particles[i].accelerationx = 0;
+    }
+    float dt = 0.1;
+    for (int i = 0; i < niters; ++i) {
+        quadtree_t qtree;
+        qtree.root = node_new(&boundary);
+        for (int ip = 0; ip < NPARTICLES; ++ip) {
+            particles[ip].velx += particles[ip].accelerationx * dt;
+            particles[ip].vely += particles[ip].accelerationy * dt;
+            particles[ip].point.x += particles[ip].velx * dt;
+            particles[ip].point.y += particles[ip].vely * dt;
+            // reflection
+            particles[ip].point.x %= boundary.x1;
+            particles[ip].point.y %= boundary.y1;
+            if (particles[ip].point.x <= boundary.x0)
+                particles[ip].point.x = boundary.x1 - 5;
+            if (particles[ip].point.y <= boundary.y0)
+                particles[ip].point.y = boundary.y1 - 5;
+
+            point_t pnew = {particles[ip].point.x , particles[ip].point.y, particles[ip].point.id};
+            node_insert(qtree.root, &pnew);
+        }
+        qtree_graph(qtree.root);
+        viz_flush();
+        usleep(33000);
+        qtree_delete(qtree.root);
     }
-    qtree_graph(qtree.root);
-    viz_flush();
-    sleep(3);
+    sleep(1);
     viz_close();
 }

src/quad.c

@@ -1,10 +1,12 @@
 #include "quad.h"
+#include "viz.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <stdbool.h>
 #include <math.h>
 #include <assert.h>
 #include <string.h>
+#include <unistd.h>
 
 #define MAX(a, b) ((a) > (b) ? (a) : (b))
 
@@ -254,6 +256,7 @@ void qtree_update_point(quadtree_t* qtree, point_t* old_point, point_t* new_poin
     node_remove_point(root, old_point);
     old_point->x = new_point->x;
     old_point->y = new_point->y;
+    old_point->id = new_point->id;
     node_insert(root, old_point);
 }
 
@@ -269,3 +272,24 @@ void qtree_delete(node_t* node) {
         free(node);
     }
 }
+
+
+void qtree_graph(node_t* node) {
+    if (node == NULL) return;
+
+    // If it's a leaf node, print the boundary and the points in it
+    if (node_is_leaf(node)) {
+        viz_write_rect(&node->boundary);
+        // Print all points in this leaf
+        for (int i = 0; i < node->count; ++i) {
+            viz_write_point(&node->points[i]);
+        }
+    } else {
+        // Recursively visit all children (if not a leaf)
+        qtree_graph(node->nw);
+        qtree_graph(node->ne);
+        qtree_graph(node->sw);
+        qtree_graph(node->se);
+    }
+}
+

src/viz.c

@@ -4,8 +4,6 @@
 #include <stdlib.h>
 #include <unistd.h> // usleep
 
-#define NOBJECTS 3000
-
 void viz_init(unsigned width, unsigned height) {
     // don't close the window upon exit
     viz_plot_pipe = popen("gnuplot -persistent", "w");