% 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;