NBody evaluation, is my model correct?

[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/

 

[tony873004]

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.

 

[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.

 

[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/

 

[Zaheer Ahmed]

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

 

[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)

#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;

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;
		}
	}
}

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);
				fx = cos(rad(theta)) * F;
        fy = sin(rad(theta)) * F;
				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;
  }
}

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)

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/

 

[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/

 

[tony873004]

Does the code work properly when TIMESTEP = 1?

 

[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/

 

[chronon]

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

 

[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/

 

[crepincdotcom]

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/

 

[DuncanM]

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/