G=2.95912208286e-4; // constante gravitationnelle
dt=200; // pas de temps
N=1000; // nombre de pas de temps
np=6; // nombre de corps

// Masse, position initiale et vitesse initiale des différents corps
// Soleil
m(1)=1.00000597682;
p(1,1)=0.;
p(1,2)=0.;
p(1,3)=0.;
v(1,1)=0.;
v(1,2)=0.;
v(1,3)=0.;
// Jupiter
m(2)=0.000954786104043;    
p(2,1)=2.2911898857058;
p(2,2)=4.4742096428394;
p(2,3)=-0.0698526085811;
v(2,1)=-0.0068140485228;
v(2,2)=0.0038022184356;
v(2,3)=0.0001366932116;
// Saturne
m(3)=0.000285583733151;
p(3,1)=-8.4077274422046;
p(3,2)=-4.9386528867475;
p(3,3)=0.4206194745212;
v(3,1)=0.0025181864149;
v(3,2)=-0.0048203311629;
v(3,3)=-0.0000164411567;
// Uranus       
m(4)=0.0000437273164546;
p(4,1)=19.9506411858858;
p(4,2)=-2.1189594536568;
p(4,3)=-0.2505524004432;
v(4,1)=-0.0004495229028;
v(4,2)=0.0037304364206;
v(4,3)=0.0000197241717;
// Neptune
m(5)=0.0000517759138449;
p(5,1)=26.3841113101207 ;
p(5,2)=-14.2637141476417;
p(5,3)=-0.3141826415529;
v(5,1)=0.0014665671695;
v(5,2)=0.0027825194375;
v(5,3)=-0.0000910521187;
// Pluton       
m(6)=1/(1.3e8);       
p(6,1)=4.7240342661351;
p(6,2)=-31.8776619468451;
p(6,3)=2.0454833905119;
v(6,1)=0.0031588952992;
v(6,2)=-0.0001620298531;
v(6,3)=-0.0008976519916;

// Création des objets graphiques qui vont afficher la position des corps et leurs trajectoires passées
drawlater();  
for i=1:np
    param3d(p(i,1),p(i,2),p(i,3));
    traj(i) = gce();
    param3d(p(i,1),p(i,2),p(i,3));
    corps(i) = gce();
    corps(i).mark_mode="on";
    corps(i).mark_style = 0;
    corps(i).mark_size = 1;
    corps(i).mark_foreground = i;
end
// Réglage des axes
axes = gca();
axes.data_bounds = [-40,-40,-40;40,40,40];

// Début de la boucle de simulation
for n=1:N
    
    // Calcul des positions et des vitesses au temps suivant par le schéma de Störmer-Verlet
    
    for i=1:np
        p(i,:)=p(i,:)+0.5*dt*v(i,:);
    end
    
    for i=1:np
        for j=1:np
            if (j<>i) then
                v(i,:)=v(i,:)-dt*G*m(j)*(p(i,:)-p(j,:))/(norm(p(i,:)-p(j,:)))^3;    
            end
        end
    end
    
    for i=1:np
        p(i,:)=p(i,:)+0.5*dt*v(i,:);
    end
    
    // Modification du contenu des objets graphiques en vue de l'affichage
    drawlater();
    for i=1:np
        corps(i).data = p(i,:);
        traj0=traj(i).data;
        traj0(n+1,:)=p(i,:);
        traj(i).data=traj0;
    end
    
    // Actualisation de l'affichage (tous les 5 pas seulement)
    if modulo(n,10)==0 then
        drawnow();
    end
    
    // Calcul de l'énergie du système
    e=0;
    for k=1:5
        for j=k+1:6
            e=e-G*m(j)*m(k)/(norm(p(j,:)-p(k,:)));
        end
    end
    for k=1:np
            e=e+0.5*m(k)*norm(v(k,:))^2;
    end
    // Affichage de l'énergie du système dans la console
    disp(e);

end
// Fin de la boucle de simulation
