% The function a_eval_deriv.m evaluates the right-hand side of a set
% of first-order differential equations which characterize the 
% intact circulation with contracting atria at a desired time step.
%
% Function arguments:
%	P - a 8x1 vector containing the current pressure values
%	  = [Pl; Pa; Pv; Pr; Ppa; Ppv; Pra; Pla]
%	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 8x1 vector containing the derivative of current
%	     pressure values
%	   = [dPl; dPa; dPv; dPr; dPpa; dPpv; dPra; dPla]
%

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

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

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

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

if (P(7) >= P(4))
	 qra = (P(7)-P(4))/th(87);
else	 
	 qra = 0;
end

if (P(3) >= 0 & P(7) >= th(91))
	qv = (P(3)-P(7))/th(17);
elseif (P(3) > th(91) & th(91) > P(7))
	qv = (P(3)-th(91))/th(17);
elseif (P(7) > th(91) & th(91) > P(3))
       qv = (th(91)-P(7))/th(17);
else
	qv = 0;
end

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

h1 = (1330/(1.06*980))*(P(6)-Pbreathe(4));
if (h1 < -10)
	h1 = -10;
elseif (h1 > 20)
	h1 = 20;
end

h2 = (1330/(1.06*980))*(P(5)-Pbreathe(4));
if (h2 < -10)
	h2 = -10;
elseif (h2 > 20)
	h2 = 20;
end

qpa = (((P(5)-P(6))/(30*th(19)))*(h1+10)) + (((P(5)-Pbreathe(4))/(30*th(19)))*(h2-h1)) - (((1.06*980)/(2660*(30*th(19))))*(h2^2-h1^2));

qpv = (P(6)-P(8))/th(20);

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

[Ela,dEla] = var_acap(th(80),th(81),th(24),sumdeltaT);
[Era,dEra] = var_acap(th(84),th(85),th(24),sumdeltaT);

% Computing dP from the above calculations and th parameter vector.
if (P(1) < Pbreathe(2))
	g = ((P(1)-Pbreathe(2)))/(El*th(9)*(2/pi));
	dPl = El*(qla-qlv)*(1+g^2) + (P(1)-Pbreathe(2))*(dEl/El) + Pbreathe(3);
else

	y = (P(1)-Pbreathe(2))/th(31);

	x = vent_vol(Qc(1),y,El);

	x0 = 1/(El+1);
	k = 50;

	mbrealscalar(-k*(x-x0));
	mbrealscalar(-k*(1-x0));
	mbrealscalar(k*x0);

	dx = (qla-qlv)/(th(26)-th(9));
	dx0 = -dEl/((El+1)^2);

	c = 1+exp(-k*(1-x0));
	d = 1+exp(k*x0);
	g = 1+exp(-k*(x-x0));

	dc = k*dx0*exp(-k*(1-x0));
	dd = k*dx0*exp(k*x0);
	dg = -k*(dx-dx0)*exp(-k*(x-x0));

	denb = (1+(1/k)*log(c/d));
	b = (1-El)/denb;
	db = (-dEl*denb - (1-El)*(1/k)*((d*dc-c*dd)/(c*d)))/(denb^2);
	dy = dEl*x + El*dx + db*(x+(1/k)*(log(g/d))) + b*(dx+(1/k)*((d*dg-g*dd)/(g*d)));

	dPl = th(31)*dy + Pbreathe(3);

end

if (P(4) < Pbreathe(2))
	g = ((P(4)-Pbreathe(2)))/(Er*th(12)*(2/pi));
	dPr = Er*(qra-qrv)*(1+g^2) + (P(4)-Pbreathe(2))*(dEr/Er) + Pbreathe(3);
else

	y = (P(4)-Pbreathe(2))/th(32);

	x = vent_vol(Qc(2),y,Er);

	x0 = 1/(Er+1);
	k = 50;

	mbrealscalar(-k*(x-x0));
	mbrealscalar(-k*(1-x0));
	mbrealscalar(k*x0);

	dx = (qra-qrv)/(th(27)-th(12));
	dx0 = -dEr/((Er+1)^2);

	c = 1+exp(-k*(1-x0));
	d = 1+exp(k*x0);
	g = 1+exp(-k*(x-x0));

	dc = k*dx0*exp(-k*(1-x0));
	dd = k*dx0*exp(k*x0);
	dg = -k*(dx-dx0)*exp(-k*(x-x0));

	denb = (1+(1/k)*log(c/d));
	b = (1-Er)/denb;
	db = (-dEr*denb - (1-Er)*(1/k)*((d*dc-c*dd)/(c*d)))/(denb^2);
	dy = dEr*x + Er*dx + db*(x+(1/k)*(log(g/d))) + b*(dx+(1/k)*((d*dg-g*dd)/(g*d)));

	dPr = th(32)*dy + Pbreathe(3);

end

dPla = (Ela*(qpv-qla))+((dEla/Ela)*(P(8)-Pbreathe(2)))+Pbreathe(3);
dPra = (Era*(qv-qra))+((dEra/Era)*(P(7)-Pbreathe(2)))+Pbreathe(3);

dP = [dPl; ((qlv-qa)/th(3)); (qa-qv)/th(4); dPr; ((qrv-qpa)/th(7))+Pbreathe(3); ((qpa-qpv)/th(8))+Pbreathe(3); dPra; dPla];