//
// probleme des trois corps
//
// reference: Gander-Hrebicek
//
function Ydot = planet(t,Y,m0,m1,m2)
//
//fonction du probleme plan a trois corps (12 composantes, 6 positions, 6 vitesses)
//
x0 = Y(1:2);x1 = Y(5:6);x2 = Y(9:10);
//
d0 = (x2-x1) / norm(x2-x1)^3;
d1 = (x0-x2) / norm(x0-x2)^3;
d2 = (x1-x0) / norm(x1-x0)^3;
Ydot(1:2) = Y(3:4);
Ydot(5:6) = Y(7:8);
Ydot(9:10) = Y(11:12);
Ydot(3:4) =   m1 * d2 - m2 * d1;
Ydot(7:8) =   m2 * d0 - m0 * d2;
Ydot(11:12) = m0 * d1 - m1 * d0;
endfunction
//
m0 = 5;m1 = 3;m2 = 4;listf=list(planet,m0,m1,m2);
x0 = [1,-1];x1 = [1, 3];x2 = [-2,-1];xp = [ 0,0];
tini=0;tfinal=20;
tpas=0.02;
Y1 = ode([x0 xp x1 xp x2 xp]',tini,tini:tpas:tfinal,listf);
xset('figure',0)
clf;
Xmin=min([Y1(1,:),Y1(5,:),Y1(9,:)]);
Ymin=min([Y1(2,:),Y1(6,:),Y1(10,:)]);
Xmax=max([Y1(1,:),Y1(5,:),Y1(9,:)]);
Ymax=max([Y1(2,:),Y1(6,:),Y1(10,:)]);
plot2d([Y1(1,:); Y1(5,:); Y1(9,:)]',[Y1(2,:); Y1(6,:); Y1(10,:)]',[1,2,3],'111',' ',[Xmin,Ymin,Xmax,Ymax]);
xtitle('trajectoires');
//
xset('figure',1)
n=prod(size(Y1(1,:)));
for i=1:(n-10)
 drawlater
 clf()
 plot2d([Y1(1,i:i+10); Y1(5,i:i+10); Y1(9,i:i+10)]',[Y1(2,i:i+10); Y1(6,i:i+10); Y1(10,i:i+10)]',[1,2,3],'111',' ',[Xmin,Ymin,Xmax,Ymax]);
  drawnow
end

