// Here there are the same comments that the program WCT_optim

function X1= infected_MAC(param)

    
t=9000;
xdata=linspace(1,9000,t)
delta_T=1

//definition of variables

X=zeros(1,t)  // the proportion of people affected
Y=zeros(1,t)  // proportion of mosquitoes affected
R_month=zeros(1,t);
q=zeros(1,t);

//definition of parameters dependant de la temperature ou des rainfalls

//defintion of exogenous parameters
//alpha=exp(-1/1.229);
//a=2500/5000; // the average number of men bitten by one mosquito in one night
b=0.3233;; // the proportion of those anophelines with sporozoites in their salivary glands which are actually infective
//m=10; // the Anopheline density in relation to man
HD=20; // the host delay (length of te interval between infection or sporozoite inoculation, and the onset of infectivity or gametocyte maturation in a host)
WN=20; // host xindow (duration of a host's infectivity to vectors, from the first to the final present of infective gametocytes)
r=1/(HD+WN); // reciprocal of the average duration of the "affected state"
p=0.95; //probability of a mosquito surviving through one whole day;
x=0.01; // x parameter
f_U=36.5; //total number of degree days needed to complete development of the ovaries
g_U=9.9; // minimum temperature to complete development of ovaries
l=5; // difference between inside and outside temperature
nuh=1.26;// length of a part of gonotrophic cycle (find a water body and find a new host)
muh=2000
n=10;
// defintion of Z and C

//Z_0=-(m*a*a*b)*p**n/(r*log(p));
//C=(r/b)*Z_0;

//Initial conditions

X(1,2)=150/50000;
Y(1,2)=0.0100;
R_month(1,1)=T(2,5);



// dynamic equations
for i=2:t-1
    
    Te(i)=T(i,4)+(1-x)*l
    
    //definition of parameters dependant de la temperature ou des rainfalls
    if 0.091678*Te(i)<=1.7982 then a(i)=0.1
    else a(i)= 0.091678*Te(i)-1.7982 ;
    end
        if i<=32 then R_month(1,i)=T(i,5)
    else R_month(1,i)=T(i,5)/30+R_month(1,i-1)-T(i-31,5)/30;
end
    q(1,i)=muh*R_month(1,i);
    m(i)=q(1,i)/50000
    U(i)=nuh + (f_U/(T(i,4)+l-g_U));
    p(i)=(param(1))**(1/U(i))
    Z_0(i)=-(m(i)*a(i)*a(i)*b)*p(i)**n/(r*log(p(i)))
    
    //adjustement of the variables
    
    
w=X(1,i)+a(i)*b*m(i)*Y(1,i)-X(1,i)*(a(i)*b*m(i)*Y(1,i)+r);
if w<=50/50000 then X(1,i+1)=50/50000;
else 
    X(1,i+1)=w;
    end
   if w>=0.03 then X(1,i+1)=0.03
end
w=Y(1,i)+a(i)*X(1,i)-Y(1,i)*(a(i)*X(1,i)-log(p(i)));
if w<=0.0001 then Y(1,i+1)=0.0001; 
else 
    Y(1,i+1)=w;
    end
     if w>=0.1 then
         Y(1,i+1)=0.1
     end



end
X1=50000*(X(1,:)')
endfunction

// different loglikelihood

function[L1]=loglikelihood(param,varargin);
    L1=0;
    In=infected_MAC(param)
    for i=2:1000 do
        L1=L1-log(1/(In(i)*(2-exp(-1)))*exp(-abs(Cases(i,2)-In(i))/(In(i))));
    end
endfunction;

function[L2]=loglikelihood2(param,varargin);
    L2=0;
    In=infected_MAC(param)
    for i=2:9000 do
        L2=L2+(Cases(i,2)-In(i))^2;
    end
    endfunction;

// here we calculate the maximum of likelihood


[cost,p_opt,g_opt]=optim(list(NDcost,loglikelihood,200),'b',[0.2],[0.9],[0.7])