- 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/createConstSymm.m
% A function that will output all of the necessary constants for the
% generalised Type II model in a structure called C. This function can
% create these constants for any arbitrary symmetric tree, with the depth
% of the tree (i.e. its order) as input.
% This function calls on "createConnMatrices.m" to construct the
% connectivity matrices.
% DATE: 18/01/18
% AUTHOR: Austin Ibarra
% VERSION: 1.2
% == LOG CHANGE == %
% v1.1: - Added in the transformation matrix T, that maps mu_hat, a vector
% containing the values of mu for order 1 airways, to mu. This will make
% the calculation more efficient, and the reason we can do this is because
% the structure of our airway doesn't change as time goes on.
% v1.2: - Added in the Horsfield constants from orders 9 to 16.
function C = createConstSymm(order)
%% Create the Connectivity Matrices based on order
[C.Cjplus,C.Cjminus,C.Ciplus,C.Ciminus1,C.Ciminus2] = createConnMatrices(order);
%% Construct the vector that defines the order of each airway
numBr = 2^order - 1;
C.order = zeros(1,numBr);
count = 0;
for i = 1:order
for j = count + 1:2^i - 1
C.order(end-count) = order - i + 1;
count = count + 1;
end
end
%% Construct the vector vtop and vbot
% vtop is 1 for the top airway position and 0 elsewhere
% vbot is 1 for the terminal airways, and 0 elsewhere
C.vtop = zeros(numBr,1);
C.vtop(end) = 1;
C.vbot = zeros(numBr,1);
for i = 1:2^(order - 1)
C.vbot(i) = 1;
end
%% v1d1: Construct the transformation T, that maps mu_hat onto mu
C.T = createSymmTMatrix(order);
%% Construct the constants based on Horsfield order
% Defined by Horsfield order. Constants taken from Politi et al (2010)
% From order 1 to 16
C.Ri = [0.058 0.065 0.073 0.083 0.096 0.113 0.132 0.156 0.185 0.222 ...
0.269 0.326 0.395 0.475 0.569 0.686];
C.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];
C.P1 = [15.728 17.342 19.475 22.747 27.140 32.305 39.429 47.104 55.704 65.407 ...
75.968 88.028 100.441 113.457 130.989 153.036];
C.n1 = [1 1 1 1 1 1 1 1 1 1 ...
1 1 1 1 1 1];
C.n2 = [7 7 7.185 7.778 8 8 8 8.148 8.741 9.333 ...
9.926 10 10 10 10 10];
C.L = [1.7 1.8778 2.0556 2.2333 2.4481 2.6852 3.0333 3.3889 3.7444 4.1667 ...
4.6407 5.0630 5.5111 6.1037 6.7556 7.4667];
C.P2 = C.P1.*C.n2.*(C.Ri.^2 - C.rimax.^2)./(C.n1.*C.Ri.^2);
% C.Ri = [0.0058 0.0065 0.0073 0.0083 0.0096 0.0113 0.0132 ...
% 0.0156 0.0185 0.0222 0.0269 0.0326 0.0395 0.0475 0.0569 0.0686];
% C.rimax = [0.0296 0.0318 0.0337 0.0358 0.0384 0.0414 0.0445...
% 0.0484 0.0539 0.0608 0.0692 0.0793 0.0913 0.1052 0.1203 0.1374];
% C.P1 = [0.160380966 0.1768391857 0.1985897325 0.231954847 0.2767509802...
% 0.3283996064 0.4020639056 0.480327125 0.5680227193 0.6669657834 0.7746580127...
% 0.897635788 1.0242131615 1.1569394238 1.3357160702 1.5605329037];
% C.n1 = [1 1 1 1 1 1 1 1 1 1 ...
% 1 1 1 1 1 1];
% C.n2 = [7 7 7.185 7.778 8 8 8 8.148 8.741 9.333 ...
% 9.926 10 10 10 10 10];
% C.L = [1.7 1.8778 2.0556 2.2333 2.4481 2.6852 3.0333 3.3889 3.7444 4.1667 ...
% 4.6407 5.0630 5.5111 6.1037 6.7556 7.4667];
% C.P2 = C.P1.*C.n2.*(C.Ri.^2 - C.rimax.^2)./(C.n1.*C.Ri.^2);
%% Convert the pressures from Pascals in cmH2O
convRate = 0.01019716213;
C.P1 = C.P1*convRate;
C.P2 = C.P2*convRate;
%% Calculate alpha for each airway
a = 10.19716213; % Conversion factor for kg*(mm)^(-1)*(s)^(-2) to cmH20
mu = 1.9008e-8; % dynamic viscosity kg(mm)^(-1)s^(-1)
C.alpha = zeros(1,numBr)';
for i = 1:numBr
C.alpha(i) = (8*a*mu*C.L(C.order(i)))/pi;%pi/(8*a*mu*C.L(C.order(i)));
end
end
