// Copyright jr Licois et Fédération Denis Poisson
//*****************************************************
//   MODELE TYPE "PHOSPHATES" suite réunion 13/11/08
//*****************************************************
// 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) : sans objet
// x(4) = nutiments A2  degradables : la proportion A1/A2 est une
// fonction µ(x) du ph.
// Le nutriment A3 directement assimilable devient la condition initiale
// pour B1. 
// x(5) = digeres pas assimilables B1
// x(6) = digeres assimilables B2
// x(7) = rayon du cylindre
// x(8) = enzymes
function f=digPhos(t,x)
  f(1)= x(2);
  f(2)= - K*x(2) + W(t-x(1)/c)*(c0+c1*x(4)*mu(x(1)))/(a+b*x(1))*(1-x(2)/c) ;
  f(3) = 0; // on calcule en fait A1 comme A1 = µ*A2.
  f(4) = (C1*ph(x(1))*x(8)/z0/rho+ dmu(x(1))*x(2)*x(4))/(1-mu(x(1)));
  f(5)= - C2*2*sqrt(Pi*l/rho)*x(5)/sqrt(abs(x(4)*(1+mu(x(1))+x(5))))..
        + C1*ph(x(1))*x(8)*x(4)*mu(x(1))/z0/rho;
  f(6)= C2*2*sqrt(Pi*l/rho)*x(5)*sqrt(abs(x(4)*(1+mu(x(1)))+x(5))) - kt*x(6);
  f(7)= - C2/rho; 
  f(8)=  - kb*x(8); 
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

// equilibre entre A1 et A2
function y = mu(x)
  y = 0.5;
endfunction

// dérivée de la précédente
function y = dmu(x)
  y = 0.0;
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