// Copyright J-R Licois, N Traoré et Fédération Denis Poisson
// Variables du problème :
// L'équation est définie dans le fichier digcylindre3.sci
// constantes du problème

entry=['K','c0','c1','a','b','c','C1','C2',..
      'kt','kb','l','rho','alpha','tau',..
      '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;
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;
B1A0=B10;
B2A0=B20;
B2B10=B20;

CI=[x0;v0;A0;z0;B10;B20;r0;B1A0;B2A0;B2B10];
y=ode("rkf",CI,0,t,digABB);
// ********* 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(3,:),style=1); // A
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(['A','B1','B2','A+B1','A+B1+B2','r'],[1,2,3,5,6,4],"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