Location: Peripheral airways matlab/CellML @ e08b9ac589b9 / Peripheral_matlab / generalised_model / coup.m

Author:

aram148 <a.rampadarath@auckland.ac.nz>

Date:

2021-12-20 14:36:36+13:00

Desc:

Generalised model

Permanent Source URI:

https://staging.physiomeproject.org/workspace/7e5/rawfile/e08b9ac589b9fffc2df026eb7ab23df1e5a75255/Peripheral_matlab/generalised_model/coup.m

function [t,NN,air,U,UVDOT,VVDOT]=coup(order,kappa)%,rinitial)
tic

numBr = 2^order - 1;     % number of airways

% 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

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

y0 = [r(:,1); uv(:,1); vv(:,1); vds(:,1)]; %ICs vector

    parm = [kappa];
    tspan = [0 80];

    % 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)); %airflow
    for j = 1:length(NN)
        [U(j,:),UVDOT(j,:),VVDOT(j,:)] = seven_airway2(t(j),NN(j,:),order,parm);
    end
toc
CM=jet(length(rinitial));
r_new = NN(:,1:numBr);

for i=1:length(rinitial)

   if r_new(end,i) <= 0.01
       air(i)=0;
   else
       air(i)=1;
   end

end
for i = 1:length(rinitial)
    figure(1)
    plot(t,r_new(:,i),'-','color',CM(i,:),'LineWidth', 1.5, 'MarkerSize', 8);
    hold on

end
% %     
% %     
legendStr = cell(1,length(rinitial));
for i = 1:length(legendStr)
    legendStr{i} = sprintf('r_{%d}',i);

end
legend(legendStr,'Location','BestOutside');