// http://www.cs.unc.edu/~andrewz/comp203/hw1/index.html /* gravitytest1Opengl.cpp 3 particles at (0,0), (5,0), (2.5,4.33) with mass 100000000 ./lg++ gravitytest1Opengl.cpp ./a.out Here's the file for lg++: #!/bin/sh g++ $1 -lm -L/usr/lib -lGLU -lGL -lglut -L/usr/X11R6/lib -lX11 -lXext -lXi -lXmu Use chmod 755 lg++ 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 _particle { 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; class NbodyParticleSystem { public: particle *P; // Dynamic array of particles double G; // Gravity constant int numParticles; // Current number of particles in array NbodyParticleSystem(int n) : numParticles(n) { P = new particle[n]; } ~NbodyParticleSystem() { delete[] P; } void UpdateSystemFast(double dt); // Delta t void UpdateSystemSlow(double dt); // Delta t void InitializeSystem(int config); // Setup the particles int CheckMomentum(); // Momentum of the System ~= 0 }; /*********************************************************************/ /** 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 NbodyParticleSystem::UpdateSystemFast(double dt) { double xdist, ydist, delta_vx, delta_vy, scalar; for (int i = 0; i < numParticles; i++) P[i].fx = P[i].fy = 0; for (int i = 0; i < numParticles; i++) { for (int 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 (int 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 NbodyParticleSystem::UpdateSystemSlow(double dt) { double xdist, ydist, delta_vx, delta_vy, scalar; for (int i = 0; i < numParticles; i++) { P[i].fx = P[i].fy = 0; for (int 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 (int 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 NbodyParticleSystem::InitializeSystem(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 NbodyParticleSystem::CheckMomentum() { double momentum_x = 0; double momentum_y = 0; for (int 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); int Pause = 0; /*********************************************************************/ /** MyInit **/ /** Input: **/ /** **/ /** Output: Creates the vertex and hull storage, random vertex **/ /*********************************************************************/ void MyInit() { p.InitializeSystem(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 < p.numParticles; i++) glVertex2f(p.P[i].x, p.P[i].y); glEnd(); glutSwapBuffers(); // p.UpdateSystemSlow(100); p.UpdateSystemSlow(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') p.InitializeSystem(1); if (Key=='2') p.InitializeSystem(2); if (Key=='3') p.InitializeSystem(3); if (Key==' ') Pause = 1-Pause; myDisplay(); } void myIdle() { glutPostRedisplay(); } /*********************************************************************/ /** Main **/ /*********************************************************************/ int main() { 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; }