clc
clear
hold on

%joint pdf van IM1 en Itop:
t=92; %1 april
mu=68+2*(t-60);
sd=mu*0.19;
X=gampdf(0:100,4.27,5.09);
Y=normpdf(0:400,mu,sd);
for j=1:101
    for h=1:301
        A(j,h)=X(j)*Y(h);
    end
end
B=A';
contour(B)

%Combinations of Itop and IM1 for start net primary production
fr=0.1;
fpa=0.45;
Efc=0.25;
Ha=13.3;
eb=0.08;
eot=0.14;
es=0.025;
fs=1.5; %veranderen voor ander jaar/zandwinscenario

for k=1:1010
    IM1(k)=0.1*k;
    
    et1(k)=eb+eot+IM1(k)*es;
    Inec1(k)=16*Efc*Ha*et1(k)/((1-fr)*fpa*(1-exp(-et1(k)*Ha)));
    
    et2(k)=eb+eot+fs*IM1(k)*es;
    Inec2(k)=16*Efc*Ha*et2(k)/((1-fr)*fpa*(1-exp(-et2(k)*Ha)));
end

plot(IM1,Inec2,'b')
plot(IM1,Inec1,'r')