course/semester2/pprog/assignment1/nbody.c

   1 // nbody.cpp : Template program for N-body
   2 // Copyright UWA, 2012
   3
   4 #define _ANSI_SOURCE
   5
   6 #include <math.h>
   7 #include <stdlib.h>
   8 #include <stdio.h>
   9 #include <string.h>
  10 #include <time.h>
  11 #include <omp.h>
  12 #include <GL/gl.h>
  13 #include <GL/glut.h>
  14
  15 #define M_PI        3.14159265358979323846264338327950288   /* pi */
  16 #define WIDTH 800
  17 #define HEIGHT 800
  18 #define POINT_SIZE 1
  19 #define POSITION_X 112
  20 #define POSITION_Y 20
  21 #define WORLD_LEFT -10000
  22 #define WORLD_RIGHT 10000
  23 #define WORLD_BOTTOM -10000
  24 #define WORLD_TOP 10000
  25 #define VIEW_ANGLE 45
  26 #define RHO 100
  27 #define WORLD_NEAR 0.1
  28 #define WORLD_FAR 1000000
  29 #define BALL_SIZE 0.5
  30 #define REFRESH_RATE 0.001
  31 #define LINE_SIZE 1000
  32 #define G 6.67428E-11
  33 #define DELTA_T 0.05
  34
  35 #define square(x) ((x)*(x))
  36
  37 typedef struct {
  38         double mass;
  39         double X;
  40         double Y;
  41         double Z;
  42         double Vx;
  43         double Vy;
  44         double Vz;
  45         double Fx;
  46         double Fy;
  47         double Fz;
  48         int color;
  49 } Particle;
  50
  51 double previousTime, eyeTheta, eyePhi, eyeRho;
  52 float look[3];
  53 int windowWidth, windowHeight, upY;
  54
  55 double SCALE = 1;
  56
  57 Particle *body;
  58 int N;
  59 int numberOfProcessors=1;
  60
  61 /*
  62  * Initialization of graphics
  63  */
  64 void Init(void) {
  65
  66         glClearColor(1.0,1.0,1.0,0.0);
  67         glColor3f(0.0f, 0.0f, 0.0f);
  68         glPointSize(POINT_SIZE);
  69         glMatrixMode(GL_PROJECTION);
  70         glLoadIdentity();
  71
  72         /*init lighting */
  73
  74         GLfloat mat_specular[] = { 1.0, 1.0, 1.0, 1.0 };
  75         GLfloat mat_shininess[] = { 50.0 };
  76         GLfloat light_position[] = { 1.0, 1.0, 0.0, 0.0 };
  77         glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
  78         glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);
  79         glLightfv(GL_LIGHT0, GL_POSITION, light_position);
  80
  81     glColorMaterial(GL_FRONT,GL_DIFFUSE);                // Set Color Capability
  82
  83         glEnable(GL_LIGHTING);
  84         glEnable(GL_LIGHT0);
  85         glEnable(GL_DEPTH_TEST);
  86
  87     glEnable(GL_COLOR_MATERIAL);                       // Enable color
  88
  89         double displayRatio = 1.0 * WIDTH / HEIGHT;
  90         windowWidth = WIDTH;
  91         windowHeight = HEIGHT;
  92         previousTime = clock();
  93         eyeTheta = 0;
  94     eyePhi = 0.5 * M_PI;
  95     eyeRho = RHO;
  96     upY = 1;
  97         look[0] = 0;
  98         look[1] = 0;
  99         look[2] = 0;
 100     gluPerspective(VIEW_ANGLE, displayRatio, WORLD_NEAR, WORLD_FAR);
 101 }
 102
 103 /*
 104  * This function redraws the screen after the positions of particles
 105  * have been updated
 106  */
 107 void Display(void) {
 108         glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
 109         glMatrixMode(GL_MODELVIEW);
 110     glLoadIdentity();
 111     gluLookAt(eyeRho * sin(eyePhi) * sin(eyeTheta), eyeRho * cos(eyePhi),
 112                     eyeRho * sin(eyePhi) * cos(eyeTheta),
 113                                         look[0], look[1], look[2], 0, upY, 0);
 114
 115         for (int i = 0; i < N; i++) {
 116         //glClearColor(1.0,1.0,1.0,0.0);
 117                 glColor3f(0.0f, body[i].mass/1e11*100, 0.0f);
 118         //glColor3f(1.0f, 0.0f, 0.0f);
 119                 glPushMatrix(); // to save the current matrix
 120                 glTranslated(SCALE*body[i].X, SCALE*body[i].Y, SCALE*body[i].Z);
 121                 glutSolidSphere (BALL_SIZE, 10, 10);
 122                 glPopMatrix(); // restore the previous matrix
 123         }
 124         glFlush();
 125 }
 126
 127 /*
 128  * Prints the body on screen
 129  */
 130 void PrintBody(int i)
 131 {
 132         printf("Body #%d M=%f X=%f Y=%f Z=%f Fx=%f Fy=%f Fz=%f Vx=%f Xy=%f Vz=%f\n",
 133            i, body[i].mass, body[i].X, body[i].Y, body[i].Z, body[i].Fx,
 134            body[i].Fy, body[i].Fz,      body[i].Vx, body[i].Vy, body[i].Vz);
 135 }
 136
 137 /*
 138  * Computing forces
 139  */
 140 void Force(int a) {
 141         double distance;
 142         double con;
 143         double gd;
 144
 145         body[a].Fx = body[a].Fy = body[a].Fz = 0;
 146
 147         for(int b=0; b<N; ++b)
 148                 if( b != a ){
 149                 distance = sqrt(square(body[b].X - body[a].X) + square(body[b].Y
 150                                 - body[a].Y) + square(body[b].Z - body[a].Z));
 151                 con = G * body[a].mass * body[b].mass / square(distance);
 152                 gd = con / distance;
 153                 body[a].Fx += gd * (body[b].X - body[a].X);
 154                 body[a].Fy += gd * (body[b].Y - body[a].Y);
 155                 body[a].Fz += gd * (body[b].Z - body[a].Z);
 156                 }
 157 }
 158
 159 /*
 160  * Compute velocities
 161  */
 162 void Velocity(int a) {
 163                 body[a].Vx += body[a].Fx/body[a].mass * DELTA_T;
 164                 body[a].Vy += body[a].Fy/body[a].mass * DELTA_T;
 165                 body[a].Vz += body[a].Fz/body[a].mass * DELTA_T;
 166 }
 167
 168 /*
 169  * Compute positions
 170  */
 171 void Position(int a) {
 172         body[a].X += body[a].Vx * DELTA_T;
 173         body[a].Y += body[a].Vy * DELTA_T;
 174         body[a].Z += body[a].Vz * DELTA_T;
 175 }
 176
 177 /*
 178  * Main compute function
 179  */
 180 void Compute() {
 181         clock_t start, finish;
 182
 183         start = clock();
 184
 185         for(int a=0; a<N; a++) {
 186                 Force(a); // compute force
 187                 Velocity(a);
 188                 Position(a);
 189         }
 190
 191         finish = clock();
 192         printf("Performance: %ld clock ticks per computation cycle\n", finish-start);
 193
 194         //      PrintBody(0);
 195 //      PrintBody(N-1);
 196 }
 197
 198 /*
 199  * This function is called repeatedly by graphics library. You can consider
 200  * it as main loop in the program.
 201  */
 202 void Animate(void) {
 203
 204         Compute(); //Compute and update new positions for the time step
 205         Display(); // Display needs to be called to redraw the screen
 206 }
 207
 208 /*
 209  * This function is to manipulate with the image
 210  */
 211 void KeyBoard(unsigned char theKey, int mouseX, int mouseY) {
 212         if (theKey == 'x' || theKey == 'X') {
 213                 free(body);
 214                 exit(EXIT_SUCCESS);
 215         }
 216                 if (theKey == 'i' || theKey == 'I') {
 217                         eyePhi -= M_PI / 20;
 218                 if (eyePhi == 0)
 219                                 eyePhi = 2 * M_PI;
 220                 } else if (theKey == 'm' || theKey == 'I') {
 221                         eyePhi += M_PI / 20;
 222                 } else if (theKey == 'j' || theKey == 'J') {
 223                         eyeTheta -= M_PI / 20;
 224                 } else if (theKey == 'k' || theKey == 'K') {
 225                         eyeTheta += M_PI / 20;
 226                 } else if (theKey == ',') {
 227                         eyeRho += 0.5;
 228                 } else if (theKey == '.' || theKey == 'I') {
 229                         eyeRho -= 0.5;
 230                 } else if (theKey == 'w' || theKey == 'W') {
 231                         look[1] += 0.5;
 232                 } else if (theKey == 'z' || theKey == 'Z') {
 233                         look[1] -= 0.5;
 234                 } else if (theKey == 'a' || theKey == 'A') {
 235                         look[0] -= 0.5;
 236                 } else if (theKey == 's' || theKey == 'S') {
 237                         look[0] += 0.5;
 238                 } else if (theKey == '+') {
 239                         SCALE *= 1.1;
 240                 } else if (theKey == '-') {
 241                         SCALE *= 0.9;
 242                 }
 243                 if (sin(eyePhi) > 0) upY = 1;
 244                 else upY = 1;
 245 }
 246
 247 void Reshape(int width, int height) {
 248         double displayRatio = 1.0 * width / height;
 249         windowWidth = width;
 250         windowHeight = height;
 251         glMatrixMode(GL_PROJECTION);
 252     glLoadIdentity();
 253     gluPerspective(VIEW_ANGLE, displayRatio, WORLD_NEAR, WORLD_FAR);
 254         glViewport(0, 0, windowWidth, windowHeight);
 255 }
 256
 257 /*
 258  * This function reads an input file. You can change it if you choose a
 259  * different file format
 260  */
 261 void readFile(char *fileName) {
 262         char line[LINE_SIZE];
 263         char *token;
 264         FILE *file;
 265
 266         file = fopen(fileName, "rt");
 267         N = atoi(fgets(line, LINE_SIZE, file));
 268         body = (Particle*) calloc((size_t)N, sizeof(Particle));
 269
 270         puts("----------------------Initial field-------------------------------");
 271
 272         for (int i=0; i<N; ++i)
 273                 if(fgets(line, LINE_SIZE, file) != NULL){
 274                 token = strtok(line, ",; ");
 275                 body[i].mass = atof(token);
 276                 token = strtok(NULL, ",; ");
 277                 body[i].X = atof(token);
 278                 token = strtok(NULL, ",; ");
 279                 body[i].Y = atof(token);
 280                 token = strtok(NULL, ",; ");
 281                 body[i].Z = atof(token);
 282                 token = strtok(NULL, ",; ");
 283                 body[i].Vx = atof(token);
 284                 token = strtok(NULL, ",; ");
 285                 body[i].Vy = atof(token);
 286                 token = strtok(NULL, ",; ");
 287                 body[i].Vz = atof(token);
 288
 289                 PrintBody(i);
 290                 }
 291         puts("--------------------------------------------------------------");
 292         puts("Use:\n X - exit\n I, J, K, M - rotate\n W, Z, A, S - move to view"
 293          " point\n ./, - zoom in/out\n +/- - scaled zoom in/out\n");
 294
 295         fclose(file);
 296 }
 297
 298
 299 // This is main function. Do not change it.
 300 int main(int argc, char** argv)
 301 {
 302         glutInit(&argc, argv);
 303
 304         if (argc < 2) {
 305                 puts("Please provide the filename, i.e. \'nbody bodiesfield.dat\'");
 306             exit(EXIT_SUCCESS);
 307         }
 308         if (argc == 2) readFile(argv[1]);
 309         if (argc == 3) numberOfProcessors = atoi(argv[2]);
 310
 311         glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
 312         glutInitWindowSize(WIDTH, HEIGHT);
 313         glutInitWindowPosition(POSITION_X, POSITION_Y);
 314         glutCreateWindow("N-Body Parallel");
 315         glutDisplayFunc(Display);
 316         glutIdleFunc(Animate);
 317         glutKeyboardFunc(KeyBoard);
 318         glutReshapeFunc(Reshape);
 319         Init();
 320         glutMainLoop();
 321 }