// 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
// constantes du problème
// constante K  atténuation de la vitesse (au pif ! )
K = 0.1;
// constante c c0 c1 a b  voir modèle vitesse
a = 2.72;
b= 1.2346;
c=7;
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 = 1.0/10; //durée impulsion = 1mn (1/10 de l'espacement).
tau = 6.0; // inverse espacement impulsions = 10 mn 
xbasc(0);xbasc(1);
// ******** conditions initiales ****************
Tmax = 1; // injection pendant une heure
alfa = input("Coefficient alpha ?"); // decroissance des injectes
Ti = [0:1/tau:Tmax]; // temps des injections
Pos = zeros(Ti);
Vit = exp(-alfa*Ti);
Nut = (Vit - [Vit(2:$) 0])*2.0; // calcul de V(Ti)-V(Ti+1)
Enz = ones(Ti);
Dig = zeros(Ti);
// matrice des conditions initiales pour l'ensemble des injectes
Minit = [Ti;Pos;Vit;Nut;Enz;Dig];
//
printf("********* debut **********\n");
// initialisation cumuls et couleur trace
nutcum=0;
digcum=0;
couleur = 1;
// Boucle (principale) sur les conditions initiales
for toto=Minit do
// resolution ODE
	init = toto(2:6);
	tzero = toto(1);
// à voir pour être sûr de sortir en 15 mètres !!
	t=tzero:0.01:6+tzero;
// methode rkf obligatoire pour éviter les artefacts !
	y=ode("rkf",init,tzero,t,dig2);   
// tracé dans la fenetre 0 de x et dx/dt
	scf(0);
	xtitle('x et dx/dt  vs t')
	plot2d(t,y(1,:),style=couleur);
	plot2d(t,y(2,:),style=couleur+1);
	legends(['x','dx/dt'],[1,2],"ur")
	// tracé dans la fenetre 1
	scf(1);
	plot2d(t,y(3,:),style=couleur);
	plot2d(t,y(5,:),style=couleur+1);
	legends(['Nutr','Digeres'],[1,2],"ur")
	xtitle('Digestion')
// instructions en fin d'iteration
	couleur = couleur+2;
	printf("Entree = %f\n",toto(4));
// affichage sortie (15 mètres)
	tsv = find(y(1,:)>15); // vecteur des temps 
	// ou on est dehors : seul compte le premier !
	if length(tsv)>0 then //on est sorti
  		ts = tsv(1)*0.01+tzero;
  		nutcum = nutcum+y(3,tsv(1));
  		digcum = digcum +y(5,tsv(1));
  	printf("Sortie a t=%f\n\tNutriments %f\n\tDigeres = %f\n\n",ts,y(3,tsv(1)),y(5,tsv(1)))
	else printf("Pas sorti\n\n")
  	end //endif
//sortie boucle
end
// Resume
printf("\nCumules\n\tNutriments=%f\n\tDigeres=%f\n\n",nutcum,digcum)
// fin.