CODEWITHSUNDEEP
cintegrationsimulation
user3722989
user3722989

Reputation: 23

Velocity Verlet implementation difficulties in C (The Sun drifts away)

I'm trying to simulate the earth-sun system via velocity verlet, but somehow the sun will not orbit around the origin (where the reduced mass is located), but drifts away. I've spent a good time looking through my algorithm, but can't find the flaw.

Would anybody have a clue of what's going wrong here?

Here's a plot of the simulation: https://i.sstatic.net/lNoSN.png

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

double xearth,yearth,vxearth,vyearth;
double xsun,ysun,vxsun,vysun;
double dt=0.5;
double fxearth,fyearth;
double fxsun,fysun;
double r;
double G;
double ms, ma;
double rx,ry;
double t;


main(){
    FILE * pFile;
    int n;
    xearth= -2.569651552438753*pow(10,-2);  /* in AU */
    yearth= -1.008909556982513;
    xsun= 2.563054664344734*pow(10,-4);
    ysun= 6.897319465467234*pow(10,-3);

    vxearth= 1.690809814669721*pow(10,-2); /* in AU per day */
    vyearth= -4.950293720310762*pow(10,-4);
    vxsun= -5.788119594348977*pow(10,-6);
    vysun= 3.335986886320253*pow(10,-6);

    G=1.488*pow(10,-34); /* G in AU, t in day */
    ms=1.9884*pow(10,30); /* kg */
    ma=5.9722*pow(10,24); /* kg */
    t=0;


    pFile = fopen ("/file.txt", "w");

    rx=xearth-xsun;
    ry=yearth-ysun;

    r=sqrt((rx*rx+ry*ry));

    fxearth= -G*ms*ma*(rx)/pow(r,3);
    fyearth= -G*ms*ma*(ry)/pow(r,3);

    fxsun= -G*ms*ma*(-rx)/pow(r,3);
    fysun= -G*ms*ma*(-ry)/pow(r,3);

    vxearth=vxearth+.5*dt/ma*fxearth;
    vyearth=vyearth+.5*dt/ma*fyearth;

    vxsun=vxsun+.5*dt/ms*fxsun;
    vysun=vysun+.5*dt/ms*fysun;


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

        xearth=xearth+dt*vxearth;
        yearth=yearth+dt*vyearth;

        xsun=xsun+dt*vxsun;
        ysun=ysun+dt*vysun;

        rx=xearth-xsun;
        ry=yearth-ysun;

        r=sqrt((rx*rx+ry*ry));


        fxearth= -G*ms*ma*(rx)/pow(r,3);
        fyearth= -G*ms*ma*(ry)/pow(r,3);

        fxsun= -G*ms*ma*(-rx)/pow(r,3);
        fysun= -G*ms*ma*(-ry)/pow(r,3);


        vxearth=vxearth+dt/ma*fxearth;
        vyearth=vyearth+dt/ma*fyearth;

        vxsun=vxsun+dt/ms*fxsun;
        vysun=vysun+dt/ms*fysun;


        t=t+dt;


        fprintf(pFile,"%f\t %f\t %f\t %f\t %f\n",xearth,yearth,xsun,ysun,t);
    }

    fclose (pFile);

    return 0;

}

Upvotes: 2

Views: 469

Answers (1)

Ilmari Karonen
Ilmari Karonen

Reputation: 50368

It happens because your initial conditions give the system a non-zero net momentum. You can fix that by calculating the initial mean velocity of the system and subtracting it from all the object velocities:

double vxavg = (vxsun*ms + vxearth*ma) / (ms + ma);
double vyavg = (vysun*ms + vyearth*ma) / (ms + ma);

vxsun -= vxavg;
vysun -= vyavg;
vxearth -= vxavg;
vyearth -= vyavg;

Upvotes: 1

Related Questions