% The function init_cond.m computes the initial values of the
% intact circulation with only linear elements. The initial
% values are computed based on conservation laws.
%
% Function arguments:
% th - vector containing the initial parameter values
%
% Function outputs:
% P - a 6x1 vector containing the six initial pressures
% = [Pl(0); Pa(0); Pv(0); Pr(0); Ppa(0); Ppv(0)]
% Q - a 6x1 vector containing the six initial volumes
% = [Ql(0); Qa(0); Qv(0); Qr(0); Qpa(0); Qpv(0)]
% q - a 6x1 vector containing the six initial flow rates
% = [qpv(0); ql(0); qa(0); qv(0); qr(0); qpa(0)]
% ve - a 2x1 vector containing the initial variable elastance values
% = [El(0); Er(0)];
%
function [P,Q,q,ve] = init_cond(th)
% Converting ventricular compliances to linear values;
th(1) = th(1)*((th(26)-th(9))/th(31));
th(2) = th(2)*((th(26)-th(9))/th(31));
th(5) = th(5)*((th(27)-th(12))/th(32));
th(6) = th(6)*((th(27)-th(12))/th(32));
% Estimating initial pressures.
Ts = .3*sqrt(1/th(22));
Td = 1/th(22) - Ts;
temp = -th(2)*th(23) + th(1)*th(23);
b = [temp+th(6)*th(23)-th(5)*th(23);
temp;
temp;
temp;
temp;
temp;
temp;
% th(21)-(th(14)+th(13)+th(12)+th(11)+th(10)+th(9))+th(2)*th(23)+(1/3)*th(3)*th(23)+th(6)*th(23)+th(7)*th(23)+th(8)*th(23)];
th(21)-(th(14)+th(13)+th(12)+th(11)+th(10)+th(9))+th(2)*th(23)+th(6)*th(23)+th(7)*th(23)+th(8)*th(23)];
A = [th(1), -th(2), 0, 0, -th(5), th(6), 0, 0;
th(1)+(Ts/th(15)), -th(2), -Ts/th(15), 0, 0, 0, 0, 0;
th(1), -th(2), 1/(th(22)*th(16)), -1/(th(22)*th(16)), 0, 0, 0, 0;
th(1), -th(2), 0, Td/th(17), 0, -Td/th(17), 0, 0;
th(1), -th(2), 0, 0, Ts/th(18), 0, -Ts/th(18), 0;
th(1), -th(2), 0, 0, 0, 0, 1/(th(22)*th(19)), -1/(th(22)*th(19));
th(1), -(th(2)+Td/th(20)), 0, 0, 0, 0, 0, Td/th(20);
0, th(2), th(3), th(4), 0, th(6), th(7), th(8)];
x = A\b;
P = [x(2:4)' x(6:8)']';
% Establishing initial flow rates.
q = zeros(6,1);
if (P(6) > P(1))
q(1) = (P(6)-P(1))/th(20);
else
q(1) = 0;
end
if (P(1) > P(2))
q(2) = (P(1)-P(2))/th(15);
else
q(2) = 0;
end
q(3) = (P(2)-P(3))/th(16);
if (P(3) > P(4))
q(4) = (P(3)-P(4))/th(17);
else
q(4) = 0;
end
if (P(4) > P(5))
q(5) = (P(4)-P(5))/th(18);
else
q(5) = 0;
end
q(6) = (P(5)-P(6))/th(19);
% Establishing initial variable ventricular elastance values.
ve = [1/th(2) 1/th(6)]';
% Establishing initial volumes.
Q = [th(2) th(3) th(4) th(6) th(7) th(8)]'.*(P(1:6)-th(23)*[1 0 0 1 1 1]') + [th(9:14)];