- 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_Crossbridge/Reduced_DI_tree/coup_DM.m
function [r1,force,air,summ,pres,q,flow_sum]=coup_DM(Tmax,n,order,kappa)
% to reduce to Austin's model, k1,k2,k7=0
% Force= dt and Rref=0.2792
delta_t=Tmax/n;
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))';
r1 = zeros(length(rinitial),n+1);
epsilon = 0;
%% Distribution Moment ICs
% M10=0.005;
% M11=0.001;
% M12=0.01;
% M20=0.005;
% M21=0.001;
% M22=0.01;
% C=1;
%
%
M10= 0.6209;
M11= 0.3106;
M12= 0.2071;
M20= 0.1749;
M21= 0.0871;
M22= 0.0579;
C= 0.222222222222222;
%
RR = [M10;M11;M12;M20;M21;M22;C];
r1(:,1)=rinitial;
%% Initialising the structure array
output.rdot = NaN;
output.R = NaN;
output.Ptm = NaN;
output.tau = NaN;
output.mu = NaN;
output.q = NaN;
output.p = NaN;
output.pbot = NaN;
output.deltap = NaN;
output.W = NaN;
output.lambda = NaN;
output.Lambda = NaN;
output.gammajplus = NaN;
output.gammajminus = NaN;
for i=1:n
% if mod(i,n) == 100
% fprintf('\niteration #: %g \n',i)
% end
% t(i)=i*delta_t;
% dt(i)=0.4*t(i);
parms = [kappa;epsilon];
[k1,F1,~,~]=Type1p2_Connect_GEN_v3d0_DM(RR,r1(:,i),order,parms);
[k2,F2,~,~]=Type1p2_Connect_GEN_v3d0_DM(RR+(delta_t/2)*F1,r1(:,i)+(delta_t/2)*k1,order,parms);
[k3,F3,~,~]=Type1p2_Connect_GEN_v3d0_DM(RR+(delta_t/2)*F2,r1(:,i)+(delta_t/2)*k2,order,parms);
[k4,F4,out,force(:,i),pbot(i),q(:,i)]=Type1p2_Connect_GEN_v3d0_DM(RR+(delta_t)*F3,r1(:,i)+(delta_t)*k3,order,parms);
Rnew = RR + (delta_t/6)*(F1+(2*F2)+(2*F3)+F4);
RR=Rnew;
r1(:,i+1) = r1(:,i) +(delta_t/6)*(k1+(2*k2)+(2*k3)+k4);
end
time = 0:delta_t:Tmax;
CM=jet(length(rinitial));
air=zeros(1,numBr);
for i=1:length(rinitial)
hold on
% subplot(8,4,i)
%
plot(time,r1(i,:));
% flow_mean(i)=mean(flow(i,:));
if r1(i,end) <= 0.01
air(i)=0;
else
air(i)=1;
end
end
numOrd1 = (numBr + 1)/2;
x = cell(1,numOrd1);
y=cell(1,numOrd1);
for k=1:length(y)
y{k}=q(k,end);
end
maxx=(max([y{:}]));
% % To normalise the flows so that it is between [0,1]
for l = 1:length(y)
y{l} = y{l}./maxx;
end
% for j=1:length(y)
% yy{j}= y{j}./max([y{:}]);
% end
flow_sum = zeros(1,length(y{1}));
for i = 1:length(y)
switch i
case 1
neighbour_diff1 = sqrt((y{i} - y{end}).^2 + (y{i} - y{i+1}).^2);
case length(y)
neighbour_diff1 = sqrt((y{i} - y{i-1}).^2 + (y{i} - y{1}).^2);
otherwise
neighbour_diff1 = sqrt((y{i} - y{i-1}).^2 + (y{i} - y{i+1}).^2);
end
flow_sum = flow_sum + neighbour_diff1;
end
flow_sum=flow_sum./numOrd1;
% flow_sum=flow_sum./(2^(order-3)*4*sqrt(2));
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));
pres=pbot(end);
