3 * @purpose Implementations for N-Body simulator
8 #include <assert.h> // To make sure I don't do anything stupid.
9 #include <math.h> //For maths and stuff
10 #include <stdio.h> //For writing data to files
11 #include "nbody.h" //Declarations of functions and structures in this file
13 // --- Body functions --- //
17 * @purpose Initialise a Body object with mass, position and velocity
18 * @param body - Target Body
19 * @param m - The mass of the Body
20 * @param x[...] - The position vector of the body
21 * @param v[...] - The velocity vector of the body
23 void Body_Init(Body * body, float m, float x[DIMENSIONS], float v[DIMENSIONS])
27 for (unsigned i = 0; i < DIMENSIONS; ++i)
33 body->x_old[i] = x[i];
34 body->v_old[i] = v[i];
38 body->collided = false;
42 * @function Body_Destroy
43 * @purpose Destroy a Body
44 * @param body - Target Body
45 * @param system - System to which the Body belongs
47 void Body_Destroy(Body * body, System * system)
56 * @purpose Perform a single step, updating target Body
57 * @param body - Target Body
58 * @param system - System that the Body is part of (used for universal constants, etc)
60 void Body_Step(Body * body, System * system)
62 assert(body != NULL && system != NULL && body->exists);
63 body->collided = false;
64 for (unsigned i = 0; i < DIMENSIONS; ++i)
67 //Each Body keeps track of position and velocity on the previous step
68 //This is used for collision handling in the overall System step.
69 body->x_old[i] = body->x[i];
70 body->v_old[i] = body->v[i];
72 body->x[i] += body->v[i] * system->dt; //Update the position vector
73 body->v[i] += (body->F[i] / body->m) * system->dt; //Update the velocity vector
74 body->F[i] = 0.00; //Reset the net force vector
80 * @function Body_Force
81 * @purpose Update the net force on Body A due to Body B
82 * @param A - Body to update; subject of force from B
83 * @param B - Body causing force on A; NOT CHANGED BY THIS FUNCTION
84 * @param system - System which Body A is a part of (used for universal constants, etc)
86 void Body_Force(Body * A, Body * B, System * system)
89 assert(A != NULL && B != NULL && system != NULL && A->exists && B->exists);
94 for (unsigned i = 0; i < DIMENSIONS; ++i)
95 r2 += powf((B->x[i] - A->x[i]), 2);
97 float k = system->G * (A->m * B->m) / powf(r2, 1.50);
98 for (unsigned i = 0; i < DIMENSIONS; ++i)
99 A->F[i] += k * (B->x[i] - A->x[i]);
104 * @function Body_Radius
105 * @purpose Calculate the radius of a body
106 * @param body - Target Body
107 * @returns The radius of body
109 float Body_Radius(Body * body)
111 assert(body != NULL && body->exists);
112 return 8.0 * sqrt(body->m);
117 * @function Body_CheckCollision
118 * @purpose Check for collision between bodies
119 * @param A - First body
120 * @param B - Second body
121 * @param system - System that the bodies belong to
122 * @param t - pointer to variable that will store the time of closest approach
123 * @returns true if collision occurred
125 bool Body_CheckCollision(Body * A, Body * B, System * system, float * t)
127 assert(A != NULL && B != NULL && system != NULL && A->exists && B->exists);
132 // Determine time and position of closest approach
133 float x_rel[DIMENSIONS];
134 float v_rel[DIMENSIONS];
137 for (unsigned i = 0; i < DIMENSIONS; ++i)
139 x_rel[i] = B->x_old[i] - A->x_old[i];
140 v_rel[i] = B->v_old[i] - A->v_old[i];
141 c += v_rel[i] * x_rel[i];
142 coeff_t += v_rel[i] * v_rel[i];
146 if (*t <= system->dt/100.0 || *t > system->dt)
150 for (unsigned i = 0; i < DIMENSIONS; ++i)
152 x_rel[i] += *t * v_rel[i];
153 r2 += pow(x_rel[i], 2);
156 return (r2 <= pow(Body_Radius(A), 2) + pow(Body_Radius(B), 2));
160 * @function Body_HandleCollision
161 * @purpose Handle a collision between two Bodies
162 * @param A - The first Body
163 * @param B - The second Body
164 * @param system - The System to which the bodies belong
165 * @param t - The time of closest approach
167 void Body_HandleCollision(Body * A, Body * B, System * system, float t)
169 fprintf(stderr, "#Collision %p %p at t %f\n", (void*)A, (void*)B, t);
170 fprintf(stderr,"# Before Collide; A velocity is %f %f, B velocity is %f %f\n", A->v[0], A->v[1], B->v[0], B->v[1]);
171 //Calculate A's new velocity and position
172 for (unsigned i = 0; i < DIMENSIONS; ++i)
174 A->v[i] = A->m * A->v_old[i] + B->m * B->v_old[i] + B->m * system->cor * (B->v_old[i] - A->v_old[i]);
175 A->v[i] = A->v[i] / (A->m + B->m);
177 A->x[i] = A->x_old[i] + A->v_old[i] * t; //Move to the collision point at the original velocity
178 A->x[i] += A->v[i] * (system->dt - t); //Move the rest of the step at the new velocity
182 //For cor = 0, bodies always merge; B is travelling with A
183 if (system->cor == 0.0)
185 A->m += B->m; //Increase A's mass
186 Body_Destroy(B, system); //Destroy B
190 //Calculate B's new velocity and position
191 for (unsigned i = 0; i < DIMENSIONS; ++i)
193 B->v[i] = B->m * B->v_old[i] + A->m * A->v_old[i] + A->m * system->cor * (A->v_old[i] - B->v_old[i]);
194 B->v[i] = B->v[i] / (B->m + A->m);
196 B->x[i] = B->x_old[i] + B->v_old[i] * t; //Move to the collision point at the original velocity
197 B->x[i] += B->v[i] * (system->dt - t); //Move the rest of the step at the new velocity
201 //Determine whether Body B has escaped Body A; if it hasn't, merge the two bodies
204 fprintf(stderr, "# After Collide; A velocity is %f %f, B velocity is %f %f\n", A->v[0], A->v[1], B->v[0], B->v[1]);
208 for (unsigned i = 0; i < DIMENSIONS; ++i)
210 v2 += pow(B->v[i], 2);
211 r2 += pow(B->x[i] - A->x[i], 2);
214 fprintf(stderr, "# Collision; v is %f, escape is %f and %f\n", v2, 2.0 * system->G * A->m / sqrt(r2),2.0 * system->G * B->m / sqrt(r2));
215 if (v2 > 2.0*system->G * A->m / sqrt(r2) && v2 > 2.0*system->G * B->m / sqrt(r2))
218 //Inelastically merge the two bodies
219 for (unsigned i = 0; i < DIMENSIONS; ++i)
221 A->x[i] = (A->m * A->x[i] + B->m * B->x[i]) / (A->m + B->m);
222 A->v[i] = (A->m * A->v[i] + B->m * B->v[i]) / (A->m + B->m);
225 Body_Destroy(B, system);
230 * @function Body_WriteData
231 * @purpose Write data about the Body to a file
232 * @param body - The target Body
233 * @param system - System to which the Body belongs
234 * @param file - File to write data to
236 void Body_WriteData(Body * body, System * system, FILE * file)
238 assert(body != NULL && system != NULL && file != NULL && body->exists);
239 fprintf(file,"%d\t%p\t%f", system->step, (void*)(body), body->m);
240 for (unsigned i = 0; i < DIMENSIONS; ++i)
241 fprintf(file,"\t%f", body->x[i]);
242 for (unsigned i = 0; i < DIMENSIONS; ++i)
243 fprintf(file,"\t%f", body->v[i]);
244 for (unsigned i = 0; i < DIMENSIONS; ++i)
245 fprintf(file,"\t%f", body->F[i]);
251 // --- System functions --- //
254 * @function System_Init
255 * @purpose Initialise a System object
256 * @param system - Target System
257 * @param x - The origin of the System (currently used for offset of graphics)
258 * @param dt - Time step
259 * @param G - Universal Gravitational Constant
260 * @param cor - The coefficient of restitution for collisions
261 * @param nMax - Number of Body objects this System will contain (DEFAULT = 10)
263 void System_Init(System * system, float x[DIMENSIONS], float dt, float G, float cor, unsigned nMax)
265 assert(system != NULL);
266 for (unsigned i = 0; i < DIMENSIONS; ++i)
276 system->bodies = (Body*)(calloc(sizeof(Body), nMax)); //Allocate space for nMax bodies
277 for (unsigned a = 0; a < nMax; ++a)
278 (system->bodies+a)->exists = false;
282 * @function System_Destroy
283 * @purpose Destroy all Bodies in a system
284 * @param system - Target system
287 void System_Destroy(System * system)
289 assert(system != NULL);
290 if (system->bodies != NULL)
291 free(system->bodies);
296 * @function System_Step
297 * @purpose Compute single step for all Bodies in System
298 * @param system - Target system
300 void System_Step(System * system)
302 assert(system != NULL);
305 for (unsigned a = 0; a < system->nMax; ++a)
307 if (!(system->bodies+a)->exists)
310 for (unsigned b = 0; b < system->nMax; ++b)
312 if (!(system->bodies+b)->exists)
315 Body_Force(system->bodies+a, system->bodies+b, system); //Update force on bodies[a] due to
321 for (unsigned a = 0; a < system->nMax; ++a)
323 if ((system->bodies+a)->exists)
324 Body_Step(system->bodies+a, system);
330 //Check for collisions
331 //Note: Only allows for one collision per step
333 for (unsigned a = 0; a < system->nMax; ++a)
335 if ((system->bodies+a)->exists == false || (system->bodies+a)->collided)
338 for (unsigned b = 0; b < system->nMax; ++b)
340 if ((system->bodies+b)->exists == false || (system->bodies+b)->collided)
344 if (Body_CheckCollision(system->bodies+a, system->bodies+b, system, &t))
346 Body_HandleCollision(system->bodies+a, system->bodies+b, system, t);
347 (system->bodies+a)->collided = true;
348 (system->bodies+b)->collided = true;
355 if ((system->bodies+a)->x[0] + system->x[0] < 0 || (system->bodies+a)->x[0] + system->x[0] > 640)
356 (system->bodies+a)->v[0] = -(system->bodies+a)->v[0];
357 if ((system->bodies+a)->x[1] + system->x[1] < 0 || (system->bodies+a)->x[1] + system->x[1] > 480)
358 (system->bodies+a)->v[1] = -(system->bodies+a)->v[1];
367 * @function System_AddBody
368 * @purpose Add a Body to a System. Resizes the System's bodies array if necessary.
369 * @param system - Target System
370 * @param m - Mass of the new Body
371 * @param x[...] - position vector for the new Body
372 * @param v[...] - velocity vector for the new Body
373 * @returns Pointer to new Body
376 Body * System_AddBody(System * system, float m, float x[DIMENSIONS], float v[DIMENSIONS])
378 assert(system != NULL);
381 for (a = 0; a < system->nMax; ++a)
383 if (!(system->bodies+a)->exists)
387 if (a >= system->nMax)
389 system->nMax = (system->nMax+1) * 2;
390 system->bodies = (Body*)(realloc(system->bodies, sizeof(Body) * system->nMax));
391 for (unsigned b = a; b < system->nMax; ++b)
392 (system->bodies+b)->exists = false;
395 Body_Init(system->bodies+a, m, x, v);
398 return (system->bodies+a);
403 * @function System_WriteData
404 * @purpose Write all information about System to a file
405 * @param system - Target System
406 * @param file - Target FILE (must be opened for writing)
408 void System_WriteData(System * system, FILE * file)
410 assert(system != NULL && file != NULL);
411 for (unsigned a = 0; a < system->nMax; ++a)
413 if ((system->bodies+a)->exists)
414 Body_WriteData(system->bodies+a, system, file);
420 // --- Graphical functions --- //
425 * @function Body_Draw
426 * @purpose Draw an individual Body
427 * @param body - Target Body
428 * @param system - System to which the Body belongs
431 void Body_Draw(Body * body, System * system)
433 assert(body != NULL && system != NULL && body->exists);
434 int pixel_pos[DIMENSIONS];
435 for (unsigned i = 0; i < DIMENSIONS; ++i)
436 pixel_pos[i] = (int)(system->x[i] - body->x[i]);
438 float rgb[3] = {0.0f, 0.0f, 0.0f};
439 if (body->m >= 100.0)
441 rgb[0] = 0.0f; rgb[1] = 0.0f; rgb[2] = 1.0f;
445 rgb[0] = 0.5f; rgb[1] = 0.5f; rgb[2] = 1.0f;
447 else if (body->m >= 1.0)
449 rgb[0] = 0.8f; rgb[1] = 0.2f; rgb[2] = 0.1f;
451 else if (body->m >= 0.1)
453 rgb[0] = 0.0f; rgb[1] = 0.5f; rgb[2] = 0.0f;
455 else if (body->m >= 0.01)
457 rgb[0] = 0.5f; rgb[1] = 0.5f; rgb[2] = 0.5f;
461 Graphics_Circle(pixel_pos,Body_Radius(body), rgb[0], rgb[1], rgb[2]);
463 //fprintf(stdout, "DEBUG - Body_Draw - Called\n");
467 * @function System_Draw
468 * @purpose Draw a System of Bodies
469 * @param system - Target System to draw
471 void System_Draw(System * system)
473 assert(system != NULL);
474 for (unsigned a = 0; a < system->nMax; ++a)
476 if ((system->bodies+a)->exists)
477 Body_Draw(system->bodies+a, system);