% The function rk4.m integrates a set of first-order differential equations
% characterizing the various pulsatile heart and circulation preparations
% usinig the fourth-order Runge-Kutta method.
%
% Function arguments:
% Pc - column vector containing the current pressure values
% Qc - a 2x1 vector containing the current ventricular volume values
% - [Ql; Qr]
% Pbreathe - 4x3 vector containing respiratory values [Qlu; Pth dPth Ppac (Palv)]
% at the current time step, at the middle of the current and next
% integration step, and at the next integration step
% th - current parameter values
% ipfm - fraction of time at in current cardiac cycle
% sumdeltaT - time surpassed in current cardiac cycle
% deltaT - numerical integration step size
% preparation - status of pulsatile heart and circulation
%
% Function outputs:
% Pn - column vector containing the pressure values at the next step
%
function Pn = rk4(Pc,Qc,Pbreathe,th,ipfm,sumdeltaT,deltaT,preparation)
% Calculating at current time step.
if (preparation == 0 | preparation == 4)
k1 = deltaT*intact_eval_deriv(Pc,Qc,Pbreathe(:,1),th,sumdeltaT);
elseif (preparation == 1)
k1 = deltaT*hlu_eval_deriv(Pc,Qc,Pbreathe(:,1),th,sumdeltaT);
elseif (preparation == 2)
k1 = deltaT*sc_eval_deriv(Pc,Qc,Pbreathe(:,1),th,sumdeltaT);
elseif (preparation == 3)
k1 = deltaT*lv_eval_deriv(Pc,Qc,Pbreathe(:,1),th,sumdeltaT);
elseif (preparation == 5)
k1 = deltaT*eval_deriv(Pc,Qc,Pbreathe(:,1),th,sumdeltaT);
elseif (preparation == 6)
k1 = deltaT*third_eval_deriv(Pc,Qc,Pbreathe(:,1),th,sumdeltaT);
elseif (preparation == 7)
k1 = deltaT*nac_eval_deriv(Pc,Qc,Pbreathe(:,1),th,sumdeltaT);
elseif (preparation == 8)
k1 = deltaT*third_nac_eval_deriv(Pc,Qc,Pbreathe(:,1),th,sumdeltaT);
elseif (preparation == 9)
k1 = deltaT*apr_eval_deriv(Pc,Qc,Pbreathe(:,1),th,sumdeltaT);
elseif (preparation == 10)
k1 = deltaT*a_eval_deriv(Pc,Qc,Pbreathe(:,1),th,sumdeltaT);
end
% Calculating in the middle of the current and next time step.
sumdeltaT1 = sumdeltaT+(deltaT/2);
% Determining if the ipfm model has fired.
nipfm1 = ipfm+(deltaT/2)*th(22);
mbrealscalar(nipfm1 >= 1);
if (nipfm1 >= 1)
sumdeltaT1 = (nipfm1 - 1)/th(22);
nipfm1 = nipfm1 - 1;
end
if (preparation == 0 | preparation == 4)
k2 = deltaT*intact_eval_deriv((Pc+(k1/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
k3 = deltaT*intact_eval_deriv((Pc+(k2/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
elseif (preparation == 1)
k2 = deltaT*hlu_eval_deriv((Pc+(k1/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
k3 = deltaT*hlu_eval_deriv((Pc+(k2/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
elseif (preparation == 2)
k2 = deltaT*sc_eval_deriv((Pc+(k1/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
k3 = deltaT*sc_eval_deriv((Pc+(k2/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
elseif (preparation == 3)
k2 = deltaT*lv_eval_deriv((Pc+(k1/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
k3 = deltaT*lv_eval_deriv((Pc+(k2/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
elseif (preparation == 5)
k2 = deltaT*eval_deriv((Pc+(k1/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
k3 = deltaT*eval_deriv((Pc+(k2/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
elseif (preparation == 6)
k2 = deltaT*third_eval_deriv((Pc+(k1/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
k3 = deltaT*third_eval_deriv((Pc+(k2/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
elseif (preparation == 7)
k2 = deltaT*nac_eval_deriv((Pc+(k1/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
k3 = deltaT*nac_eval_deriv((Pc+(k2/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
elseif (preparation == 8)
k2 = deltaT*third_nac_eval_deriv((Pc+(k1/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
k3 = deltaT*third_nac_eval_deriv((Pc+(k2/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
elseif (preparation == 9)
k2 = deltaT*apr_eval_deriv((Pc+(k1/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
k3 = deltaT*apr_eval_deriv((Pc+(k2/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
elseif (preparation == 10)
k2 = deltaT*a_eval_deriv((Pc+(k1/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
k3 = deltaT*a_eval_deriv((Pc+(k2/2)),Qc,Pbreathe(:,2),th,sumdeltaT1);
end
% Calculating at the next time step.
sumdeltaT2 = sumdeltaT+deltaT;
% Determining if the ipfm model has fired.
nipfm2 = ipfm + deltaT*th(22);
mbrealscalar(nipfm2 >= 1);
if (nipfm2 >= 1)
sumdeltaT2 = (nipfm2 - 1)/th(22);
nipfm2 = nipfm2 - 1;
end
if (preparation == 0 | preparation == 4)
k4 = deltaT*intact_eval_deriv((Pc+k3),Qc,Pbreathe(:,3),th,sumdeltaT2);
elseif (preparation == 1)
k4 = deltaT*hlu_eval_deriv((Pc+k3),Qc,Pbreathe(:,3),th,sumdeltaT2);
elseif (preparation == 2)
k4 = deltaT*sc_eval_deriv((Pc+k3),Qc,Pbreathe(:,3),th,sumdeltaT2);
elseif (preparation == 3)
k4 = deltaT*lv_eval_deriv((Pc+k3),Qc,Pbreathe(:,3),th,sumdeltaT2);
elseif (preparation == 5)
k4 = deltaT*eval_deriv((Pc+k3),Qc,Pbreathe(:,3),th,sumdeltaT2);
elseif (preparation == 6)
k4 = deltaT*third_eval_deriv((Pc+k3),Qc,Pbreathe(:,3),th,sumdeltaT2);
elseif (preparation == 7)
k4 = deltaT*nac_eval_deriv((Pc+k3),Qc,Pbreathe(:,3),th,sumdeltaT2);
elseif (preparation == 8)
k4 = deltaT*third_nac_eval_deriv((Pc+k3),Qc,Pbreathe(:,3),th,sumdeltaT2);
elseif (preparation == 9)
k4 = deltaT*apr_eval_deriv((Pc+k3),Qc,Pbreathe(:,3),th,sumdeltaT2);
elseif (preparation == 10)
k4 = deltaT*a_eval_deriv((Pc+k3),Qc,Pbreathe(:,3),th,sumdeltaT2);
end
Pn = Pc + k1/6 + k2/3 + k3/3 + k4/6;