# A function that solves the n coupled airway BG model

function seven_airway(t,y,order,parm)

numBr = 2^order - 1;
r = y[1:numBr]'
uv = y[numBr+1: 2*numBr]'
vv = y[2*numBr+1:3*numBr]'
vds = y[3*numBr+1:(
return [NN]
end-4)]'
# RR = y[end-6:end]'
RR = y[end-3:end]'

M = RR[1]; Mp = RR[2]; AMp = RR[3]; AM = RR[4]
# r1 = [y[1] y[2] y[3] y[4] y[5] y[6] y[7]]
# uv1 = [y[8] y[9] y[10] y[11] y[12] y[13] y[14]]
# vv1 = [y[15] y[16] y[17] y[18] y[19] y[20] y[21]]
# vds1 =[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)
    @sprintf("Contrusting the structure array of constants...\n")
    C = createConstSymm[order]
    @sprintf("Finished.\n")
end

## Parameters for MoC model
rho=2
## Define the constants
f = 0.2
A = 5e5; # Used to be 5e5
# Pref = 10+ 10*sin(2*pi*f*t).*(t<=60)+ 120*sin(2*pi*f*t).*(t>60 && t<62)+ 10*sin(2*pi*f*t).*(t>=62)
# Rref = 0.2792
P0 = 10

Pref = 10 + P0*sin(2*pi*f*t)
Rref=0.4140

u_top = Pref


# kappa = 0.25*parm[1,:].*(t<=60) + 0.5*parm[1,:].*(t>60 && t<=120) + 0.75*parm[1,:].*(t>120 && t<=180) + (parm[1,:]).*(t>=180);
kappa = parm[1,:]
# epsilon = 10^(parm[2,:]);  # A parameter that controls sliding between Type II & Type I + II
epsilon = 0
qhat = -50

N = 2^order - 1; # number of branches for symmetric trees
numOrd1 = (N+1)/2
NT = (N-1)/2
E1 = 2.5; #cmH20/ml


## Setting up BG parameters
E = 25;#kPa
rho1 = 1.225
k11 = 0.0057
k22 = 0.2096
k33 = 0.0904
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)
u_bot = zeros(N,1)
u=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)
    WT[i] = (k11.*(2*r[i])^2)/4 + (k22.*(2.*r[i]))/2 +k33;#wall thickness
    CC[i] = (C.L[air_ord]*(2*r[i])^3)/(4*E*WT[i]);# taken from Yoon et al 2020
    I[i] = rho1*C.L[air_ord]/(pi.*r[i]^2);#airway Inertance
    R[i] = C.alpha(air_ord).*r[i]^(-4); #airway resistance
    Rv[i] = 0.01/CC[i];#air viscoelastic effect
    pa_bar[i] = (P0*E1)./sqrt(E1^2 + (2*pi*f*R[i])^2)
    alph[i] = atan((2*pi*f*R[i])./E1)
    u_bot[i] = 10+pa_bar[i]*sin(2*pi*f*t-alph[i]);#acinar pressure
end



##
for a = numOrd1+1:N
    uv[a] = (vv[a]-vv[2*(a-numOrd1)]-vv[2*(a-numOrd1)-1])/CC[a]
    u[a] = uv[a] +Rv[a]*((2*(a-numOrd1))-vv[2*(a-numOrd1)-1])
end

for b = 1:numOrd1
    uv[b] = (vv[b]-vds[b])/CC[b]
    u[b] = uv[b] + Rv[b]*(vv[b]-vds[b])
    vds[b] = (uv[b]-u_bot[b]-(R[b]/2)*vds[b])/(I[b]/2)
end


for d = 1:N-1
    a[d] =((2*d) - 1)
    aa = a[a<N-1];
    b[d] = 2*d
    bb = b[b<=N-1]
    z = [aa bb]
end

# for i = 1: length(aa)
# #     if (aa[i]) <=numBr && bb[i] <=numBr
#       vv[aa[i]] = (uv[numOrd1+i] - u[aa[i]] - (R[aa[i]]*vv[aa[i]])/2)./(I[aa[i]]/2)
#       vv[bb[i]] = (uv[numOrd1+i] - u[bb[i]] - (R[bb[i]]*vv[bb[i]])/2)./(I[bb[i]]/2)
# #     elseif [aa[i]] > numBr && bb[i] .> numBr
# #         vv[aa[i]] = (uv[numOrd1+i] - u[aa[i]] - (R[aa[i]]*vv[aa[i]]))./(I[aa[i]])
# #         vv[bb[i]] = (uv[numOrd1+i] - u[bb[i]] - (R[bb[i]]*vv[bb[i]]))./(I[bb[i]]);
# #     end
# end
# vv[N] = (Pref - u[N] - R[N]*vv[N])/I[N]

        ## 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)

        for i = 1: (numOrd1/2)
          vv[aa[i]] = (u[numOrd1+i] - u[aa[i]] - (R[aa[i]]*vv[aa[i]])/2)./(I[aa[i]]/2)
          vv[bb[i]] = (u[numOrd1+i] - u[bb[i]] - (R[bb[i]]*vv[bb[i]])/2)./(I[bb[i]]/2)
        end

        for i = (numOrd1/2)+1 : length(aa)
        vv[aa[i]] = (u[numOrd1+i] - u[aa[i]] - (vv[aa[i]]*R[aa[i]]))/(I[aa[i]])
        vv[bb[i]] = (u[numOrd1+i] - u[bb[i]] - (vv[bb[i]]*R[bb[i]]))/(I[bb[i]])
        end
       # vv[7] = (ptop - u[7] - (vv[7]*R[7]))/(I[7])
        vv[N] = (Pref - u[N] - R[N]*vv[N])/I[N]


## 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 & flows

Ci = C.Ciplus - C.Ciminus1 - C.Ciminus2
# size(C.alpha^(-1))
# size(r')
Dalphar = diag(C.alpha^(-1).*r'^4)
# size(Dalphar)
# Dalphar1 = diag(R)
# return

W = Ci*Dalphar*(C.Cjplus - C.Cjminus)
# size(W)
# return
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*(-u_top*C.vtop - 1/lambda*qhat*C.vbot))); #optimising code

gammajplus = C.Cjplus*p + u_top.*C.vtop
gammajminus = C.Cjminus*p + u_bot.*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)




k2 = 0.5; k5 = 0.5; k3 = 0.4; k4 = 0.1; k7 = 0.01; k1 = 0.55.*(t<=5)+0.06.*(t>5); k6 =k1

# dt=0.5*t
# force=kappa*(rr[2]+rr[5])
force = kappa*(AMp+AM)
#
pmid = 0.5*(gammajplus + gammajminus)

# Ptm = zeros(N,1)
Ptm = pmid - (force*Rref./r')+ tau
for i=1:N
    air_ord = C.order[i];  # Determine the order of the current airway

# Ptm[i] = pmid[i] - (force*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

F = [-k1*M+k2*Mp+k7*AM k4*AMp+k1*M-(k2+k3)*Mp k3*Mp+k6*AM-(k5+k4)*AMp k5*AMp-(k6+k7)*AM]'
RRDOT=RDOT
# size(RRDOT)
UVDOT = uv'
# size(UVDOT)
VVDOT = vv'
# size(VVDOT)
VDSDOT = vds'
# size(VDSDOT)
NN = [RRDOT UVDOT VVDOT VDSDOT F]'

end