%% JANE WEINSTOCK, 12.823 PSET 01
% Date created: 03.05.2021
% Date updated: 03.17.2021
% Collaborators: Talked with Evie in starting and re: error calculation

load Riley_dat.mat

% input Riley's parameter values
P_t0 = 3.4; % g C m^-2
p = 2.5; % day^-1 (g cal cm^-2 min^-1)^-1
R_0 = 0.0175; % day^-1
r = 0.069; % (deg C)^-1
g = 0.0075; % day^-1 (g C m^-2)^-1

params = [p,R_0,r,g]';
P_est = zeros(24,1);
P_est(1) = P_t0;
for i = 1:23
    G_t = params(4)*Z(i);
    R_t = params(2) * exp(params(3)*T(i));
    P_ht = ((params(1)*I(i))/(k(i)*z1(i))) * (1 - exp(-k(i)*z1(i))) * OneMinusN(i)*OneMinusV(i);
    lnPt15_Minus_lnPt = (P_ht - R_t - G_t)*15;
    lnPt15 = lnPt15_Minus_lnPt + log(P_est(i));
    P_est(i+1) = exp(lnPt15);
end

t = [0:15:345]'; % time in # days after starting value

figure(1)
xlabel('t (days)')
ylim([0 45]) 
ylabel('Phytoplankton (g C m^{-2})')
title('Riley NPZ Model')
hold on
plot(t,P_est,'-o')
plot(Obs_time,Obs_P,'o')
hold off

%% Sensitivity Analysis

P_est_SA_up = [];
par_toggle = 1.2;

for j = 1:4
    params = [p,R_0,r,g]';
    params(j) = params(j)*par_toggle;
    P_est_SA = zeros(24,1);
    P_est_SA(1) = P_t0;
    for i = 1:23
        G_t = params(4)*Z(i);
        R_t = params(2) * exp(params(3)*T(i));
        P_ht = ((params(1)*I(i))/(k(i)*z1(i))) * (1 - exp(-k(i)*z1(i))) * OneMinusN(i)*OneMinusV(i);
        lnPt15_Minus_lnPt = (P_ht - R_t - G_t)*15;
        lnPt15 = lnPt15_Minus_lnPt + log(P_est_SA(i));
        P_est_SA(i+1) = exp(lnPt15);
    end
    P_est_SA_up = [P_est_SA_up P_est_SA];
end

params = [p,R_0,r,g]';
P_est_SA = zeros(24,1);
P_est_SA(1) = P_t0;
OneMinusN_SA = OneMinusN * par_toggle;
for i = 1:23
    G_t = params(4)*Z(i);
    R_t = params(2) * exp(params(3)*T(i));
    P_ht = ((params(1)*I(i))/(k(i)*z1(i))) * (1 - exp(-k(i)*z1(i))) * OneMinusN_SA(i)*OneMinusV(i);
    lnPt15_Minus_lnPt = (P_ht - R_t - G_t)*15;
    lnPt15 = lnPt15_Minus_lnPt + log(P_est_SA(i));
    P_est_SA(i+1) = exp(lnPt15);
end
P_est_SA_up = [P_est_SA_up P_est_SA];

P_est_SA_low = [];
par_togglee = 0.8;

for j = 1:4
    params = [p,R_0,r,g]';
    params(j) = params(j)*par_togglee;
    P_est_SA = zeros(24,1);
    P_est_SA(1) = P_t0;
    for i = 1:23
        G_t = params(4)*Z(i);
        R_t = params(2) * exp(params(3)*T(i));
        P_ht = ((params(1)*I(i))/(k(i)*z1(i))) * (1 - exp(-k(i)*z1(i))) * OneMinusN(i)*OneMinusV(i);
        lnPt15_Minus_lnPt = (P_ht - R_t - G_t)*15;
        lnPt15 = lnPt15_Minus_lnPt + log(P_est_SA(i));
        P_est_SA(i+1) = exp(lnPt15);
    end
    P_est_SA_low = [P_est_SA_low P_est_SA];
end

params = [p,R_0,r,g]';
P_est_SA = zeros(24,1);
P_est_SA(1) = P_t0;
OneMinusN_SA = OneMinusN * par_togglee;
for i = 1:23
    G_t = params(4)*Z(i);
    R_t = params(2) * exp(params(3)*T(i));
    P_ht = ((params(1)*I(i))/(k(i)*z1(i))) * (1 - exp(-k(i)*z1(i))) * OneMinusN_SA(i)*OneMinusV(i);
    lnPt15_Minus_lnPt = (P_ht - R_t - G_t)*15;
    lnPt15 = lnPt15_Minus_lnPt + log(P_est_SA(i));
    P_est_SA(i+1) = exp(lnPt15);
end
P_est_SA_low = [P_est_SA_low P_est_SA];

%% Plots

figure(2)
subplot(2,2,1);
xlabel('t (days)')
ylabel('Phytoplankton (g C m^{-2})')
title('Sensitivity Analysis: p, (1-N)')
hold on
plot(t,P_est,'-')
plot(Obs_time,Obs_P,'ko')
plot(t,P_est_SA_up(:,1),'r--')
plot(t,P_est_SA_low(:,1),'r-')
hold off

subplot(2,2,2);
xlabel('t (days)')
ylabel('Phytoplankton (g C m^{-2})')
title('Sensitivity Analysis: R_0')
hold on
plot(t,P_est,'-')
plot(Obs_time,Obs_P,'ko')
plot(t,P_est_SA_up(:,2),'r--')
plot(t,P_est_SA_low(:,2),'r-')
hold off

subplot(2,2,3);
xlabel('t (days)') 
ylabel('Phytoplankton (g C m^{-2})')
title('Sensitivity Analysis: r')
hold on
plot(t,P_est,'-')
plot(Obs_time,Obs_P,'ko')
plot(t,P_est_SA_up(:,3),'r--')
plot(t,P_est_SA_low(:,3),'r-')
hold off

subplot(2,2,4);
xlabel('t (days)')
ylabel('Phytoplankton (g C m^{-2})')
title('Sensitivity Analysis: g')
hold on
plot(t,P_est,'-')
plot(Obs_time,Obs_P,'ko')
plot(t,P_est_SA_up(:,4),'r--')
plot(t,P_est_SA_low(:,4),'r-')
hold off

% subplot(2,3,5);
% xlabel('t (days)')
% ylim([0 45]) 
% ylabel('Phytoplankton (g C m^{-2})')
% title('Sensitivity Analysis: 1-N')
% hold on
% plot(t,P_est,'-')
% plot(t,P_est_SA_up(:,5),'r--')
% plot(t,P_est_SA_low(:,5),'r-')
% hold off

%% Estimating Error

t_steps = [1:1:345]';
Riley_interp = [interp1(t,P_est,t_steps)];

obs = [Obs_time' Obs_P'];
Riley_comp = [Riley_interp(10) Riley_interp(90) Riley_interp(110) Riley_interp(130) Riley_interp(170) Riley_interp(255)]';

Err1 = sqrt((Obs_P(1)-Riley_interp(10))^2)/Obs_P(1);
Err2 = sqrt((Obs_P(2)-Riley_interp(90))^2)/Obs_P(2);
Err3 = sqrt((Obs_P(3)-Riley_interp(110))^2)/Obs_P(3);
Err4 = sqrt((Obs_P(4)-Riley_interp(130))^2)/Obs_P(4);
Err5 = sqrt((Obs_P(5)-Riley_interp(170))^2)/Obs_P(5);
Err6 = sqrt((Obs_P(6)-Riley_interp(255))^2)/Obs_P(6);

Riley_error = mean([Err1,Err2,Err3,Err4,Err5,Err6])*100;

Riley_SA_interp_up = interp1(t,P_est_SA_up,t_steps);
Riley_SA_interp_low = interp1(t,P_est_SA_low,t_steps);

Riley_error_SA_up = zeros(4,1);
Riley_error_SA_low = zeros(4,1);
for i = 1:4
    Err1u = sqrt((Obs_P(1)-Riley_SA_interp_up(10,i))^2)/Obs_P(1);
    Err2u = sqrt((Obs_P(2)-Riley_SA_interp_up(90,i))^2)/Obs_P(2);
    Err3u = sqrt((Obs_P(3)-Riley_SA_interp_up(110,i))^2)/Obs_P(3);
    Err4u = sqrt((Obs_P(4)-Riley_SA_interp_up(130,i))^2)/Obs_P(4);
    Err5u = sqrt((Obs_P(5)-Riley_SA_interp_up(170,i))^2)/Obs_P(5);
    Err6u = sqrt((Obs_P(6)-Riley_SA_interp_up(255,i))^2)/Obs_P(6);
    Riley_error_SA_up(i) = mean([Err1u,Err2u,Err3u,Err4u,Err5u,Err6u])*100;
    
    Err1 = sqrt((Obs_P(1)-Riley_SA_interp_low(10,i))^2)/Obs_P(1);
    Err2 = sqrt((Obs_P(2)-Riley_SA_interp_low(90,i))^2)/Obs_P(2);
    Err3 = sqrt((Obs_P(3)-Riley_SA_interp_low(110,i))^2)/Obs_P(3);
    Err4 = sqrt((Obs_P(4)-Riley_SA_interp_low(130,i))^2)/Obs_P(4);
    Err5 = sqrt((Obs_P(5)-Riley_SA_interp_low(170,i))^2)/Obs_P(5);
    Err6 = sqrt((Obs_P(6)-Riley_SA_interp_low(255,i))^2)/Obs_P(6);
    Riley_error_SA_low(i) = mean([Err1,Err2,Err3,Err4,Err5,Err6])*100;
end

%% Periodic Phytos

g_adj = 0.007435;
params = [p,R_0,r,g_adj]';
P_estt = zeros(240,1);
P_estt(1) = P_t0;
P_estt_orig = zeros(240,1);
P_estt_orig(1) = P_t0;

Z_rep = repmat(Z,10,1);
T_rep = repmat(T,10,1);
I_rep = repmat(I,10,1);
k_rep = repmat(k,10,1);
z1_rep = repmat(z1,10,1);
OneMinusN_rep = repmat(OneMinusN,10,1);
OneMinusV_rep = repmat(OneMinusV,10,1);

for i = 1:239
    Z_var = Z_rep * (0.8 + 0.4*rand);
    % Z_var = Z_rep; %toggle this on and the line above off to plot without
        % random variation in Z
    G_t = params(4)*Z_var(i);
    G_t_orig = g*Z_rep(i);
    R_t = params(2) * exp(params(3)*T_rep(i));
    P_ht = ((params(1)*I_rep(i))/(k_rep(i)*z1_rep(i))) * (1 - exp(-k_rep(i)*z1_rep(i))) * OneMinusN_rep(i)*OneMinusV_rep(i);
    lnPt15_Minus_lnPt = (P_ht - R_t - G_t)*15;
    lnPt15_Minus_lnPt_orig = (P_ht - R_t - G_t_orig)*15;
    lnPt15 = lnPt15_Minus_lnPt + log(P_estt(i));
    lnPt15_orig = lnPt15_Minus_lnPt_orig + log(P_estt_orig(i));
    P_estt(i+1) = exp(lnPt15);
    P_estt_orig(i+1) = exp(lnPt15_orig);
end

t = [0:15:3585]'; % time in # days after starting value

figure(3)
xlabel('t (days)') 
ylabel('Phytoplankton (g C m^{-2})')
title('Random Variation in Z')
hold on
plot(t,P_estt_orig,'-')
plot(t,P_estt,'-','Linewidth',2)
hold off


