% Problem Set 1
%Kharis Schrage


%Ph photosynthetic rate = radiation 
%z1 depth
%k extinction coefficient
%I radiation 
%G rate of grazing 
%g rate of reduction of phyto by unit
%Z quantity of zooplankton 
%R is R.0e^rT temperature effect 
%R.T is respiratory rate at any temp T
%Ro rate at 0 C
%r is constant expressing rate of change of respriatory rate with temp 
%constants
Ro = 0.0175;
r = 0.069;
g = 0.0075;
p = 2.5;
P = zeros(24,1)
P(1)= 3.4; %Po
for t=1:24
    RT(t) = Ro*exp(r*T(t));
    G(t)= g*Z(t);
    Ph(t) = ((p*I(t))/(k(t)*z1(t)))*(1-exp(-k(t)*z1(t)))*OneMinusN(t)*OneMinusV(t);
    P(t+1)=P(t)*exp(15*(Ph(t)- RT(t) - G(t)));
end


%calc error 
Obs_time_adj = (Obs_time./15)+1
Pmod = interp1(P,Obs_time_adj,"spline")
OG_error = (sum((abs(Obs_P-Pmod))./Obs_P)/6)*100
%41.6% error w/ linear
%with spline, 44.5% error

%plot that with observed p and time
time = [1:15:375]
figure(1)
scatter(Obs_time,Obs_P)

hold on
plot(time,P)
title('P vs observed')

hold on
scatter(Obs_time, Pmod)
hold off




%sensitivity analysis 
%p 
p = 2.5;
pU =1.2*p
PpU = zeros(24,1)
PpU(1)= 3.4; %Po
for t=1:24
    RT(t) = Ro*exp(r*T(t));
    G(t)= g*Z(t);
    Ph(t) = ((pU*I(t))/(k(t)*z1(t)))*(1-exp(-k(t)*z1(t)))*OneMinusN(t)*OneMinusV(t);
    PpU(t+1)=PpU(t)*exp(15*(Ph(t)- RT(t) - G(t)));
end

pD =0.8*p
PpD = zeros(24,1)
PpD(1)= 3.4; %Po
for t=1:24
    RT(t) = Ro*exp(r*T(t));
    G(t)= g*Z(t);
    Ph(t) = ((pD*I(t))/(k(t)*z1(t)))*(1-exp(-k(t)*z1(t)))*OneMinusN(t)*OneMinusV(t);
    PpD(t+1)=PpD(t)*exp(15*(Ph(t)- RT(t) - G(t)));
end
figure(2)
pvals = tiledlayout(1,2);
nexttile
plot(time,PpU)
title('p increased')
hold on
plot(time, P)
hold off
nexttile
plot(time, PpD)
hold on 
title('p decreased')
plot(time, P)
hold off


%Ro 
Ro = 0.0175;
RoU =1.2*Ro
PRoU = zeros(24,1)
PRoU(1)= 3.4; %Po
for t=1:24
    RT(t) = RoU*exp(r*T(t));
    G(t)= g*Z(t);
    Ph(t) = ((p*I(t))/(k(t)*z1(t)))*(1-exp(-k(t)*z1(t)))*OneMinusN(t)*OneMinusV(t);
    PRoU(t+1)=PRoU(t)*exp(15*(Ph(t)- RT(t) - G(t)));
end

RoD =0.8*Ro
PRoD = zeros(24,1)
PRoD(1)= 3.4; %Po
for t=1:24
    RT(t) = RoD*exp(r*T(t));
    G(t)= g*Z(t);
    Ph(t) = ((p*I(t))/(k(t)*z1(t)))*(1-exp(-k(t)*z1(t)))*OneMinusN(t)*OneMinusV(t);
    PRoD(t+1)=PRoD(t)*exp(15*(Ph(t)- RT(t) - G(t)));
end
figure(3)
Rovals = tiledlayout(1,2);
nexttile
plot(time,PRoU)
hold on
plot(time, P)
title('Ro increased')
hold off
nexttile
plot(time, PRoD)
hold on 
plot(time, P)
title('Ro decreased')
hold off


%r
r = 0.069;
rU =1.2*r
PrU = zeros(24,1)
PrU(1)= 3.4; %Po
for t=1:24
    RT(t) = R*exp(rU*T(t));
    G(t)= g*Z(t);
    Ph(t) = ((p*I(t))/(k(t)*z1(t)))*(1-exp(-k(t)*z1(t)))*OneMinusN(t)*OneMinusV(t);
    PrU(t+1)=PrU(t)*exp(15*(Ph(t)- RT(t) - G(t)));
end

rD =0.8*r
PrD = zeros(24,1)
PrD(1)= 3.4; %Po
for t=1:24
    RT(t) = Ro*exp(rD*T(t));
    G(t)= g*Z(t);
    Ph(t) = ((p*I(t))/(k(t)*z1(t)))*(1-exp(-k(t)*z1(t)))*OneMinusN(t)*OneMinusV(t);
    PrD(t+1)=PrD(t)*exp(15*(Ph(t)- RT(t) - G(t)));
end
figure(4)
rvals = tiledlayout(1,2);
nexttile
plot(time,PrU)
hold on
plot(time, P)
title('r increased')
hold off
nexttile
plot(time, PrD)
hold on 
plot(time, P)
title('r decreased')
hold off


%g
g = 0.0075;
gU =1.2*g
PgU = zeros(24,1)
PgU(1)= 3.4; %Po
for t=1:24
    RT(t) = Ro*exp(r*T(t));
    G(t)= gU*Z(t);
    Ph(t) = ((p*I(t))/(k(t)*z1(t)))*(1-exp(-k(t)*z1(t)))*OneMinusN(t)*OneMinusV(t);
    PgU(t+1)=PgU(t)*exp(15*(Ph(t)- RT(t) - G(t)));
end

gD =0.8*g
PgD = zeros(24,1)
PgD(1)= 3.4; %Po
for t=1:24
    RT(t) = Ro*exp(r*T(t));
    G(t)= gD*Z(t);
    Ph(t) = ((p*I(t))/(k(t)*z1(t)))*(1-exp(-k(t)*z1(t)))*OneMinusN(t)*OneMinusV(t);
    PgD(t+1)=PgD(t)*exp(15*(Ph(t)- RT(t) - G(t)));
end
figure(5)
gvals = tiledlayout(1,2);
nexttile
plot(time,PgU)
hold on
plot(time, P)
title('g increased')
hold off
nexttile
plot(time, PgD)
hold on 
plot(time, P)
title('g decreased')
hold off

%(1-N)

PNU = zeros(24,1)
PNU(1)= 3.4; %Po
for t=1:24
    NU(t)=1.2*OneMinusN(t);
    RT(t) = Ro*exp(r*T(t));
    G(t)= g*Z(t);
    Ph(t) = ((p*I(t))/(k(t)*z1(t)))*(1-exp(-k(t)*z1(t)))*NU(t)*OneMinusV(t);
    PNU(t+1)=PNU(t)*exp(15*(Ph(t)- RT(t) - G(t)));
end

PND = zeros(24,1);
PND(1)= 3.4; %Po
for t=1:24
    ND(t)=0.8*OneMinusN(t);
    RT(t) = Ro*exp(r*T(t));
    G(t)= g*Z(t);
    Ph(t) = ((p*I(t))/(k(t)*z1(t)))*(1-exp(-k(t)*z1(t)))*ND(t)*OneMinusV(t);
    PND(t+1)=PND(t)*exp(15*(Ph(t)- RT(t) - G(t)));
end
figure(6)
Nvals = tiledlayout(1,2);
nexttile
plot(time,PNU)
hold on
plot(time, P)
title('(1-N) increased')
hold off
nexttile
plot(time, PND)
hold on 
plot(time, P)
title('(1-N) decreased')
hold off


%error based on model
SE_pU =(sum((PpU-P)./P)/25)*100
SE_pD =(sum((PpD-P)./P)/25)*100
SE_RoU =(sum((PRoU-P)./P)/25)*100
SE_RoD =(sum((PRoD-P)./P)/25)*100
SE_rU =(sum((PrU-P)./P)/25)*100
SE_rD =(sum((PrD-P)./P)/25)*100
SE_gU =(sum((PgU-P)./P)/25)*100
SE_gD =(sum((PgD-P)./P)/25)*100
SE_NU =(sum((PNU-P)./P)/25)*100
SE_ND =(sum((PND-P)./P)/25)*100

%error based on Observed values
PpU_a = interp1(PpU,Obs_time_adj);
PpU_err = (sum((abs(Obs_P-PpU_a))./Obs_P)/6)*100
PpD_a = interp1(PpD,Obs_time_adj);
PpD_err = (sum((abs(Obs_P-PpD_a))./Obs_P)/6)*100

PRoU_a = interp1(PRoU,Obs_time_adj);
PRoU_err = (sum((abs(Obs_P-PRoU_a))./Obs_P)/6)*100
PRoD_a = interp1(PRoD,Obs_time_adj);
PRoD_err = (sum((abs(Obs_P-PRoD_a))./Obs_P)/6)*100

PrU_a = interp1(PrU,Obs_time_adj);
PrU_err = (sum((abs(Obs_P-PrU_a))./Obs_P)/6)*100
PrD_a = interp1(PrD,Obs_time_adj);
PrD_err = (sum((abs(Obs_P-PrD_a))./Obs_P)/6)*100

PgU_a = interp1(PgU,Obs_time_adj);
PgU_err = (sum((abs(Obs_P-PgU_a))./Obs_P)/6)*100
PgD_a = interp1(PgD,Obs_time_adj);
PgD_err = (sum((abs(Obs_P-PgD_a))./Obs_P)/6)*100

PNU_a = interp1(PNU,Obs_time_adj);
PNU_err = (sum((abs(Obs_P-PNU_a))./Obs_P)/6)*100
PND_a = interp1(PND,Obs_time_adj);
PND_err = (sum((abs(Obs_P-PND_a))./Obs_P)/6)*100


% #2 Periodicity of P

P_equal = zeros(24,1);
P_equal(1)= 3.4; %Po
g_P25 = zeros(5001,1);
gs =0.0073:.0000001:0.0078;
for g = 1:length(gs);
    gs(g);
    for t=1:24
        RT(t) = Ro*exp(r*T(t));
        G(t)= gs(g)*Z(t);
        Ph(t) = ((p*I(t))/(k(t)*z1(t)))*(1-exp(-k(t)*z1(t)))*OneMinusN(t)*OneMinusV(t);
        P_equal(t+1)=P_equal(t)*exp(15*(Ph(t)- RT(t) - G(t)));
    end
    g_P25(g)= P_equal(25);
end
format longE
for a = 1:length(gs);
    if abs(g_P25(a)- 3.4) <=0.001
        g_P25(a)
        gs(a)  
    end
end
%To be periodic, P(0) must equal P(25) - this occurs when g is .0074348

%vary other parameters
r_a = 0.8*r;
P_equal_a = zeros(24,1);
P_equal_a(1)= 3.4; %Po
g_P25_a = zeros(7001,1);
gs_a =0.003:.000001:0.01;
for g = 1:length(gs_a);
    gs_a(g);
    for t=1:24
        RT(t) = Ro*exp(r_a*T(t));
        G(t)= gs_a(g)*Z(t);
        Ph(t) = ((p*I(t))/(k(t)*z1(t)))*(1-exp(-k(t)*z1(t)))*OneMinusN(t)*OneMinusV(t);
        P_equal_a(t+1)=P_equal_a(t)*exp(15*(Ph(t)- RT(t) - G(t)));
    end
    g_P25_a(g)= P_equal_a(25);
end
format longE
for a = 1:length(gs_a);
    if abs(g_P25_a(a)- 3.4) <=0.01
        g_P25_a(a)
        g_new = gs_a(a)
        g_new
    end
end



format longE

%3 10 year time series 
%constants
Ro = 0.0175;
r = 0.069;
g = 0.0074348;
p = 2.5;

T_years =repmat(T,10,1);
Z_years =repmat(Z,10,1);
I_years =repmat(I,10,1);
k_years =repmat(k,10,1);
z1_years =repmat(z1,10,1);
N_years =repmat(OneMinusN,10,1);
V_years =repmat(OneMinusV,10,1);


P_years = zeros(240,1);
P_years(1)= 3.4; %Po
    for t=1:240
        RT(t) = Ro*exp(r*T_years(t));
        Z_rand(t) = Z_years(t)*(0.8 + 0.4*rand());
        G(t)= g*Z_rand(t);
        Ph(t) = ((p*I_years(t))/(k_years(t)*z1_years(t)))*(1-exp(-k_years(t)*z1_years(t)))*N_years(t)*V_years(t);
        P_years(t+1)=P_years(t)*exp(15*(Ph(t)- RT(t) - G(t)));
    end

figure(8)
time = [1:15:3601];
plot(time,P_years)
title('10 year time series with Z varying randomly')
xlabel('Days')


