clear all
% set up global variables for use in the DE routine.
global q un uc Kn Kc Vn Vc r
global yv nv uv cv
global icount

t0  = 0;     %initial time
t1  = 40;    % final time
c0  = 0;     % initial carbon conc
n0  = 0;     % initial nitrogen conc
y0  = 4.9;   % initial yeast concentration

% initial state for use in DE solver.

S0  = [c0;n0;y0];

q  = 0.15; % dilution rate

%
% terSchure data from Figure 1 A.
%
uv = [29 44	61	66	78	90	96	114	118];
yv = [4.9	7	8.1	8.1	7.5	8.1	7.9	8.1	8.2];
nv = [0	0	0	8	21	30	40	55	60];

% terSchure paper observation on carbon.
cv = ones(size(uv))*0.1;

%
% Initial guesses of possible parameter values
%
r   = 0.14;
Kn  = 10
Kc  = 10;
Vn  = 2;
uc  = 100;
Vc  = 2;

x0 = [Vn;Vc;Kn;Kc;r];

icount = 0;

x1 = fmincon('steady_state_fit_fn',x0,[],[],[],[],[0.1;0.1;0.1;0.1;q],[1e2;1e2;1e2;1e2;1e2])

Vn = x1(1);
Vc = x1(2);
Kn = x1(3);
Kc = x1(4);
r  = x1(5);

% % % V = 8;
% % % Kn = 5;
% % % Kc = 10;
% % % V2 = 20;

sym = '+o*.xsd^v><ph+o*.xsd^v><ph+o*.xsd^v><ph+o*.xsd^v><ph';

nm = zeros(size(nv));
cm = zeros(size(nv));
ym = zeros(size(nv));

figure(1)

for ii = 1:length(uv)

    un = uv(ii);
    
    [t,St] = ode45('chemostat_dynamics',[t0,t1],S0);
    
    subplot(311)
    plot(t,St(:,1),t,cv(ii)*ones(size(t))),title('carbon')
    
    subplot(312)
    plot(t,St(:,2),t,nv(ii)*ones(size(t))),title('nitrogen')
    
    subplot(313)
    plot(t,St(:,3),t,yv(ii)*ones(size(t))),title('biomass')
    
    pause
    
    cm(ii) = St(end,1);
    nm(ii) = St(end,2);
    ym(ii) = St(end,3);
    
end
figure(2)
plot(uv,nv,uv,nm),title('Nitrogen residual')
figure(3)
plot(uv,cv,uv,cm),title('Carbon residual')
figure(4)
plot(uv,yv,uv,ym),title('Biomass')