clear
clc

nsim=100000;

for n=1:nsim
    
    %Growth factor
    Gr=[10.8 
        3.36 
        1.86 
        1.37 
        1.17 
        1.09];

    %Biomass of recruits
    BmRcr=0.12;
    
    %Aantal jaren na start zandwinning:
    nyears=13;

    %Mismatch
    for y=1:8
        
        %H(y)=normrnd(139,10); %18 mei, grote std
        H(y)=normrnd(139,8); %18 mei
        %H(y)=normrnd(146,8); %25 mei
        
        %Aref(y)=normrnd(108,15); %17 april, std=15
        Aref(y)=normrnd(106,15); %15 april, std=15
        %Aref(y)=normrnd(100,14); %9 april, std=14
        %Aref(y)=normrnd(96,13); %5 april, std=13
    
        Dfa=[1.41 1.91 2.37 2.78 3.15];
        Dfb=[1.43 1.64 1.78 1.97 2.08];
        
        Da=[Dfa(2) Dfa(5) Dfa(5) Dfa(3) Dfa(2) Dfa(2) Dfa(1) Dfa(1)];
        Db=[Dfb(2) Dfb(5) Dfb(5) Dfb(3) Dfb(2) Dfb(2) Dfb(1) Dfb(1)];
        
        D(y)=gamrnd(Da(y),Db(y));
        ZR(y)=-H(y)+Aref(y);
        ZS(y)=-H(y)+Aref(y)+D(y);
    end

    %Helemaal geen mismatch, door naar volgende simulatie
    ZSmax=max(ZS);
    if ZSmax<0
        EffectSandMining(n,1:nyears)=0;
    else
        for y=1:8;
            
            %Mismatch in jaar y voor ref?
            if ZR(y)<0
                MMR(y)=0;
            else
                if ZR(y)<21
                    MMR(y)=ZR(y);
                else
                    MMR(y)=21; 
                end
            end
        
            %Mismatch in jaar y voor sand?
            if ZS(y)<0
                MMS(y)=0;
            else
                if ZS(y)<21
                    MMS(y)=ZS(y);
                else
                    MMS(y)=21; 
                end
            end
        
            fR(y)=Calg_MM(MMR(y));
            fS(y)=Calg_MM(MMS(y));
            
            if fS(y)>fR(y)
                fS(y)=fR(y);
            else
                fS(y)=fS(y);
            end
            
            gR(y)=4.8/(1+exp(-0.0025*(fR(y)-1850)));
            gS(y)=4.8/(1+exp(-0.0025*(fS(y)-1850)));
            
            GLR(y)=1-( (154+MMR(y)*gR(y)+(21-MMR(y))*4.8)/(154+4.8*21) )^2.8;
            GLS(y)=1-( (154+MMS(y)*gS(y)+(21-MMS(y))*4.8)/(154+4.8*21) )^2.8;
        
            CF(y)=corrGL;
        
            GLCR(y)=GLR(y)*CF(y);
            GLCS(y)=GLS(y)*CF(y);
        
            ER(y)=1-GLCR(y);
            ES(y)=1-GLCS(y);
        end
    
        %Vector EMR en EMS vullen voor jaren 1 t/m 13
        for y=1:8
            EMR(y)=ER(y);
            EMS(y)=ES(y);
        end

        for y=9:nyears
            EMR(y)=1;
            EMS(y)=1;
        end
    
        %biomass at start
        for yp=1:6
            BmP(1,yp)=recrrnd*BmRcr;
            BmP(2,yp)=BmP(1,yp)*Gr(1);
    
            MrtWP(yp)=mrtwrnd_35;
    
            MrtSP(yp)=mrtsrnd;
    
        end

        for ap=3:7
            BmP(ap,1)=1; %dummy value
            for yp=2:6
                BmP(ap,yp)=BmP(ap-1,yp-1)*Gr(ap-1)*MrtSP(yp)*MrtWP(yp);
            end
        end
        
        %Voor jaren 1 t/m 13
        for Y=1:nyears
            BmRecr(Y)=recrrnd*BmRcr;
            MrtS(Y)=mrtsrnd;
            MrtW(Y)=mrtwrnd_35;
        end
    
        for Y=1:nyears
            BmR(1,Y)=BmRecr(Y); %Biomass reference scenario
            BmS(1,Y)=BmRecr(Y); %Biomass sand mining scenario
        
            BmR(2,Y)=BmR(1,Y)*Gr(1)*EMR(Y);
            BmS(2,Y)=BmS(1,Y)*Gr(1)*EMS(Y);
        end
    
        for A=3:7
            BmR(A,1)=BmP(A,6);
            BmS(A,1)=BmP(A,6);
        
            for Y=2:nyears
                BmR(A,Y)=BmR(A-1,Y-1)*Gr(A-1)*MrtS(Y)*MrtW(Y);
                BmS(A,Y)=BmS(A-1,Y-1)*Gr(A-1)*MrtS(Y)*MrtW(Y);
            end
        end
    
        for Y=1:nyears
            BmREat(Y)=BmR(3,Y)+BmR(4,Y)+BmR(5,Y)+BmR(6,Y)+BmR(7,Y);
            BmSEat(Y)=BmS(3,Y)+BmS(4,Y)+BmS(5,Y)+BmS(6,Y)+BmS(7,Y);
        end
    
        %vertaling naar eidereenden
    
        %lineair:
%         for Y=1:nyears
%             EffectSandMining(n,Y)=1-BmSEat(Y)/BmREat(Y);
%         end
        
        %exponential, SY=0.35:
        for Y=1:nyears
            EiderR(Y)=2.3*(1-exp(-0.102*BmREat(Y)));
            EiderS(Y)=2.3*(1-exp(-0.102*BmSEat(Y)));
            EffectSandMining(n,Y)=(1-EiderS(Y)/EiderR(Y));
            EffectEiders(n,Y)=(EiderR(Y)-EiderS(Y))/1.35*2500;
            EffectEiders((nsim+1),Y)=0;
            %Eidereenden(n,Y)=EiderR(Y);
        end
        
        %exponential, SY=0.35, 40% van alle eiders:
%         for Y=1:nyears
%             EiderR(Y)=2.5*2.3*(1-exp(-0.102*BmREat(Y)));
%             EiderS(Y)=2.3*(1-exp(-0.102*BmSEat(Y)))+1.5*2.3*(1-exp(-0.102*BmREat(Y)));
%             EffectSandMining(n,Y)=(1-EiderS(Y)/EiderR(Y));
%         end
        
        %exponential, SY=0.5:
%         for Y=1:nyears
%             EiderR(Y)=2.3*(1-exp(-0.0529*BmREat(Y)));
%             EiderS(Y)=2.3*(1-exp(-0.0529*BmSEat(Y)));
%             EffectSandMining(n,Y)=1-EiderS(Y)/EiderR(Y);
%         end
        
    end
end

%maximale effect in periode van jaren 1 t/m 13
for t=1:nsim
     %Neg(t)=max(EffectSandMining(t,[1:13]));
     Neg(t)=max(EffectEiders(t,[1:13]));
end

Neg2=sort(Neg);
x=0:(1/nsim):(1-(1/nsim));
x2=1-x;
plot(Neg2,x2,'k','LineWidth',2)
hold on

% Imp2=sort(EffectSandMining(1:nsim,2));
% plot(Imp2,x2,'y','LineWidth',2)
% 
% Imp3=sort(EffectSandMining(1:nsim,3));
% plot(Imp3,x2,'m','LineWidth',2)
% 
% Imp4=sort(EffectSandMining(1:nsim,4));
% plot(Imp4,x2,'r','LineWidth',2);

% Imp5=sort(EffectSandMining(1:nsim,5));
% plot(Imp5,x2,'r','LineWidth',2)
% 
% Imp6=sort(EffectSandMining(1:nsim,6));
% plot(Imp6,x2,'g','LineWidth',2)
% 
% Imp7=sort(EffectSandMining(1:nsim,7));
% plot(Imp7,x2,'c','LineWidth',2)

% Imp8=sort(EffectSandMining(1:nsim,8));
% plot(Imp8,x2,'--y')

% Imp9=sort(EffectSandMining(1:nsim,9));
% plot(Imp9,x2,'--m')
% 
% Imp10=sort(EffectSandMining(1:nsim,10));
% plot(Imp10,x2,'--b')

% Imp11=sort(EffectSandMining(1:nsim,11));
% plot(Imp11,x2,'--r')

% Imp12=sort(EffectSandMining(1:nsim,12));
% plot(Imp12,x2,'--g')
% 
% Imp13=sort(EffectSandMining(1:nsim,13));
% plot(Imp13,x2,'--c')



% hist(EffectSandMining(1:nsim,5),50)
    


