Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

NBody evaluation, is my model correct?

  1. May 25, 2005 #1
    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. jcsd
  3. May 25, 2005 #2

    tony873004

    User Avatar
    Science Advisor
    Gold Member

    Can't tell without seeing the code, especially the lines that contain the time step variable.
     
  4. May 25, 2005 #3
    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.
     
  5. May 25, 2005 #4
    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? :bugeye:

    -Jack Carrozzo
    http://www.crepinc.com/
     
  6. May 26, 2005 #5
    Hello,
    Try to work with different DT. Best of Luck. Zaheer A
     
  7. May 26, 2005 #6
    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 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 (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);
                    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;
      }
    }
     
    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
  8. May 26, 2005 #7
    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
  9. May 26, 2005 #8

    tony873004

    User Avatar
    Science Advisor
    Gold Member

    Does the code work properly when TIMESTEP = 1?
     
  10. May 26, 2005 #9
    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/
     
  11. May 26, 2005 #10
    In fx = cos(rad(theta)) * F; theta is already in radians, so multiplying it by pi/180 is an error
     
  12. May 26, 2005 #11
    Aha! I'm an idiot :redface: , 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/
     
  13. May 27, 2005 #12
  14. Jun 10, 2005 #13
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: NBody evaluation, is my model correct?
  1. Is this correct? (Replies: 3)

Loading...