// Copyright JR Licois, N Traoré et Fédération Denis Poisson
// 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 pour dégradation A vers B1, et C2 pour A et B1 vers B2
// x(1) = position
// x(2) = derivee position (ie vitesse)
// x(3) = nutriments (A)
// x(4) = enzymes
// x(5) = digeres pas assimilables B1
// x(6) = digeres assimilables B2
// x(7) = rayon du cylindre
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)= - C2*2*sqrt(Pi*l/rho)*x(3)/sqrt(abs(x(3)+x(5))) - C1/rho*ph(x(1))*x(4)/z0*x(3)..
    +log(1.25)/alpha*x(2)*(x(3)+x(5))*indic((x(1)-0.85)/alpha);
  f(4)= - kb*x(4)+log(1.25)/alpha*x(2)*x(4)*indic((x(1)-0.85)/alpha);//saut enzyme
  f(5)= - C2*2*sqrt(Pi*l/rho)*x(5)/sqrt(abs(x(3)+x(5)))..
        + C1/rho*ph(x(1))*x(4)/z0*x(3);
  f(6)= C2*2*sqrt(Pi*l/rho)*sqrt(abs(x(3)+x(5))) - kt*x(6);
  f(7)= - C2/rho + log(1.25)/alpha/2*x(2)*indic((x(1)-0.85)/alpha)*x(7);
  f(8)= C1/rho*ph(x(1))*x(4)/z0*x(3);
  f(9)= C2*2*sqrt(Pi*l/rho)*x(3)/sqrt(abs(x(3)+x(5)));
  f(10)=C2*2*sqrt(Pi*l/rho)*x(5)/sqrt(abs(x(3)+x(5)));
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