// Tracé de l'ellipse
e=0.6;
for i=1:101
  theta=2*%pi*i/100;
  r=(1-e^2)/(1+e*cos(theta));
  x(i)=r*cos(theta);
  y(i)=r*sin(theta);
end
plot2d(x,y,rect=[-2,-1.5,1,1.5]);

// Résolution par la méthode d'Euler
function kepler_euler(T,N)
  h=T/N;
  e=0.6;
  pold=[1-e,0];
  vold=[0,sqrt((1+e)/(1-e))];
  for i=1:N
    pnew=pold+h*vold;
    vnew=vold-h*pold/(norm(pold))^3;
    pold=pnew;
    vold=vnew;
    // stockage de la position pour l'affichage de la trajectoire
    px(i)=pold(1); 
    py(i)=pold(2);
  end
  plot2d(px,py);
endfunction

// Résolution par la méthode d'Euler symplectique
function kepler_euler_symp(T,N)
  h=T/N;
  e=0.6;
  pold=[1-e,0];
  vold=[0,sqrt((1+e)/(1-e))];
  for i=1:N
    pnew=pold+h*vold;
    vnew=vold-h*pnew/(norm(pnew))^3;
    pold=pnew;
    vold=vnew;
    // stockage de la position pour l'affichage de la trajectoire
    px(i)=pold(1); 
    py(i)=pold(2);
  end
  plot2d(px,py);
endfunction

// Résolution par la méthode de Störmer-Verlet
function kepler_stormer_verlet(T,N)
  h=T/N;
  e=0.6;
  pold=[1-e,0];
  vold=[0,sqrt((1+e)/(1-e))];
  for i=1:N
    pint=pold+0.5*h*vold;
    vnew=vold-h*pint/(norm(pint))^3;
    pnew=pint+0.5*h*vnew;
    pold=pnew;
    vold=vnew;
    // stockage de la position pour l'affichage de la trajectoire
    px(i)=pold(1); 
    py(i)=pold(2);
  end
  plot2d(px,py);
endfunction

// Résolution par la méthode de Störmer-Verlet avec animation
function kepler_anime(T,N)
  
  h=T/N;
  e=0.6;
  pold=[1-e,0];
  vold=[0,sqrt((1+e)/(1-e))];
  
  drawlater();    
  plot(0,0,'o');
  corps1 = gce();
  corps1.children.mark_size = 10;
  corps1.children.mark_background = 3;
  plot(pold(1),pold(2),'o');
  corps2 = gce();
  corps2.children.mark_size = 10;
  corps2.children.mark_background = 2;
  axes = gca();
  axes.data_bounds = [-2,-2;1.5,1.5];
  
  for i=1:N
    pint=pold+0.5*h*vold;
    vnew=vold-h*pint/(norm(pint))^3;
    pnew=pint+0.5*h*vnew;
    pold=pnew;
    vold=vnew;
    corps2.children.data = pold;
    drawnow();
  end
  
endfunction
