/* simulation of N-body problem (gravitational forces) - Chuck Liang compile with (on deepdish cluster): mpecc grav3d.c -o grav3d Please run only on 1 processor until I figure out how to fix the X11 forwarding (os/network problem). */ #include #include #include #include #include "mpe_graphics.h" #include "mpe.h" #include "threed.c" /* Physics stuff: Given two objects with masses M1 and M2, the gravitationl force between these objects is F = (G*M1*M2)/(d*d), where d is the distance between the center of gravity of the two objects. Newton also says that F=MA, so A1 = G*M2/(d*d) is the accelartion caused by the pull of M2 on M1, and A2 = G*M1/(d*d) is that caused by M1 on M2. The gravitational constant G is 6.673e-11 meters^3/(kg*sec^2) in our universe (yours may be different :) Acceleration is given in meters/sec^2 Distances are in meters and mass in kg. */ /* universal constants */ // screen coordinates will be in units of 10e6 meters //#define UDIST 1 #define UDIST 1.0e8 // time unit between frames - in seconds //#define UTIME 1 //#define UTIME 60 #define UTIME 360 //#define UTIME 3600 //#define UTIME 86400 //#define UTIME 2629800 //#define UTIME 31557600 //#define UTIME 31557600.0e3 // wait for mouse click between frames? #define SYNCH 0 #define TRACE 0 #define ITERS 10000000 // window width, height: #define SW 1024 #define SH 768 #define SD 512 // wrapround? #define WRAPAROUND 1 // number of particles #define PSIZE 200 // gravitational constant #define G 6.673e-11 // mass of earth in kg, radius in meters, mass of sun #define MEARTH 5.9742e24 #define REARTH 6.4e6 #define MSUN 1.99e30 // body will be a vector of 6 doubles: xcord, ycord, mass, xvel, yvel, radius. // radius is used only for drawing (unit in pixels) #define XCORD(AV) AV[0] #define YCORD(AV) AV[1] #define ZCORD(AV) AV[2] #define MASS(AV) AV[3] #define XVEL(AV) AV[4] #define YVEL(AV) AV[5] #define ZVEL(AV) AV[6] #define RADIUS(AV) AV[7] // not used for now: if radius is <0, particle has been destroyed/merged #define DESTROYED(AV) (AV[7]<0) double BODY[PSIZE][8]; MPE_XGraph graph; int minb, maxb; // partition process is responsible for MPE_Color mycolor; /**** grid ****/ #define DROW 24 #define DCOL 32 #define DSPN 12 int grid[DCOL][DROW][DSPN]; // contains indices into BODY array // calculate effect of force on A and B bodies. // returns if bodies have collided (not used now) //(i is local, j maybe remote) int calcforce(int i, int j) { double *A, *B; double d2, dx, dy, dist, sx, sy, sm, sz, dz; double accAX, accAY, accAZ, deeper; // accelerations A = BODY[i]; B = BODY[j]; dx = XCORD(B)-XCORD(A); dy = YCORD(B)-YCORD(A); dz = ZCORD(B)-ZCORD(A); d2 = (dx*dx) + (dy*dy) + (dz*dz); if (d2==0) d2=1e-100; dist = sqrt(d2); // now determine if bodies have collide: if (ZCORD(A)>ZCORD(B)) deeper = ZCORD(A); else deeper=ZCORD(B); if (dist < UDIST*(dscale(RADIUS(A) + RADIUS(B), deeper/UDIST))) { // destroy smaller i // printf("collusion!\n"); /* if (i1) srandom(rank+atoi(argv[1])); // seed random number generator if (argc>1) tiltR = toRadians(atoi(argv[1])); // MPE_Open_graphics(&graph,MPI_COMM_WORLD,(char*)0,0,0,1024,768, // MPE_GRAPH_INDEPENDENT); minb = rank * (PSIZE/size); // assuming no remainder maxb = minb + (PSIZE/size) - 1; tiltR=toRadians(220); ZR = 0.5; //SR = 0.8; XDIM = 800; YDIM = 640; ZDIM = 400; xshift = (int)(ZDIM * cos(tiltR)); yshift = (int)(ZDIM * sin(tiltR)); sprintf(display,"%s","10.1.0.2:10.0"); MPE_Open_graphics(&graph,MPI_COMM_WORLD,(char*)0,1,1, XDIM+(int)abs(xshift)+10,YDIM+(int)abs(yshift)+10, MPE_GRAPH_INDEPENDENT); drawbox(graph,1,1,1,XDIM,YDIM,ZDIM,MPE_GREEN); initialize(rank); simulate(rank); MPI_Finalize(); exit(0); }