// Copyright jr Licois et Fédération Denis Poisson
// Variables du problème :
// f(1) = position
// f(2) = derivee position (ie vitesse)
// f(3) = nutriments
// f(4) = enzymes
// f(5) = digérés
// L'équation est définie dans le fichier dig2.sci
// vitesse avec 3 impulsions d'amplitude décroissantes.
// constantes du problème
// constante K  atténuation de la vitesse (au pif ! )
K = 1.5;
// constante c c0 c1 a b  voir modèle vitesse
a = 2.72;
b= 1.2346;
c= 7.0;
c0 = 1.0; // au pif
c1 = 2.0; // au pif
// constante  kb  disparition des enzymes
kb=log(4/3); 
// constante kt  absorption des nutriments.
kt=1/3;
// constante C  correspond à 36\pi/\rho^2 .... (au pif!)
C = 1.5;
alpha = 0.1; //durée impulsion = 1mn
tau = 6.0; // espacement impulsions = 10 mn 
// tracé de la fonction W 
scf(2);
xbasc(2);
fplot2d([0:0.001:1],W)
//xtitle('ph vs x')
// resolution ODE
init=[0;1;1;1;0];
t=0:0.005:5;
y=ode("rkf",init,0,t,dig2);
// tracé dans la fenetre 0 de x et dx/dt
scf(0);
xbasc(0);
xtitle('x et dx/dt  vs t')
plot2d(t,y(1,:),style=1);
plot2d(t,y(2,:),style=2);
legends(['x','dx/dt'],[1,2],"ur")
// tracé dans la fenetre 1
scf(1);
xbasc(1);
plot2d(t,y(3,:),style=1);
plot2d(t,y(5,:),style=2);
legends(['Nutr','Digeres'],[1,2],"ur")
xtitle('Digestion')
// retour à la fenetre 0 pour tracer A en fonction de x
scf(0);
// pour extraire tous les 10 pas de temps les valeurs de x et a
//tt=10*[1:35];
//xx=y(1,tt);
//yy=y(2,tt);
// diagramme en bâtons
//plot2d3(xx,yy,style=2)