% runmodel.m

% global
global I_n f s I_p r q m k alpha;

% set default parameter values;
I_n = 0.2; % external N inputs (ug N/L/d)
I_p = 0.0008;  % external P inputs (ug P/L/d)
s=0.01; % sinking/loss rate for both P and N (per d)
r=0.15; % max rate of P recycling from sediments (ug P/L/d)
q=2; % exponent on terms of P recycling 
m=10; % half-saturation rate for P recycling from sediments (ug P/L)
f=0.88; % max N fixation rate (ug N/L per d)
k=30; % half-saturation rate for N fixation (N to P ratio);
alpha=0.2; % control on the slope of the logistic function

% simulation conditions
t0=0; % initial time
tfinal=1000; % final time (varied, but often run for 20000 iterations)

%  set up initial conditions - TN then TP (these varied considerably among runs)
initial=[40 2]'; 

% solve the equations for the request times and initial conditions
[t y]=ode45('derivsub_ratio', [t0 tfinal], initial);

% plot output to visualize the simulation
subplot(2,2,1);
plot(t,y(:,1),'b',t,y(:,2),'r');
title('Simulated Time Series'); ylabel('Numbers'); xlabel('Time'); legend('N','P');

subplot(2,2,2);
plot(t,y(:,1)./y(:,2),'k-');
ylabel('TN/TP'); xlabel('time');

subplot(2,2,3);
plot(y(:,1),y(:,2),'rx-');
title('Phase Plane'); ylabel('P (ug/L)'); xlabel('N (ug/L)');
