% Kevin Archibald
% 12.823 Problem Set 1

%% Part 1: Run Riley (1946) model)

load Riley_dat

p = 2.5;
R0 = 0.0175;
r = 0.069;
g = 0.0075;

P = riley(p, R0, r, g, I, k, OneMinusN, OneMinusV, T, Z, z1);

% plot model results vs. observations

mod_time = 1:15:360;
figure;
plot(mod_time, P, '-k'); hold on;
scatter(Obs_time, Obs_P, 'filled'); hold off;
xlabel('Time (days)');
ylabel('Phytoplankton Biomass (g C m^{-2})');
legend('Theoretical Curve', 'Observed Population');

%% Part 2: Sensitivity Analysis

% sensitivity to photosynthetic constant (p)

lowp = 0.8*p;
highp = 1.2*p;
P_lowp = riley(lowp, R0, r, g, I, k, OneMinusN, OneMinusV, T, Z, z1);
P_highp = riley(highp, R0, r, g, I, k, OneMinusN, OneMinusV, T, Z, z1);

figure;
subplot(3,2,1)
plot(mod_time, P_lowp, '--b', mod_time, P, '-k', mod_time, P_highp, '--r');
xlabel('Time (days)');
ylabel('Phytoplankton Biomass (g C m^{-2})');
legend('p = 2', 'p = 2.5', 'p = 3');

% sensitivity to respiration rate (R0)

lowR0 = 0.8*R0;
highR0 = 1.2*R0;
P_lowR0 = riley(p, lowR0, r, g, I, k, OneMinusN, OneMinusV, T, Z, z1);
P_highR0 = riley(p, highR0, r, g, I, k, OneMinusN, OneMinusV, T, Z, z1);

subplot(3,2,2)
plot(mod_time, P_lowR0, '--b', mod_time, P, '-k', mod_time, P_highR0, '--r');
xlabel('Time (days)');
ylabel('Phytoplankton Biomass (g C m^{-2})');
legend('R0 = 0.014', 'R0 = 0.0175', 'R0 = 0.021');

% sensitivity to temperature dependence (r)

lowr = 0.8*r;
highr = 1.2*r;
P_lowr = riley(p, R0, lowr, g, I, k, OneMinusN, OneMinusV, T, Z, z1);
P_highr = riley(p, R0, highr, g, I, k, OneMinusN, OneMinusV, T, Z, z1);

subplot(3,2,3)
plot(mod_time, P_lowr, '--b', mod_time, P, '-k', mod_time, P_highr, '--r');
xlabel('Time (days)');
ylabel('Phytoplankton Biomass (g C m^{-2})');
legend('r = 0.0552', 'r = 0.069', 'r = 0.0828');

% sensitivity to grazing rate (g)

lowg = 0.8*g;
highg = 1.2*g;
P_lowg = riley(p, R0, r, lowg, I, k, OneMinusN, OneMinusV, T, Z, z1);
P_highg = riley(p, R0, r, highg, I, k, OneMinusN, OneMinusV, T, Z, z1);

subplot(3,2,4)
plot(mod_time, P_lowg, '--b', mod_time, P, '-k', mod_time, P_highg, '--r');
xlabel('Time (days)');
ylabel('Phytoplankton Biomass (g C m^{-2})');
legend('g = 0.006', 'g = 0.0075', 'g = 0.009');

% sensitivity to nutrient depletion(1-N)

lowOneMinusN = 0.8*OneMinusN;
highOneMinusN = 1.2*OneMinusN;
P_lowOneMinusN = riley(p, R0, r, g, I, k, lowOneMinusN, OneMinusV, T, Z, z1);
P_highOneMinusN = riley(p, R0, r, g, I, k, highOneMinusN, OneMinusV, T, Z, z1);

subplot(3,2,5)
plot(mod_time, P_lowOneMinusN, '--b', mod_time, P, '-k', mod_time, P_highOneMinusN, '--r');
xlabel('Time (days)');
ylabel('Phytoplankton Biomass (g C m^{-2})');
legend('(1-N) 20% smaller', '(1-N)', '(1-N) 20% larger');

% compute the average percent error for each sensitivity experiment

sens_exp = [P P_lowp P_highp P_lowR0 P_highR0 P_lowr P_highr P_lowOneMinusN P_highOneMinusN];
error = mean(abs(interp1(mod_time, sens_exp, Obs_time) - repmat(Obs_P', 1, 9))./interp1(mod_time, sens_exp, Obs_time));

%% Part 3: Periodic Solutions

g = 0.005;
P_periodic = riley(p, R0, r, g, I, k, OneMinusN, OneMinusV, T, Z, z1);

% increase value of g until the solution has periodic seasonality

while P_periodic(end) > P_periodic(1)
    g = g + 1e-6;
    P_periodic = riley(p, R0, r, g, I, k, OneMinusN, OneMinusV, T, Z, z1);
end

g_periodic = g;

% repeat this calculation for values of p, R0, r, and 1-N increased by 20%

g = 0.005;
P_periodic = riley(highp, R0, r, g, I, k, OneMinusN, OneMinusV, T, Z, z1);
while P_periodic(end) > P_periodic(1)
    g = g + 1e-6;
    P_periodic = riley(highp, R0, r, g, I, k, OneMinusN, OneMinusV, T, Z, z1);
end
g_periodic(2) = g;

g = 0.005;
P_periodic = riley(p, highR0, r, g, I, k, OneMinusN, OneMinusV, T, Z, z1);
while P_periodic(end) > P_periodic(1)
    g = g + 1e-6;
    P_periodic = riley(p, highR0, r, g, I, k, OneMinusN, OneMinusV, T, Z, z1);
end
g_periodic(3) = g;

g = 0.005;
P_periodic = riley(p, R0, highr, g, I, k, OneMinusN, OneMinusV, T, Z, z1);
while P_periodic(end) > P_periodic(1)
    g = g + 1e-6;
    P_periodic = riley(p, R0, highr, g, I, k, OneMinusN, OneMinusV, T, Z, z1);
end
g_periodic(4) = g;

g = 0.005;
P_periodic = riley(p, R0, r, g, I, k, highOneMinusN, OneMinusV, T, Z, z1);
while P_periodic(end) > P_periodic(1)
    g = g + 1e-6;
    P_periodic = riley(p, R0, r, g, I, k, highOneMinusN, OneMinusV, T, Z, z1);
end
g_periodic(5) = g;

%% Part 4: Random perturbations to Z

g = g_periodic(1);

I_10yr = repmat(I, 10, 1);  % construct 10 year time series for drivers
k_10yr = repmat(k, 10, 1);
OneMinusN_10yr = repmat(OneMinusN, 10, 1);
OneMinusV_10yr = repmat(OneMinusV, 10, 1);
T_10yr = repmat(T, 10, 1);
z1_10yr = repmat(z1, 10, 1);
Z_10yr = repmat(Z, 10, 1);

Z_10yr = Z_10yr.*(0.8 + 0.4*rand(240,1));

P_10yr = riley(p, R0, r, g, I_10yr, k_10yr, OneMinusN_10yr, OneMinusV_10yr, T_10yr, Z_10yr, z1_10yr);

figure;
plot(P_10yr);
xlabel('Time (days)');
ylabel('Phytoplankton Biomass')




