// http://www.cs.unc.edu/~andrewz/comp203/hw1/index.html /* gravitytestOpengl.c 3 particles at (0,0), (5,0), (2.5,4.33) with mass 100000000 ./lgcc gravitytest1Opengl.c ./a.out Here's the file for lgcc: #!/bin/sh gcc $1 -lm -L/usr/lib -lGLU -lGL -lglut -L/usr/X11R6/lib -lX11 -lXext -lXi -lXmu Use chmod 755 lgcc inorder to make this executable */ #include #include #include #include #include //#define NUM_PARTICLES 1000 //#define NUM_ITERATIONS 300 #define NUM_PARTICLES 3 #define NUM_ITERATIONS 100 typedef struct ParticleStruct { double m; // particle mass double x; // x-coordinate double y; // y-coordinate double vx; // x-velocity double vy; // y-velocity double fx; // force accumulation on x double fy; // force accumulation on y } Particle; /*********************************************************************/ /** Andrew Zaferakis **/ /** UNC Chapel Hill **/ /** COMP 203 Spring 2000 **/ /** Programming Homework 1 **/ /** N-body particle system **/ /** File: nbody_particle.cpp **/ /*********************************************************************/ #define PI 3.1415927 #define G 6.673E-11 #define EPSILON 0.0001 #define cube(a) a*a*a #define square(a) a*a void UpdateSystemFast(Particle P[], int numParticles, double dt) { double xdist, ydist, delta_vx, delta_vy, scalar; int i,j; for (i = 0; i < numParticles; i++) P[i].fx = P[i].fy = 0; for (i = 0; i < numParticles; i++) { for (j = i+1; j < numParticles; j++) { xdist = P[j].x - P[i].x; ydist = P[j].y - P[i].y; scalar = G * P[i].m * P[j].m / (cube(sqrt(square(xdist) + square(ydist))) + EPSILON); P[i].fx += scalar * xdist; P[i].fy += scalar * ydist; P[j].fx -= scalar * xdist; // must reverse the subtraction for jth P[j].fy -= scalar * ydist; } } for (i = 0; i < numParticles; i++) { delta_vx = P[i].fx * dt / P[i].m; // Finalize the ith body, not the jth delta_vy = P[i].fy * dt / P[i].m; P[i].x += (P[i].vx + delta_vx / 2) * dt; P[i].y += (P[i].vy + delta_vy / 2) * dt; P[i].vx += delta_vx; P[i].vy += delta_vy; } } void UpdateSystemSlow(Particle P[], int numParticles, double dt) { double xdist, ydist, delta_vx, delta_vy, scalar; int i, j; for (i = 0; i < numParticles; i++) { P[i].fx = P[i].fy = 0; for (j = 0; j < numParticles; j++) { xdist = P[j].x - P[i].x; ydist = P[j].y - P[i].y; scalar = G * P[i].m * P[j].m / (cube(sqrt(square(xdist) + square(ydist))) + EPSILON); P[i].fx += scalar * xdist; P[i].fy += scalar * ydist; } } for (i = 0; i < numParticles; i++) { delta_vx = P[i].fx * dt / P[i].m; delta_vy = P[i].fy * dt / P[i].m; P[i].x += (P[i].vx + delta_vx / 2) * dt; P[i].y += (P[i].vy + delta_vy / 2) * dt; P[i].vx += delta_vx; P[i].vy += delta_vy; printf("\t%f\t%f", P[i].x, P[i].y); } printf("\n"); } void InitializeSystem(Particle P[], int numParticles, int config) { int i; double tmp; if (config == 1) for (i = 0; i < numParticles; i++) { P[i].m = 1.0; P[i].x = P[i].y = ((double)i) / ((double)numParticles); P[i].vx = P[i].vy = P[i].fx = P[i].fy = 0; } else if (config == 2) for (i = 0; i < numParticles; i++) { tmp = ((double)i) / ((double)numParticles); P[i].m = ((double)(numParticles - i))/ ((double)numParticles); P[i].x = 0.5 * (1 + tmp) * cos(2*PI*tmp); P[i].y = 0.5 * (1 + tmp) * sin(2*PI*tmp); P[i].vx = P[i].vy = P[i].fx = P[i].fy = 0; } P[0].m = 100000000; P[0].x = 0.0; P[0].y = 0.0; P[0].vx = 0.0; P[0].vx = 0.0; P[1].m = 100000000; P[1].x = 5.0; P[1].y = 0.0; P[1].vx = 0.0; P[1].vx = 0.0; P[2].m = 100000000; P[2].x = 2.5; P[2].y = 4.33; P[2].vx = 0.0; P[2].vx = 0.0; } int CheckMomentum(Particle P[], int numParticles) { double momentum_x = 0; double momentum_y = 0; int i, j; for (i = 0; i < numParticles; i++) { momentum_x += P[i].m * P[i].vx; momentum_y += P[i].m * P[i].vy; } return ((momentum_x <= 0.01) && (momentum_x >= -0.01) && \ (momentum_y <= 0.01) && (momentum_y >= -0.01)) ? 1 : 0; } /*********************************************************************/ /** Andrew Zaferakis **/ /** UNC Chapel Hill **/ /** COMP 203 Spring 2000 **/ /** Programming Homework 1 **/ /** N-body particle system **/ /** File: draw_system.cpp **/ /*********************************************************************/ //NbodyParticleSystem p(1000); // NbodyParticleSystem p(3); C++ version Particle P[NUM_PARTICLES]; int Pause = 0; /*********************************************************************/ /** MyInit **/ /** Input: **/ /** **/ /** Output: Creates the vertex and hull storage, random vertex **/ /*********************************************************************/ void MyInit() { InitializeSystem(P, NUM_PARTICLES, 2); glPointSize(4.0f); } /*********************************************************************/ /** Main Display Function, calls to draw the scene **/ /*********************************************************************/ void myDisplay() { if (Pause) return; int i; glClear(GL_COLOR_BUFFER_BIT); glColor3f(1.0f, 1.0f, 1.0f); glBegin(GL_POINTS); for (i = 0; i < NUM_PARTICLES; i++) glVertex2f(P[i].x, P[i].y); glEnd(); glutSwapBuffers(); // p.UpdateSystemSlow(100); UpdateSystemSlow(P, NUM_PARTICLES, 1); } /*********************************************************************/ /** OpenGL Callback Functions for resize, keys, mouse, etc. **/ /*********************************************************************/ void myReshape(int w, int h) { glViewport(0,0,w,h); glMatrixMode (GL_PROJECTION); glLoadIdentity (); glOrtho(-10,10, -10,10, -1,1); } void KeyboardFunc(unsigned char Key, int x, int y) { if (Key==27) exit(0); // ESC if (Key=='1') InitializeSystem(P, NUM_PARTICLES, 1); if (Key=='2') InitializeSystem(P, NUM_PARTICLES, 2); if (Key=='3') InitializeSystem(P, NUM_PARTICLES, 3); if (Key==' ') Pause = 1-Pause; myDisplay(); } void myIdle() { glutPostRedisplay(); } /*********************************************************************/ /** Main **/ /*********************************************************************/ int main(int argc, char** argv) { glutInit(&argc, argv); glClearColor(0.0, 0.0, 0.0, 0.0); glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB); glutInitWindowSize(400, 400); glutInitWindowPosition(180,100); glutCreateWindow("Andrew Zaferakis - Nbody Particle Simulation"); MyInit(); glutDisplayFunc(myDisplay); glutReshapeFunc(myReshape); glutIdleFunc(myIdle); glutKeyboardFunc(KeyboardFunc); glutMainLoop(); return 0; }