// Copyright jr Licois et Fédération Denis Poisson
//  Modèle avec vitesse de type impulsion
// et deux types de digeres B1 et B2
// constante K  atténuation de la vitesse
// constante c c0 c1 a b  voir modèle vitesse
// constante  kb  disparition des enzymes
// constante kt  absorption des nutriments.
// constantes C1 et C2 pour dégradation A vers B1 et B2
// x(1) = position
// x(2) = derivee position (ie vitesse)
// x(3) = nutriments
// x(4) = enzymes
// x(5) = digeres pas assimilables B1
// x(6) = digeres assimilables B2
function f=digABB(t,x)
  f(1)=x(2);
  f(2)=-K*x(2) + W(t-x(1)/c)*(c0+c1*x(3))/(a+b*x(1))*(1-x(2)/c) ;
  // pour eviter les problemes on prend la valeur absolue
  // avec ajout du terme de boost
  f(3)=-C1*ph(x(1))*x(4)*x(3) - C2*ph(x(1))*x(3)/abs(x(3)+x(5))^(1/3)..
    +1.25*x(3)*indic(x(1)-0.85);
  f(4)=-kb*x(4);
  f(5)=C1*ph(x(1))*x(3)*x(4)- C2*ph(x(1))*x(5)/abs(x(3)+x(5))^(1/3);
  f(6)= C2*ph(x(1))*abs(x(3)+x(5))^(2/3)-kt*x(6);
endfunction

// activité fonction du ph
function kph=ph(x)
  kph = 2.0*abs(x)/3-5*abs(x-3/4)/3+67*abs(x-1)/68+abs(x-18)/68;
endfunction

// echelon de longueur alpha entre 0 et 1
function y=ech(x)
  y=(alpha+abs(x)-abs(x-alpha) )/2/alpha;
endfunction

// derivee de la precedente
function z=indic(x)
  if x  <0 then z=0; 
    elseif   x<alpha then z=1/alpha; else z= 0; end;
endfunction


// somme d'echelons espaces de tau
//function y=W(x)
//  if x < 0 then y=0; 
//  else y=indic(6*x-int(6*x));end;
//endfunction

function yw=W(x)
  //y= indic(2*x)+indic(2*x-1/2)+indic(2*x-1)+indic(2*x-3/2);
  yw = tau*indic(tau*x-fix(tau*x));
endfunction