function xdot = suspens(t,x,u)
  // etat [z,v];
  xdot=[ x(2);(1/m)*( (k*(profil(t)+h - x(1)) -m*g  + G(u(t),profildot(t)-x(2))))]
endfunction 

function g = G(u,v) 
  g = mu*v + (2*(a0*u+b0)/%pi)*atan(v*(a1*u+b1)/(a0*u+b0))
  g = g* ( c + (d/%pi)*2*atan(v/beta))
endfunction

function h = wsin(t) 
  // profil de route sinusoidal 
  h = Hp*sin(2*%pi*V*t/T)
endfunction 

function h = wsindot(t) 
  // derivee du profil de route sinusoidal 
  h = Hp*2*%pi*V*cos(2*%pi*V*t/T)/T
endfunction 


function h = obstm(t) 
  // profil de route avec montée d'obstacle 
  h = H*(1/2 + (1/%pi)*atan(10*(t-5)))
endfunction 

function h = obstmdot(t) 
  // dérivée du profil de route avec montée d'obstacle 
  h = H*( (10/%pi)/( 1+ (10*(t-5))^2))
endfunction 

function h = obstd(t) 
  // profil de route avec descente d'obstacle 
  h = - H*(1/2 + (1/%pi)*atan(10*(t-5)))
endfunction 

function h = obstddot(t) 
  // dérivée du profil de route avec descente d'obstacle 
  h = - H*( (10/%pi)/( 1+ (10*(t-5))^2))
endfunction 

function x = u(t)
  // coeficient de raidissment de l'amortisseur
  // entre -1 et 1.5 
  x= -1
endfunction 

// recherche d'un linéarisé tangent 
// pour un profil constant 

function xdot = suspensu(t,x,u1)
  // etat [z,v];
  xdot=[ x(2);(1/m)*( (k*(h - x(1)) -m*g  + G(u1,-x(2))))]
endfunction 

function xe = equilu(u1) 
  xe=[h-m*g/k-G(u1,0);0];
endfunction

m= 215; // kg 
g= 9.81; // N/kg 
mu= 2000/3;// Ns/m 
a0= 600;
b0= 700;
a1 = 7000;
b1 = 26000;
alpha = 155/170;
beta= 0.05; 
k= 50000; // N/m: raideur du ressort 
h= 0.8 ; // m : hauteur du ressort a l'équilibre 
V= 10; //m/s  vitesse horizontale 
// constantes 
c  = (1+ alpha)/2 ;
d  = (1- alpha)/2;

// Ction initiales 
z0= h - m*g/k //m  hauteur initiale de la caisse 
v0= 0.0 //m/s  vitesse initiale verticale 

// Paramètre de route 
H = 0.05 // 5 cm  hauteur de l'obstacle 
Hp = 0.02 // amplitude de la sinusoide (1/2 hauteur des pavés)
T = 0.10 // s période de la route sinusoidale 

// Pour un profil constant =0 
// et ue fixé on calcule un point d'équilibre 
// et le linéarisé tangent 

ue=0;
xe=equilu(ue);
// vérification 
if norm(suspensu(0,xe,ue)) > 1.e-8 then pause;end 


//XXXXX a corriger 
//[F,G]=tangent('suspensu',xe,ue)


choice = 1; 



if choice == 1 
  // exemple profil sinusoidal 
  profil = wsin 
  profildot = wsindot 
  
  pas = 1.e-4 ; 
  x0=[z0;v0];
  t=0:pas:3;
  x= ode(x0,0,t,suspens);
  xset('window',0)
  xbasc();
  plot(t,x(1,:)) ;
  xtitle('Hauteur de caisse','temps(s)','Hauteur(m)')
  // dessin de la route 
  xset('window',1)
  xbasc();
  xx= feval(t,profil);
  plot2d(t,xx);
  xtitle('Profil de route','temps(s)','Hauteur(m)')

elseif choice == 2 
   // exemple de montée d'obstacle 
   profil = obstm
   profildot = obstmdot
   
   pas = 1.e-4 ; 
   x0=[z0;v0];
   t=0:pas:10;
   x= ode(x0,0,t,suspens);
   xset('window',0);xbasc();
   plot(t,x(1,:)) ;
   xtitle('Hauteur de caisse','temps(s)','Hauteur(m)')
   // dessin de la route 
   xset('window',1);xbasc();
   xx= feval(t,profil);
   plot2d(t,xx);
   xtitle('Profil de route','temps(s)','Hauteur(m)')

elseif choice == 3 

   //exemple de descente d'obstacle 
   profil = obstd
   profildot = obstddot
   
   pas = 1.e-4 ; 
   x0=[z0;v0];
   t=0:pas:10;
   x= ode(x0,0,t,suspens);
   xset('window',0);xbasc();
   plot(t,x(1,:)) ;
   xtitle('Hauteur de caisse','temps(s)','Hauteur(m)')
   // dessin de la route 
   xx= feval(t,profil);
   xset('window',1);xbasc();
   plot2d(t,xx);
   xtitle('Profil de route','temps(s)','Hauteur(m)')

end