%12.823 Andrew Hirzel Homework 1

load Riley_dat.mat;

%Replicate Figure 21:
p=2.5;%Photosynthetic constant
R0=0.0175;%Respiratory rate at 0 deg C
r=0.069;%rate of change constant of respiratory rate with temperature
g=0.0075;%rate of reduction of phytoplankton by the quantity of zooplankton
Ph=p.*I./(k.*z1).*(1-exp(-k.*z1)).*OneMinusN.*OneMinusV;%Photosynthetic rate
RT=R0*exp(r*T);%Phytoplankton Respiration Rate
G=g*Z;

t=0:15:12*30;%Time step of 15 days, 12 months of 30 days
P0=3.4;
P=Riley(t,P0,Ph,RT,G);
figure(1)
plot(t,P,Obs_time,Obs_P,'*')
xlabel('Time (days)');ylabel('Phytoplankton P (g C m-2)')
title('Original Model, Replicating Riley Figure 21')
legend('Model Results','Observations','Location','NorthEast')

figure(2)
%Sensitivity Analysis for p:
Phlower=(p.*.8).*I./(k.*z1).*(1-exp(-k.*z1)).*OneMinusN.*OneMinusV;% 20% Decrease in p
Phhigher=(p.*1.2).*I./(k.*z1).*(1-exp(-k.*z1)).*OneMinusN.*OneMinusV;% 20% Increase in p
Plower=Riley(t,P0,Phlower,RT,G);
Phigher=Riley(t,P0,Phhigher,RT,G);
subplot(231)
plot(t,P,t,Plower,'r--',t,Phigher,'g--');
xlabel('Time (days)');ylabel('Phytoplankton P (g C m-2)')
title('Sensitivity Analysis for p')
legend('Original','-20% in p','+20% in p','Location','southoutside','Orientation','horizontal')
clear Phlower Phhigher Plower Phigher;

%Sensitivity Analysis for R0:
RTlower=(R0*0.8)*exp(r*T);% 20% Decrease in R0
RThigher=(R0*1.2)*exp(r*T);% 20% Increase in R0
Plower=Riley(t,P0,Ph,RTlower,G);
Phigher=Riley(t,P0,Ph,RThigher,G);
subplot(232)
plot(t,P,t,Plower,'r--',t,Phigher,'g--');
xlabel('Time (days)');ylabel('Phytoplankton P (g C m-2)')
title('Sensitivity Analysis for R_0')
legend('Original','-20% in R_0','+20% in R_0','Location','southoutside','Orientation','horizontal')
clear RTlower RThigher Plower Phigher;

%Sensitivity Analysis for r:
RTlower=R0*exp((r*0.8)*T);% 20% Decrease in r
RThigher=R0*exp((r*1.2)*T);% 20% Increase in r
Plower=Riley(t,P0,Ph,RTlower,G);
Phigher=Riley(t,P0,Ph,RThigher,G);
subplot(233)
plot(t,P,t,Plower,'r--',t,Phigher,'g--');
xlabel('Time (days)');ylabel('Phytoplankton P (g C m-2)')
title('Sensitivity Analysis for r')
legend('Original','-20% in r','+20% in r','Location','southoutside','Orientation','horizontal')
clear RTlower RThigher Plower Phigher;

%Sensitivity Analysis for g:
Glower=(g*0.8)*Z;% 20% Decrease in g
Ghigher=(g*1.2)*Z;% 20% Increase in g
Plower=Riley(t,P0,Ph,RT,Glower);
Phigher=Riley(t,P0,Ph,RT,Ghigher);
subplot(234)
plot(t,P,t,Plower,'r--',t,Phigher,'g--');
xlabel('Time (days)');ylabel('Phytoplankton P (g C m-2)')
title('Sensitivity Analysis for g')
legend('Original','-20% in g','+20% in g','Location','southoutside','Orientation','horizontal')
clear Glower Ghigher Plower Phigher;

%Sensitivity Analysis for 1-N:
Phlower=p.*I./(k.*z1).*(1-exp(-k.*z1)).*(OneMinusN.*.8).*OneMinusV;% 20% Decrease in 1-N
Phhigher=p.*I./(k.*z1).*(1-exp(-k.*z1)).*(OneMinusN.*1.2).*OneMinusV;% 20% Increase in 1-N
OneNnew=OneMinusN*1.2;OneNnew(OneNnew>1)=1;
Phnew=p.*I./(k.*z1).*(1-exp(-k.*z1)).*(OneNnew).*OneMinusV;% 20% Increase Caps at 1
Plower=Riley(t,P0,Phlower,RT,G);
Phigher=Riley(t,P0,Phhigher,RT,G);
Pnew=Riley(t,P0,Phnew,RT,G);
subplot(235)
plot(t,P,t,Plower,'r--',t,Phigher,'g--',t,Pnew,'k--');
xlabel('Time (days)');ylabel('Phytoplankton P (g C m-2)')
title('Sensitivity Analysis for 1-N')
legend('Original','-20% in 1-N','+20% in 1-N','+20%, Limit at 1','Location','southoutside','Orientation','horizontal')
clear OneNnew Phlower Phhigher Phnew Plower Phigher Pnew;

%27% Error in P:
subplot(236)
plot(t,P,Obs_time,Obs_P,'*',t,P*.73,'k--',t,P*1.23,'k--');
xlabel('Time (days)');ylabel('Phytoplankton P (g C m-2)')
title('Error Around Original Model')
legend('Original','Observations','\pm 27% Error in P','Location','southoutside','Orientation','horizontal')


%Find g where P is periodic
gnew=0.0075;%Initialize
for i=1:10000
    Pnew=Riley(t,P0,Ph,RT,gnew*Z);%Calculate P
    %Multiply the value of g by a tenth of the difference between the end
    %value and P0. This is to guide the end P towards the initial value.
    %Once they are the same, the Riley model will be completely periodic.
    gnew=gnew*((Pnew(end)-P0)/P0/10+1);
end
%10 Year Periodic
t10=0:15:12*30*10;%10 years
Ph10=[Ph;Ph;Ph;Ph;Ph;Ph;Ph;Ph;Ph;Ph];
RT10=[RT;RT;RT;RT;RT;RT;RT;RT;RT;RT];
G1=g*Z;G10=[G1;G1;G1;G1;G1;G1;G1;G1;G1;G1];
Pold=Riley(t10,P0,Ph10,RT10,G10);
G1=gnew*Z;G10=[G1;G1;G1;G1;G1;G1;G1;G1;G1;G1];
Pnew=Riley(t10,P0,Ph10,RT10,G10);
figure(3)
plot(t10,Pold,t10,Pnew);axis([0 4000 0 40])
xlabel('Time (days)');ylabel('Phytoplankton P (g C m-2)')
title('10 year Periodic g Model')
legend('Original Model','Periodic Model','Location','northeast')
%Various periodic g values (to different conditions)
gtemp=0.0075;%Initialize
for i=1:10000
    Pnew=Riley(t,P0,Ph*0.8,RT,gtemp*Z);%Calculate P
    %Multiply the value of g by a tenth of the difference between the end
    %value and P0. This is to guide the end P towards the initial value.
    %Once they are the same, the Riley model will be completely periodic.
    gtemp=gtemp*((Pnew(end)-P0)/P0/10+1);
end
disp(['Ph*0.8: g=',num2str(gtemp)])
gtemp=0.0075;%Initialize
for i=1:10000
    Pnew=Riley(t,P0,Ph*1.2,RT,gtemp*Z);%Calculate P
    %Multiply the value of g by a tenth of the difference between the end
    %value and P0. This is to guide the end P towards the initial value.
    %Once they are the same, the Riley model will be completely periodic.
    gtemp=gtemp*((Pnew(end)-P0)/P0/10+1);
end
disp(['Ph*1.2: g=',num2str(gtemp)])

gtemp=0.0075;%Initialize
for i=1:10000
    Pnew=Riley(t,P0,Ph,RT*0.8,gtemp*Z);%Calculate P
    %Multiply the value of g by a tenth of the difference between the end
    %value and P0. This is to guide the end P towards the initial value.
    %Once they are the same, the Riley model will be completely periodic.
    gtemp=gtemp*((Pnew(end)-P0)/P0/10+1);
end
disp(['RT*0.8: g=',num2str(gtemp)])
gtemp=0.0075;%Initialize
for i=1:10000
    Pnew=Riley(t,P0,Ph,RT*1.2,gtemp*Z);%Calculate P
    %Multiply the value of g by a tenth of the difference between the end
    %value and P0. This is to guide the end P towards the initial value.
    %Once they are the same, the Riley model will be completely periodic.
    gtemp=gtemp*((Pnew(end)-P0)/P0/10+1);
end
disp(['RT*1.2: g=',num2str(gtemp)])

%10 Year Random, use 10 year values from periodic
P101=Riley2(t10,P0,Ph10,RT10,G10);%Run 1
P102=Riley2(t10,P0,Ph10,RT10,G10);%Run 2
P103=Riley2(t10,P0,Ph10,RT10,G10);%Run 3
P10all=NaN(100,length(t10));
for i=1:100 %100 replicates of 10 years
    P10all(i,:)=Riley2(t10,P0,Ph10,RT10,G10);
end
P10mean=mean(P10all);%Mean of 100 runs of 10 years
P10std=std(P10all);%STD of 100 runs of 10 years
P10=Riley(t10,P0,Ph10,RT10,G10);%Original Run, 10 years
figure(4)
subplot(121)
plot(t10,P101,t10,P102,t10,P103);axis([0 4000 0 100])
xlabel('Time (days)');ylabel('Phytoplankton P (g C m-2)')
title('10 year Random Zooplankton Grazing Model')
legend('Run 1','Run 2','Run 3','Location','northwest')
subplot(122)
plot(t10,P10,t10,P10mean,t10,P10mean-P10std,'r--',t10,P10mean+P10std,'r--');
axis([0 4000 0 100])
xlabel('Time (days)');ylabel('Phytoplankton P (g C m-2)')
title('Mean and Standard Deviation of 10 Year Random Model (100 replicates)')
legend('Periodic Model','Mean of Random Model','\pm 1 Standard Deviation From Mean','Location','northwest')

function P=Riley(t,P0,Ph,RT,G)
P=P0;
for i=1:length(t)-1
    P(i+1)=P(i)*exp(15*(Ph(i)-RT(i)-G(i)));
end
end

function P=Riley2(t,P0,Ph,RT,G)
P=P0;
for i=1:length(t)-1
    P(i+1)=P(i)*exp(15*(Ph(i)-RT(i)-G(i)*(0.8+0.4*rand)));
end
end