%
% 生体システム工学の基礎
% 図4.70 tolic.mから呼び出される関数
%
function [Ip Ii G t] = gi4()
clear all

global Vp Vi Vg E tp ti td Rm a1 C1 Ub C2 C3 U0 Um beta C4 Rg alpha C5

    function a = f1(G)
        a = Rm / (1 + exp((C1 - G/Vg)/a1));
    end

    function a = f2(G)
        a = Ub * (1 - exp(-G / (C2 * Vg)));
    end

    function a = f3(G)
        a = G / (C3 * Vg);
    end

    function a = f4(Ii)
        a = U0 + (Um - U0) / (1 + exp(-beta * log(Ii/C4*(1/Vi + 1/(E * ti)))));
    end

    function a = f5(x3)
        a = Rg / (1 + exp(alpha * (x3/Vp - C5)));
    end

    function a = Gin(t)
        a = 216;
    end

Vp = 3;
Vi = 11;
Vg = 10;
E = 0.2;
tp= 6;
ti = 100;
td = 36;
Rm = 210;
a1 = 300;
C1 = 2000;
Ub = 72;
C2 = 144;
C3 = 1000;
U0 = 40;
Um = 940;
beta = 1.77;
C4 = 80;
Rg = 180;
alpha = 0.29;
C5 = 26;

Ip(1) = 30 * Vp ; 
Ii(1) = 20 * Vi;
G(1) = 12000; % dLで10倍，体積10Lで10倍含まれるので100倍した値
x1(1) = 0;
x2(1) = 0;
x3(1) = 0;

dt = 0.1;
step = 10000;
t = 0:step-1;
t = dt * t;
for i = 2:step
    Ip(i) = Ip(i-1) + dt * (f1(G(i-1)) - E * (Ip(i-1)/Vp - Ii(i-1)/Vi) - Ip(i-1)/tp);
    Ii(i) = Ii(i-1) + dt * (E * (Ip(i-1)/Vp - Ii(i-1)/Vi) - Ii(i-1)/ti);
    G(i) = G(i-1) + dt * (Gin(i) - f2(G(i-1)) - f3(G(i-1))*f4(Ii(i-1)) + f5(x3(i-1)));
    x1(i) = x1(i-1) + dt * (3 / td * (Ip(i-1) - x1(i-1)));
    x2(i) = x2(i-1) + dt * (3 / td * (x1(i-1) - x2(i-1)));
    x3(i) = x3(i-1) + dt * (3 / td * (x2(i-1) - x3(i-1)));
end

G = G/(10 * Vg); %mg/dLになる
Ip = Ip/(Vp);
end


