using PyPlot

## This model incorporates the Hai-Murphy crossbridge ODEs to the BG-AW-heterogeneity
## model being developed.

function coup(order,kappa)#,rinitial)
tic()

numBr = 2^order - 1;     # number of airways
numOrd1 = (numBr+1)/2
# initial values for airways [mm]
rimax = [0.296 0.318 0.337 0.358 0.384 0.414 0.445 0.484 0.539 0.608 ...
    0.692 0.793 0.913 1.052 1.203 1.374]
rinitial=[]
for k=1:order
    rinitial=[rinitial rimax[k]*ones(2^(order-k),1)]'
end
rinitial=rinitial+0.01*rand(1,length(rinitial))'; # peturbation of rinitial

M = 1; Mp = 0; AMp = 0; AM = 0

RR0= [M Mp AMp AM]'
r[:,1]=rinitial; #airways, order 1 : order n
uv[:,1] = zeros(numBr,1); # pressure in airway
vv[:,1] = zeros(numBr,1); # airflow within the airway
vds[:,1] = zeros((numBr+1)/2,1); # airflow entering alveoli
RR[:,1] = RR0
y0 = [r[:,1] uv[:,1] vv[:,1] vds[:,1] RR[:,1]]'; #ICs vector

#Stepwise kappa value
    tspan = range(0,100,length=1000)


    parm = kappa


    # Set numerical accuracy options for ODE solver
    options = odeset("RelTol', 1e-08, 'AbsTol', 1e-08, 'MaxStep", 0.1)

    [t,NN]=ode15s(@(t,y) seven_airway[t,y,order,parm], tspan, y0, options)
    U = zeros(length(NN),length(r)); # pressure including viscous damping
    UVDOT = zeros(length(NN),length(r))
    VVDOT = zeros(length(NN),length(r)); #airflow
    force = zeros(length(NN),length(r))
    phosphorylation = zeros(length(NN),length(r))
    pmid = zeros(length(NN),length(r))
    for j = 1:length(NN)
        [U[j,:],UVDOT[j,:],VVDOT[j,:],force[j,:],phosphorylation[j,:],pmid[j,:]] = seven_airway2[t[j],NN[j,:],order,parm]
    end
toc()
# CM=jet[length(rinitial)]
r_new = NN[:,1:numBr]
air = zeros(length(NN),length(r))
for i=1:length(rinitial)
    for j = 1:length(t)
     if r_new[j,i] <= 0.01
       air[j,i]=0
      else()
       air[j,i]=1
     end
    end

end
for i = 1:(numOrd1)
#     figure(1)
#     subplot(order,order,i)
#     plot(t,r_new[:,i])
#     title(sprintf("Radius of airway [#d]",i))
# #
# #     figure(2)
# #     subplot(order,10,i)
# # #     M[i] = findpeaks[r_new[:,i]]
# #     DataInv[:,i] = 1.01*max(r_new[:,i]) - r_new[:,i]
# #     [Minima_P(i},MinIIdx{i)] = findpeaks[DataInv[:,i]]
# #     abc = r_new[:,i];
# #     Minima_P[i] = abc[MinIIdx[i]];
# #     plot(t[MinIIdx(i}],Minima_P{i),"-o")
#
    figure(2)
    subplot(order,order,i)
    plot(t,VVDOT[:,i])
    title(sprintf("Flow in airway [#d]",i))
# # #
#     figure(3)
#     subplot(order,order,i)
#     plot(t,U[:,i])
#     title(sprintf("Pressure in airway [#d]",i))
#
#     figure(5)
#     subplot(order,order,i)
#     plot(t,pmid[:,i])
end
# #
# #
# legendStr = cell(1,length(rinitial))
# for i = 1:length(legendStr)
#     legendStr[i] = sprintf("r_[#d]",i)
#
# end
# legend(legendStr,"Location','BestOutside")
# ## Measure for state of heterogeniety
x = cell(1,numOrd1)
for i = 1:length(x)
    x[i] = air[:,i]
end
summ = zeros(1,length(x[1]))
for i = 1:length(x)
    switch i
        case 1
            neighbour_diff = sqrt((x(:,i) - x(:,end))^2 + (x(:,i) - x(:,i+1))^2)
        case length(x)
            neighbour_diff = sqrt((x(:,i) - x(:,i-1))^2 + (x(:,i) - x(:,1))^2)
        otherwise()
            neighbour_diff = sqrt((x(:,i) - x(:,i-1))^2 + (x(:,i) - x(:,i+1))^2)
    end
    summ = summ + neighbour_diff'


end
summ=summ/(2^(order-3)*4*sqrt(2))

y=cell(1,numOrd1)
 for k=1:length(y)
     y[k]=VVDOT[:,k]
 end
    flow_sum = zeros(1,length(y[1]))

    for ii = 1:length(y)
      switch ii
          case 1
              neighbour_diff1 = sqrt((y(:,ii) - y(:,end))^2 + (y(:,ii) - y(:,ii+1))^2)
          case length(y)
              neighbour_diff1 = sqrt((y(:,ii) - y(:,ii-1))^2 + (y(:,ii) - y(:,1))^2)
          otherwise()
              neighbour_diff1 = sqrt((y(:,ii) - y(:,ii-1))^2 + (y(:,ii) - y(:,ii+1))^2)
      end
     flow_sum = flow_sum + neighbour_diff1'


return [t,NN,air,U,UVDOT,VVDOT,force,phosphorylation,summ,flow_sum]
end
flow_sum=flow_sum./numOrd1


# ## graphing airways
node1 = zeros(1,numBr)
    node2 = zeros(1,numBr)

    # e.g. for an order 4: node1 = [1 2 2 3 3 4 4 5 5 6 6 7 7 8 8]
    node1[1] = 1
    for i = 1:(numBr-1)/2
        node1[2*i:2*i+1] = [i+1 i+1]
    end

    # e.g. for an order 4: node2 = [2 3 4 5 6 7 8 9 10 11 12 13 14 15 16]
    for i = 1:numBr
        node2[i] = i+1
    end

    ind = zeros(1,numBr)
    for i = 1:numBr
        ind[i] = i
    end
#     figure()
    A=reverse(air[end,:], dims = 2)

    figure(4)
    G = graph[node1,node2,A]
    h=plot(G,"LineWidth",2)
    h.EdgeCData=A
    labelnode[h,[node1 node2],'']
    labeledge[h,ind,G.Edges.Weight[ind]]
    colormap jet



# for i = 1:8
# subplot(4,2,i)
# plot(time,VVDOT[:,i])
# title(sprintf("Flow of airway [#d]",i))
# end

end