function BatchBR_CellCult
% 厭氧間歇發酵動態模擬
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/06/01 $

clear all
clc
global umax KS pmax m YX2S YX2Shat MS

% 已知參數
umax = 0.5;
pmax = 0.109;
KS = 0.001;
m = 1;
YX2Shat = 2;
YX2S = 1;
MS = 0.008;

% 已知初始濃度
X0 = 0.0005;
S0 = 0.04;
p0 = 0

% 求解ODE-IVPs
tspan = [0  15];
y0 = [X0  S0  p0];
[t y] = ode45(@ModelEqs, tspan, y0)

% 繪圖
plot(t,y)
xlabel('反應時間，h')
ylabel('無量綱濃度，w/w')
axis([0,t(end),0,max(max(y))+0.005])
str = {['活細胞'],['葡萄糖'],['乙醇']};
for i = 1:3
    gtext(str{i})
end


% ------------------------------------------------------------------
function dydt = ModelEqs(t,y,k)     % 模型方程
global umax KS pmax m YX2S YX2Shat MS
X = y(1);
S = y(2);
p = y(3);
RX = umax*X * S/(KS+S) * (1-p/pmax)^m;
RS = -RX/YX2S - MS*X;
Rp = RX/YX2S - RX/YX2Shat + MS*X;
dydt = [RX; RS; Rp];