% The function eval_deriv.m evaluates the right-hand side of a set
% of first-order differential equations which characterize the 
% intact circulation with only linear elements at a desired time
% step.
%
% Function arguments:
%	P - a 6x1 vector containing the current pressure values
%	  = [Pl; Pa; Pv; Pr; Ppa; Ppv]
%	Qc - a 2x1 vector containing the current ventricular volume values
%	   - [Ql; Qr]
%	Pbreathe - 4x1 vector containing current respiratory values
%		   [Qlu; Pth dPth Ppac (Palv)]
%	th - current parameter values
%	sumdeltaT - time surpassed in current cardiac cycle
%
% Function outputs:
%	dP - a 6x1 vector containing the derivative of current
%	     pressure values
%	   = [dPl; dPa; dPv; dPr; dPpa; dPpv]
%

function dP = eval_deriv(P,Qc,Pbreathe,th,sumdeltaT);

% Computing flow rates from the current pressure values.
if (P(6) >= P(1))
	qli = (P(6)-P(1))/th(20);
else
	qli = 0;
end

if (P(1) >= P(2))
	qlo = (P(1)-P(2))/th(15);
else
	qlo = 0;
end

qa = (P(2)-P(3))/th(16);

if (P(3) >= P(4))
	qri = (P(3)-P(4))/th(17);
else
	qri = 0;
end

if (P(4) >= P(5))
	qro = (P(4)-P(5))/th(18);
else
	qro = 0;
end

qpa = (P(5)-P(6))/th(19);

% Computing ventricular elastances.
[El,dEl] = var_cap(th(1),th(2),th(24),sumdeltaT);
[Er,dEr] = var_cap(th(5),th(6),th(24),sumdeltaT);

% Computing dP from the above calculations and th parameter vector.
dP = [(El*(qli-qlo))+((dEl/El)*(P(1)-Pbreathe(2)))+Pbreathe(3);
      ((qlo-qa)/th(3))+(Pbreathe(3)/3);
      (qa-qri)/th(4);
      (Er*(qri-qro))+((dEr/Er)*(P(4)-Pbreathe(2)))+Pbreathe(3);
      ((qro-qpa)/th(7))+Pbreathe(3);
      ((qpa-qli)/th(8))+Pbreathe(3);
      0;
      0;];