interstellar interference (2008)
[maximus:dr1.git] / interstellar-interference / nbody.c
1 #include <math.h>
2 #include <stdlib.h>
3
4 #include <m_pd.h>
5
6 t_class *nbody_class;
7
8 struct vec3 {
9   double x;
10   double y;
11   double z;
12 };
13
14 struct body {
15   double m;
16   struct vec3 p;
17   struct vec3 v;
18   struct vec3 a;
19 };
20
21 struct nbody {
22   t_object pd;
23   t_outlet *o1;
24   t_outlet *o2;
25   unsigned int n;
26   double g;
27   double h;
28   struct body *bs;
29 };
30
31 void nbody_g(struct nbody *self, t_float g) {
32   self->g = g;
33 }
34
35 void nbody_step(struct nbody *self, t_float sf) {
36   int s = floor(sf);
37   if (s < 1) s = 1;
38   for (int k = 0; k < s; ++k) {
39     for (int i = 0; i < self->n; ++i) {
40       struct vec3 *pi = &self->bs[i].p;
41       struct vec3 f = { 0, 0, 0 };
42       for (int j = 0; j < self->n; ++j) {
43          if (i != j) {
44            struct vec3 *pj = &self->bs[j].p;
45            struct vec3 d = { pj->x - pi->x, pj->y - pi->y, pj->z - pi->z };
46            double d2 = d.x * d.x + d.y * d.y + d.z * d.z;
47            double k = self->g * self->bs[j].m / d2;
48            f.x += k * d.x; f.y += k * d.y; f.z += k * d.z;
49          }
50       }
51       self->bs[i].a.x = f.x; self->bs[i].a.y = f.y; self->bs[i].a.z = f.z;
52     }
53     double h = self->h;
54     for (int i = 0; i < self->n; ++i) {
55       struct body *b = &self->bs[i];
56       b->v.x += b->a.x * h; b->v.y += b->a.y * h; b->v.z += b->a.z * h;
57       b->p.x += b->v.x * h; b->p.y += b->v.y * h; b->p.z += b->v.z * h;
58     }
59   }
60 }
61
62 void nbody_randomize(struct nbody *self) {
63   double m = 0;
64   struct vec3 c = { 0, 0, 0 };
65   for (int i = 0; i < self->n; ++i) {
66     struct body *b = &self->bs[i];
67     b->m = 1;
68     b->p.x = 20.0 * ((double) rand() / (double) RAND_MAX - 0.5);
69     b->p.y = 20.0 * ((double) rand() / (double) RAND_MAX - 0.5);
70     b->p.z = 20.0 * ((double) rand() / (double) RAND_MAX - 0.5);
71     b->v.x = 0; b->v.y = 0; b->v.z = 0;
72     b->a.x = 0; b->a.y = 0; b->a.z = 0;
73     m += b->m;
74     c.x += b->m * b->p.x; c.y += b->m * b->p.y; c.z += b->m * b->p.z;
75   }
76   c.x /= m; c.y /= m; c.z /= m;
77   for (int i = 0; i < self->n; ++i) {
78     struct body *b = &self->bs[i];
79     b->p.x -= c.x; b->p.y -= c.y; b->p.z -= c.z;
80   }  
81 }
82
83 void nbody_positions(struct nbody *self) {
84   t_atom atoms[4];
85   for (int i = 0; i < self->n; ++i) {
86     struct vec3 *p = &self->bs[i].p;
87     SETFLOAT(&atoms[0], i+1);
88     SETFLOAT(&atoms[1], p->x);
89     SETFLOAT(&atoms[2], p->y);
90     SETFLOAT(&atoms[3], p->z);
91     outlet_list(self->o1, &s_list, 4, atoms);
92   }
93 }
94
95 void nbody_distances(struct nbody *self) {
96   t_atom atoms[3];
97   for (int i = 0; i < self->n; ++i) {
98     struct vec3 *pi = &self->bs[i].p;
99     for (int j = 0; j < self->n; ++j) {
100       if (i < j) {
101         struct vec3 *pj = &self->bs[j].p;
102         struct vec3 d = { pj->x - pi->x, pj->y - pi->y, pj->z - pi->z };
103         SETFLOAT(&atoms[0], i+1);
104         SETFLOAT(&atoms[1], j+1);
105         SETFLOAT(&atoms[2], sqrt(d.x * d.x + d.y * d.y + d.z * d.z));
106         outlet_list(self->o2, &s_list, 3, atoms);
107       }
108     }
109   }
110 }
111
112 struct nbody *nbody_new(t_floatarg n) {
113   struct nbody *self = pd_new(nbody_class);
114   self->o1 = outlet_new(&self->pd, &s_list);
115   self->o2 = outlet_new(&self->pd, &s_list);
116   self->n = floor(n);
117   if (self->n < 1) self->n = 1;
118   self->g = 1;
119   self->h = 0.001;
120   self->bs = malloc(self->n * sizeof(struct body));
121   nbody_randomize(self);
122   return self;
123 }
124
125 void nbody_free(struct nbody *self) {
126   if (self && self->bs) {
127     free(self->bs);
128   }
129 }
130
131 extern void nbody_setup() {
132   nbody_class = class_new(gensym("nbody"), nbody_new, nbody_free, sizeof(struct nbody), 0, A_DEFFLOAT, 0);
133   if (!nbody_class) return;
134   class_addmethod(nbody_class, nbody_distances, gensym("distances"), 0);
135   class_addmethod(nbody_class, nbody_positions, gensym("positions"), 0);
136   class_addmethod(nbody_class, nbody_randomize, gensym("randomize"), 0);
137   class_addmethod(nbody_class, nbody_step, gensym("step"), A_FLOAT, 0);
138   class_addmethod(nbody_class, nbody_g, gensym("g"), A_FLOAT, 0);
139 }