restricted 3 body gravity simulation c

Question

I'm currently trying to do an assignment where I have to write a simulation for the restricted 3 body gravitational problem, with two fixed masses and one test mass. I have been given the equations I need to use in order to do so, but either I'm not understanding them correctly or I'm not implementing correctly. I would be very grateful is someone could help to push me in the right direction.

The information I've been given is as follows:

http://www.flickr.com/photos/91029993@N07/8269806430/in/photostream

basically I've been trying to test for a single mass,but my program is giving me a straight line for the movement of the test mass, over small times it looks as though the program is functioning correctly but then it doesn't work as you get higher. (see images, inputs for these were 1 0.4 0.5 0 -1)

http://www.flickr.com/photos/91029993@N07/8268764897/in/photostream

I originally wrote my program quite faithfully to the equations given as follows:

#include<stdlib.h>
#include<stdio.h>
#include <math.h>

int main (int argc, char* argv[])
{
    double dt=0.005, x[20000],y[20000],xv,yv,ax[20000],ay[20000],mneg,mpos,time;
    int n;
    FILE* output=fopen("proj1.out", "w");

    sscanf(argv[1], "%lf", &mneg);
    sscanf(argv[2], "%lf", &mpos);
    sscanf(argv[3], "%lf", &x[0]);
    sscanf(argv[4], "%lf", &y[0]);
    sscanf(argv[5], "%lf", &xv);
    sscanf(argv[6], "%lf", &yv);

    x[1]=x[0]+(xv*dt);
    y[1]=y[0]+(yv*dt);



    for(n=1;n<150;n++)
    {
        ax[n]= (-mneg*(x[n]+1)/(pow((sqrt(pow((x[n]+1),2))),3))) -(mpos*(x[n]-1)/(pow((sqrt(pow((x[n]-1),2))),3)));
        ay[n]= (-mneg*(y[n])/(pow(y[n],3))) -(mpos*(y[n])/(pow(y[n],3)));



        x[n+1]=((2*x[n])-x[n-1] +(dt*dt*ax[n]));
        y[n+1]=((2*y[n])-y[n-1]+(dt*dt*ay[n]));


        fprintf(output, "%lf %lf\n",x[n-1], y[n-1]); 
    }
return(0);
}

I then tried a new method following the advice given here :Simulate the gravitational pull of a star? , so that i now have:

#include<stdlib.h>
#include<stdio.h>
#include <math.h>

int main (int argc, char* argv[])
{
    double dt=0.005, x[20000],y[20000],xv,yv,ax[20000],ay[20000],mneg,time,r,a;
    int n;
    FILE* output=fopen("proj1.out", "w");

    sscanf(argv[1], "%lf", &mneg);
    sscanf(argv[2], "%lf", &x[0]);
    sscanf(argv[3], "%lf", &y[0]);
    sscanf(argv[4], "%lf", &xv);
    sscanf(argv[5], "%lf", &yv);

    x[1]=x[0]+(xv*dt);
    y[1]=y[0]+(yv*dt);



for(n=1;n<150;n++)
    {

        r=sqrt(pow((x[n]+1),2)+pow(y[n],2));
        a=mneg/(r*r);



        ax[n]=a*((x[n]+1)/r);
        ay[n]=a*((y[n])/r);


        x[n+1]=((2*x[n])-x[n-1] +(dt*dt*ax[n]));
        y[n+1]=((2*y[n])-y[n-1]+(dt*dt*ay[n]));


        fprintf(output, "%lf %lf\n",x[n-1], y[n-1]); 
    }
return(0);
}

however, this doesn't give me any better results.

I really don't know where to go from here, I think I'm using the method right but I cant see any issues in the actual programming so I really don't know whats going on, so any advice or pointers in the right direction would be very very much appreciated!

score 1 · Accepted Answer · answered Dec 13 '12 at 13:42

1

You did not translate the equations properly into code, and in particular the expression for the distance between the test mass and the fixed masses. The proper way to compute the acceleration is:

double dxm = x[n] + 1.0;
double dxp = x[n] - 1.0;
double dy = y[n];
double denom_minus_inv = pow(dxm*dxm + dy*dy, -1.5);
double denom_plus_inv = pow(dxp*dxp + dy*dy, -1.5);
ax[n] = -mneg*dxn*denom_minus_inv - mpos*dxp*denom_plus_inv;
ay[n] = -mneg*dy*denom_minus_inv - mpos*dy*denom_plus_inv;

Please, use temporary variables to store intermediate expressions - do not put it all in one single complicated expression. Modern compilers are very good at eliminating redundant temporary expressions when they can optimise the code. The code above uses the fact that multiplication is usually a bit faster than division and also that 1.0/pow(sqrt(x), 3.0) == pow(x, -1.5).

I would recommend that you replace the Verlet integrator with velocity Verlet instead. It stores particle velocities explicitly which allows you to compute the total energy of the system (sum of the kinetic and the potential energy) at each step. It should remain almost the same through the entire run (give or take the rounding and discretisation errors). If it deviates wildly, then you known that you do not compute your forces correctly.

answered Dec 13 '12 at 13:42

Hristo Iliev

72,659
12
135
186

thank you very much, i dont have time to properly look through at the moment but quikly putting in you expressions for acceleration is giving me a nice orbit for a single mass :) ill have a proper look this evening. – user1831711 Dec 13 '12 at 14:35
I think `a=dxm*dxm+dy*dy; 1.0/sqrt(a*a*a);` is an order of magnitude faster than pow(a,-1.5). – Aki Suihkonen Dec 13 '12 at 15:05
@AkiSuihkonen, most compilers are smart enough to transform `pow(a, -1.5)` into `1/(a*sqrt(a))`, but you are right - in general one should not rely on compiler smartness. The idea here was to keep the code as short as possible. – Hristo Iliev Dec 13 '12 at 15:34
Don't know about most, but gcc 4.6.1. seems to convert pow(a,2) and pow(a,3) to faster code, but not pow(a, +-1.5); – Aki Suihkonen Dec 13 '12 at 15:38
@AkiSuihkonen, GCC 4.6.1, `-ffast-math -msse42` compiles `pow(a,-1.5)` into `sqrtsd %xmm0,%xmm1; mulsd %xmm0,%xmm1; movsd .LC1(%rip),%xmm0; divsd %xmm1,%xmm0` where `.LC1` is `1.0`. `icc` compiles to the equivalent of `1.0/sqrt(a*a*a)`. – Hristo Iliev Dec 13 '12 at 15:41
It now seems to be working fine for one and two masses, ill try out the velocity verlet tommorow. i dont suppose anyone knows of a simulator i could download to test my results against? – user1831711 Dec 13 '12 at 16:28
@user1831711, it is a bit of an overkill, but the molecular dynamics simulation package [LAMMPS](http://lammps.sandia.gov/) includes long-range Coulombic interactions. – Hristo Iliev Dec 13 '12 at 16:48

restricted 3 body gravity simulation c

1 Answers1