clear
clc

%Growth factor
Gr=[10.8 
    3.36 
    1.86 
    1.37 
    1.17 
    1.09];

BmRcr=0.12;

%-------------------------------------------------------------------------
% Biomass of shellfish 
%-------------------------------------------------------------------------
ny=10005;

for Y=1:ny
    BmP(1,Y)=recrrnd*BmRcr;
    MrtS(Y)=mrtsrnd;
    MrtW(Y)=mrtwrnd_35;
    BmP(2,Y)=BmP(1,Y)*Gr(1);
end
for A=3:7
    BmP(A,1)=1; %dummy
    for Y=2:ny
        BmP(A,Y)=BmP(A-1,Y-1)*Gr(A-1)*MrtS(Y)*MrtW(Y);
    end
end

for Y=6:ny
    BmTot(Y-5)=BmP(2,Y)+BmP(3,Y)+BmP(4,Y)+BmP(5,Y)+BmP(6,Y)+BmP(7,Y);
    BmEat(Y-5)=BmP(3,Y)+BmP(4,Y)+BmP(5,Y)+BmP(6,Y)+BmP(7,Y);
end

%MeanTot=mean(BmTot)
MeanEat=mean(BmEat)
StdEat=std(BmEat)

for n=1:10000
    if BmEat(n)>MeanEat
        N(n)=1;
    else
        N(n)=0;
    end
end

Perc=sum(N)/10000
Max=max(BmEat)
Min=min(BmEat)
BmEat2=BmEat';
BmEat3=BmEat2(100:199);
plot(BmEat3)

Popb=BmP';
Pb=Popb(100:199,3:7);
