% 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