clear;

//chikungunya

gamm=7/7; //incubation moustique
mu=7/30; //par semaine
delta=7/4; //incubation homme
alpha=7/7; //par semaine
bet=7/4; //par semaine
phi=2*%pi/12; 
omega=2*%pi/52;
P=785000;

r=0.06;
epsi=(1-r)/(1+r); 
pmax=1.2*P/bet; //par semaine
pmin=r*pmax;
p0=0.5*(pmin+pmax);

T=2*%pi/omega;

//order 2
//disp("R0 order 2 = "+string(a*b*(1-alpha*mu*e^2/2/(omega^2+(alpha+mu)^2))/alpha/mu));

//dichotomie
deff("[dqdt]=floquet(t,q)",["e=q(1)";
"i=q(2)";...
"E=q(3)";...
"I=q(4)";...
"dqdt(1)=-(gamm+mu)*e+(bet*p0/R0/P)*(1+epsi*cos(omega*t-phi))*I";...
"dqdt(2)=gamm*e-mu*i";...
"dqdt(3)=bet*i-delta*E";...
"dqdt(4)=delta*E-alpha*I"]);
R0min=0.5;
R0max=100;
while R0max-R0min>1e-7,
  R0=(R0min+R0max)/2;
  c1=ode([1;0;0;0],0,[0 T],floquet);
  colonne1=c1(:,2);
  c2=ode([0;1;0;0],0,[0 T],floquet);
  colonne2=c2(:,2);
  c3=ode([0;0;1;0],0,[0 T],floquet);
  colonne3=c3(:,2);
  c4=ode([0;0;0;1],0,[0 T],floquet);
  colonne4=c4(:,2);
  sp=spec([colonne1 colonne2 colonne3 colonne4]);
  rho=max(real(sp));
  if rho>1 then
    R0min=R0;
  else
    R0max=R0;
  end;
end;
disp("R0 dichotomy = "+string(R0));
  
//eigenvalue  
if 0==0 then
  N=100;
  matrice=zeros(N,N);
  dt=T/N;
  deff("[y]=Ghat(x)","y=(bet^2*p0*gamm*delta/P)*(...
  exp(-alpha*x)/(1-exp(-alpha*T))/(gamm+mu-alpha)/(mu-alpha)/(delta-alpha)...
  +exp(-mu*x)/(1-exp(-mu*T))/(delta-mu)/(alpha-mu)/gamm...
  +exp(-delta*x)/(1-exp(-delta*T))/(alpha-delta)/(gamm+mu-delta)/(mu-delta)...
  +exp(-(mu+gamm)*x)/(1-exp(-(mu+gamm)*T))/(-gamm)/(delta-mu-gamm)/(alpha-mu-gamm))");
  for i=1:N,
    ti=(i-1)*dt;
    f(i)=1+epsi*cos(omega*ti-phi);
  end;
  for i=1:N,
    ti=(i-1)*dt;
    for j=1:i-1,
      tj=(j-1)*dt;
      matrice(i,j)=f(i)*Ghat(ti-tj);
    end;
    for j=i:N,
      tj=(j-1)*dt;
      matrice(i,j)=f(i)*Ghat(ti-tj+T);
    end;
  end; 
  matrice=matrice*dt;
  sp2=spec(matrice);
  R0=max(real(sp2));
  disp("R0 eigenvalue = "+string(R0));
end;

//Fourier
if 0==0 then
deff("[y]=G(n)","y=(bet^2*p0*gamm*delta/P)/(mu+n*%i*omega)/(alpha+n*%i*omega)/(delta+n*%i*omega)/(gamm+mu+n*%i*omega)");
kappa=12;
Lambda=zeros(kappa+1,1);
c=zeros(kappa+2,kappa+1);
Lambda(1)=G(0);
c(1,1)=1;
for k=1:kappa,
  Lambda(k+1)=real(G(1)*c(2,k));
  for n=1:k,
    c(n+1,k+1)=(0.5*G(n-1)*c(n,k)+0.5*G(n+1)*c(n+2,k)-c(n+1,k:-1:2)*Lambda(2:k))/(G(0)-G(n));
  end;
end;
R0=(epsi.^[0:kappa])*Lambda(1:kappa+1);
disp("R0 Fourier = "+string(R0));
end;

//Dye
if 0==0 then
N=4;
M=zeros(2*N+1,2*N+1);
deff("[y]=f(n)","y=1*(n==0)+(epsi/2)*(n==1)+(epsi/2)*(n==-1)");
deff("[y]=G(n)","y=(bet^2*p0*gamm*delta/P)/(mu+n*%i*omega)/(alpha+n*%i*omega)/(delta+n*%i*omega)/(gamm+mu+n*%i*omega)");

fourier=zeros(2*N+1,2*N+1);
for n=-N:N,
  fourier(N+n+1,N+n+1)=f(0)*G(n);
  if n>-N then
    fourier(N+n+1,N+n)=f(1)*G(n-1);
  end;
  if n<N then
    fourier(N+n+1,N+n+2)=f(-1)*G(n+1);
  end;
end;
R0=max(abs(spec(fourier)));
disp("R0 Williams = "+string(R0));
end;