# A function that solves the n coupled airway BG model

function seven_airway2(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 [U,UVDOT,VVDOT,force,phosphorylation,pmid]
end-4)]'
# RR = y[end-6:end]'
RR = y[end-3:end]'
M = RR[1]; Mp = RR[2]; AMp = RR[3]; AM = RR[4]

## 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("Constructing 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

ptop = 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
E1 = 2.5; #cmH20/ml

## Setting up BG parameters
E = 25
h1 = 0.9732
k11 = 0.0057
k22 = 0.2096
k33 = 0.0904
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)
u=zeros(N,1)
pmid = 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]);#
    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] = (P0*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



##
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]-pbot[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]
end



        ## Set up uv1 - uv[n]

        ## Set up v1 - v[n]

        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)
# return
# Dalphar = diag(R)
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)


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)
# stiffness =rr[1]+rr[4]
phosphorylation = AMp+Mp
#
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

#
#    drdt[i]=rho*(R[i]-(r[i]));
#    V[i]=-((gamma)/(2*pi*Rref))*drdt[i];
   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]'

U=u'
end