% 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.L = [1.966675000000000,2.888875000000000,4.403700000000000,6.459275000000000];

C.n1 = [1,1,1,1];

C.n2 = [7.240749999999999,8.036999999999999,9.500000000000000,10];

C.P1 = [18.823000000000000,36.494500000000000,71.276750000000000,1.244807500000000e+02];



C.Ri = [0.069750000000000,0.124250000000000,0.250500000000000,0.531250000000000];

C.rimax = [0.327250000000000,0.431750000000000,0.658000000000000,1.135500000000000];

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