- Author:
- aram148 <42922407+aram148@users.noreply.github.com>
- Date:
- 2022-07-26 11:49:42+12:00
- Desc:
- updated documentation and renamed model
- Permanent Source URI:
- https://models.fieldml.org/workspace/7e5/rawfile/30e6db7fb82a9bb715b0121cdfc131c19e5cb2ff/Peripheral_matlab/Order3/coup.asv
function [t,NN,air,U]=coup(order,kappa)%,rinitial)
tic
numBr = 2^order - 1; % number of airways
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))';
r(:,1)=rinitial;
uv(:,1) = zeros(numBr,1);
vv(:,1) = zeros(numBr,1);
vds(:,1) = zeros((numBr+1)/2,1);
y0 = [r(:,1); uv(:,1); vv(:,1); vds(:,1)];
parm = [kappa];
tspan = [0 80];
% Set numerical accuracy options for ODE solver
options = odeset('RelTol', 1e-06, 'AbsTol', 1e-06, 'MaxStep', 0.01);
[t,NN]=ode15s(@(t,y) seven_airway(t,y,order,parm), tspan, y0, options);
U = zeros(length(NN),length(r));
for j = 1:length(NN)
[U(j,:),UVDOT(J)] = seven_airway2(t(j),NN(j,:),order,parm);
end
toc
CM=jet(length(rinitial));
r_new = NN(:,1:7);
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');