function [out]=Fooo(t,Y,Yprim)
// Function generated by Maple to Scilab interface
out=[
Yprim(1)-Y(7);
Yprim(2)-Y(8);
Yprim(3)-Y(9);
Yprim(4)-Y(10);
Yprim(5)-Y(11);
Yprim(6)-Y(12);
M*Yprim(7)+m1*Yprim(7)+m1*Yprim(10)*sin(Y(3))+2.0*m1*Y(10)*cos(Y(3))*Y..
(9)-m1*Y(4)*sin(Y(3))*Y(9)**2+m1*Y(4)*cos(Y(3))*Yprim(9)+m2*Yprim(7)+m..
2*Yprim(10)*sin(Y(3))+2.0*m2*Y(10)*cos(Y(3))*Y(9)-m2*Y(4)*sin(Y(3))*Y(..
9)**2+m2*Y(4)*cos(Y(3))*Yprim(9)+m2*Yprim(12)*sin(Y(5))+2.0*m2*Y(12)*c..
os(Y(5))*Y(11)-m2*Y(6)*sin(Y(5))*Y(11)**2+m2*Y(6)*cos(Y(5))*Yprim(11)-..
2.0*Y(13)*Y(1)+0.11D1*Y(13)*sin(0.33D1*Y(1))+k*Y(7);
M*Yprim(8)+m1*Yprim(8)-m1*Yprim(10)*cos(Y(3))+2*m1*Y(10)*sin(Y(3))*Y(9..
)+m1*Y(4)*cos(Y(3))*Y(9)**2+m1*Y(4)*sin(Y(3))*Yprim(9)+m2*Yprim(8)-m2*..
Yprim(10)*cos(Y(3))+2*m2*Y(10)*sin(Y(3))*Y(9)+m2*Y(4)*cos(Y(3))*Y(9)**..
2+m2*Y(4)*sin(Y(3))*Yprim(9)-m2*Yprim(12)*cos(Y(5))+2*m2*Y(12)*sin(Y(5..
))*Y(11)+m2*Y(6)*cos(Y(5))*Y(11)**2+m2*Y(6)*sin(Y(5))*Yprim(11)+M*g+m1..
*g+m2*g+Y(13)+k*Y(8);
Y(4)*(m1*g*sin(Y(3))+m2*g*sin(Y(3))+m1*cos(Y(3))*Yprim(7)+m1*sin(Y(3))..
*Yprim(8)+m2*cos(Y(3))*Yprim(7)+m2*sin(Y(3))*Yprim(8)+m2*cos(Y(3))*Ypr..
im(12)*sin(Y(5))-m2*sin(Y(3))*Yprim(12)*cos(Y(5))+2*m2*cos(Y(3))*Y(12)..
*cos(Y(5))*Y(11)-m2*cos(Y(3))*Y(6)*sin(Y(5))*Y(11)**2+m2*cos(Y(3))*Y(6..
)*cos(Y(5))*Yprim(11)+2*m2*sin(Y(3))*Y(12)*sin(Y(5))*Y(11)+m2*sin(Y(3)..
)*Y(6)*cos(Y(5))*Y(11)**2+m2*sin(Y(3))*Y(6)*sin(Y(5))*Yprim(11)+m1*Y(4..
)*Yprim(9)+2*m1*Y(9)*Y(10)+m2*Y(4)*Yprim(9)+2*m2*Y(9)*Y(10));
-m2*g*cos(Y(3))+m1*Yprim(10)-m2*cos(Y(3))*Yprim(8)-m1*g*cos(Y(3))+m1*s..
in(Y(3))*Yprim(7)-m1*cos(Y(3))*Yprim(8)+m2*sin(Y(3))*Yprim(12)*sin(Y(5..
))+m2*cos(Y(3))*Yprim(12)*cos(Y(5))+m2*sin(Y(3))*Yprim(7)+k1*Y(4)-k1*l..
1bar+2*m2*sin(Y(3))*Y(12)*cos(Y(5))*Y(11)-m2*sin(Y(3))*Y(6)*sin(Y(5))*..
Y(11)**2+m2*sin(Y(3))*Y(6)*cos(Y(5))*Yprim(11)-2*m2*cos(Y(3))*Y(12)*si..
n(Y(5))*Y(11)-m2*cos(Y(3))*Y(6)*cos(Y(5))*Y(11)**2-m2*cos(Y(3))*Y(6)*s..
in(Y(5))*Yprim(11)-m1*Y(4)*Y(9)**2-m2*Y(4)*Y(9)**2+m2*Yprim(10);
m2*Y(6)*(-cos(Y(5))*Y(4)*sin(Y(3))*Y(9)**2+cos(Y(5))*Y(4)*cos(Y(3))*Yp..
rim(9)+2*cos(Y(5))*Y(10)*cos(Y(3))*Y(9)-sin(Y(5))*Yprim(10)*cos(Y(3))+..
Y(6)*Yprim(11)+g*sin(Y(5))+cos(Y(5))*Yprim(10)*sin(Y(3))+2*Y(11)*Y(12)..
+cos(Y(5))*Yprim(7)+sin(Y(5))*Yprim(8)+2*sin(Y(5))*Y(10)*sin(Y(3))*Y(9..
)+sin(Y(5))*Y(4)*cos(Y(3))*Y(9)**2+sin(Y(5))*Y(4)*sin(Y(3))*Yprim(9));
-m2*g*cos(Y(5))+m2*Yprim(12)-m2*cos(Y(5))*Yprim(8)+m2*sin(Y(5))*Yprim(..
10)*sin(Y(3))+m2*cos(Y(5))*Yprim(10)*cos(Y(3))-m2*Y(6)*Y(11)**2+2*m2*s..
in(Y(5))*Y(10)*cos(Y(3))*Y(9)+m2*sin(Y(5))*Yprim(7)-m2*sin(Y(5))*Y(4)*..
sin(Y(3))*Y(9)**2+m2*sin(Y(5))*Y(4)*cos(Y(3))*Yprim(9)-2*m2*cos(Y(5))*..
Y(10)*sin(Y(3))*Y(9)-m2*cos(Y(5))*Y(4)*cos(Y(3))*Y(9)**2-m2*cos(Y(5))*..
Y(4)*sin(Y(3))*Yprim(9)+k2*Y(6)-k2*l2bar;
Y(7)*(-2*Y(1)+0.11D1*sin(0.33D1*Y(1)))+Y(8);
]

endfunction


function [r,ir]=FFFF(t,Y,Yprim)
r=Fooo(t,Y,Yprim);
ir=0;
endfunction

// pendule vertical de longueur l
// pendu en (x,y) avec un angle theta
function [xv,yv]=pend_vert(x,y,theta,l)
xv=[0;0;(-1).^(1:nspire)'*ep;0;0];
yv=[0;-ef;-(1:nspire)'*(l-2*ef)/(nspire+1)-ef;-l+ef;-l]
// translation
xv=xv+x; yv=yv+y;
// rotation
rot=[cos(theta) -sin(theta); sin(theta) cos(theta)];
N=nspire+4;
res=[x*ones(1,N);y*ones(1,N)]+rot*[xv'-x;yv'-y];
xv=res(1,:)'; yv=res(2,:)';
endfunction

function  [B,P1,P2]=contruit_pendule()
// construit la figure du pendule avec ses entités graphiques
// retourne les handles sur la boule et les pendules
  clf();
  // définition de la figure
  f=gcf();a=gca();
  toolbar(f.figure_id,"off");
  a.isoview="on";
  f.pixmap="on";
  a.box="off";
  
  drawlater();
  xmin=-1.25; xmax=1.25; ymin=-3; ymax=1;
  a.data_bounds=[xmin ymin;xmax ymax]
  // le cadre
//  xrect(xmin,ymax,xmax-xmin,ymax-ymin)
  // la courbe
  vx=[xmin:0.01:xmax]'; vy=vx.*vx+cos(3.3*vx)/3;
  xpoly(vx,vy,"lines")
  c=gce();c.foreground=color("red"); c.thickness=2;
  // la boule
  r=0.05
  x=0; y=0; theta1=0; l1=l10;
  xfarc(x-r,y+r,2*r,2*r,0,360*64)
  B=gce();
  // la tige du pendule 1
  x1=x+l1*sin(theta1); y1=y-l1*cos(theta1);
  r=0.05 // le rayon de la boule du pendule 1
  [xv,yv]=pend_vert(x,y,theta1,l1)
  xpoly(xv,yv,"lines")
  p=gce();p.thickness=2;
  // la boule du pendule 1
  xfarc(x1-r,y1+r,2*r,2*r,0,360*64)
  b=gce()
  P1=glue([p,b]) //retourne le handle sur la tige et la boule 1
  // la tige du pendule 2
  theta2=0; l2=l20;
  x2=x1+l1*sin(theta1); y2=y1-l1*cos(theta1);
  r=0.05 // le rayon de la boule du pendule 2
  [xv,yv]=pend_vert(x1,y1,theta2,l2)
  xpoly(xv,yv,"lines")
  p=gce();p.thickness=2;
  // la boule du pendule 2
  xfarc(x2-r,y2+r,2*r,2*r,0,360*64)
  b=gce()
  P2=glue([p,b]) //retourne le handle sur la tige et la boule 1
endfunction

function dessine_pendule(B,P1,P2,position)
// dessine une position du pendule
  drawlater();
  x=position(1); y=position(2); theta1=position(3); l1=position(4);
  theta2=position(5); l2=position(6);
// boule
  r=B.data(3)/2;
  B.data=[x-r,y+r,2*r,2*r,0,360*64];
// premier pendule
  b = P1.children(1);r=b.data(3)/2
  x1=x+l1*sin(theta1); y1=y-l1*cos(theta1);
  p = P1.children(2);
  [xv,yv]=pend_vert(x,y,theta1,l1)
  p.data=[xv yv];
  b = P1.children(1); b.data=[x1-r,y1+r,2*r,2*r,0,360*64];
// deuxieme pendule
  b = P2.children(1);r=b.data(3)/2
  x2=x1+l2*sin(theta2); y2=y1-l2*cos(theta2);
  p = P2.children(2);
  [xv,yv]=pend_vert(x1,y1,theta2,l2)
  p.data=[xv yv];
  b = P2.children(1); b.data=[x2-r,y2+r,2*r,2*r,0,360*64];
  drawnow();
  show_pixmap()
endfunction

function demo_sliding_pendulum()
//
// Sliding pendulum
// Claude Gomez
// Copyright INRIA

  demo_help  demo_sliding_pendulum
// données
g=10; M=1; k=0.35; l1bar=1; m1=1; k1=40; l2bar=1; m2=1; k2=40;
// conditions initiales données
x0=1; y0=1+cos(3.3)/3; theta10=0; l10=l1bar; theta20=0; l20=l2bar; lambda0=0;
xprim0=0; yprim0=0; theta1prim0=0; l1prim0=0; theta2prim0=0; l2prim0=0; lambdaprim0=0;
xprimprim0=0; yprimprim0=-g; theta1primprim0=0; l1primprim0=0; theta2primprim0=0; l2primprim0=0;

Y0=[x0;y0;theta10;l10;;theta20;l20;xprim0;yprim0;theta1prim0;l1prim0;theta2prim0;l2prim0;lambda0];
Yprim0=[xprim0;yprim0;theta1prim0;l1prim0;theta2prim0;l2prim0;
        xprimprim0;yprimprim0;theta1primprim0;l1primprim0;theta2primprim0;l2primprim0;
	      lambdaprim0];

t0=0; t=t0:0.05:20;
atol=[0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.001;0.001;0.001;0.001;0.001;0.001;0.01];
rtol=atol;
Y=dassl([Y0,Yprim0],t0,t,rtol,atol,FFFF);

nspire=10;
ep=0.1;
ef=0.1;  
// dessine le pendule dans sa position initiale
[B,P1,P2]=contruit_pendule();
dessine_pendule(B,P1,P2,Y0(1:8));

// animation du pendule
realtimeinit(0.1);realtime(0) 
for i=1:size(Y,2)
  realtime(i);
  dessine_pendule(B,P1,P2,Y(2:7,i))
end
endfunction

