% This model incorporates crossbridge PDEs (Mijailovich et al) to the BG-AW-heterogeneity 
% model being developed.
% The model produces heterogeneity for both static and dynamic pressure oscillations
% as well as the re-narrowing and ASM fluidization effects post-DI together with the
% the expected changes in airway compliance and resistance.

function [t,NN,air,U,UVDOT,VVDOT,summ,force,ptop,Rr,rad_change,Ptm]=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];
  M10= 0.6209;
  M11=  0.3106;
  M12= 0.2071;
  M20= 0.1749;
  M21= 0.0871;
  M22= 0.0579;
  C= 0.222222222222222;
RR0 = [M10;M11;M12;M20;M21;M22;C];
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

    parm = [kappa];
    tspan = linspace(0,80,2000);

    % Set numerical accuracy options for ODE solver
    options = odeset('RelTol', 1e-08, 'AbsTol', 1e-08, 'MaxStep', 0.01);
    
    [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)); %volumetric airflow
    for j = 1:length(NN)
        [U(j,:),UVDOT(j,:),VVDOT(j,:),force(j,:), ptop(j,:),Rr(j,:),Ptm(j,:)] = seven_airway2(t(j),NN(j,:),order,parm);
    end
toc
% CM=jet(length(rinitial));
r_new = NN(:,1:numBr);
r_end = r_new(end,:);
air = zeros(length(NN),length(r));

% for i=1:length(rinitial)
% 
%    if r_new(end,i) <= 0.01
%        air(i)=0;
%    else
%        air(i)=1;
%    end
% 
% end
rad_change = (r_end-rinitial')./rinitial';
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(4,2,i);
    plot(t,r_new(:,i))
    title(sprintf('Radius in airway {%d}',i))
    
    figure(2);
    subplot(4,2,i);
    plot(t,VVDOT(:,i))
    title(sprintf('Flow in airway {%d}',i))
    ylim([-0.5, 0.5])
% %     
    figure(3);
    subplot(4,2,i);
    plot(t,U(:,i))
    title(sprintf('Pressure of airway {%d}',i))
    
    figure(4);
    subplot(4,2,i);
    plot(t,Rr(:,i))
    title(sprintf('Resistance of airway {%d}',i))
    
%     figure(4);
%     subplot(4,2,i)
%     plot(r_new(:,i),force)
end

Lung_volume = sum(VVDOT(:,1:end),2);
Lung_resistance = sum(Rr(:,1:end),2);
% % %     [
% % %     
% legendStr = cell(1,length(rinitial));
% for i = 1:length(legendStr)
%     legendStr{i} = sprintf('r_{%d}',i);
% 
% end
% figure,
% plot(t,force/max(force))
% 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));
% 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));
% %% 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=[fliplr(air(end,:))];
%     
% %     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
