course/semester2/pprog/assignment0/nbody.c

   1 /**
   2  * @file nbody.c
   3  * @purpose Implementations for N-Body simulator
   4  * @author Sam Moore
   5  * @date August 2012
   6  */
   7
   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
  12
  13 // --- Body functions --- //
  14
  15 /**
  16  * @function Body_Init
  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
  22  */
  23 void Body_Init(Body * body, float m, float x[DIMENSIONS], float v[DIMENSIONS])
  24 {
  25         assert(body != NULL);
  26         body->m = m;
  27         for (unsigned i = 0; i < DIMENSIONS; ++i)
  28         {
  29                 body->x[i] = x[i];
  30                 body->v[i] = v[i];
  31                 body->F[i] = 0.00;
  32
  33                 body->x_old[i] = x[i];
  34                 body->v_old[i] = v[i];
  35         }
  36
  37         body->exists = true;
  38         body->collided = false;
  39 }
  40
  41 /**
  42  * @function Body_Destroy
  43  * @purpose Destroy a Body
  44  * @param body - Target Body
  45  * @param system - System to which the Body belongs
  46  */
  47 void Body_Destroy(Body * body, System * system)
  48 {
  49         assert(body != NULL);
  50         body->exists = false;
  51         system->n -= 1;
  52 }
  53
  54 /**
  55  * @function Body_Step
  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)
  59  */
  60 void Body_Step(Body * body, System * system)
  61 {
  62         assert(body != NULL && system != NULL && body->exists);
  63         body->collided = false;
  64         for (unsigned i = 0; i < DIMENSIONS; ++i)
  65         {
  66
  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];
  71
  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
  75         }
  76
  77 }
  78
  79 /**
  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)
  85  */
  86 void Body_Force(Body * A, Body * B, System * system)
  87 {
  88
  89         assert(A != NULL && B != NULL && system != NULL && A->exists && B->exists);
  90         if (A == B)
  91                 return;
  92
  93         float r2 = 0.00;
  94         for (unsigned i = 0; i < DIMENSIONS; ++i)
  95                 r2 += powf((B->x[i] - A->x[i]), 2);
  96
  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]);
 100
 101 }
 102
 103 /**
 104  * @function Body_Radius
 105  * @purpose Calculate the radius of a body
 106  * @param body - Target Body
 107  * @returns The radius of body
 108  */
 109 float Body_Radius(Body * body)
 110 {
 111         assert(body != NULL && body->exists);
 112         return 8.0 * sqrt(body->m);
 113 }
 114
 115
 116 /**
 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
 124  */
 125 bool Body_CheckCollision(Body * A, Body * B, System * system, float * t)
 126 {
 127         assert(A != NULL && B != NULL && system != NULL && A->exists && B->exists);
 128
 129         if (A == B)
 130                 return false;
 131
 132         // Determine time and position of closest approach
 133         float x_rel[DIMENSIONS];
 134         float v_rel[DIMENSIONS];
 135         float coeff_t = 0.0;
 136         float c = 0.0;
 137         for (unsigned i = 0; i < DIMENSIONS; ++i)
 138         {
 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];
 143         }
 144
 145         *t = -c / coeff_t;
 146         if (*t <= system->dt/100.0 || *t > system->dt)
 147                 return false;
 148
 149         float r2 = 0.0;
 150         for (unsigned i = 0; i < DIMENSIONS; ++i)
 151         {
 152                 x_rel[i] += *t * v_rel[i];
 153                 r2 += pow(x_rel[i], 2);
 154         }
 155
 156         return (r2 <= pow(Body_Radius(A), 2) + pow(Body_Radius(B), 2));
 157 }
 158
 159 /**
 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
 166  */
 167 void Body_HandleCollision(Body * A, Body * B, System * system, float t)
 168 {
 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)
 173         {
 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);
 176
 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
 179
 180         }
 181
 182         //For cor = 0, bodies always merge; B is travelling with A
 183         if (system->cor == 0.0)
 184         {
 185                 A->m += B->m; //Increase A's mass
 186                 Body_Destroy(B, system); //Destroy B
 187                 return;
 188         }
 189
 190         //Calculate B's new velocity and position
 191         for (unsigned i = 0; i < DIMENSIONS; ++i)
 192         {
 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);
 195
 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
 198
 199         }
 200
 201         //Determine whether Body B has escaped Body A; if it hasn't, merge the two bodies
 202
 203
 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]);
 205
 206         float v2 = 0.0;
 207         float r2 = 0.0;
 208         for (unsigned i = 0; i < DIMENSIONS; ++i)
 209         {
 210                 v2 += pow(B->v[i], 2);
 211                 r2 += pow(B->x[i] - A->x[i], 2);
 212         }
 213
 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))
 216                 return;
 217
 218         //Inelastically merge the two bodies
 219         for (unsigned i = 0; i < DIMENSIONS; ++i)
 220         {
 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);
 223         }
 224         A->m += B->m;
 225         Body_Destroy(B, system);
 226
 227 }
 228
 229 /**
 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
 235  */
 236 void Body_WriteData(Body * body, System * system, FILE * file)
 237 {
 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]);
 246
 247         fprintf(file,"\n");
 248
 249 }
 250
 251 // --- System functions --- //
 252
 253 /**
 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)
 262  */
 263 void System_Init(System * system, float x[DIMENSIONS], float dt, float G, float cor, unsigned nMax)
 264 {
 265         assert(system != NULL);
 266         for (unsigned i = 0; i < DIMENSIONS; ++i)
 267                 system->x[i] = x[i];
 268
 269
 270         system->dt = dt;
 271         system->G = G;
 272         system->cor = cor;
 273         system->nMax = nMax;
 274         system->n = 0;
 275
 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;
 279 }
 280
 281 /**
 282  * @function System_Destroy
 283  * @purpose Destroy all Bodies in a system
 284  * @param system - Target system
 285  */
 286
 287 void System_Destroy(System * system)
 288 {
 289         assert(system != NULL);
 290         if (system->bodies != NULL)
 291                 free(system->bodies);
 292         system->nMax = 0;
 293 }
 294
 295 /**
 296  * @function System_Step
 297  * @purpose Compute single step for all Bodies in System
 298  * @param system - Target system
 299  */
 300 void System_Step(System * system)
 301 {
 302         assert(system != NULL);
 303
 304         // Update forces
 305         for (unsigned a = 0; a < system->nMax; ++a)
 306         {
 307                 if (!(system->bodies+a)->exists)
 308                         continue;
 309
 310                 for (unsigned b = 0; b < system->nMax; ++b)
 311                 {
 312                         if (!(system->bodies+b)->exists)
 313                                 continue;
 314
 315                         Body_Force(system->bodies+a, system->bodies+b, system); //Update force on bodies[a] due to
 316                 }
 317
 318         }
 319
 320         // Perform steps
 321         for (unsigned a = 0; a < system->nMax; ++a)
 322         {
 323                 if ((system->bodies+a)->exists)
 324                         Body_Step(system->bodies+a, system);
 325
 326
 327
 328         }
 329
 330         //Check for collisions
 331         //Note: Only allows for one collision per step
 332
 333         for (unsigned a = 0; a < system->nMax; ++a)
 334         {
 335                 if ((system->bodies+a)->exists == false || (system->bodies+a)->collided)
 336                         continue;
 337
 338                 for (unsigned b = 0; b < system->nMax; ++b)
 339                 {
 340                         if ((system->bodies+b)->exists == false || (system->bodies+b)->collided)
 341                                 continue;
 342
 343                         float t = 0.0;
 344                         if (Body_CheckCollision(system->bodies+a, system->bodies+b, system, &t))
 345                         {
 346                                 Body_HandleCollision(system->bodies+a, system->bodies+b, system, t);
 347                                 (system->bodies+a)->collided = true;
 348                                 (system->bodies+b)->collided = true;
 349                                 break;
 350                         }
 351
 352                 }
 353
 354                 if ((system->bodies+a)->x[0] + system->x[0] < 0 || (system->bodies+a)->x[0] + system->x[0] > 640)
 355                         (system->bodies+a)->v[0] = -(system->bodies+a)->v[0];
 356                 if ((system->bodies+a)->x[1] + system->x[1] < 0 || (system->bodies+a)->x[1] + system->x[1] > 480)
 357                         (system->bodies+a)->v[1] = -(system->bodies+a)->v[1];
 358
 359         }
 360
 361
 362 }
 363
 364 /**
 365  * @function System_AddBody
 366  * @purpose Add a Body to a System. Resizes the System's bodies array if necessary.
 367  * @param system - Target System
 368  * @param m - Mass of the new Body
 369  * @param x[...] - position vector for the new Body
 370  * @param v[...] - velocity vector for the new Body
 371  * @returns Pointer to new Body
 372  */
 373
 374 Body * System_AddBody(System * system, float m, float x[DIMENSIONS], float v[DIMENSIONS])
 375 {
 376         assert(system != NULL);
 377
 378         unsigned a = 0;
 379         for (a = 0; a < system->nMax; ++a)
 380         {
 381                 if (!(system->bodies+a)->exists)
 382                         break;
 383         }
 384
 385         if (a >= system->nMax)
 386         {
 387                 system->nMax = (system->nMax+1) * 2;
 388                 system->bodies = (Body*)(realloc(system->bodies, sizeof(Body) * system->nMax));
 389                 for (unsigned b = a; b < system->nMax; ++b)
 390                         (system->bodies+b)->exists = false;
 391         }
 392
 393         Body_Init(system->bodies+a, m, x, v);
 394         system->n += 1;
 395
 396         return (system->bodies+a);
 397
 398 }
 399
 400 /**
 401  * @function System_WriteData
 402  * @purpose Write all information about System to a file
 403  * @param system - Target System
 404  * @param file - Target FILE (must be opened for writing)
 405  */
 406 void System_WriteData(System * system, FILE * file)
 407 {
 408         assert(system != NULL && file != NULL);
 409         for (unsigned a = 0; a < system->nMax; ++a)
 410         {
 411                 if ((system->bodies+a)->exists)
 412                         Body_WriteData(system->bodies+a, system, file);
 413         }
 414         fprintf(file, "\n");
 415
 416 }
 417
 418 // --- Graphical functions --- //
 419
 420 #ifdef _GRAPHICS_H
 421
 422 /**
 423  * @function Body_Draw
 424  * @purpose Draw an individual Body
 425  * @param body - Target Body
 426  * @param system - System to which the Body belongs
 427  */
 428
 429 void Body_Draw(Body * body, System * system)
 430 {
 431         assert(body != NULL && system != NULL && body->exists);
 432         int pixel_pos[DIMENSIONS];
 433         for (unsigned i = 0; i < DIMENSIONS; ++i)
 434                 pixel_pos[i] = (int)(system->x[i] - body->x[i]);
 435
 436         float rgb[3] = {0.0f, 0.0f, 0.0f};
 437         if (body->m >= 100.0)
 438         {
 439                 rgb[0] = 0.0f; rgb[1] = 0.0f; rgb[2] = 1.0f;
 440         }
 441         if (body->m >= 10.0)
 442         {
 443                 rgb[0] = 0.5f; rgb[1] = 0.5f; rgb[2] = 1.0f;
 444         }
 445         else if (body->m >= 1.0)
 446         {
 447                 rgb[0] = 0.8f; rgb[1] = 0.2f; rgb[2] = 0.1f;
 448         }
 449         else if (body->m >= 0.1)
 450         {
 451                 rgb[0] = 0.0f; rgb[1] = 0.5f; rgb[2] = 0.0f;
 452         }
 453         else if (body->m >= 0.01)
 454         {
 455                 rgb[0] = 0.5f; rgb[1] = 0.5f; rgb[2] = 0.5f;
 456         }
 457
 458
 459         Graphics_Circle(pixel_pos,Body_Radius(body), rgb[0], rgb[1], rgb[2]);
 460
 461         //fprintf(stdout, "DEBUG - Body_Draw - Called\n");
 462 }
 463
 464 /**
 465  * @function System_Draw
 466  * @purpose Draw a System of Bodies
 467  * @param system - Target System to draw
 468  */
 469 void System_Draw(System * system)
 470 {
 471         assert(system != NULL);
 472         for (unsigned a = 0; a < system->nMax; ++a)
 473         {
 474                 if ((system->bodies+a)->exists)
 475                         Body_Draw(system->bodies+a, system);
 476         }
 477 }
 478
 479 /**
 480  * @function Process_Events
 481  * @purpose Handle any SDL events recieved.
 482  */
 483 void Process_Events()
 484 {
 485         SDL_Event event;
 486         while (SDL_PollEvent(&event))
 487         {
 488                 switch(event.type)
 489                 {
 490                         case SDL_QUIT:
 491                                 exit(EXIT_SUCCESS);
 492                                 break;
 493                 }
 494         }
 495         //SDL_Delay(1);
 496 }
 497
 498 #endif //_GRAPHICS_H
 499
 500 //EOF
 501
 502