% The function a_init_cond.m computes the initial values of the 
% intact circulation with contracting atria.  The initial values
% are computed based on conservation laws of a linearized model.
%
% Function arguments:
%	th - vector containing the initial parameter values
%
% Function outputs:
%	P - a 8x1 vector containing the six initial pressures
%	  = [Pl(0); Pa(0); Pv(0); Pr(0); Ppa(0); Ppv(0); Pra(0); Pla(0)]
%	Q - a 8x1 vector containing the six initial volumes
%	  = [Ql(0); Qa(0); Qv(0); Qr(0); Qpa(0); Qpv(0); Qra(0); Qla(0)]
%	q - a 8x1 vector containing the six initial flow rates
%	  = [qpv(0); ql(0); qa(0); qv(0); qr(0); qpa(0); qv(0); qpv(0)]
%	ve - a 4x1 vector containing the initial variable elastance values
%	   = [El(0); Er(0); Ela(0); Era(0)];
%

function [P,Q,q,ve] = a_init_cond(th)

% Storing nonlinear ventricular compliance values for 
% subsequent initial ventricular volume calculation.
Cld = th(2);
Crd = th(6);

% Converting ventricular compliances to linear values which
% is necessary for estimating initial pressures.
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;
     temp;
     temp;
     th(21)-(th(14)+th(13)+th(12)+th(11)+th(10)+th(9)+th(82)+th(86))+th(2)*th(23)+th(6)*th(23)+th(7)*th(23)+th(8)*th(23)+th(81)*th(23)+th(85)*th(23)];

A = [-th(2), th(1), th(6), th(5), 0, 0, 0, 0, 0, 0;
     -th(2), (th(1)+Ts/th(15)), 0, 0, -Ts/th(15), 0, 0, 0, 0, 0;
     -th(2), th(1), 0, 0, 1/(th(22)*th(16)), -1/(th(22)*th(16)), 0, 0, 0, 0;
     -th(2), th(1), 0, 0, 0, 1/(th(22)*th(17)), -1/(th(22)*th(17)), 0, 0, 0;
     -th(2), th(1), 0, 0, 0, 0, Td/th(87), -Td/th(87), 0, 0;
     -th(2), th(1), 0, Ts/th(18), 0, 0, 0, -Ts/th(18), 0, 0;
     -th(2), th(1), 0, 0, 0, 0, 0, 1/(th(22)*th(19)), -1/(th(22)*th(19)), 0;
     -th(2), th(1), 0, 0, 0, 0, 0, 0, 1/(th(22)*th(20)), -1/(th(22)*th(20));
     (-th(2)-Td/th(83)), th(1), 0, 0, 0, 0, 0, 0, 0, Td/th(83); 
     th(2), 0, th(6), 0, th(3), th(4), th(85), th(7), th(8), th(81)];

x = A\b;

P = [x(1) x(5) x(6) x(3) x(8) x(9) x(7) x(10)]';

% Establishing initial flow rates.
q = zeros(8,1);

if (P(8) >= P(1))
	q(1) = (P(8)-P(1))/th(83);
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(7) > P(4))
	q(4) = (P(7)-P(4))/th(87);
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);

q(7) = (P(3)-P(7))/th(87);

q(8) = (P(8)-P(1))/th(83);

% Establishing initial variable elastance values.
ve = [1/th(2) 1/th(6) 1/th(81) 1/th(85)]';

% Establishing initial volumes.
Q = [th(2) th(3) th(4) th(6) th(7) th(8) th(85) th(81)]'.*(P-th(23)*[1 0 0 1 1 1 1 1]')+[th(9:14); th(86); th(82)];

% Integrating a Fermi function approach

if (P(1) < th(23))
	Q(1) = th(9)*(2/pi)*atan(((P(1)-th(23))/((2/pi)*th(9)*ve(1)))) + th(9);
else
	yl = (P(1)-th(23))/th(31);
	xl = vent_vol(0.5,yl,1/Cld);
	Q(1) = (th(26)-th(9))*xl+th(9);
end

if (P(4) < th(23))
	Q(4) = th(12)*(2/pi)*atan(((P(4)-th(23))/((2/pi)*th(12)*ve(2)))) + th(12);
else
	yr = (P(4)-th(23))/th(32);
	xr = vent_vol(0.5,yr,1/Crd);
	Q(4) = (th(27)-th(12))*xr+th(12);
end