//lancement affichage 
xset("color",1) 
xset("pixmap",1)//ouverture du buffer video
xtitle("(click droit pour démarrer l''animation, click gauche pour arrêter)") 
titlepage(["résolution numérique d''une EDO  le pendule à ressort : ";" ";"rd²a/dt²+2dr/dt * da/dt=g*sin(a)";" ";"d²r/dt²-(k/m)*(r-r0)=r*(da/dt)²+g*cos(a)";" ";"avec a(0)=?  et da(0)/dt=0   avec r(0)=r0=?  et dr(0)/dt=0"])
show_pixmap() 
rect=[-1,-1,1,0]; 
deff('M=rot(a)','M=[cos(a),sin(a);-sin(a),cos(a)]');
n=40;//spires du ressort
T=10;//durée de la simulation 
[c_i,c_x,c_y,c_w]=xclick();
while (c_i<>2)&(c_i<>5) //tant qu'on a pas fait un click droit 
      xbasc() 
      clear_pixmap()
      xtitle("(click droit pour initialiser la position du pendule)") 
      winnum=winsid() ;
      toolbar(winnum(1),'off');
      plot2d(0,0,-1,frameflag=4,axesflag=5,rect=rect)
      show_pixmap() 
      [c_i,x,y,c_w]=xclick();//coordonnées de la position initiale du pendule
      a=sign(x)*abs(atan(x/y));a0=a;da=0;//calcul de l'angle a(0)   
      l=sqrt(x^2+y^2);r=l;,dr=0;//calcul de  r(0)
      rect=[-1,-max(4*l,3),1,0];//fenêtre
      g=9.81;k=10;dt=0.01;//initialisation 
      //le ressort  
      ressortr=linspace(0,l,n)';
      ressorta=[0;(-1).^[0:n-3]';0]*(l/10);
      for t=0:dt:T
           //mise à jour du calcul de a(t) par méthode incrémentale
           dda=-(g*sin(a)+2*dr*da)/r;
           ddr=r*(da)^2-k*(r-l)+g*cos(a);
           da=da+dt*dda; 
           dr=dr+dt*ddr;
           a=a+dt*da;
           r=r+dt*dr;
           //ressort
           ressortr=linspace(0,r,n)';
           //changement de coordonnées angle--->(x,y)
           x=[x;r*sin(a)];
           y=[y;-r*cos(a)];
           M=-rot(-a);
           N=[ressortr,ressorta]*M;
           ressorty=N(:,1);ressortx=N(:,2);
           //affichage
           xbasc()
           clear_pixmap()
           plot2d(ressortx,ressorty,[1],frameflag=3,axesflag=5,rect=rect)//le ressort  du pendule
           xstring(0.8,-0.5,["t=" string(t)]) // temps écoulé
           xset("color",1)
           xfarc(r*sin(a)-0.05,-r*cos(a)+0.05,0.1,0.1,0,360*64) // la boule du prendule
           show_pixmap()
     end 
      plot2d(x,y,[5],frameflag=4,axesflag=5,rect=rect)
      //xs2gif(0,'pendule_elastique.gif')//pour sauver l'image
      show_pixmap()
      [c_i,x,y,c_w]=xclick();
      rect=[-1,-1,1,0]; 
      xbasc(); 
      xtitle("(click droit pour initialiser la position du pendule)")
      plot2d(0,0,-1,frameflag=4,axesflag=5,rect=rect) 
      show_pixmap() 
      [c_i,x,y,c_w]=xclick();//arrêter ou  continuer 
end
xset("pixmap",0)//fermeture du buffer video 
a=winsid();
xdel(a); //destruction de la fenêtre graphique.