% 12.823 Problem Set 1 - Evie Fachon, 3/15/21
% based on Riley et al. 1946
%% Set up model

clear; clc;
load('Riley_dat.mat');

%fixed conditions
p = 2.5; %day^-1 (g cal cm^-2 min^-1)^-1 : photosynthetic constant
R0 = 0.0175; %day^-1 : respiration rate at 0 degrees C
r = 0.069; %degC^-1 : constant for rate of change of respiration
g = 0.0075; %day^-1 (g C m^-2)^-1 : rate of phytoplankton reduction per unit grazers
tstep = 15; %days

%initial conditions
P = 3.4; %g C m^-2 (Initial plankton population)
t = 0;
i= 1;

%set up and run simulation, save outputs in vector format for plotting

time = [t]
population = [P]

while i<=24
    Pht = (p*I(i))/(k(i)*z1(i))*(1-exp(-k(i)*z1(i)))*OneMinusN(i)*OneMinusV(i);
    Rt = (R0*exp(r*T(i))); %Rate of respiration
    Gt = (g*Z(i)); %Grazing rate
    Pnew = exp((tstep*(Pht-Rt-Gt)))*P
    P = Pnew
    t = t + tstep
    i = i + 1
    time = [time t]
    population = [population P]
end

figure
plot(time,population,'b-o')
xlabel('Day')
ylabel('Phytoplankton (g C m^-2)')
hold on;
scatter(Obs_time, Obs_P, 'k', 'filled')

%% Sensitivity Analyses

%p
P = 3.4; t = 0; i= 1;
population_p80 = [P]

while i<=24
    Pht = ((0.8*p)*I(i))/(k(i)*z1(i))*(1-exp(-k(i)*z1(i)))*OneMinusN(i)*OneMinusV(i);
    Rt = (R0*exp(r*T(i))); %Rate of respiration
    Gt = (g*Z(i)); %Grazing rate
    Pnew = exp((tstep*(Pht-Rt-Gt)))*P
    P = Pnew
    t = t + tstep
    i = i + 1
    population_p80 = [population_p80 P]
end

P = 3.4; t = 0; i= 1;
population_p120 = [P]

while i<=24
    Pht = ((1.2*p)*I(i))/(k(i)*z1(i))*(1-exp(-k(i)*z1(i)))*OneMinusN(i)*OneMinusV(i);
    Rt = (R0*exp(r*T(i))); %Rate of respiration
    Gt = (g*Z(i)); %Grazing rate
    Pnew = exp((tstep*(Pht-Rt-Gt)))*P
    P = Pnew
    t = t + tstep
    i = i + 1
    population_p120 = [population_p120 P]
end

%R0
P = 3.4; t = 0; i= 1;
population_R080 = [P]

while i<=24
    Pht = (p*I(i))/(k(i)*z1(i))*(1-exp(-k(i)*z1(i)))*OneMinusN(i)*OneMinusV(i);
    Rt = ((R0*0.8)*exp(r*T(i))); %Rate of respiration
    Gt = (g*Z(i)); %Grazing rate
    Pnew = exp((tstep*(Pht-Rt-Gt)))*P
    P = Pnew
    t = t + tstep
    i = i + 1
    population_R080 = [population_R080 P]
end

P = 3.4; t = 0; i= 1;
population_R0120 = [P]

while i<=24
    Pht = (p*I(i))/(k(i)*z1(i))*(1-exp(-k(i)*z1(i)))*OneMinusN(i)*OneMinusV(i);
    Rt = ((R0*1.2)*exp(r*T(i))); %Rate of respiration
    Gt = (g*Z(i)); %Grazing rate
    Pnew = exp((tstep*(Pht-Rt-Gt)))*P
    P = Pnew
    t = t + tstep
    i = i + 1
    population_R0120 = [population_R0120 P]
end

%r
P = 3.4; t = 0; i= 1;
population_r80 = [P]

while i<=24
    Pht = (p*I(i))/(k(i)*z1(i))*(1-exp(-k(i)*z1(i)))*OneMinusN(i)*OneMinusV(i);
    Rt = (R0*exp((r*.8)*T(i))); %Rate of respiration
    Gt = (g*Z(i)); %Grazing rate
    Pnew = exp((tstep*(Pht-Rt-Gt)))*P
    P = Pnew
    t = t + tstep
    i = i + 1
    population_r80 = [population_r80 P]
end

P = 3.4; t = 0; i= 1;
population_r120 = [P]

while i<=24
    Pht = (p*I(i))/(k(i)*z1(i))*(1-exp(-k(i)*z1(i)))*OneMinusN(i)*OneMinusV(i);
    Rt = (R0*exp((1.2*r)*T(i))); %Rate of respiration
    Gt = (g*Z(i)); %Grazing rate
    Pnew = exp((tstep*(Pht-Rt-Gt)))*P
    P = Pnew
    t = t + tstep
    i = i + 1
    population_r120 = [population_r120 P]
end

%g

P = 3.4; t = 0; i= 1;
population_g80 = [P]

while i<=24
    Pht = (p*I(i))/(k(i)*z1(i))*(1-exp(-k(i)*z1(i)))*OneMinusN(i)*OneMinusV(i);
    Rt = (R0*exp(r*T(i))); %Rate of respiration
    Gt = (0.8*g*Z(i)); %Grazing rate
    Pnew = exp((tstep*(Pht-Rt-Gt)))*P
    P = Pnew
    t = t + tstep
    i = i + 1
    population_g80 = [population_g80 P]
end

P = 3.4; t = 0; i= 1;
population_g120 = [P]

while i<=24
    Pht = (p*I(i))/(k(i)*z1(i))*(1-exp(-k(i)*z1(i)))*OneMinusN(i)*OneMinusV(i);
    Rt = (R0*exp(r*T(i))); %Rate of respiration
    Gt = (1.2*g*Z(i)); %Grazing rate
    Pnew = exp((tstep*(Pht-Rt-Gt)))*P
    P = Pnew
    t = t + tstep
    i = i + 1
    population_g120 = [population_g120 P]
end

%1-N

P = 3.4; t = 0; i= 1;
population_N80 = [P]

while i<=24
    Pht = (p*I(i))/(k(i)*z1(i))*(1-exp(-k(i)*z1(i)))*(0.8*OneMinusN(i))*OneMinusV(i);
    Rt = (R0*exp(r*T(i))); %Rate of respiration
    Gt = (g*Z(i)); %Grazing rate
    Pnew = exp((tstep*(Pht-Rt-Gt)))*P
    P = Pnew
    t = t + tstep
    i = i + 1
    population_N80 = [population_N80 P]
end

P = 3.4; t = 0; i= 1;
population_N120 = [P]

while i<=24
    Pht = (p*I(i))/(k(i)*z1(i))*(1-exp(-k(i)*z1(i)))*(1.2*OneMinusN(i))*OneMinusV(i);
    Rt = (R0*exp(r*T(i))); %Rate of respiration
    Gt = (g*Z(i)); %Grazing rate
    Pnew = exp((tstep*(Pht-Rt-Gt)))*P
    P = Pnew
    t = t + tstep
    i = i + 1
    population_N120 = [population_N120 P]
end

%% Multipanel figure
subplot(5,1,1)
plot(time,population,'b-o', time,population_p80,'r-o', time,population_p120,'g-o')
legend('100%','80%','120%')
title('ST: Photosynthetic Constant (p)')
subplot(5,1,2)
plot(time,population,'b-o', time,population_R080,'r-o', time,population_R0120,'g-o')
title('ST: R0')
subplot(5,1,3)
plot(time,population,'b-o', time,population_r80,'r-o', time,population_r120,'g-o')
title('ST: r')
ylabel('Phytoplankton (g C m^-2)')
subplot(5,1,4)
plot(time,population,'b-o', time,population_g80,'r-o', time,population_g120,'g-o')
title('ST: g')
subplot(5,1,5)
plot(time,population,'b-o', time,population_N80,'r-o', time,population_N120,'g-o')
title('ST: 1-N')
xlabel('Day')

%% interpolate to calculate error 

dayvec = [1:1:360]
population_interp = interp1(time, population, dayvec);

error_vec = [];
for i=1:length(Obs_time);
    error_vec = [error_vec (abs((Obs_P(i) - population_interp(Obs_time(i))))/(Obs_P(i)))];
end
mean(error_vec)

%calculate errors for different parameters in the sensitivity tests

SensitivityTable = table(population_p80, population_p120, population_R080, population_R0120, population_r80, population_r120, population_g80, population_g120, population_N80, population_N120)
sens_vec = [];

for col = 1 : width(SensitivityTable)
  senstest = table2array(SensitivityTable(:, col));
  sens_interp = interp1(time, senstest, dayvec);
  error_vec = [];
  for i = 1:length(Obs_time);
      error = ((Obs_P(i) - sens_interp(Obs_time(i)))/(Obs_P(i)))
      error_vec = [error_vec error];
  end
  sens_vec = [sens_vec mean(error_vec)];
end

%% g for periodicity

Pi = 3.4
Pf = 50
gper = 0.00750
while abs(Pf - Pi) > 0.002
    gper = gper - 0.000001
    i = 1
    P = Pi
    while i<=24
        Pht = (p*I(i))/(k(i)*z1(i))*(1-exp(-k(i)*z1(i)))*OneMinusN(i)*OneMinusV(i);
        Rt = (R0*exp(r*T(i))); %Rate of respiration
        Gt = (gper*Z(i)); %Grazing rate
        Pf = exp((tstep*(Pht-Rt-Gt)))*P
        P = Pf
        t = t + tstep
        i = i + 1
    end
end

%changing other parameters (manipulated various parameters within the 
%equations and re-ran same section of code)
Pi = 3.4
Pf = 50
gper = 0.00750
diff = 10000000000
while abs(Pf - Pi) > 0.1
    gper = gper + 0.00001
    i = 1;
    P = Pi;
    while i<=24;
        Pht = (p*I(i))/(k(i)*z1(i))*(1-exp(-k(i)*z1(i)))*1.2*OneMinusN(i)*OneMinusV(i);
        Rt = (R0*exp(r*T(i))); %Rate of respiration
        Gt = (gper*Z(i)); %Grazing rate
        Pf = exp((tstep*(Pht-Rt-Gt)))*P
        P = Pf
        t = t + tstep
        i = i + 1
    end
    if abs(Pf-Pi) > diff
        break
    end
    diff = abs(Pf - Pi)
end



%% 10 year time series
gper = 0.00743

rand1 = (randi([0 400],1,240)+800)/1000
rand2 = (randi([0 400],1,240)+800)/1000
rand3 = (randi([0 400],1,240)+800)/1000
P=3.4
j=1
i=1
t=0
population_10yr_2 = [P]
time_10yr = [t]

while j < 240
    Pht = (p*I(i))/(k(i)*z1(i))*(1-exp(-k(i)*z1(i)))*OneMinusN(i)*OneMinusV(i);
    Rt = (R0*exp(r*T(i))); %Rate of respiration
    Gt = (gper*(Z(i)*rand2(j))); %Grazing rate
    Pnew = exp((tstep*(Pht-Rt-Gt)))*P
    P = Pnew
    t = t + tstep
    i = i + 1
    time_10yr = [time_10yr t]
    population_10yr_2 = [population_10yr_2 P]
    if i > 24
       i=i-24
    end
    j = j + 1
end

% Figure 3
figure
subplot(3,1,1)
plot(time_10yr,population_10yr_1)
yyaxis right
scatter(time_10yr, rand1)
ylim([0.8 1.2])
subplot(3,1,2)
plot(time_10yr,population_10yr_2)
ylabel('Phytoplankton (g C m^-2)')
yyaxis right
ylabel('Random g (day^-1 (g C m^-2)^-1)')
scatter(time_10yr, rand2)
subplot(3,1,3)
plot(time_10yr,population_10yr_3)
yyaxis right
scatter(time_10yr, rand3)