write(%io(2),'etude du pendule pesant (amorti ou non, entretenu ou non)');
val=evstr(x_dialog('type de pendule (0:non amorti/1:amorti/2:amorti entretenu)','1'));
X=evstr(x_dialog('valeurs initiales de vitesse angulaire (2)?','[1.99,3]'));
rho=zeros(1,2);
if (val==1)
  rho=evstr(x_dialog('valeurs de rho (frottement)?','[0.2,0.2]'));
elseif (val==2)
  rho=evstr(x_dialog('valeurs de mu (acceleration a l''origine)?','[0.1,0.01]'));
end
Tmax=evstr(x_dialog('temps de simulation?','50'));
//
//
deff('Y=pendule(t,X,rho,val)',[
'Y(1)=X(2)';
'if (val==2)';
   'Y(2)=-sin(X(1))-(X(1)^2-rho)*X(2)';
'else';
   'Y(2)=-sin(X(1))-rho*X(2)';
'end']);
//
y0=[0,X(1)]';t0=0;t=0:0.1:Tmax;
//solution donnee 1
[yr1]=ode(y0,t0,t,list(pendule,rho(1)));
//
y1=[0,X(2)]';
//solution donnee 2
//
[yr2]=ode(y1,t0,t,list(pendule,rho(2)));
//
// trace des resultats
//
xset('window',0)
xbasc(0)
xsetech([0.,0.,0.5,1])
plot2d(t,yr1(1,:))
xtitle('donnee initiale 1','t','theta')
xsetech([0.5,0.,0.5,1])
plot2d(t,yr2(1,:))
xtitle('donnee initiale 2','t','theta')
xset('window',1)
xbasc(1)
Xmin1=min(yr1(1,:));Xmax1=max(yr1(1,:));
Ymin1=min(yr1(2,:));Ymax1=max(yr1(2,:));
Xmin2=min(yr2(1,:));Xmax2=max(yr2(1,:));
Ymin2=min(yr2(2,:));Ymax2=max(yr2(2,:));
xsetech([0.,0.,0.5,1])
plot2d(yr1(1,:),yr1(2,:),1,'011',' ',[Xmin1,Ymin1,Xmax1,Ymax1])
xtitle('donnee initiale 1','theta','v')
xsetech([0.5,0.,0.5,1])
plot2d(yr2(1,:),yr2(2,:),1,'011',' ',[Xmin2,Ymin2,Xmax2,Ymax2])
xtitle('donnee initiale 2','theta','v')
pause
xbasc(1)
nc=prod(size(yr1(1,:)));
for i=1:nc-10
xbasc()
xsetech([0.,0.,0.5,1])
plot2d(yr1(1,i:i+10),yr1(2,i:i+10),1,'011',' ',[Xmin1,Ymin1,Xmax1,Ymax1])
xsetech([0.5,0.,0.5,1])
plot2d(yr2(1,i:i+10),yr2(2,i:i+10),1,'011',' ',[Xmin2,Ymin2,Xmax2,Ymax2])
end
//
xset('window',2)
xbasc(2)
xsetech([0.,0.,0.5,1])
param3d(cos(yr1(1,:)),sin(yr1(1,:)),yr1(2,:))
xtitle('donnee initiale 1','theta mod 2pi','v')
xsetech([0.5,0.,0.5,1])
param3d(cos(yr2(1,:)),sin(yr2(1,:)),yr2(2,:))
xtitle('donnee initiale 2','theta mod 2pi','v')
// xset('window',3)
pause
for i=1:int(Tmax/0.1)
xbasc()
xsetech([0.,0.,0.5,1])
plot2d([0;sin(yr1(1,i))],[0;-cos(yr1(1,i))],1,'011',' ',[-1,-1,1,1])
xsetech([0.5,0.,0.5,1])
plot2d([0;sin(yr2(1,i))],[0;-cos(yr2(1,i))],1,'011',' ',[-1,-1,1,1])
end


