course/semester2/pprog/assignment1/nbody.c

   1 /**
   2  * @file nbody.c
   3  * @author Sam Moore (20503628) 2012
   4  * @purpose N-Body simulator - Definition of simulation functions; single threaded version
   5  */
   6
   7 #include <math.h>
   8 #include <stdlib.h>
   9 #include <stdio.h>
  10 #include <string.h>
  11 #include <time.h>
  12 #include <omp.h>
  13 #include <GL/gl.h>
  14 #include <GL/glut.h>
  15
  16 double previousTime, eyeTheta, eyePhi, eyeRho;
  17 float look[3];
  18 int windowWidth, windowHeight, upY;
  19
  20 double SCALE = 1;
  21
  22 Particle *body;
  23 int N;
  24 int numberOfProcessors=1;
  25
  26
  27
  28 /*
  29  * Prints the body on screen
  30  */
  31 void PrintBody(int i)
  32 {
  33         printf("Body #%d M=%f X=%f Y=%f Z=%f Fx=%f Fy=%f Fz=%f Vx=%f Xy=%f Vz=%f\n",
  34            i, body[i].mass, body[i].X, body[i].Y, body[i].Z, body[i].Fx,
  35            body[i].Fy, body[i].Fz,      body[i].Vx, body[i].Vy, body[i].Vz);
  36 }
  37
  38 /*
  39  * Computing forces
  40  */
  41 void Force(int a) {
  42         double distance;
  43         double con;
  44         double gd;
  45
  46         body[a].Fx = body[a].Fy = body[a].Fz = 0;
  47
  48         for(int b=0; b<N; ++b)
  49                 if( b != a ){
  50                 distance = sqrt(square(body[b].X - body[a].X) + square(body[b].Y
  51                                 - body[a].Y) + square(body[b].Z - body[a].Z));
  52                 con = G * body[a].mass * body[b].mass / square(distance);
  53                 gd = con / distance;
  54                 body[a].Fx += gd * (body[b].X - body[a].X);
  55                 body[a].Fy += gd * (body[b].Y - body[a].Y);
  56                 body[a].Fz += gd * (body[b].Z - body[a].Z);
  57                 }
  58 }
  59
  60 /*
  61  * Compute velocities
  62  */
  63 void Velocity(int a) {
  64                 body[a].Vx += body[a].Fx/body[a].mass * DELTA_T;
  65                 body[a].Vy += body[a].Fy/body[a].mass * DELTA_T;
  66                 body[a].Vz += body[a].Fz/body[a].mass * DELTA_T;
  67 }
  68
  69 /*
  70  * Compute positions
  71  */
  72 void Position(int a) {
  73         body[a].X += body[a].Vx * DELTA_T;
  74         body[a].Y += body[a].Vy * DELTA_T;
  75         body[a].Z += body[a].Vz * DELTA_T;
  76 }
  77
  78 /*
  79  * Main compute function
  80  */
  81 void Compute() {
  82         clock_t start, finish;
  83
  84         start = clock();
  85
  86         for(int a=0; a<N; a++) {
  87                 Force(a); // compute force
  88                 Velocity(a);
  89                 Position(a);
  90         }
  91
  92         finish = clock();
  93         printf("Performance: %ld clock ticks per computation cycle\n", finish-start);
  94
  95         //      PrintBody(0);
  96 //      PrintBody(N-1);
  97 }
  98
  99 /*
 100  * This function is called repeatedly by graphics library. You can consider
 101  * it as main loop in the program.
 102  */
 103 void Animate(void) {
 104
 105         Compute(); //Compute and update new positions for the time step
 106         Display(); // Display needs to be called to redraw the screen
 107 }
 108
 109 /*
 110  * This function is to manipulate with the image
 111  */
 112 void KeyBoard(unsigned char theKey, int mouseX, int mouseY) {
 113         if (theKey == 'x' || theKey == 'X') {
 114                 free(body);
 115                 exit(EXIT_SUCCESS);
 116         }
 117                 if (theKey == 'i' || theKey == 'I') {
 118                         eyePhi -= M_PI / 20;
 119                 if (eyePhi == 0)
 120                                 eyePhi = 2 * M_PI;
 121                 } else if (theKey == 'm' || theKey == 'I') {
 122                         eyePhi += M_PI / 20;
 123                 } else if (theKey == 'j' || theKey == 'J') {
 124                         eyeTheta -= M_PI / 20;
 125                 } else if (theKey == 'k' || theKey == 'K') {
 126                         eyeTheta += M_PI / 20;
 127                 } else if (theKey == ',') {
 128                         eyeRho += 0.5;
 129                 } else if (theKey == '.' || theKey == 'I') {
 130                         eyeRho -= 0.5;
 131                 } else if (theKey == 'w' || theKey == 'W') {
 132                         look[1] += 0.5;
 133                 } else if (theKey == 'z' || theKey == 'Z') {
 134                         look[1] -= 0.5;
 135                 } else if (theKey == 'a' || theKey == 'A') {
 136                         look[0] -= 0.5;
 137                 } else if (theKey == 's' || theKey == 'S') {
 138                         look[0] += 0.5;
 139                 } else if (theKey == '+') {
 140                         SCALE *= 1.1;
 141                 } else if (theKey == '-') {
 142                         SCALE *= 0.9;
 143                 }
 144                 if (sin(eyePhi) > 0) upY = 1;
 145                 else upY = 1;
 146 }
 147
 148 void Reshape(int width, int height) {
 149         double displayRatio = 1.0 * width / height;
 150         windowWidth = width;
 151         windowHeight = height;
 152         glMatrixMode(GL_PROJECTION);
 153     glLoadIdentity();
 154     gluPerspective(VIEW_ANGLE, displayRatio, WORLD_NEAR, WORLD_FAR);
 155         glViewport(0, 0, windowWidth, windowHeight);
 156 }
 157
 158 /*
 159  * This function reads an input file. You can change it if you choose a
 160  * different file format
 161  */
 162 void readFile(char *fileName) {
 163         char line[LINE_SIZE];
 164         char *token;
 165         FILE *file;
 166
 167         file = fopen(fileName, "rt");
 168         N = atoi(fgets(line, LINE_SIZE, file));
 169         body = (Particle*) calloc((size_t)N, sizeof(Particle));
 170
 171         puts("----------------------Initial field-------------------------------");
 172
 173         for (int i=0; i<N; ++i)
 174                 if(fgets(line, LINE_SIZE, file) != NULL){
 175                 token = strtok(line, ",; ");
 176                 body[i].mass = atof(token);
 177                 token = strtok(NULL, ",; ");
 178                 body[i].X = atof(token);
 179                 token = strtok(NULL, ",; ");
 180                 body[i].Y = atof(token);
 181                 token = strtok(NULL, ",; ");
 182                 body[i].Z = atof(token);
 183                 token = strtok(NULL, ",; ");
 184                 body[i].Vx = atof(token);
 185                 token = strtok(NULL, ",; ");
 186                 body[i].Vy = atof(token);
 187                 token = strtok(NULL, ",; ");
 188                 body[i].Vz = atof(token);
 189
 190                 PrintBody(i);
 191                 }
 192         puts("--------------------------------------------------------------");
 193         puts("Use:\n X - exit\n I, J, K, M - rotate\n W, Z, A, S - move to view"
 194          " point\n ./, - zoom in/out\n +/- - scaled zoom in/out\n");
 195
 196         fclose(file);
 197 }
 198
 199
 200 // This is main function. Do not change it.
 201 int main(int argc, char** argv)
 202 {
 203         glutInit(&argc, argv);
 204
 205         if (argc < 2) {
 206                 puts("Please provide the filename, i.e. \'nbody bodiesfield.dat\'");
 207             exit(EXIT_SUCCESS);
 208         }
 209         if (argc == 2) readFile(argv[1]);
 210         if (argc == 3) numberOfProcessors = atoi(argv[2]);
 211
 212         glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
 213         glutInitWindowSize(WIDTH, HEIGHT);
 214         glutInitWindowPosition(POSITION_X, POSITION_Y);
 215         glutCreateWindow("N-Body Parallel");
 216         glutDisplayFunc(Display);
 217         glutIdleFunc(Animate);
 218         glutKeyboardFunc(KeyBoard);
 219         glutReshapeFunc(Reshape);
 220         Init();
 221         glutMainLoop();
 222 }