1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
| double f(double x,double y){
double G=6.674*pow(10,-11);
double Ms=1.989*pow(10,30);
double Aso=13.34*pow(10,19);
double F;
double r=sqrt(x*x+y*y);
F=-((Aso*x)/(r*r*r));
return F;
}
double g(double x,double y){
double G=6.674*pow(10,-11);
double Ms=1.989*pow(10,30);
double Aso=13.34*pow(10,19);
double r=sqrt(x*x+y*y);
double g;
g=-((Aso*y)/(r*r*r));
return g;
}
double RK4_PLANETE(double x, double y, double u, double v, double n){
double h=0.1;
double kx1,kx2,kx3,kx4,A,B;
double ky1,ky2,ky3,ky4;
double ku1,ku2,ku3,ku4;
double kv1,kv2,kv3,kv4;
ofstream traj("traj.txt");
for (float t=0 ; t<1 ; t+=h){
kx1=u*h;
ky1=v*h;
ku1=f(x,y)*h;
kv1=g(x,y)*h;
kx2=(u+ku1/2.0)*h;
ky2=(v+kv1/2.0)*h;
ku2=(f(x+kx1/2.0,y+ky1/2.0))*h;
kv2=(g(x+kx1/2.0,y+ky1/2.0))*h;
kx3=(u+ku2/2.0)*h;
ky3=(v+kv2/2.0)*h;
ku3=(f(x+kx2/2.0,y+ky2/2))*h;
kv3=(g(x+kx2/2.0,y+ky2/2))*h;
kx4=(u+ku3)*h;
ky4=(v+kv3)*h;
ku4=(f(x+kx3,y+ky3))*h;
ku4=(g(x+kx3,y+ky3))*h;
x=x+((1/6.0)*(kx1+2*kx2+2*kx3+kx4));
y=y+((1/6.0)*(ky1+ky2*2+ky3*2+ky4));
u=u+((1/6.0)*(ku1+ku2*2+ku3*2+ku4));
v=v+((1/6.0)*(kv1+kv2*2+kv3*2+kv4));
traj<<x<<" "<<y<<endl;
};
} |
Partager