% Zoe Aarons
% March 15, 2021
% 12.823, Problem St #1

clear all
close all

%% Run Riley (1946) plankton model

load('Riley_dat.mat'); %observations and forcing functions

%calculate Ph, mean photosynthetic rate in euphotic zone
p = 2.5; %photosynthetic constant
Ph = ((p*I)./(k.*z1)) .* (1-exp(-k.*z1)) .* OneMinusN .* OneMinusV;

%calculate R, respitory rate
R_0 = 0.0175; %rate at 0 deg C
r = 0.069; %constant, rate of change with temp
R = R_0 * exp(r*T);

%Calculate G, rate of zooplankton grazing
g = 0.0075; %rate of reduction of phytoplankton
G = g*Z;

%Use model to solve for seasonal cycle of P
P = zeros(1,25); %build array to store [phytoplankton] values
P(1) = 3.4; %P at t=0
t = 0;
days = zeros(1,25);

for i=1:24
    ln_newP = 15 * (Ph(i) - R(i) - G(i)) + log(P(i)); %calc ln(P(t+15))
    P(i+1) = exp(ln_newP); %calc P(t+15) and store in array
    t = t+15; %timestep
    days(i+1) = t;
end

%plot P over year
figure(1)
plot(days,P)
%hold on
%scatter(Obs_time, Obs_P)
ylabel('Phytoplankton (g C m^{-2})');
ylim([0 40]);
xlim([0 360]);
xlabel('Day')

%Compare to observations
error = zeros(1,length(Obs_P));
P_day = interp1(days,P,0:5:360);
for i = 1:length(Obs_P)
    idx = Obs_time(i)/5 + 1;
    error(i) = abs(P_day(idx) - Obs_P(i))/Obs_P(i);
end
mean_error_baseline = mean(error) * 100;

%% Sensitivity Analysis

%factors (will change all variables to matrix w/ 3 columns)
sensitivity = [0.8, 1, 1.2];
baseline = [1,1,1];

%place variables for sensitivity analysis in cells
vars = cell(5,5);
for i=1:5
    vars{1,i} = p*baseline; vars{2,i} = R_0*baseline;
    vars{3,i} = r*baseline; vars{4,i} = g*baseline;
    vars{5,i} = OneMinusN*baseline;
end
%multiply by sensitivities
vars{1,1} = p*sensitivity; vars{2,2} = R_0*sensitivity; 
vars{3,3} = r*sensitivity; vars{4,4} = g*sensitivity;
vars{5,5} = OneMinusN*sensitivity;

I_base = I*baseline; k_base = k*baseline; z1_base = z1*baseline; 
oneMinusV_base = OneMinusV*baseline;
T_base = T*baseline;
Z_base = Z*baseline;

name = {'p', 'R_0', 'r', 'g', '1-N'}; %names for figure
%table for sensitivity analysis
varNames = {'-20%','Baseline','+20%'};
varTypes = {'double','double','double'};
sens = table('Size',[5 3],'VariableTypes',varTypes,'VariableNames',...
    varNames, 'RowNames', name);

figure(2)
for j = 1:5
    %calculate inputs for model
    Ph = ((vars{1,j}.*I_base)./(k_base.*z1_base)) .* (1-exp(-k_base.*z1_base)) .* vars{5,j} .* oneMinusV_base;
    R = vars{2,j} .* exp(vars{3,j}.*T_base);
    G = vars{4,j}.*Z_base;
    P = zeros(3,25); %build array to store [phytoplankton] values
    P(:,1) = 3.4; %P at t=0
    t = 0;
    days = zeros(1,25);

    %run timestepping
    for i=1:24
        ln_newP = 15 * (Ph(i,:) - R(i,:) - G(i,:)) + log(P(:,i))'; %calc ln(P(t+15))
        P(:,i+1) = exp(ln_newP); %calc P(t+15) and store in array
        t = t+15; %timestep
        days(i+1) = t;
    end
    
    %plot sensitivity analysis
    subplot(3,2,j)
    plot(days,P(1,:), 'LineWidth', 1)
    hold on
    plot(days,P(2,:),'--', 'LineWidth', 1)
    hold on
    plot(days,P(3,:), 'LineWidth', 1)
    legend('-20%', 'Baseline', '+20%', 'Location', 'northwest');
    ylabel('Phytoplankton (g C m^{-2})');
    xlim([0 360]);
    xlabel('Day')
    title(['Sensitivity Analysis: ', name{j}], 'FontSize', 12);
    
    %Compare to observations (error analysis)
    mean_error = zeros(1,3);
    for m = 1:length(mean_error)
        error = zeros(1,length(Obs_P));
        P_mod = interp1(days,P(m,:),0:5:360);
        for l = 1:length(Obs_P)
            idx = Obs_time(l)/5 + 1;
            error(l) = abs(P_mod(idx) - Obs_P(l))/Obs_P(l);
        end
        mean_error(m) = mean(error) * 100;
    end
    
    %put percent error in table
    sens(j,:) = {mean_error(1),mean_error(2),mean_error(3)};
end

%% Question 2 (find g for periodic)
%calculate Ph, mean photosynthetic rate in euphotic zone
p = 2.5; %photosynthetic constant
Ph = ((p*I)./(k.*z1)) .* (1-exp(-k.*z1)) .* OneMinusN .* OneMinusV;

%calculate R, respitory rate
R_0 = 0.0175; %rate at 0 deg C
r = 0.069; %constant, rate of change with temp
R = R_0 * exp(r*T);

%variable to save g value that is closest to periodic
min_diff = 1; min_g = 0;
for g = 0.0074:0.0000001:0.0076
    G = g*Z;

    %Use model to solve for seasonal cycle of P
    P = zeros(1,25); %build array to store [phytoplankton] values
    P(1) = 3.4; %P at t=0
    t = 0;
    days = zeros(1,25);

    for i=1:24
        ln_newP = 15 * (Ph(i) - R(i) - G(i)) + log(P(i)); %calc ln(P(t+15))
        P(i+1) = exp(ln_newP); %calc P(t+15) and store in array
        t = t+15; %timestep
        days(i+1) = t;
    end
    
    %calc difference between start and end of year
    diff = P(25)-P(1);
    if abs(diff)<abs(min_diff) %if new g is closer to exact periodic, save
        min_diff = diff;
        min_g = g;
    end
end

%% Question 3

%allow Z to vary by 20%
Z_rep = repmat(Z,10);
Z_var(:) = Z_rep(:,1) .* (0.8 + 0.4*rand(240,1));
G_rep = min_g*Z_var; %use min g value calculated above (question 2)

Ph_rep = repmat(Ph,10); Ph_rep = Ph_rep(:,1);
R_rep = repmat(R,10); R_rep = R_rep(:,1);

%timestep over 10 years using varying Z value
P_rep = zeros(1,241); %build array to store [phytoplankton] values
P_rep(1) = 3.4; %P at t=0
t = 0;
days = zeros(1,241);

for i=1:240
    ln_newP = 15 * (Ph_rep(i) - R_rep(i) - G_rep(i)) + log(P_rep(i)); %calc ln(P(t+15))
    P_rep(i+1) = exp(ln_newP); %calc P(t+15) and store in array
    t = t+15; %timestep
    days(i+1) = t;
end

%plot P over year
figure(3)
plot(days,P_rep)
ylabel('Phytoplankton (g C m^{-2})');
%ylim([0 40]);
xlim([0 360*10]);
xlabel('Day')
