% The function vent_vol.m inverts the ventricular pressure-volume functions 
% which are established by integrating a hyperbolic tangent function.  
% That is, the function determines the volume associated with a given 
% pressure using Newton's method which converges quadratically to the root.
%
% Function arguments:
%	xi - initial guess of normalized volume (must be between 0 and 1)
%	yv - normalized pressure for which the corresponding normalized
%	     volume is desired (must be between 0 and 1)
%	Ev - the current ventricular elastance OR differential
%            ventricular elastance at unstressed volume (mmHg/ml)
%
% Function output:
%	x - desired normalized volume
%

function x = vent_vol(xi,yv,Ev);

% Assigning variables.
TOL = 10e-6;
k = 50;

% Assigning initial guess.
x = xi;

% Calculating normalized ventricular pressure-volume function
% and its derivative at the initial guess.
f = Ev*x + (1-Ev)/(1+(1/k)*(log((1+exp(-k*(1-(1/(Ev+1)))))/(1+exp(k*(1/(Ev+1)))))))*(x+(1/k)*(log((1+exp(-k*(x-(1/(Ev+1)))))/(1+exp(k*(1/(Ev+1))))))) - yv;

df = Ev + (1-Ev)/(1+(1/k)*(log((1+exp(-k*(1-(1/(Ev+1)))))/(1+exp(k*(1/(Ev+1)))))))*(1-(exp(-k*(x-(1/(Ev+1))))/(1+exp(-k*(x-(1/(Ev+1)))))));

% Calculating next normalized volume based on this function
% and its derivative via a single Newton's step.
xn = x - (f/df);

% Repeating Newton's iteration until convergence of normalized
% volume is achieved.
q = 0;
zz = ((abs(xn-x) > TOL));
mbintscalar(zz)
while (zz)

	x = xn;

	f = Ev*x + (1-Ev)/(1+(1/k)*(log((1+exp(-k*(1-(1/(Ev+1)))))/(1+exp(k*(1/(Ev+1)))))))*(x+(1/k)*(log((1+exp(-k*(x-(1/(Ev+1)))))/(1+exp(k*(1/(Ev+1))))))) - yv;

	df = Ev + (1-Ev)/(1+(1/k)*(log((1+exp(-k*(1-(1/(Ev+1)))))/(1+exp(k*(1/(Ev+1)))))))*(1-(exp(-k*(x-(1/(Ev+1))))/(1+exp(-k*(x-(1/(Ev+1)))))));

	xn = x - (f/df);
	
	zz = ((abs(xn-x) > TOL));
	mbintscalar(zz);

end

x = xn;

% Exiting program with error if normalized volume is 
% less than zero or greater than one.
zz = (x < -TOL | x > (1+TOL));
mbintscalar(zz);
if (zz)

	error('No root is found');

end