// suspension en boucle fermée 

function xdot = suspens(t,x)
  // etat [z,v];
  xdot=[ x(2);(1/m)*( (k*(profil(t)+h - x(1)) -m*g  + U(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 obstacle 
  h = H*(1/2 + (1/%pi)*atan(10*(t-5)))
endfunction 

function h = obstmdot(t) 
  // profil de route avec obstacle 
  h = H*( (10/%pi)/( 1+ (10*(t-5))^2))
endfunction 


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

function h = obstddot(t) 
  // profil de route avec obstacle 
  h = - H*( (10/%pi)/( 1+ (10*(t-5))^2))
endfunction 


function rep = U(t,vz)
  // coeficient de raidissment de l'amortisseur
  // entre -1 et 1.5 
  // commande en boucle fermée 
  umin=-1 
  umax=1.5 
  vr = profildot(t)-vz;
  a=(-r*vz -G(umin,vr))*vr 
  b=(r*vz + G(umax,vr))*vr
  if a >= 0 & b >=0 
    rep = -r*vz; 
  elseif a < 0 
     rep = G(umin,vr); 
  else 
     rep = G(umax,vr); 
  end
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 

// paramètre de controle (gain)

r = evstr(x_dialog('valeur du gain','0.0'))

choice = 1; 

// exemple profil sinusoidal 

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

elseif choice == 2 
   // exemple profil sinusoidal 

   profil = obstm
   profildot = obstmdot
   
   pas = 1.e-4 ; 
   x0=[z0;v0];
   t=0:pas:10;
   x= ode(x0,0,t,suspens);
   xset('window',0)
   plot(t,x(1,:)) ;
   // dessin de la route 
   xset('window',1)
   xx= feval(t,profil);
   plot2d(t,xx);

elseif choice == 3 
   // exemple profil sinusoidal 

   profil = obstd
   profildot = obstddot
   
   pas = 1.e-4 ; 
   x0=[z0;v0];
   t=0:pas:10;
   x= ode(x0,0,t,suspens);
   xset('window',0)
   plot(t,x(1,:)) ;
   // dessin de la route 
   xset('window',1)
   xx= feval(t,profil);
   plot2d(t,xx);

end