% The function var_vcap.m computes the ventricular elastance and 
% derivative of the elastance for a desired time step.  This
% function is valid for the intact circulation with contracting
% atria.
%
% Function arguments:
%	Cs - end-systolic OR differential end-systolic compliance
%	     value at the unstressed volume (ml/mmHg)
%	Cd - diastolic OR differential diastolic compliance value
%	     at the dead unstressed (ml/mmHg)
%	Tp - previous cardiac cycle length (s)
%	sumdeltaT -  time surpassed in current cardiac cycle (s)
%
% Function outputs:
%	E - elastance or differential elastance value at the unstressed volume
%	    for the desired time step (mmHg/ml)
%	dE - derivative of E at the desired time ste[ (mmHg/ml-s)
%

function [E,dE] = var_vcap(Cs,Cd,Tp,sumdeltaT)

% Making the ventricle a passive element, if Cd is less 
% than Cs.
if (Cd <= Cs)

	E = 1/Cd;
	dE = 0;

else

	% Assigning variables.
	Tinit = 0.125*sqrt(Tp);
	Ts = .3*sqrt(Tp);
	Tir = Ts/2;

	x = (sumdeltaT < Tinit);
	y = (sumdeltaT < Tinit+Ts); 
	z = (sumdeltaT < Tinit+Ts+Tir);
	mbintscalar(x);
	mbintscalar(y);
	mbintscalar(z);

	% Setting values during diastolic filling phase.
	if (x)

	        E = 1/Cd;
		dE = 0;

	% Calculating values during isovolumic contraction phase.
	elseif (y)

		E = (((1/Cs)-(1/Cd))/2) * (1-cos((pi/Ts)*(sumdeltaT-Tinit))) + (1/Cd);
		dE = (pi/Ts)*(((1/Cs)-(1/Cd))/2) * sin((pi/Ts)*(sumdeltaT-Tinit));


	% Calculating values during isovolumic relaxation phase.
	elseif (z)

		E = (((1/Cs)-(1/Cd))/2) * (1+cos((pi/Tir)*(sumdeltaT-Ts-Tinit))) + (1/Cd);
		dE = -(pi/Tir)*(((1/Cs)-(1/Cd))/2) * sin((pi/Tir)*(sumdeltaT-Ts-Tinit));

	% Setting values during diastolic filling phase.
	else

		E = 1/Cd;
		dE = 0;

	end

end