Ciao. Ho un problema questo codice. L'ho scritto per simulare il problema a tre corpi (completo) con il metodo di Euler-Cromer. Purtroppo il programma non da il risultato sperato... E da ore che ci lavoro ma senza risultato...

Qualcuno saprebbe aiutarmi? Grazie!

codice:
#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;
}