// Copyright jr Licois et Fédération Denis Poisson
//*****************************************************
//   MODELE TYPE "PHOSPHATES" suite réunion 13/11/08
//*****************************************************
// L'équation est définie dans le fichier digphos.sci
// constantes du problème
entry=['K attenuation vitesse','c0 cte vitesse','c1 facteur de A vitesse',..
'a attenuation v','b attenuation v','c ','C1 A->B1','C2 A -> B2',..
      'kt abs B2','kb disp enzymes','l longueur cylindre','rho ','alpha impulsion','tau impulsion',..
      'x0','v0','z0','B10','B20','r0'];
init=['15','15','1000','2.72','1.2346','8','1200','1',..
      '.5','log(4/3)','.1','1000','.1','5',..
      '0','1','.001','0','0','.032'];

Pi=4*atan(1);


t=0:0.005:7;
x_message_modeless(['HELP'
'K atténuation vitesse'
'c0 c1 facteur (c0 + c1x) de la vitesse'
'a b atténuation vitesse facteur 1/(a+bx)'
'c  ????'
'kb disparition des enzymes'
'kt absorption nutriments B2'
'C1 A->B1'
'C2 B1->B2'
'alpha tau : impulsions'
'rho l ???'
'x0,v0,z0,B10,B20,r0  conditions initiales'
])
while %t do
donnee= x_mdialog('Parametres',entry,init);
if donnee == [] then break; end
init=donnee;

K=evstr(donnee(1));
c0=evstr(donnee(2));
c1=evstr(donnee(3));
a=evstr(donnee(4));
b=evstr(donnee(5));
c=evstr(donnee(6));
C1=evstr(donnee(7));
C2=evstr(donnee(8));
kt=evstr(donnee(9));
kb=evstr(donnee(10));
l=evstr(donnee(11));
rho=evstr(donnee(12));
alpha=evstr(donnee(13));
tau=evstr(donnee(14));

//conditions initiales
x0=evstr(donnee(15));
v0=evstr(donnee(16));
z0=evstr(donnee(17));
B10=evstr(donnee(18));
B20=evstr(donnee(19));
r0=evstr(donnee(20));

A0=rho*Pi*l*(r0)^2; //
A10 = 0.55*A0;
A20 = 0.4*A0;
B10 = 0.05*A0;
CI=[x0;v0;A10;A20;B10;B20;r0;z0];
y=ode("rkf",CI,0,t,digPhos);
// ********* affichage des constantes ************
ctes_nom = ["K =", msprintf('%4.2f',K)
                     "a =",  msprintf('%4.2f',a) 
                     "b =",  msprintf('%4.2f',b)
                     "c =",  msprintf('%4.2f',c)
                     "c0=",  msprintf('%4.2f',c0)
                     "c1=",  msprintf('%4.2f',c1)
                     "kb=", msprintf('%4.2f',kb)
                     "C1=", msprintf('%4.2f',C1)
                     "C2=", msprintf('%4.2f',C2)
                     "kt=", msprintf('%4.2f',kt)];
// ************************************************     
// tracé dans la fenetre 0
scf(0);
xbasc(0);
plot2d(y(1,:),y(4,:),style=1); // A2
plot2d(y(1,:),y(5,:),style=2); // B1
plot2d(y(1,:),y(6,:),style=3); // B2
//plot2d(y(1,:),y(3,:)+y(5,:),style=5);//A+B1
//plot2d(y(1,:),rho*Pi*l*(y(7,:))^2,style=4);
//plot2d(y(1,:),y(7,:),style=4);  //r
//plot2d(y(1,:),y(3,:)+y(5,:)+y(6,:),style=6);//A+B1+B2
legends(['A2','B1','B2'],[1,2,3],"ur");
xtitle('Digestion vs distance en m');
//[w1,w2,w3,w4]=xgetech();
//xstring(9*w2(3)/10,w2(4)/4,ctes_nom);

scf(1);
xbasc(1);
plot2d(t,y(2,:),style=2);
plot2d(t,y(1,:),style=3);
xtitle('Vitesse et position vs temps (en m/h)');
legends(['Vitesse','Position'],[2,3],"ul");
[w1,w2,w3,w4]=xgetech();
xstring(9*w2(3)/10,w2(4)/4,ctes_nom);

//scf(2);
//xbasc(2);
//plot2d(y(1,:),y(8,:),style=2);//quantite de B1
//plot2d(y(1,:),y(9,:),style=5);//quantite de B2 direct
//plot2d(y(1,:),y(10,:),style=6);//quantite de B2 qui vient de B1
//legends(['B1 crée','B2 par A','B2 par B1'],[2,5,6]);
//xtitle('Cummul B');

end