I wrote the following program to solve the three body problem (Sun, Earth, Jupiter), initializing the system such that the total angular momentum is 0. The result of the simulation is really bad (and the Sun tend to infinity). I looked at the code for more than an hour trying to find the problem. Does anybody see what is wrong with my program?
Thank you very much.
#include <iostream>
using std::cout;
using std::endl;
#include <fstream>
using std::ofstream;
#include <cmath>
using std::sqrt;
int main()
{
double time(0), dt(0.05);
double G(1.9838e-29);
double m1(6e24), m2(1.9e27), ms(1.99e30);
double x1(1), y1(0), x2(5.2), y2(0), xs(0), ys(0);
double T1(1), T2(11.8618);
double vx1(0),vy1(2*M_PI*x1/T1), vx2(0), vy2(2*M_PI*x2/T2);
double Rx( (m1*x1+m2*x2+ms*xs)/(m1+m2+ms) ), Ry( (m1*y1+m2*y2+ms*ys)/(m1+m2+ms) );
x1 = x1 - Rx;
x2 = x2 - Rx;
xs = xs - Rx;
y1 = y1 - Ry;
y2 = y2 - Ry;
ys = ys - Ry;
double Lx1(x1*m1*vx1), Lx2(x2*m2*vx2);
double vxs(0), vys( -(Lx1+Lx2)/(ms*xs) );
double r12(0), r1s(0), r2s(0);
ofstream output("TB_jupiter-earth.dat");
output << time << ' ' << x1 << ' ' << y1 << ' ' << x2 << ' ' << y2 << ' ' << xs << ' ' << ys <<endl;
while(time < 12)
{
r12 = sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) );
r1s = sqrt( (x1-xs)*(x1-xs) + (y1-ys)*(y1-ys) );
r2s = sqrt( (x2-xs)*(x2-xs) + (y2-ys)*(y2-ys) );
vx1 = vx1 - G*ms*dt*(x1-xs)/(r1s*r1s*r1s) - G*m2*(x1-x2)*dt/(r12*r12*r12);
vy1 = vy1 - G*ms*dt*(y1-ys)/(r1s*r1s*r1s) - G*m2*(y1-y2)*dt/(r12*r12*r12);
vx2 = vx2 - G*ms*dt*x2/(r2s*r2s*r2s) - G*m1*(x2-x1)*dt/(r12*r12*r12);
vy2 = vy2 - G*ms*dt*y2/(r2s*r2s*r2s) - G*m1*(y2-y1)*dt/(r12*r12*r12);
vxs = vxs - G*m1*(xs-x1)*dt*(r1s*r1s*r1s) - G*m2*(xs-x2)*dt*(r2s*r2s*r2s);
vys = vys - G*m1*(ys-y1)*dt*(r1s*r1s*r1s) - G*m2*(ys-y2)*dt*(r2s*r2s*r2s);
x1 = x1 + vx1*dt;
y1 = y1 + vy1*dt;
x2 = x2 + vx2*dt;
y2 = y2 + vy2*dt;
xs = xs + vxs*dt;
ys = ys + vys*dt;
time += dt;
output << time << ' ' << x1 << ' ' << y1 << ' ' << x2 << ' ' << y2 << ' ' << xs << ' ' << ys <<endl;
}
output.close();
return 0;
}
EDIT:
The correction of the first answer were good but now I changed the calculations of the angula momentum (they were wrong) and the program doesn't work properly anymore. The code is the following:
#include <iostream>
using std::cout;
using std::endl;
#include <fstream>
using std::ofstream;
#include <cmath>
using std::sqrt;
int main()
{
double time(0), dt(0.005);
double G(1.9838e-29);
double m1(6e24), m2(500*1.9e27), ms(1.99e30);
double x1(1), y1(0), x2(5.2), y2(0), xs(0), ys(0);
double T1(1), T2(11.8618);
double vx1(0),vy1(2*M_PI*x1/T1), vx2(0), vy2(2*M_PI*x2/T2);
double Rx( (m1*x1+m2*x2+ms*xs)/(m1+m2+ms) ), Ry( (m1*y1+m2*y2+ms*ys)/(m1+m2+ms) );
x1 = x1 - Rx;
x2 = x2 - Rx;
xs = xs - Rx;
y1 = y1 - Ry;
y2 = y2 - Ry;
ys = ys - Ry;
double Ly1(x1*m1*vy1), Ly2(x2*m2*vy2);
double vxs(0), vys( -(Ly1+Ly2)/(ms*xs) );
double r12(0), r1s(0), r2s(0);
ofstream output("TB_jupiter-earth.dat");
output << time << ' ' << x1 << ' ' << y1 << ' ' << x2 << ' ' << y2 << ' ' << xs << ' ' << ys <<endl;
cout << x1 << ' ' << y1 << ' '<< vx1 << ' ' << vy1 << endl;
cout << x2 << ' ' << y2 << ' '<< vx2 << ' ' << vy2 << endl;
cout << xs << ' ' << ys << ' '<< vxs << ' ' << vys << endl;
while(time < 4)
{
r12 = sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) );
r1s = sqrt( (x1-xs)*(x1-xs) + (y1-ys)*(y1-ys) );
r2s = sqrt( (x2-xs)*(x2-xs) + (y2-ys)*(y2-ys) );
vx1 = vx1 - G*ms*dt*(x1-xs)/(r1s*r1s*r1s) - G*m2*(x1-x2)*dt/(r12*r12*r12);
vy1 = vy1 - G*ms*dt*(y1-ys)/(r1s*r1s*r1s) - G*m2*(y1-y2)*dt/(r12*r12*r12);
vx2 = vx2 - G*ms*dt*(x2-xs)/(r2s*r2s*r2s) - G*m1*(x2-x1)*dt/(r12*r12*r12);
vy2 = vy2 - G*ms*dt*(y2-ys)/(r2s*r2s*r2s) - G*m1*(y2-y1)*dt/(r12*r12*r12);
vxs = vxs - G*m1*(xs-x1)*dt/(r1s*r1s*r1s) - G*m2*(xs-x2)*dt/(r2s*r2s*r2s);
vys = vys - G*m1*(ys-y1)*dt/(r1s*r1s*r1s) - G*m2*(ys-y2)*dt/(r2s*r2s*r2s);
x1 = x1 + vx1*dt;
y1 = y1 + vy1*dt;
x2 = x2 + vx2*dt;
y2 = y2 + vy2*dt;
xs = xs + vxs*dt;
ys = ys + vys*dt;
time += dt;
output << time << ' ' << x1 << ' ' << y1 << ' ' << x2 << ' ' << y2 << ' ' << xs << ' ' << ys <<endl;
}
output.close();
return 0;
}