NBody evaluation, is my model correct?

1. May 25, 2005

crepincdotcom

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/

2. May 25, 2005

tony873004

Can't tell without seeing the code, especially the lines that contain the time step variable.

3. May 25, 2005

chronon

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.

4. May 25, 2005

crepincdotcom

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/

5. May 26, 2005

Zaheer Ahmed

Hello,
Try to work with different DT. Best of Luck. Zaheer A

6. May 26, 2005

crepincdotcom

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 (Text):

#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 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 (Text):

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<NUMPARTS;i++) //starts at 1, not 0
{
Particles[i].x = ((double)rand()/((double)pow(2,31))) * GRID_X;
Particles[i].y = ((double)rand()/((double)pow(2,31))) * GRID_Y;
Particles[i].m = MASS;
Particles[i].vi = 0;
Particles[i].vj = 0;
}

for (i=0;i<GRID_X;i++)
{
for (j=0;j<GRID_Y;j++)
{
ImageGrid[i][j]=0;
}
}
}

Code (Text):

void ComputeForces( void )
{
int i,j,fxsign,fysign; //signs: 0 +, 1 -
double r,F,theta,dx,dy,fx,fy;
double G=0.667; //has 10 zeros taken off

for (i=0;i<NUMPARTS;i++)
{
ForceSum.i=0;
ForceSum.j=0;
for (j=0;j<NUMPARTS;j++)
{
if (j!=i)
{
fxsign=0;
fysign=0;
if ((Particles[j].x - Particles[i].x) < 0) fxsign++;
if ((Particles[j].y - Particles[i].y) < 0) fysign++;
dx = abs(Particles[i].x - Particles[j].x);
dy = abs(Particles[i].y - Particles[j].y);
theta = atan(dy/dx); //theta in radians
r = sqrt((dx * dx) + (dy * dy)); //distance
F = ((G * (double)(Particles[i].m) * (double)(Particles[j].m)) / (r * r));
if (DEBUG) printf("F= %e, G= %e m1= %d, m2= %d, r= %f\n",F,G,Particles[i].m,Particles[j].m,r);
if (DEBUG) if (i==2) printf("compareing to %d, r - %f, force x %f, force y %f, theta %f\n",j,r,fx,fy,theta);
if (fxsign==0)
ForceSum.i += fx;
else
ForceSum.i -= fx;

if (fysign==0)
ForceSum.j += fy;
else
ForceSum.j -= fy;
}
}
MoveForce[i].i = ForceSum.i;
MoveForce[i].j = ForceSum.j;
}
}

Code (Text):

void MoveParticles( void )
{
int i;
double A; //acceleration

for (i=0;i<NUMPARTS;i++)
{
//F=MA, A=F/M
A = (MoveForce[i].i / (Particles[i].m));
Particles[i].vi += (A * TIMESTEP);
A = (MoveForce[i].j / (Particles[i].m));
A *= 25; //don't know why, but need this
Particles[i].vj += (A * TIMESTEP);

Particles[i].x += (int)(Particles[i].vi);
Particles[i].y += (int)(Particles[i].vj);

if (DEBUG) if (i==2) printf("part 2, y %f, f %f, %f.\n",Particles[i].vj,MoveForce[i].i,MoveForce[i].j);
}
Steps++;
}

("DataOut and bitmap functions left out)
Code (Text):

int main( int argc, char *argv[] )
{
int i, frames;

printf("Simple NBODY solver by Jack Carrozzo.\n");

if (argc!=2)
{
printf("Usage: ./nbody <frames>\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<NUMPARTS;i++)
{
printf("----------------------\n");
printf("Particle Num     : %d\n",i);
printf("Particle Position: (%d,%d)\n",Particles[i].x,Particles[i].y);
printf("Particle Velocity: (%f,%f)\n",Particles[i].vi,Particles[i].vj);
printf("----------------------\n");
}
}
DataOut();

for (i=0;i<(frames*IMAGESTEP);i++)
{
ComputeForces();
MoveParticles();
DataOut();
}

if (DEBUG)
{
printf("\n   [+] Time: %d [+]\n",Steps);
for (i=0;i<NUMPARTS;i++)
{
printf("----------------------\n");
printf("Particle Num     : %d\n",i);
printf("Particle Position: (%d,%d)\n",Particles[i].x,Particles[i].y);
printf("Particle Velocity: (%f,%f)\n",Particles[i].vi,Particles[i].vj);
printf("----------------------\n");
}
}
printf("Done.\n");

return 0;
}

Thanks a lot,

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

Last edited: May 26, 2005
7. May 26, 2005

crepincdotcom

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/

Last edited: May 26, 2005
8. May 26, 2005

tony873004

Does the code work properly when TIMESTEP = 1?

9. May 26, 2005

crepincdotcom

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/

10. May 26, 2005

chronon

In fx = cos(rad(theta)) * F; theta is already in radians, so multiplying it by pi/180 is an error

11. May 26, 2005

crepincdotcom

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/

12. May 27, 2005

crepincdotcom

13. Jun 10, 2005