clc
clear

%fs en waarden voor lognormale verdeling van Dstart:
fs=1.1; D1=0.029; D2=0.55; mD=1;
%fs=1.2; D1=0.49; D2=0.55; mD=2;
%fs=1.3; D1=0.8; D2=0.55; mD=3;
%fs=1.4; D1=1.07; D2=0.55; mD=3;
%fs=1.5; D1=1.23; D2=0.55; mD=4;

Ha=13.3;
eb=0.08;
sal=29.9;
eot=0.97*(1-sal/34.92);
es=0.025;
fr_fpa=(1-0.1)*0.45;
T=10.5;
CalgC=0.4; %critical algal concentration [gC/m3]

%Gamma-verdeling IM1:
GIM11=4.27;
GIM12=5.09;

%Normale verdeling Itop:
%m=a*(t-b)+c, s=d*m,
a=2;
b=60;
c=68;
d=0.19;

%Phaeocystis:
%kgp=0.084*(T+3.25);krsp=0.06*1.066^T;kmrt=0.07*1.072^T;eat=0.45; 

%Marine flagellates:
kgp=0.09*(T+1);krsp=0.06*1.066^T;kmrt=0.07*1.072^T;eat=0.25; 

Efc=(krsp+kmrt)/kgp;

%Maximale efficiency factor 1e periode van 3 dagen, zandwinscenario
EfmP1=0.33;

%Grenzen uniforme verdeling voor Calg op t=0
Cmin=0.03;
Cmax=0.3;

%Grenzen uniforme verdeling voor factor Calg_sand/Calg_ref op t0+Dstart
Rmin=0.6;
Rmax=1;

%Startdatum: 
t0=106; %15 april
%t0=121; %30 april
t0Ds=t0+mD; %20 april

%Aantal simulaties:
nsimulations=10000;

for j=1:nsimulations
    
    %IM1 en Itop voor eerste periode van 3 dagen na t0+Dstart
    for t=t0Ds
        Ef(t)=0.1; %dummy
        Efm(t)=0.4; %dummy
        while Ef(t)<Efc
            IM1d(t)=gamrnd(GIM11,GIM12);
            m(t)=a*(t-b)+c;
            s(t)=d*m(t);
            Itopd(t)=normrnd(m(t),s(t));
            et(t)=IM1d(t)*fs*es+eb+eot;
            Efm(t)=1/16*fr_fpa*Itopd(t)*(1-exp(-et(t)*Ha))/(et(t)*Ha);
            if Efm(t)<EfmP1
                Ef(t)=Efm(t);
            else
                Ef(t)=0.1;
            end
        end
    end
    
    %IM1 en Itop voor 2e t/m laatste periode van 3 dagen
    for t=(t0Ds+3):3:(t0Ds+48)
        Ef(t)=0.1; %dummy
        while Ef(t)<Efc
            IM1d(t)=gamrnd(GIM11,GIM12);
            m(t)=a*(t-b)+c;
            s(t)=d*m(t);
            Itopd(t)=normrnd(m(t),s(t));
            et(t)=IM1d(t)*fs*es+eb+eot;
            Ef(t)=1/16*fr_fpa*Itopd(t)*(1-exp(-et(t)*Ha))/(et(t)*Ha);
        end
    end
    
    %Opvullen periodes:
    for t=t0Ds:3:(t0Ds+48)
        IM1(t)=IM1d(t);
        Itop(t)=Itopd(t);
    end

    for t=(t0Ds+1):3:(t0Ds+49)
        IM1(t)=IM1d(t-1);
        Itop(t)=Itopd(t-1);
    end

    for t=(t0Ds+2):3:(t0Ds+50)
        IM1(t)=IM1d(t-2);
        Itop(t)=Itopd(t-2);
    end

    %Dstart
    D=lognrnd(D1,D2);
    
    %Calg op t0:
    Calgs=unifrnd(Cmin,Cmax);
    
    %Itop en IM1 voor periode t0 tot t0+Dstart:
    for t=t0
        Efr(t)=0.1; %dummy
        Efsand(t)=0.4; %dummy
        while Efr(t)<Efc
            IM1s(t)=gamrnd(GIM11,GIM12);
            m(t)=a*(t-b)+c;
            s(t)=d*m(t);
            Itops(t)=normrnd(m(t),s(t));
            
            etref(t)=IM1s(t)*es+eb+eot;
            Efref(t)=1/16*fr_fpa*Itops(t)*...
                    (1-exp(-etref(t)*Ha))/(etref(t)*Ha);
            
            et(t)=IM1s(t)*fs*es+eb+eot;
            Efsand(t)=1/16*fr_fpa*Itops(t)*(1-exp(-et(t)*Ha))/(et(t)*Ha);
         
            if Efsand(t)<Efc
                Efr(t)=Efref(t);
            else
                Efr(t)=0.1;
            end
        end
    end
    
    %Groei van CalgRef in periode van t0 tot t0+Dstart:
    Efs=Efr(t0);
    Calg1(t0Ds-1)=Calgs*(exp((kgp*Efs-krsp-kmrt)))^D;
    
    %Kritische Calg al bereikt voor t0+Dstart?
    Gn=CalgC/Calgs;
    Dn=log(Gn)/(kgp*Efs-krsp-kmrt);

    if D<Dn
        tel1(t0Ds-1)=400; %dummy
    else
        tel1(t0Ds-1)=110-D+Dn;
    end

    %Calg op t0+Dstart voor zandwinscenario
    R=unifrnd(Rmin,Rmax);
    Calg2(t0Ds-1)=R*Calgs;
    tel2(t0Ds-1)=500; %dummy

    %Groei van Calg
    for t=t0Ds:(t0Ds+50)
        et1(t)=eb+eot+eat*Calg1(t-1)+es*IM1(t);
        Ef1(t)=1/16*fr_fpa*Itop(t)*(1-exp(-et1(t)*Ha))/(et1(t)*Ha);
    
        et2(t)=eb+eot+eat*Calg2(t-1)+es*fs*IM1(t);
        Ef2(t)=1/16*fr_fpa*Itop(t)*(1-exp(-et2(t)*Ha))/(et2(t)*Ha);

        Calg1(t)=Calg1(t-1)*exp((kgp*Ef1(t)-krsp-kmrt));
        Calg2(t)=Calg2(t-1)*exp((kgp*Ef2(t)-krsp-kmrt));
        gf1(t)=exp((kgp*Ef1(t)-krsp-kmrt));
        
        %CalgC overschreden?
        if Calg1(t)>CalgC
            tel1(t)=t;
        else
            tel1(t)=400;
        end
    
        if Calg2(t)>CalgC
            tel2(t)=t;
        else
            tel2(t)=500;
        end
    end

    Ex1(j)=min(tel1((t0Ds-1):(t0Ds+50)));
    Ex2(j)=min(tel2((t0Ds-1):(t0Ds+50)));
    Delta(j)=Ex2(j)-Ex1(j);

end
hbins=0:40;
hist(Delta,hbins)
mean(Delta)

%hist(gf1(111:160))

% for t=(t0Ds-1):(t0Ds+50)
%     X(t)=t-(t0Ds-1);
% end
% 
% plot(X((t0Ds-1):(t0Ds+50)),Calg1((t0Ds-1):(t0Ds+50)),'b')
% hold on
% plot(X((t0Ds-1):(t0Ds+50)),Calg2((t0Ds-1):(t0Ds+50)),'r')
