// dynamique 
function qdot=Cart(t,q)
// q(1) q(2) le vecteur position
// q(3) q(4) le vecteur vitesse
global  alpha ;
// la constante de gravité redimensionée

    qdot(1)=q(3);
    qdot(2)=q(4);
    mod=sqrt(q(1)**2+q(2)**2)**3;
    if (mod>0) then
        qdot(3)=-alpha**2 *q(1)/  mod;
        qdot(4)=-alpha**2 *q(2)/  mod;
     else
         disp("Pb")
     end    
    
endfunction
// programme principal

//unités
// temps : h
// masse : kg
// distance : Mm

global alpha ;

pf=42.165;// orbite geo
p0=11.625;// orbite basse


// les constantes
mu=398600.47*10^(-9)*(3600^2);
alpha=sqrt(mu/(pf^3));


//integration directe en coordonnées q1 q2 q1dot q2dot

// cas geostationnaire orbite haute
// condition initiale
q1_0=0
q2_0=1
q1dot_0=-alpha
q2dot_0=0
// resolution de l'equation differentielle
t=0:1:24;
z=ode("rk",[q1_0;q2_0;q1dot_0;q2dot_0],0,t,Cart);
// visualisation trajectoire
plot(z(1,1:size(t,2)),z(2,1:size(t,2)),'b.')
// on recalcule l ode sur 24h
tE=0:0.1:24;
z=ode("rk",[q1_0;q2_0;q1dot_0;q2dot_0],0,tE,Cart);
// calcul d'energie
// $$$ for i=1:1:size(tE,2)
// $$$  E(1,i)=0.5*(z(3,i)**2+z(4,i)**2)
// $$$  E(2,i)=-alpha**2/ sqrt(z(1,i)**2+z(2,i)**2) 
// $$$  E(3,i)=E(1,i)+E(2,i)
// $$$ end


// cas geostationnaire orbite basse
// condition initiale
q1_0=p0/pf/(1+0.75)
q2_0=0
q1dot_0=0
q2dot_0=alpha*sqrt(pf*(1+0.75)/p0)
// resolution de l'equation differentielle
//t2=0:0.1:24
t2=0:0.1:5;
z2=ode("rk",[q1_0;q2_0;q1dot_0;q2dot_0],0,t2,Cart);
// visualisation trajectoire
plot(z2(1,1:size(t2,2)),z2(2,1:size(t2,2)),'g.')
// on recalcule l ode sur 24h
z=ode("rk",[q1_0;q2_0;q1dot_0;q2dot_0],0,tE,Cart);
// calcul d'energie
// $$$ for i=1:1:size(tE,2)
// $$$  E2(1,i)=0.5*(z(3,i)**2+z(4,i)**2)
// $$$  E2(2,i)=-alpha**2/ sqrt(z(1,i)**2+z(2,i)**2) 
// $$$  E2(3,i)=E2(1,i)+E2(2,i)
// $$$ end

// cas ellipse
// condition initiale
q1_0=-p0/pf/(1-0.75)
q2_0=0
q1dot_0=0
q2dot_0=-(1-0.75)*alpha*sqrt(pf/p0)
// resolution de l'equation differentielle
t1=0:0.2:12;
//t1=0:0.1:24
z1=ode("rk",[q1_0;q2_0;q1dot_0;q2dot_0],0,t1,Cart);
// visualisation trajectoire
plot(z1(1,1:size(t1,2)),z1(2,1:size(t1,2)),"r+")
// on recalcule l ode sur 24h
z=ode("rk",[q1_0;q2_0;q1dot_0;q2dot_0],0,tE,Cart);
// Calcul d'energie
// $$$ for i=1:1:size(tE,2)
// $$$  E1(1,i)=0.5*(z(3,i)**2+z(4,i)**2);
// $$$  E1(2,i)=-alpha**2/ sqrt(z(1,i)**2+z(2,i)**2); 
// $$$  E1(3,i)=E1(1,i)+E1(2,i);
// $$$ end

// cas ellipse
// condition initiale
q1_0=p0/pf/(1+0.75)
q2_0=0
q1dot_0=0
q2dot_0=1.25*alpha*sqrt(pf*(1+0.75)/p0)
// resolution de l'equation differentielle
t3=0:0.2:12;
//t1=0:0.1:24
z3=ode("rk",[q1_0;q2_0;q1dot_0;q2dot_0],0,t3,Cart);
// visualisation trajectoire
plot(z3(1,1:size(t3,2)),z3(2,1:size(t3,2)),"y+")
// on recalcule l ode sur 24h
z=ode("rk",[q1_0;q2_0;q1dot_0;q2dot_0],0,tE,Cart);


// *** decommenter si on veut regarder les energies ***

// Calcul d'energie
// $$$ for i=1:1:size(tE,2)
// $$$  E3(1,i)=0.5*(z(3,i)**2+z(4,i)**2);
// $$$  E3(2,i)=-alpha**2/ sqrt(z(1,i)**2+z(2,i)**2);
// $$$  E3(3,i)=E3(1,i)+E3(2,i);
// $$$ end
// scf
// $$$ plot(tE(1,1:size(tE,2)), E(1,1:size(tE,2)),'*b',tE(1,1:size(tE,2)), E(2,1:size(tE,2)),'.b',tE(1,1:size(tE,2)), E(3,1:size(tE ,2)),'b')
// $$$ plot(tE(1,1:size(tE,2)),E1(1,1:size(tE,2)),'*r',tE(1,1:size(tE,2)),E1(2,1:size(tE,2)),'.r',tE(1,1:size(tE,2)),E1(3,1:size(tE,2)),'r')
// $$$ plot(tE(1,1:size(tE,2)),E2(1,1:size(tE,2)),'*g',tE(1,1:size(tE,2)),E2(2,1:size(tE,2)),'.g',tE(1,1:size(tE,2)),E2(3,1:size(tE,2)),'g')
// $$$ plot(tE(1,1:size(tE,2)),E3(1,1:size(tE,2)),'*y',tE(1,1:size(tE,2)),E3(2,1:size(tE,2)),'.y',tE(1,1:size(tE,2)),E3(3,1:size(tE,2)),'y')