% A function that solves the 7 coupled airway BG model

function [NN] = seven_airway(t,y,order,parm)

r = [y(1) y(2) y(3) y(4) y(5) y(6) y(7)];
uv = [y(8) y(9) y(10) y(11) y(12) y(13) y(14)];
vv = [y(15) y(16) y(17) y(18) y(19) y(20) y(21)];
vds =[y(22) y(23) y(24) y(25)];
%% Construct the constant array as a persistent variable

% So that repeated calls to this function doesn't result in repeated
% calculations of C
persistent C   

if isempty(C)
    fprintf('Contrusting the structure array of constants...\n');
    C = createConstSymm(order);
    fprintf('Finished.\n');
end

%% Parameters for MoC model
rho=2;
%% Define the constants
f = 0.33;
A = 5e5; % Used to be 5e5;
Pref = 10+ 10*sin(2*pi*f*t);
% Rref = 0.2792;
Rref=0.4140;

ptop = Pref;

kappa = parm(1,:);
% epsilon = 10^(parm(2,:));  % A parameter that controls sliding between Type II and Type I + II
epsilon = 0;
qhat = -50;

N = 2^order - 1; % number of branches for symmetric trees
numOrd1 = (N+1)/2;
E1 = 0.25;


%% Setting up BG parameters
E = 25;
h1 = 0.9732;
rho1 = 1.225;
R = zeros(N,1);
Rv = zeros(N,1);
CC = zeros(N,1);
I = zeros(N,1);
pa_bar = zeros(N,1);
alph = zeros(N,1);
pbot = zeros(N,1);

for i=1:N
    air_ord = C.order(i);
    CC(i) = (2*pi*C.L(air_ord).*r(i).^3)/(E*h1);
    I(i) = rho1*C.L(air_ord)/(pi.*r(i).^2);
    R(i) = C.alpha(air_ord).*r(i).^(-4);
    Rv(i) = 0.01/CC(i);
    pa_bar(i) = (10*E1)./sqrt(E1.^2 + (2*pi*f*R(i)).^2);
    alph(i) = atan((2*pi*f*R(i))./E1);
    pbot(i) = 10+pa_bar(i)*sin(2*pi*f*t-alph(i));
end


%% Set up uv1 - uv7
uv(1) = (vv(1)-vds(1))/CC(1);
    u(1) = uv(1)+Rv(1)*(vv(1)-vds(1));
uv(2) = (vv(2)-vds(2))/CC(2);
    u(2) = uv(2)+Rv(2)*(vv(2)-vds(2));
uv(3) = (vv(3)-vds(3))/CC(3);
    u(3) = uv(3)+Rv(3)*(vv(3)-vds(3));
uv(4) = (vv(4)-vds(4))/CC(4);
    u(4) = uv(4)+Rv(4)*(vv(4)-vds(4));
uv(5) = (vv(5)-vv(2)-vv(1))/CC(5);
    u(5) = uv(5)+Rv(5)*(vv(5)-vv(2)-vv(1));
uv(6) = (vv(6)-vv(4)-vv(3))/CC(6);
    u(6) = uv(6)+Rv(6)*(vv(6)-vv(4)-vv(3));
uv(7) = (vv(7)-vv(6)-vv(5))/CC(5);
    u(7) = uv(7)+Rv(7)*(vv(7)-vv(6)-vv(5));
%% Set up v1 - v7
vv(1) = (u(5) - u(1) - (vv(1)*R(1))/2)/(I(1)/2);
vv(2) = (u(5) - u(2) - (vv(2)*R(2))/2)/(I(2)/2);
vv(3) = (u(6) - u(3) - (vv(3)*R(3))/2)/(I(3)/2);
vv(4) = (u(6) - u(4) - (vv(4)*R(4))/2)/(I(4)/2);
vv(5) = (u(7) - u(5) - (vv(5)*R(5)))/(I(5));
vv(6) = (u(7) - u(6) - (vv(6)*R(6)))/(I(6));
vv(7) = (ptop - u(7) - (vv(7)*R(7)))/(I(7));


    
%% Set up vds1 - vds4
vds(1) = (u(1)-pbot(1)-(R(1)/2)*vds(1))/(I(1)/2);
vds(2) = (u(2)-pbot(2)-(R(2)/2)*vds(2))/(I(2)/2);
vds(3) = (u(3)-pbot(3)-(R(3)/2)*vds(3))/(I(3)/2);
vds(4) = (u(4)-pbot(4)-(R(4)/2)*vds(4))/(I(4)/2);
%% Determine what R is (airway radius)

RDOT=zeros(N,1);

%% Find expressions for the pressure and flows

Ci = C.Ciplus - C.Ciminus1 - C.Ciminus2;

Dalphar = diag(C.alpha.^(-1).*r'.^4);

W = Ci*Dalphar*(C.Cjplus - C.Cjminus);

lambda = dot(C.alpha.^(-1).*r'.^4, C.vbot);

temp = (W\(Ci*Dalphar*(C.vbot*C.vbot')*Dalphar*C.Cjplus))/lambda; %optimising code
Lambda = eye(size(temp)) - temp;
p = Lambda\(W\(Ci*Dalphar*(-ptop*C.vtop - 1/lambda*qhat*C.vbot))); %optimising code

gammajplus = C.Cjplus*p + ptop.*C.vtop;
gammajminus = C.Cjminus*p + pbot.*C.vbot;
deltap = gammajplus - gammajminus;

% 
T = C.T;    % The transformation matrix that will map mu_hat to mu
% 
% %% Determine what mu is (parenchymal shear modulus)
% 
mu = zeros(N,1);
% 
for i = 1:numOrd1  % Wanting to loop through only the order 1 airways
        
    % Add the Type 1 coupling by putting in the dependence on the
    % neighbours. 
    if i == 1
        mu(i) = abs(A/3*(deltap(i)*r(i)^4 + epsilon*(deltap(numOrd1)*r(numOrd1)^4 + deltap(i+1)*r(i+1)^4)))/2;
    elseif i == numOrd1
        mu(i) = abs(A/3*(deltap(i)*r(i)^4 + epsilon*(deltap(i-1)*r(i-1)^4 + deltap(1)*r(1)^4)))/2;
    else
        mu(i) = abs(A/3*(deltap(i)*r(i)^4 + epsilon*(deltap(i-1)*r(i-1)^4 + deltap(i+1)*r(i+1)^4)))/2;
    end

end
% 
% % Use the transformation matrix, T, to create the full mu vector
mu = T*mu;


% 
% %% Determine what tau is (parenchymal tethering pressure)
tau = 2*mu.*(((Rref - r')/Rref) + 1.5*((Rref - r')/Rref).^2);


% 
% %% Determine what Ptm is (transmural pressure)
% 
% 
pmid = 0.5*(gammajplus + gammajminus);
Ptm = zeros(N,1);
for i=1:N
    air_ord = C.order(i);  % Determine the order of the current airway

Ptm(i) = pmid(i) - (kappa*Rref/r(i))+ tau(i);


    if Ptm(i) <= 0
        R(i) = sqrt((C.Ri(air_ord).^2)*(1 - Ptm(i)./C.P1(air_ord)).^(-C.n1(air_ord)));
    elseif Ptm(i) >= 0
        R(i) = sqrt(C.rimax(air_ord).^2 - (C.rimax(air_ord).^2 - C.Ri(air_ord).^2)*(1 - Ptm(i)/C.P2(air_ord)).^(-C.n2(air_ord)));
    end
    
   RDOT=rho*((R)-r');
end

RRDOT=RDOT;
% size(RRDOT)
UVDOT = uv';
% size(UVDOT)
VVDOT = vv';
% size(VVDOT)
VDSDOT = vds';
% size(VDSDOT)
NN = [RRDOT;UVDOT;VVDOT;VDSDOT];