write(%io(2),'production de testosterone: recherche d''un cycle periodique ([Mur], p. 166)')
A=evstr(x_dialog('valeurs de b1 b2 b3 g1 et g2','[1,1,1,1,1]'));
tauca=evstr(x_dialog('valeur de tau (temps de retard)','3'));
dt=evstr(x_dialog('pas de temps (<tau calcul)','0.05'));n=int(tauca/dt);dt=tauca/n;
N=evstr(x_dialog('nombre de pas de temps de calcul','5000'));
b1=A(1);b2=A(2);b3=A(3);g1=A(4);g2=A(5);
// 
// fonction de feeback
//
m=8;deff('y=f(x)','y=1./(1+x.^m)');
//
deff('y=fdot(x)','y=-m.*x.^(m-1)./(1+x.^m).^2');
//
deff('y=fonct(x)','y=b1*b2*b3/(g1*g2)*x-f(x)');
//
T0=fsolve(0,fonct);// point fixe 
a=b1+b2+b3;
b=b1*b2+b2*b3+b3*b1;
c=b1*b2*b3;
q=fdot(T0);
d=-fdot(T0)*g1*g2;
//
deff('y=h(x,tau)','if (x<sqrt(b))&(0<x)&(x*tau<%pi),y=1/tan(tau*x)-(a*x^2-c)/(x*(b-x^2)),else y=1,end;');
deff('y=nu(tau)','y=fsolve(0.0001,list(h,tau))');
deff('y=g(tau)','y=(b-nu(tau).^2).*nu(tau)-d*sin(nu(tau).*tau)');
//
tau=[0.005,0.01,0.05:0.05:5,10,15];
xbasc()
xsetech([0,0,0.5,0.5]);
  fplot2d(tau,nu);
  xtitle('evolution de nu(tau)','tau');
xsetech([0.5,0,0.5,0.5]);
     fplot2d(tau,g)
  xtitle('evolution de g(tau)','tau');
//
// resolution et trace du systeme differentiel avec retard
//
x=0.1*ones(N+n+1,1);
y=0.1*ones(N+n+1,1);
z=0.1*ones(N+n+1,1);
for i=n+2:N+n
  x(i+1)=x(i)+dt*(q*z(i)-b1*x(i));
  y(i+1)=y(i)+dt*(g1*x(i)-b2*y(i));
  z(i+1)=z(i)+dt*(g2*y(i-n-1)-b3*z(i));
end
t=0:dt:N*dt;
xsetech([0,0.5,1,0.5])
plot2d([t',t',t'],[x(n+1:N+n+1),y(n+1:N+n+1),z(n+1:N+n+1)],[1,2,3],'121','hormone LHRH@hormone LH@testosterone');
xtitle('concentration d''hormones','t');
//