% A function that will compute the necessary connectivity matrices of a
% symmetric tree given its order, and will output all 5 connectivity
% matrices.

% This function will be structured as a recursive function, with the base
% case being the trivial 3-branch case (order 2)

% DATE: 18/01/18
% AUTHOR: Austin Ibarra
% VERSION: 1.1

% == LOG CHANGE == %

% v1.1: - Updated the code because we now have Horsfield constants up to
% order 16

function [Cjplus,Cjminus,Ciplus,Ciminus1,Ciminus2] = createConnMatrices(order)
%% Create the Connectivity Matrices based on order

if order < 2      
    %% Throw an error if the depth of the tree is less than 2
    error('A tree must be at least of order 2!');
    
elseif order > 16
    %% Throw an error if the depth of the tree is greater than 16
    error('Currently, a tree cannot be greater than order 16, because the Horsfield constants for higher orders have yet to be coded in.');
elseif order == 2   
    %% The base case
    Cjplus = [1 1 0]';
    Cjminus = [0 0 1]';
    Ciplus = [0 0 1];
    Ciminus1 = [1 0 0];
    Ciminus2 = [0 1 0];
    
else    
    %% Call the function again recursively
    M = 2^(order-1) - 1;   % Number of junctions
    N = 2^order - 1;   % Number of airways
    
    [A,~,~,F,~] = createConnMatrices(order - 1);
    
    B = zeros((N+1)/2,(M-1)/2);
    C = zeros(N,1);
    C((N-1)/2) = 1;
    C((N+1)/2) = 1;
    D = zeros((N+1)/2,M);
    E = eye((N-1)/2);
    G = zeros(1,N);
    G((N-1)/2) = 1;
    
    %% Concatenating the appropriate matrices
    ApB = cat(1,A,B);
    BpA = cat(1,B,A);
    Cjplus = cat(2,ApB,C,BpA);
    
    Cjminus = cat(1,D,E);
    Ciplus = Cjminus';
    
    FpBt = cat(2,F,B');
    BtpF = cat(2,B',F);
    Ciminus1 = cat(1,FpBt,G,BtpF);
    
    Ciminus2 = circshift(Ciminus1,1,2);
end