## NBody evaluation, is my model correct?

Hello all,

I've written a very simple nbody solver in c. Basically, I just went through an array of particles and solved F = G(m1)(m2) / r ^2

The output is here: (sorry for the windows format)
http://www.crepinc.com/files/galax.wmv

As you can see, some of the particles don't do what one would expect....

So far the only possible reason I can come up for is that I'm solving the equations at a set dT. Thus if the dT is too large, "stuff happens" between the two different solves. But if dT is too small, the whole thing basically blows up.

Am I overlooking something, or is my math ok and my code is just bad?

Thanks a lot,

-Jack Carrozzo
http://www.crepinc.com/

 PhysOrg.com science news on PhysOrg.com >> Front-row seats to climate change>> Attacking MRSA with metals from antibacterial clays>> New formula invented for microscope viewing, substitutes for federally controlled drug

Recognitions:
Gold Member
 Quote by crepincdotcom Am I overlooking something, or is my math ok and my code is just bad?
Can't tell without seeing the code, especially the lines that contain the time step variable.

 Recognitions: Gold Member The time step can intrdouce such problems - you tend to find that energy isn't conserved. I would tend to put in energy conservation 'by hand'. However, there does seem to be another problem with your animation - it takes some time for all of the particles to start moving, and then some of them stop moving. I would guess that this is due to a problem with the code.

## NBody evaluation, is my model correct?

OK, I will post the code at 8:00, sorry I can't get it before then.

As to the particles that stop, It usually doesn't do that. I think I've decreased dT too much or something.... I didn't make any hardcore code changes.

But, I don't want to bother you guys with code issues. In terms of theory, does evaluating F and then making a sum of all forces do the trick? Or do I need to find a way to use the differentia l equations and such?

-Jack Carrozzo
http://www.crepinc.com/

 Hello, Try to work with different DT. Best of Luck. Zaheer A
 OK, at the bottom of this post is the code. I got the model to look FAIRLY accurate by playing with the dT. But if I go higher or lower, the whole system no longer works accuratly. IS there a logical explanation for this? The other issue is that of the Y axis. You can see in the "MoveParticles" function that when taking the sum of all y-forces and adding them to the velocity, I multiply the acceleration by 25. I have no idea why this needs to happen. All I know is that, whithout that line, the simulation seems to "favor" the X axis, with the particles hardly moving in the Y. Finally, one should notice that I removed some 0's from the constant G. This is because multiplying something X 10^ -11 with a number (Mass) like 10 ^ 6, there are a lot of extra spaces in there. I was told it was easier and more accurate in the evaluation of that code if I simply lopped off some zero's. So, any mass you see is actually large by a factor of 10^5. Here's the code: (it will only let me post little bits, the whole thing is at crepinc.com/files/nbody.c) Code: #define NUMPARTS 100 //how many particles #define GRID_X 640 #define GRID_Y 480 //known working values for these are 10, 100000, 0.004, 2 #define MASS 10 //plus 10 zeros #define CENTERMASS 100000 //plus 10 zeros #define TIMESTEP 0.004 //seconds per step #define IMAGESTEP 2 //make bmp every n steps #define rad(d) (3.1415926 * d) / 180 //radians from degress #define DEBUG 0 //0 off, 1 on struct PartType { int x; //x coord for body int y; //y coord unsigned long int m; double vi;//i-hat of velocity vector double vj;//j-hat " " }; struct ForceType { double i; double j; }; struct PartType Particles[NUMPARTS]; struct ForceType ForceSum; struct ForceType MoveForce[NUMPARTS]; //where forces are stored until move() int ImageGrid[GRID_X][GRID_Y]; int Steps,ERR; int ImgSteps=0; Code: `void InitParticles( void ) { int i,j; srand(time()*getpid()); Particles[0].x = (int)(GRID_X / 2); Particles[0].y = (int)(GRID_Y / 2); Particles[0].m = CENTERMASS; Particles[0].vi = 0; Particles[0].vj = 0; for (i=1;i\n"); exit(1); } frames = atoi(argv[1]); printf("Will go for %d iterations, writing %d frames.\n",(int)(frames*IMAGESTEP),frames); InitParticles(); if (DEBUG) { printf("\n [+] Time: %d [+]\n",Steps); for (i=0;i
 OK the source above can by found fully at: crepinc.com/files/nbody.c The movie it put out is at: crepinc.com/files/galax2.wmv (sorry for WMV) Notice how it's better than the first movie (crepinc.com/files/galax.wmv), but how some particles still stop. I can find no reason for this. All the particles are evaluated exactly the same way, why are some treated differently? FYI to compile the code: gcc -O3 -lm -onbody nbody.c Thanks -Jack Carrozzo http://www.crepinc.com/
 Recognitions: Gold Member Science Advisor Does the code work properly when TIMESTEP = 1?
 It DID work, though everything went so fast it was a little hard to tell. That's a good point, I'll check it out tomorrow when I can get to my other comp. Thanks -Jack Carrozzo http://www.crepinc.com/
 Recognitions: Gold Member In fx = cos(rad(theta)) * F; theta is already in radians, so multiplying it by pi/180 is an error
 Aha! I'm an idiot , quite sorry for that. I suppose that would account for the X-favouring issue. I'm not sure why it would have worked before and not not, but tomorrow I'll fix this and see what I get. Thanks a lot, -Jack Carrozzo http://www.crepinc.com/
 OK, I fixed the Radians problem. It seems to be working now. The only issue seems to be the singularity issue. I would imagine that this is an inherant issue in every code of this sort, am I right? The pretty-much ok video is at: http://creponc.com/files/galax3.wmv A funny mess up: http://crepinc.com/files/blowup.wmv Thanks. -Jack Carrozzo http://www.crepinc.com/
 Check out the NEMO site at the University of Maryland: http://bima.astro.umd.edu/nemo/ In addition to posting information about their own work, they also list several sites/resources/software packages relating to N-Body work. Regards, DuncanM http://www.rocketscientists.ca/