- Author:
- aram148 <42922407+aram148@users.noreply.github.com>
- Date:
- 2022-07-22 15:40:16+12:00
- Desc:
- Added type up of muscle model
- Permanent Source URI:
- https://models.physiomeproject.org/workspace/6b0/rawfile/d39d197151fcd8e3935a13c0c200987bbb41cd15/peripheral_airway/Peripheral_matlab/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);
%% 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
