% This file calls the built-in ode45, 
% to solve the fermentation problem in section 4.3.2;
%
init = [70 220 40 8];
[t, variable] = ode45('model432', 0, 150, init(:),1e-9,1);

G = variable(:,1);
M = variable(:,2);
N = variable(:,3);
T = variable(:,4);

% calculate the control action
Kp = 1*10^6;            % J/(hr.C^2.m^3)
Tset = 12;              % C
umin = 0;               % J/(hr.C.m^3)
umax = 40000;

u = Kp*(T-Tset);
for I = 1:length(u),
  if u(I) < umin
    u(I) = umin;
  end
  if u(I) > umax
    u(I) = umax;
  end
end

% calculate Ethanol and Biomass
G0 = 70;
M0 = 220;
N0 = 40;
X0 = 175;
REG = 1.92;     RXG = 0.134;  
REM = 3.84;     RXM = 0.268; 
REN = 5.76;     RXN = 0.402;

X = X0 + RXG*(G0-G) + RXM*(M0-M) + RXN*(N0-N);
E = REG*(G0-G) + REM*(M0-M) + REN*(N0-N);

plot(t, G, '-', t, M, '--', t, N, ':');
axis([0 150 0 250]);
title('Sugar Transient Response (Example 4.3)');
xlabel('Time (hours)');
ylabel('Concentration (mol/m^3)');
text(70, 40, 'Maltotriose');
text(10, 60, 'Glucose');
text(80, 140, 'Maltose');
print -deps fig421.eps;
pause;
plot(t, T);
title('Temperature Response');
xlabel('Time (hours)');
ylabel('Temperature (degree C)');
print -deps fig422.eps;
pause;
plot(t, u)
title('Control Response');
xlabel('Time (hours)');
ylabel('Cool Control (J/(hr-C-m^3))');
print -deps fig423.eps;
pause;
plot(t, X, '-', t, E, '--');
title('Production Response');
xlabel('Time (hours)');
ylabel('Concentration (mol/m^3)');
text(80, 180, 'Biomass');
text(45, 500, 'Ethanol');
print -deps fig424.eps;
