%% This program corresponds to Sections 7.4 through 7.6 of the thesis.
% The capacity of a multistorey frame subjected to variable gravity loading
%   is investigated. Two examples - a 2x2 frame and a 3x1 frame can be run
%   using this program by changing the "example" field.
%   The 2x2 frame corresponds to the elastic analysis example.
%   The 3x1 frame corresponds to the inelastic analysis example.
clc, clear, close all
example = '2x2';    % 2x2 or 3x1; corresponding to the various examples being analyzed
proportional = 1;   % For Section 7.4, choose 1 to conduct proportional analysis in Code Block #2, 0 to skip it
variable = 0;       % For Section 7.6, choose 1 for variable gravity loading in Code Block #3, 0 to skip
parametric = 0;     % For Section 7.5, choose 1 for parametric analysis in Code Block #4, 0 to toggle off

% Toggles for considering shear deformations
bool_shear = 0;         % 1 to consider shear deformations
shear_assumption = 1;   % 0 for Engesser, 1 for Haringx
bool_buckle = 1;        % For shear with inelasticity - toggle buckling on in Becque (2016) analysis since we only care about the critical loads
% Currently this is only implemented for the following:
%  ENABLED - Proportional Loading (§2) - Appendix 7.4
%  ENABLED - Variable Loading (§3) - Appendix 7.6
%  HARD_CODE DISABLED - Parametric Study (§4) - not reported in Section 7.5

% Toggles for considering beam axial Deformations
bool_axial = 0;
% Currently this is only implemented for the following:
%  ENABLED - Proportional Loading (§2A - DECOMPOSITION METHOD ONLY (MATRIX METHOD IGNORES AXIAL DEFORMATIONS)) - Appendix 7.4
%  ENABLED - Variable Loading (§3) - Appendix 7.6
%  HARD_CODE DISABLED - Parametric Study (§4) - not reported in Section 7.5

%% 1A. Problem definition - 2-storey elastic example (Section 7.4.1)
if strcmp(example,'2x2')
    % Global properties
    E0 = 2e11; Fy = 350e6;
    kappa = 0.44;
    v0 = 0.3;
    Hk = [3.5296; 3.048];
    H = sum(Hk);
    bool_inel = 0; % No tangent modulus for this example (elastic analysis)

    % Beam properties
    zL = [1 1; 1 1];
    zR = [1 1; 1 1];
    Ib = 8.3246e-5*[10.51 10.51; 8.41 8.41];
    Lb = [6.0168 4.4928; 6.0168 4.4928];
    Ab = 7420e-6*[1.95 1.95; 1.72 1.72];
    
    % Column properties
    Rl = 1e50*[1 1 1]; % rigid
    Lc = [Hk Hk Hk];
    Ic = 8.3246e-5*[4.60 3.13 2.36; 2.35 2.35 1];
    Ac = 7420e-6*[1.40 1.29 1.23; 1.48 1.48 1];
    lambda = [1.1009 1.8731 0.7722; 0.4671 0.8158 0.3487];
    N_beams = [0 0;0 0];
    zl = [1 1 1;1 1 1]; 
    zu = [1 1 1;1 1 1]; 
    n_storeys = size(Ic,1);
    n_cols = size(Ic,2);
    Lr = Lc./(Ic./Ac).^0.5;
    
    % Bracing
    B = zeros(2,3);
    sumB = zeros(2,2);

    % Calibration
    % Without shear deformations - validation in Section 7.4.1
    if bool_shear == 0 && bool_axial == 0
        vb = [-0.0220405/-0.130514 -0.0793291/-0.0220405;-0.0135264/-0.00194714 -0.000244452/-0.0135264]; %all storeys fixed, no bracing
        vc = [1/0 1/0 1/0; -0.00193903/-0.130349 -0.0135094/-0.0219637 -0.000235727/-0.0792165]; %upper over lower, no bracing
        wul = [1/0 1/0 1/0; 0.4544 2.6908 0.7476]; % based on FEA results (note: CW +ve, wFN)
        wlu = [1.7784 10.5308 2.9257; 30.455 4.3846 242.614]; % based on FEA results
    end
    
    % With shear deformations only - validation in Appendix 7.4.1
    if bool_shear == 1 && bool_axial == 0
        vc = [	Inf	Inf	Inf	;	0.030448082	0.492583107	0.035451185	];
        vb = [	0.193938152	2.911618531	;	3.137493472	0.209549058	];		
        wul = [	Inf	Inf	Inf	;	0.457785923	2.360473778	0.810708461	];
        wlu = [	1.926786375	9.935055861	3.412210685	;	15.03496762	4.792031523	22.86830382	];
    end
    
    % With axial deformations only - validation in Appendix 7.4.3
    if bool_shear == 0 && bool_axial == 1
        vc = [	Inf	Inf	Inf	;	0.014252441	0.673614459	0.012893809	];
        vb = [	0.142183192	3.589378502	;	6.720017528	0.068705119	];		
        wul = [	Inf	Inf	Inf	;	0.448638501	2.951416701	0.830352767	];
        wlu = [	1.770468079	12.59424308	3.501552403	;	31.47801124	4.38146281	64.39933729	];
    end
    
    % With both shear and axial deformations - validation in Appendix 7.4.3
    if bool_shear == 1 && bool_axial == 1
        % Original data missing - but can be obtained from the FEA model if desired
    end
    
    wb = [0 0;0 0]; % right to left
    wc_uncal = 1*[1 1 1; 1 1 1]; % used later for uncalibrated results
    
    % Uncalibrated analysis - Section 7.4.1
    % wul = wc_uncal; wlu = wc_uncal; % uncalibrated analysis
    % vb = [1 1;1 1];         % uncalibrated analysis
    % vc = [1 1 1; 1 1 1];    % uncalibrated analysis
    % disp('Warning: uncalibrated analysis is in effect')
    
    % Note - the last uncommented and applicable instance of 
    %  vc, vb, wul, wlu and wb will be used in the analysis
    
%% 1B. Problem definition - 3-storey inelastic example (Section 7.4.2)
elseif strcmp(example,'3x1')
    % Global properties
    E0 = 2e11; Fy = 350e6;
    kappa = 0.44;
    v0 = 0.3;
    Hk = 3*[3.5296; 3.048; 3.048];
    H = sum(Hk);
    bool_inel = 1; % With tangent modulus (inelastic analysis)

    % Beam properties
    Ib = 8.3246e-5*[10.51; 8.41; 8.41];
    Lb = [6.0168; 6.0168; 6.0168];
    Ab = 7420e-6*[1.95; 1.72; 1.72];
    R_con = 25e6;
    zL = (1+3*E0*Ib./Lb./R_con).^-1;
    zR = zL;
    
    % Column propert0ies
    Rl = 1e50*[1 1]; % rigid
    Lc = [Hk Hk];
    Ic = 8.3246e-5*[4.60 3.13; 2.36 2.35; 2.35 1];
    Ac = 7420e-6*[1.40 1.29; 1.23 1.48; 1.48 1];
    lambda = [1.1009 1.8731; 0.7722 0.4671; 0.8158 0.3487];
    N_beams = [0;0;0];
    zl = [1 1;1 1;1 1];
    zu = [1 1;1 1;1 1];
    n_storeys = size(lambda,1);
    n_cols = size(lambda,2);
    
    % Bracing
    sway_dir = "RIGHT"; % Direction of sway chosen for analysis 
    BL = 1e5*[1 0; 1 0; 1 0]; % Left sway bracing stiffness
    BR = 1e5*[0 1; 0 1; 0 1]; % Right sway bracing stiffness
    B = zeros(n_storeys,n_cols);
    if sway_dir == "LEFT"
        B = BL; 
    elseif sway_dir == "RIGHT"
        B = BR;
    else %sway_dir == "NONE"
    end
    sumB = Assemble_Bracing_Matrix(B,sway_dir);

    % Calibration    
    % Without shear deformations - validation in Section 7.4.2
    if bool_shear == 0
        vc = [	Inf	Inf	;	0.099394661	0.349081366	;	0.318587933	-0.003625689	];
        vb = [	0.603416617	;	2.119243578	;	-0.024118045	];			
        wul = [	Inf	Inf	;	0.540697371	0.896059797	;	0.79276503	0.374079242	];
        wlu = [	1.792708428	2.970929833	;	5.439903577	2.566908133	;	-2.488371178	103.1746619	];
    end

    % With shear deformations - validation in Appendix 7.4.2
    if bool_shear == 1
        vc = [	Inf	Inf	;	0.103888915	0.346019616	;	0.308077371	0.000785069	];
        vb = [	0.613128598	;	2.042128577	;	0.005203927	];			
        wul = [	Inf	Inf	;	0.546015919	0.890540615	;	0.778972141	0.381451075	];
        wlu = [	1.785136365	2.911520308	;	5.25576689	2.573670899	;	-2.528495159	-485.8821288	];
    end
    
    wb = [0;0;0];
    wc_uncal = [0 0; 0 0;0 0];  % used later for uncalibrated analysis
    
    % Uncalibrated analysis
    % wul = [1 1;1 1;1 1];
    % wlu = [1 1;1 1;1 1];
    % vc = [1 1;1 1;1 1];
    % vb = [1;1;1];
    % disp('Warning: uncalibrated analysis is in effect')

    % Note - the last uncommented and applicable instance of 
    %  vc, vb, wul, wlu and wb will be used in the analysis
end

%% 1C. Pre-Processing
G0 = E0/(2*(1+v0));
if bool_shear == 1
    nb = E0.*Ib./(Lb.^2)./(kappa.*Ab*G0);
else
    nb = 0.*Ab; % make zero array with the same dimensions
end
RL = 6*E0.*Ib.*zL./(Lb.*(4-zL.*zR + 12.*nb.*(zL+zR+zL.*zR))).*(2+6*nb.*zR+   vb.*zR.*(1-6.*nb) );
RR = 6*E0.*Ib.*zR./(Lb.*(4-zL.*zR + 12.*nb.*(zL+zR+zL.*zR))).*(2+6*nb.*zL+1./vb.*zL.*(1-6.*nb) );
kAGc_0 = kappa*Ac*G0;
kAGb_0 = kappa*Ab*G0;

%% 2. Proportional Loading (Section 7.4)
if strcmp(example,'2x2')
    P_vec = 1e5*[170:1:225]; % Loading vector - data points
else %strcump(example,'3x1')
    if bool_inel == 1
        if bool_shear == 0
            % For no shear portion of the example (Section 7.4.2)
            % For plotting decomposition curve
            P_vec_A = 1e5*[0.01:0.01:9.0];
            P_vec_B = 1e5*[9.01:0.001:9.075];
            P_vec_C = 1e5*[9.076:0.000002:9.087];
            P_vec_D = 1e5*[9.088:0.01:9.7];
            P_vec = [P_vec_A P_vec_B P_vec_C P_vec_D];
            
            % For plotting matrix curve (Section 7.4.2)
            P_vec_A = 1e5*[0.01:0.1:9.01]; % Loading vector - data points
            P_vec_B = 1e5*[9.07:0.01:9.1]; % Loading vector - data points
            P_vec_C = 1e5*[9.2:0.1:10]; % Loading vector - data points
            P_vec = [P_vec_A P_vec_B P_vec_C];
            
            % For uncalibrated analysis (Section 7.4.2)
            %P_vec = 1e5*[9:0.001:10.5];
        else %%bool_shear == 1
            % For with shear portion of the example (Appendix 7.4.2)
            P_vec_A = 1e5*[8.5:0.02:9.03];
            P_vec_B = 1e5*[9.04:0.001:9.07];
            P_vec_C = 1e5*[9.08:0.01:9.2]; % Loading vector - data points
            P_vec = [P_vec_A P_vec_B P_vec_C];
        end
    else
    end
end

if proportional == 1 
%% 2A. Proportional Loading via Decomposition Method (Section 7.4)
% analysis = 1: Analysis with Taylor series, no column contribution to R - obsolete
% analysis = 2: Analysis without Taylor series, no column contribution to R - obsolete
% analysis = 3: Analysis with Taylor series, with column contribution to R - obsolete
% analysis = 4: Analysis without Taylor series, with column contribution to R (this should give an exact result to the matrix method)
XW = zeros(4,length(P_vec)); % Results of lateral stiffness versus gravity loads of each case
Pcr = zeros(1,4); % Critical load for each case
ru_data = zeros(4,length(P_vec),n_storeys,n_cols);
rl_data = zeros(4,length(P_vec),n_storeys,n_cols);
storey_data = zeros(length(P_vec),n_storeys); % for exact method (Case 4), track the stiffness of each storey
for analysis = 1:4
    lcr_storey = zeros(1,n_storeys); % Critical load multiplier for Taylor Series analysis
    for i = 1:length(P_vec)
        P = P_vec(i).*lambda;
        K_storey_based = 1;
        
        for storey = 1:n_storeys
            S_storey = 0;
            lcr_num = 0;
            lcr_den = 0;
            for column = 1:n_cols

                % I. Calculate terms for axial load coefficient
                N = 0;
                for ss = storey:1:n_storeys
                    N = N + P(ss,column);
                end
                if bool_inel == 1
                    E = Get_E_IE(N,Ac(storey,column),E0,Fy);
                    if bool_buckle == 1
                        e = 0.5*(E0/E-1);
                    else
                        e = 0;
                    end
                    G = E0/(2*(1+v0) + 3*e);
                else
                    E = E0;
                    G = E0/(2*(1+v0));
                end
                EI = E*Ic(storey,column);
                kAG = G*kappa*Ac(storey,column);
                
                % II. Calculate end fixities
                alternative = 1; % Choose 1 for decomposition method, 2 for Liu & Xu (2005) method
                col_first_order = 1; % 1 = make column stiffness first order when comparing against Liu and Xu (2005)
                override = 1;        % 1 = use the R' formula directly from Liu and Xu (2005), which is incorrect (for comparison only)
                rl = -109; ru = -109; % Error 109 - if rl or ru don't get updated
                
                % II-1. Use the exact storey decomposition method to calculate the end fixities and phi of the columns
                if alternative == 1
                    if analysis == 3 || analysis == 4
                        [ru,rl,phi,omega] = Get_End_Fixities(...
                            E0,Ic,Hk,P,RL,RR,Rl,vc,wul,wlu,zu,zl,1,bool_inel,Ac,Fy,0,...
                            bool_shear,shear_assumption,kappa,v0,bool_buckle); % Include column contributions
                    else
                        [ru,rl,phi,omega] = Get_End_Fixities(...
                            E0,Ic,Hk,P,RL,RR,Rl,vc,wul,wlu,zu,zl,0,bool_inel,Ac,Fy,0,...
                            bool_shear,shear_assumption,kappa,v0,bool_buckle); % Exclude column contributions
                    end
                    ru = ru(storey,column); ru_data(analysis,i,storey,column) = ru;
                    rl = rl(storey,column); rl_data(analysis,i,storey,column) = rl;
                    phi = phi(storey,column);
                    omega = omega(storey,column);
                    
                % II-2. Alternative method (obsolete) - Distribute Rb via Liu & Xu (2005) for comparison purposes (Section 7.4.1(a))
                elseif alternative == 2
                    % Obsolete - only used for comparison purposes
                    disp('Warning: Liu & Xu method active')
                    RTu = 0; RTl = 0; % Tracks the total rotational restraint at each end of the column in question

                    % a) Beam contributions
                    Rbu = 0; Rbl = 0; % Tracks the total available rotational restraint of the beams at the end of the column in question
                    if column > 1 % Right end of beams
                        Rbu = Rbu + RR(storey,column-1);
                        if storey>1
                            Rbl = Rbl + RR(storey-1,column-1);
                        end
                    end
                    if column < n_cols % Left end of beams
                        Rbu = Rbu + RL(storey,column);
                        if storey > 1
                            Rbl = Rbl + RL(storey-1,column);
                        end
                    end

                    % b) Column contributions
                    if storey > 1 % Lower end of current column
                        if bool_inel == 1
                            E_below = Get_E_IE(N+P(storey-1,column),Ac(storey-1,column),E0,Fy);
                        else
                            E_below = E0;
                        end
                        EI_below = E_below*Ic(storey-1,column);
                        RCu_below = R_prime(zu(storey-1,column),zl(storey-1,column),N+P(storey-1,column),EI_below,Hk(storey-1),1/vc(storey-1,column),wlu(storey-1,column),0,0,0);
                        RCl_curr = R_prime(zl(storey,column),zu(storey,column),N,EI,Hk(storey),vc(storey,column),wul(storey,column),0,0,0);
                        if col_first_order == 1 % R based on first order and also neglects w in Xu and Liu (2007)
                            RCu_below = R_prime(zu(storey-1,column),zl(storey-1,column),0,EI_below,Hk(storey-1),1/vc(storey-1,column),0,0,0,0);
                            RCl_curr = R_prime(zl(storey,column),zu(storey,column),0,EI,Hk(storey),vc(storey,column),0,0,0,0);
                        end
                        
                        % Override - use the formula directly from Liu & Xu (2005)
                        if override == 1
                            Rb_below = 0;
                            if column > 1 % Right end of beams
                                Rb_below = Rb_below + RR(storey-1,column-1);
                            end
                            if column < n_cols % Left end of beams
                                Rb_below = Rb_below + RL(storey-1,column);
                            end
                            phi_below = sqrt(N+P(storey-1,column) / EI_below)*Hk(storey-1);
                            phi_curr = sqrt(N/EI)*Hk(storey);
                            RCu_below = EI_below/Hk(storey-1) *( Rb_below*Hk(storey-1)/EI_below - phi_below*tan(phi_below) ) /...
                                        ( 1+ Rb_below*Hk(storey-1)/EI_below*tan(phi_below)/phi_below );
                            RCl_curr = EI/Hk(storey) * (Rbl*Hk(storey)/EI - phi_curr*tan(phi_curr) ) /...
                                        ( 1+ Rbl*Hk(storey)/EI*tan(phi_curr)/phi_curr);
                        end
                        % End override - comment this block out to use the above expressions
                        
                        mu_l = RCl_curr/(RCl_curr + RCu_below);
                        RTl = mu_l*Rbl;
                        if analysis == 3 || analysis == 4
                            RTl = RTl + RCu_below; % Add column contribution to R
                        end
                    else % Lower end to ground
                        RTl = Rl(column);
                    end
                    if storey < 2 % Upper end of current column
                        if bool_inel == 1
                            E_above = Get_E_IE(N-P(storey,column),Ac(storey+1,column),E0,Fy);
                        else
                            E_above = E0;
                        end
                        EI_above = E_above*Ic(storey+1,column);
                        RCl_above = R_prime(zl(storey+1,column),zu(storey+1,column),N-P(storey,column),EI_above,Hk(storey+1),vc(storey+1,column),wul(storey+1,column),0,0,0);
                        RCu_curr = R_prime(zu(storey,column),zl(storey,column),N,EI,Hk(storey),1/vc(storey,column),wlu(storey,column),0,0,0);
                        if col_first_order == 1 % R based on first order and also neglects w in Xu and Liu (2007)
                            RCl_above = R_prime(zl(storey+1,column),zu(storey+1,column),0,EI_above,Hk(storey+1),vc(storey+1,column),0,0,0,0);
                            RCu_curr = R_prime(zu(storey,column),zl(storey,column),0,EI,Hk(storey),1/vc(storey,column),0,0,0,0);
                        end
                        
                        % Override - use the formula in Liu & Xu (2005) 
                        if override == 1
                            Rb_above = 0;
                            if column > 1 % Right end of beams
                                Rb_above = Rb_above + RR(storey+1,column-1);
                            end
                            if column < n_cols % Left end of beams
                                Rb_above = Rb_above + RL(storey+1,column);
                            end
                            phi_above = sqrt(N-P(storey,column) / EI_above)*Hk(storey+1);
                            phi_curr = sqrt(N/EI)*Hk(storey);
                            RCl_above = EI_above/Hk(storey+1) *( Rb_above*Hk(storey+1)/EI_above - phi_above*tan(phi_above) ) /...
                                        ( 1+ Rb_above*Hk(storey+1)/EI_above*tan(phi_above)/phi_above );
                            RCu_curr = EI/Hk(storey) * (Rbu*Hk(storey)/EI - phi_curr*tan(phi_curr) ) /...
                                        ( 1+ Rbu*Hk(storey)/EI*tan(phi_curr)/phi_curr);
                        end
                        % End override - comment this block out to use the above expressions
                        
                        mu_u = RCu_curr/(RCu_curr + RCl_above);
                        RTu = mu_u*Rbu;
                        if analysis == 3 || analysis == 4
                            RTu = RTu + RCl_above; % Add column contribution to R
                        end
                    else % Highest storey
                        RTu = Rbu;
                    end

                    % c) Calculate resulting end fixity factors
                    ru = 1/(1+3*EI/(RTu*Hk(storey))); ru_data(analysis,i,storey,column) = ru;
                    rl = 1/(1+3*EI/(RTl*Hk(storey))); rl_data(analysis,i,storey,column) = rl;
                    phi = sqrt(N/EI)*Hk(storey);
                    omega = 1;
                end
                
                % For the 2x2 example, the following overrides gives Pcr = 108.2 MN and is
                % obtained from directly substituting the end fixity factors reported in Liu & Xu (2005)
                % if column == 1
                %    if storey == 1
                %        ru = 0.635; rl = 1;
                %    else
                %        ru = 0.784; rl = 0.614;
                %    end
                % elseif column ==2
                %    if storey == 1
                %        ru = 0.834; rl = 1;
                %    else
                %        ru = 0.895; rl = 0.828;
                %    end
                % else
                %    if storey == 1
                %        ru = 0.826; rl = 1;
                %    else
                %        ru = 0.919; rl = 0.822;
                %    end
                % end

                % III. Calculate Stiffness
                a1 = 3*(rl*(1-ru) + ru*(1-rl));
                a2 = 9*rl*ru - (1-rl)*(1-ru)*(omega*phi)^2;
                a3 = 18*rl*ru + a1*(omega*phi)^2;
                beta_0 = (rl+ru+rl*ru)/(4-rl*ru);
                beta_1 = (8*(5+ru^2)-(34-ru)*ru*rl + (8+ru+3*ru^2)*rl^2)/(30*(4-rl*ru)^2); % obsolete - Taylor series is garbage, will probably give wrong answer for shear
                beta = 1/12*omega*phi^3 * ...
                        (a1*omega*phi*cos(omega*phi) + a2*sin(omega*phi)) /...
                        (18*rl*ru - a3*cos(omega*phi) + (a1-a2)*omega*phi*sin(omega*phi));
                zeta = (omega^2-1)*( 18*rl*ru*(1-cos(omega*phi)) + a1*omega*phi*sin(omega*phi) )/...
                        (18*rl*ru - a3*cos(omega*phi) + (a1-a2)*omega*phi*sin(omega*phi));
                zeta_0 = 12*EI/Hk(storey)^2/kAG*beta_0;
                
                if phi == 0 % no load
                    beta = beta_0;
                    zeta = zeta_0;
                end
                if analysis == 1 || analysis == 3
                    Si = 12*(EI/Hk(storey)^3*beta_0 - N/Hk(storey)*beta_1); % obsolete - Taylor series is garbage, will probably give wrong answer for shear
                    if bool_shear == 1
                        disp('Warning - Taylor series not implemented for shear deformations in proportional loading analysis')
                    end
                else
                    Si = 12*EI/Hk(storey)^3*beta/(1+zeta);
                end
                S_storey = S_storey + Si + B(storey,column);
                if bool_axial == 1 && column < n_cols % if approximating the effect of axial deformations
                    kb = E0*Ab(storey,column)/Lb(storey,column);
                	S_storey = 1/(1/S_storey + 1/kb);
                end
                
                % IV. Track the critical load multiplier (comparison with Liu & Xu (2005) Taylor series results)
                lcr_num = lcr_num + EI*beta_0/Hk(storey)^3;            
                cumulative_lambda = N/P_vec(i);                         
                lcr_den = lcr_den + cumulative_lambda*beta_1/Hk(storey);
            end
            
            % Product sum of storey lateral stiffness
            K_storey_based = K_storey_based * S_storey;
            storey_data(i,storey) = S_storey;
            lcr_storey(storey) = lcr_num/lcr_den;
        end
        
        % Estimate critical load multiplier (comparison with Liu & Xu (2005) Taylor series results)
        if analysis == 1 || analysis == 3 % 
            Pcr(analysis) = min(lcr_storey)*sum(sum(lambda)); 
        else
            Pcr(analysis) = 1/0;
        end
        XW(analysis,i) = K_storey_based;
    end
end

% Plot the equivalent lateral stiffness of each storey
figure('Renderer', 'painters', 'Position', [10 10 550 330])
set(gca, 'FontName', 'Times New Roman');
legend_titles = {};
for storey=1:n_storeys
    if storey == 1
        str_style = 'k-';
    elseif storey == 2
        str_style = 'k--';
    else
        str_style = 'k-.';
    end
    plot(P_vec*sum(sum(lambda)),storey_data(:,storey),str_style)
    hold on
    grid on
    legend_titles(storey) = {strcat('Storey ',num2str(storey))};
end
lgnd = legend(legend_titles,'location','southwest')
xlabel('Total Applied Load [N]')
ylabel('Storey Lateral Stiffness, \SigmaS [N/m]')
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');
set(gca, 'FontName', 'Times New Roman')
if strcmp(example,'3x1')
    axis([3e6 5e6 -1e7 1e7]) 
end

%% 2B. Proportional Loading: Matrix Method Results
RL_1O = 6*E0.*Ib.*zL./(Lb.*(4-zL.*zR)).*(2+vb.*zR); % first-order no shear only - fyi
RR_1O = 6*E0.*Ib.*zR./(Lb.*(4-zL.*zR)).*(2+1./vb.*zL); % first-order no shear only - fyi
data_full = zeros(1,length(P_vec));
data_dia = zeros(1,length(P_vec));
for i = 1:length(P_vec)
    Load = P_vec(i)
    
    % Calculate second-order rotational stiffness of beam ends
    Nb = Load*N_beams;
    for ii = 1:size(Lb,1)
        for jj = 1:size(Lb,2)
            RL(ii,jj) = R_prime(zL(ii,jj),zR(ii,jj),Nb(ii,jj),E0*Ib(ii,jj),Lb(ii,jj),vb(ii,jj),wb(ii,jj),...
                bool_shear,kappa*Ab(ii,jj)*E0/(2*(1+v0)),shear_assumption);
            RR(ii,jj) = R_prime(zR(ii,jj),zL(ii,jj),Nb(ii,jj),E0*Ib(ii,jj),Lb(ii,jj),1/vb(ii,jj),wb(ii,jj),...
                bool_shear,kappa*Ab(ii,jj)*E0/(2*(1+v0)),shear_assumption);
        end
    end
    %RL_pct = RL./RL_1O
    %RR_pct = RR./RR_1O
    
    K_cols = zeros(n_storeys,n_storeys,n_cols);
    K_tot = zeros(n_storeys,n_storeys);
    for col = 1:n_cols % Determine the K matrix contribution of each column
        Pk = Load*lambda(:,col);
        
        % Calculate axial forces
        N = zeros(1,n_storeys);
        EI = zeros(1,n_storeys);
        kAG = zeros(1,n_storeys);
        phi = zeros(1,n_storeys);
        for storey = 1:n_storeys
            N(storey) = 0;
            for ss = storey:n_storeys
                N(storey) = N(storey) + Pk(ss);
            end
            if bool_inel == 1
                E = Get_E_IE(N(storey),Ac(storey,col),E0,Fy);
                if bool_buckle == 1
                    e = 0.5*(E0/E-1);
                else
                    e = 0;
                end
                G = E0/(2*(1+v0) + 3*e);
            else
                E = E0;
                G = E0/(2*(1+v0));
            end
            EI(storey) = E*Ic(storey,col);
            kAG(storey) = G*kappa*Ac(storey,col);
            phi(storey) = H*(N(storey)/EI(storey))^0.5;
        end
        
        % Calculate R at each level
        Rk(1) = Rl(col); % The lowest level is connected to the ground (R0)
        for storey = 1:n_storeys % Each level goes to Rk(storey+1) since MATLAB sucks and doesn't have zero indexing
            Rk(storey+1) = 0;
            if col < n_cols % Left end contributions
                Rk(storey+1) = Rk(storey+1) + RL(storey,col);
            end % Right end contributions
            if col > 1
                Rk(storey+1) = Rk(storey+1) + RR(storey,col-1);
            end
        end
        
        % Solve for the K matrix for the current column
        K_cols(:,:,col) = Subs_K_matrix_m_storeys(Hk, Rk, Pk, EI, bool_shear, kAG, shear_assumption);
        %K_cols(:,:,col) = Subs_Continuous_column_2_storey(Hk, Rk, Pk, EI)
        K_tot = K_tot + K_cols(:,:,col);
    end
    K_tot = K_tot + sumB; % Add the global bracing stiffness matrix
    K_tot
    det_K = det(K_tot)
    K_dia = 1;
    for ss = 1:n_storeys
        K_dia = K_dia*K_tot(ss,ss);
    end
    data_full(i) = det_K;
    data_dia(i) = K_dia;
   
end

%% 2C. Proportional Loading: Comparison
figure('Renderer', 'painters', 'Position', [10 10 500 325])
plot(P_vec,data_full) % Full matrix method
grid on
xlabel('P'); ylabel('|K_{eq}|');
hold on
plot(P_vec,data_dia) % Matrix method, diagonal only
for i=1:4
    plot(P_vec,XW(i,:),'--');
end
legend('Full matrix', 'Diagonal matrix',...
    'Xu & Wang - Case 1','Xu & Wang - Case 2','Xu & Wang - Case 3','Xu & Wang - Case 4')
figure('Renderer', 'painters', 'Position', [10 10 500 200])
set(gca, 'FontName', 'Times New Roman');
plot(P_vec*sum(sum(lambda)),data_full,'k-') % Full matrix method
grid on
xlabel('Total Applied Load [N]');
if strcmp(example,'3x1')
    ylabel('|K_{eq}| [N/m]^3');
    %axis([4.85e6 5.375e6 -2e18 2e18]); %10^6
    axis([3e6 5e6 -1e17 2e17]); %10^5
elseif strcmp(example,'2x2')
    ylabel('|K| [N/m]^2');
end
set(gca, 'FontName', 'Times New Roman');
end

%% 3. Variable Loading Analysis
if variable == 1

% 3.1   Shape parameters for iterating the variable loading analysis (runs the last block of vc,vb,wul,wlu)
vc(1,[1:n_cols]) = 1; % was 1/0
vc([2:n_storeys],[1:n_cols]) = 1;   % Uncalibrated - assume vc since it is not known 
vb = ones(n_storeys,n_cols-1);      % Uncalibrated - assume vb since it is not known
wlu = ones(n_storeys,n_cols);       % Uncalibrated - assume wlu since it is not known
wul = ones(n_storeys,n_cols);       % Uncalibrated - assume wul since it is not known
wb = zeros(n_storeys,n_cols-1);
Pl = 300e3*ones(n_storeys,n_cols);

% 3.2	Get Pu based on rotational buckling (obsolete - not very useful)
if strcmp(example,'3x1')  % Calibration for determining Pu (assuming vb = 1 and using the calibrated rotationally buckled shape -- not very useful)
    vc_buckle(:,:,1) = [1/0 1/0; -0.175107/0.781289 -0.190059/0.67832; 0.0619285/-0.175107 0.0482274/-0.190059]; % Col 1
    vc_buckle(:,:,2) = [1/0 1/0; 0.618458/-0.68923 0.723726/-0.817626; -0.218853/0.618458 -0.183645/0.723726]; % Col 2
    vc_buckle(:,:,3) = [1/0 1/0; -0.853888/0.364554 -0.477887/0.16874; 0.859668/-0.853888 0.732261/-0.477887]; % Col 3
    
    vc_buckle(:,:,1) = vc;
    vc_buckle(:,:,2) = vc;
    vc_buckle(:,:,3) = vc;
else %strcmp(example,'2x2')
    vc_buckle(:,:,1) = vc; % I didn't bother finding these - they're not very useful
    vc_buckle(:,:,2) = vc; % I didn't bother finding these - they're not very useful
end
RL = 6*E0.*Ib.*zL./(Lb.*(4-zL.*zR)).*(2+vb.*zR);
RR = 6*E0.*Ib.*zR./(Lb.*(4-zL.*zR)).*(2+1./vb.*zL);
Nu = Iterate_Nu_Multi(E0,Ic,Hk,RL,RR,Rl,vc_buckle,... % should really use vc during rotational buckling, which would require separate analysis -- but doesn't matter since Pu is a dummy input and can be anything
    ones(n_storeys,n_cols),ones(n_storeys,n_cols),1,...
    bool_inel,Ac,Fy,0,[],[],[],[],0,wul,wlu,...
    bool_shear,shear_assumption,kappa,v0);
Pu = Nu;
for i=1:n_storeys
    for j=1:n_cols
        for ss = 1:i
            if Nu(ss,j) < Pu(i,j)
                Pu(i,j) = Nu(ss,j);
            end
        end
    end
end
Pu = 2*Pu; % Needs to be a positive input but doesn't actually change the problem - doesn't matter what Pu is as long as it's reasonable

% 3.3. Variable Loading - Decomposition Method
[P_opt, P_vec, err_flag] = var_analysis_decomp(...
    "MIN",E0,Hk,Lb,Ib,Ic,zL,zR,Rl,vb,vc,wb,wul,wlu,B,Pl,Pu,bool_inel,Ac,Ab,Fy,...
    bool_shear,shear_assumption,kappa,v0,bool_axial)
[P_opt, P_vec, err_flag] = var_analysis_decomp(...
    "MAX",E0,Hk,Lb,Ib,Ic,zL,zR,Rl,vb,vc,wb,wul,wlu,B,Pl,Pu,bool_inel,Ac,Ab,Fy,...
    bool_shear,shear_assumption,kappa,v0,bool_axial)
end

%% 4. Parametric Analyses (2x2 frame) - Section 7.5
bool_shear = 0; % hard-code disabled
bool_axial = 0; % hard-code disabled
if parametric == 1
if strcmp(example,'2x2')	
    fprintf('\nBeginning parametric study for 2x2 frame example')
else
    fprintf('\nError - 3x1 frame example not supported in parametric analysis')
end

%% 4.1   Effect of shape parameters on critical load (Section 7.5.1)
%  Case 1 - Asymmetric Buckling (v = 1, w = w0)
%  Case 2 - *OBSOLETE* Symmetric Buckling Mode I (v = -1, w = +/- w0) - gives bad results
%  Case 3 - *OBSOLETE* Symmetric Buckling Mode II (v = -1, w = -/+ w0) - gives bad results
%  Case 4 - Actual calibrated solution
n_mode_cases = 4;
lambda = [1.1009 1.8731 0.7722; 0.4671 0.8158 0.3487];
P_vec = 1e5*[0:0.1:400];
w0_vec = [0 0.2 0.4 0.6 0.8 1];
K_data = zeros(length(P_vec),n_storeys,length(w0_vec),n_mode_cases);
K_prod_data= zeros(length(P_vec),length(w0_vec),n_mode_cases);
K_min_data = zeros(length(P_vec),length(w0_vec),n_mode_cases);
for w0_case = 1:length(w0_vec) 
    % Value of w
    w0 = w0_vec(w0_case);
    
    % Mode cases
    for mode_case = 1:n_mode_cases 
        wb = [0 0;0 0];
       if mode_case == 1 % Asymmetrical buckling (in thesis)
            vc = [1 1 1 ;1 1 1]; 
            vb = [1 1;1 1];		
            wul = [	w0 w0 w0 ; w0 w0 w0];
            wlu = [w0 w0 w0; w0 w0 w0];
       elseif mode_case == 2 % Symmetric buckling (not in thesis)
            vc = [-1 -1 -1 ;-1 -1 -1];
            vb = [-1 -1;-1 -1];		
            wul = [	-w0 w0 -w0	; w0 -w0 w0];
            wlu = [w0 -w0 w0; -w0 w0 -w0];
       elseif mode_case == 3 % Symmetric buckling (not in thesis)
            vc = [-1 -1 -1 ;-1 -1 -1]; 
            vb = [-1 -1;-1 -1];		
            wul = [	w0 -w0 w0	; -w0 w0 -w0];
            wlu = [-w0 w0 -w0; w0 -w0 w0];
       elseif mode_case == 4 % Calibrated results
            vb = [-0.0220405/-0.130514 -0.0793291/-0.0220405;-0.0135264/-0.00194714 -0.000244452/-0.0135264]; 
            vc = [1/0 1/0 1/0; -0.00193903/-0.130349 -0.0135094/-0.0219637 -0.000235727/-0.0792165];
            wul = [1/0 1/0 1/0; 0.4544 2.6908 0.7476]; 
            wlu = [1.7784 10.5308 2.9257; 30.455 4.3846 242.614];
       else % invalid
       end
       RL = 6*E0.*Ib.*zL./(Lb.*(4-zL.*zR)).*(2+vb.*zR);
       RR = 6*E0.*Ib.*zR./(Lb.*(4-zL.*zR)).*(2+1./vb.*zL);
       
       % Run K vs P proportional loading
       for P_case = 1:length(P_vec) 
            P = P_vec(P_case)*lambda;
            K_storey_based = 1;
            for storey = 1:n_storeys
                S_storey = 0;
                lcr_num = 0;
                lcr_den = 0;
                for column = 1:n_cols
                    N = 0;
                    for ss = storey:1:n_storeys
                        N = N + P(ss,column);
                    end
                    if bool_inel == 1
                        E = Get_E_IE(N,Ac(storey,column),E0,Fy);
                    else
                        E = E0;
                    end
                    EI = E*Ic(storey,column);
                    [ru,rl,phi] = Get_End_Fixities(...
                        E0,Ic,Hk,P,RL,RR,Rl,vc,wul,wlu,zu,zl,1,bool_inel,Ac,Fy,0,...
                        bool_shear,shear_assumption,kappa,v0,bool_buckle); % Include column contributions
                    ru = ru(storey,column);
                    rl = rl(storey,column);
                    phi = phi(storey,column);
                    a1 = 3*(rl*(1-ru) + ru*(1-rl));
                    a2 = 9*rl*ru - (1-rl)*(1-ru)*phi^2;
                    a3 = 18*rl*ru + a1*phi^2;
                    beta_0 = (rl+ru+rl*ru)/(4-rl*ru);
                    beta_1 = (8*(5+ru^2)-(34-ru)*ru*rl + (8+ru+3*ru^2)*rl^2)/(30*(4-rl*ru)^2);
                    beta = 1/12*phi^3 * ...
                            (a1*phi*cos(phi) + a2*sin(phi)) /...
                            (18*rl*ru - a3*cos(phi) + (a1-a2)*phi*sin(phi));
                    if phi == 0 % no load
                        beta = beta_0;
                    end
                    Si = 12*EI/Hk(storey)^3*beta;
                    S_storey = S_storey + Si;
                end
                S_storey = S_storey + sum(B(storey,:)); 
                K_data(P_case,storey,w0_case,mode_case) = S_storey;
                K_storey_based = K_storey_based * S_storey;
            end
            K_min_data(P_case,w0_case,mode_case) = min(K_data(P_case,:,w0_case,mode_case));
            K_prod_data(P_case,w0_case,mode_case) = K_storey_based;
       end
    end
end
for mode_case = 1:n_mode_cases
    for storey = 1:n_storeys
        figure % Plot for each storey, each mode case
        for w0_case = 1:length(w0_vec)
            colour = w0_case/length(w0_vec)*0.75*[1 1 1];
            plot(sum(sum(lambda))*P_vec,K_data(:,storey,w0_case,mode_case),'-','Color',colour)
            legend_title(w0_case) = strcat("w0 = ", num2str(w0_vec(w0_case)));
            hold on
        end
        grid on
        title_string = strcat("Mode case ",num2str(mode_case),", Storey ",num2str(storey));
        title(title_string);
        lgnd = legend(legend_title);
        temp = [lgnd; lgnd.ItemText];
        set(temp, 'FontSize', 10);
        set(temp, 'FontName', 'Times New Roman');
        axis([0, sum(sum(lambda))*P_vec(length(P_vec)), -1e8, 1e8]);
        xlabel('Total Load [N]')
        ylabel('Storey lateral stiffness, {\Sigma}S [N/m]')
        set(gca, 'FontName', 'Times New Roman');
    end
    % Since Storey 1 governs instability for all of the above plots, plot prodK
    figure
    for w0_case = 1:length(w0_vec)
        colour = w0_case/length(w0_vec)*0.75*[1 1 1];
        plot(sum(sum(lambda))*P_vec,K_prod_data(:,w0_case,mode_case),'-','Color',colour)
        legend_title(w0_case) = strcat("w_0 = ", num2str(w0_vec(w0_case)));
        hold on
    end
    grid on
    title_string = strcat("Mode case ",num2str(mode_case),", Product K"); 
    title_string = (''); % Include this line for writing the manuscript
    title(title_string);
    lgnd = legend(legend_title);
    temp = [lgnd; lgnd.ItemText];
    set(temp, 'FontSize', 10);
    set(temp, 'FontName', 'Times New Roman');
    axis([0, sum(sum(lambda))*P_vec(length(P_vec)), -5e14, 1.5e15]);
    xlabel('Total Load [N]')
    ylabel('Product of storey lateral stiffness, {\Pi\Sigma}S [N/m]^2')
    set(gca, 'FontName', 'Times New Roman');
    
end

%% 4.2   Accuracy of asymmetrical buckling assumption on critical load
% Generates 1000 test cases to be run in ABAQUS
w0 = 1;
vc = [1 1 1;1 1 1];
vb = [1 1;1 1];		
wul = [w0 w0 w0; w0 w0 w0];
wlu = [w0 w0 w0; w0 w0 w0];
num_sims = 1000;
rng('default');
Pcr_data = zeros(num_sims,1);
filename = strcat('para_cases_output_new_',string(datetime('now','TimeZone','local','Format','yMMddHHmmss')),'.txt');
fileID = fopen(filename,'w');
fprintf(fileID,'\nCase\t...X(1,1)\tX(1,2)\tX(1,3)\tX(2,1)\tX(2,2)\tX(2,3)\tIc(1,1)\tIc(1,2)\tIc(1,3)\tIc(2,1)\tIc(2,2)\tIc(2,3)\tIb(1,1)\tIb(1,2)\tIb(2,1)\tIb(2,2)\tRl(1)\tRl(2)\tRl(3)\tR_conn\tPcr');
for sim_case = 1:num_sims
    Ic = 83.246e-6*[ 1.0 + 10*rand(n_storeys,n_cols)]; % From 1 to 11 I
    Ib = 83.246e-6*[ 1.0 + 10*rand(n_storeys,n_cols-1)]; % From 1 to 11 I
    R_conn = rand()*100e6; % Beam onnection fixity - vary from zero to 100e6
    Rl = 10.^(rand(1,n_cols)*12); % Vary from 10^0 (virtually pinned) to 10^12 (virtually fixed)
    zL = (1+3*E0.*Ib./Lb./R_conn).^(-1); % Left end same as right end
    zR = (1+3*E0.*Ib./Lb./R_conn).^(-1); % Left end same as right end
    RL = 6*E0.*Ib.*zL./(Lb.*(4-zL.*zR)).*(2+vb.*zR);
    RR = 6*E0.*Ib.*zR./(Lb.*(4-zL.*zR)).*(2+1./vb.*zL);
    X = rand(n_storeys,n_cols);
    X = X./sum(sum(X)); % normalize so that the total load is Pcr
    
    Pl = zeros(n_storeys,n_cols);
    Pu = 1e8*ones(n_storeys,n_cols); % make it around Pu (it doesn't really matter what it is)
    [soln,prodS_val,num_evals,err_flag] = get_crit_loads_decomp_Brent(...
        X,[],Pl,Pu,E0,Ic,Rl,RL,RR,vc,wul,wlu,Hk,B,bool_inel,Ac,Fy,...
                        bool_shear,shear_assumption,kappa,v0,bool_axial,Ab,Lb);
    Pcr_data(sim_case) = sum(sum(soln)); % Store the data point
    fprintf(fileID,'\n%.0d\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e\t%.15e',...
        sim_case,X(1,1),X(1,2),X(1,3),X(2,1),X(2,2),X(2,3),...
        Ic(1,1),Ic(1,2),Ic(1,3),Ic(2,1),Ic(2,2),Ic(2,3),Ib(1,1),Ib(1,2),Ib(2,1),Ib(2,2),...
        Rl(1),Rl(2),Rl(3),R_conn,Pcr_data(sim_case));
end
fclose(fileID);
% Now run the following files:
% 1. case_creator.m (make sure the input file path is updated)
% 2. Run all of the generated ABAQUS files via Command prompt (NOTE: THE VERSION OF ABAQUS USED TO COMPLETE THIS WAS 6.14 - IT MIGHT NOT WORK FOR OTHER VERSIONS)
% 3. Run case_analyzer.m

end

%% 5. Helper functions 

%% Variable Loading Algorithm: Pseudo-Proportional (Section 7.6 of the thesis)
% MIN_OR_MAX = string, "MIN" for minimization, "MAX" for maximization
% E0 = scalar, elastic modulus of all members
% Hk = [n_storeys] vector, height of each storey from bottom to top
% Lb = [n_storey × n_bays] matrix, lengths of beams
% Ib = [n_storey × n_bays] matrix, inertias of beams
% Ic = [n_storey × n_cols] matrix, inertias of columns
% rL = [n_storey × n_bays] matrix, left end fixity factors of beams
% rR = [n_storey × n_bays] matrix, right end fixity factors of beams
% Rl = [n_cols] vector, ground connection rotational rigidities of columns
% vb = [n_storey × n_bays] matrix, end connection rotation ratios of beams (assumed) -- theta_right / theta_left
% vc = [n_storey * n_cols] matrix, end connection rotation ratios of columns (assumed) -- theta_upper / theta_lower
% wb = [n_storey × n_bays] matrix, ratio of chord rotation to left end rotation for beams
% wul = [n_storey × n_cols] matrix, ratio of chord rotation to lower end rotation for columns
% wlu = [n_storey × n_cols] matrix, ratio of chord rotation to upper end rotation for columns
% B = [n_storey × n_cols] matrix, bracing stiffnesses at each node
% Pl = [n_storey × n_cols] matrix, lower bound applied load at each column
function [P_opt, P_vec, err_flag] = var_analysis_decomp(MIN_OR_MAX,E0,Hk,Lb,Ib,Ic,rL,rR,Rl,vb,vc,wb,wul,wlu,B,Pl,Pu,bool_inel,Ac,Ab,Fy,...
                        bool_shear,shear_assumption,kappa,v0,bool_axial)

    %% A. Setup the problem
    num_evals = 0;      % Number of times the objective function has been called (note: not the same as the number of root finds)
    num_iters = 1;      % Number of iterations of PSO that have been completed
    max_iters = 15;     % Maximum number of iterations of PSO before quitting
    n_multistart = 1;   % Number of starting points in Newton's method for finding roots (doesn't work for inelastic yet)
    n_PSO = 3;          % Total number of times that the entire PSO algorithm will run before finalizing (will be different if the mask changes)
    swarm_size = 50;    % Swarm size for PSO
    max_newton = 100;   % Maximum number of Newton's method iterations before an error is returned  
    opt_tol = 10;       % Convergence tolerance of the total load for the optimal solution, in Newtons
    rng default;
    K_tol = (1E3)^(size(Ic,1)); % Value of prod(K) below which it can sufficiently be considered to be zero (unit: [N/m]^(# storeys))
    G0 = E0/(2*(1+v0));
    if bool_shear == 1
        nb = E0.*Ib./(Lb.^2)./(kappa.*Ab*G0);
    else
        nb = 0.*Ab; % make zero array with the same dimensions
    end
    RL = 6*E0.*Ib.*rL./(Lb.*(4-rL.*rR + 12.*nb.*(rL+rR+rL.*rR))).*(2+6*nb.*rR+   vb.*rR.*(1-6.*nb) );
    RR = 6*E0.*Ib.*rR./(Lb.*(4-rL.*rR + 12.*nb.*(rL+rR+rL.*rR))).*(2+6*nb.*rL+1./vb.*rL.*(1-6.*nb) );
    n_storeys = size(Ic,1);
    n_cols = size(Ic,2);
    n_bays = n_cols-1;
    %filename = strcat('var_opt_output_v2_',string(datetime('now','TimeZone','local','Format','yMMddHHmmss')),'.txt');
    %fileID = fopen(filename,'w');
    
    
    %% B. Verify input
    % i. Confirm that the input variables have the correct dimensions
    if (length(E0)~=1 || length(Hk) ~= n_storeys || size(Lb,1) ~= n_storeys || size(Lb,2) ~= n_bays ||...
            size(Ib,1) ~= n_storeys || size(Ib,2) ~= n_bays || size(Ic,1) ~= n_storeys || size(Ic,2) ~= n_cols ||...
            size(rL,1) ~= n_storeys || size(rL,2) ~= n_bays || size(rR,1) ~= n_storeys || size(rR,2) ~= n_bays ||...
            size(vb,1) ~= n_storeys || size(vb,2) ~= n_bays || size(vc,1) ~= n_storeys || size(vc,2) ~= n_cols ||...
            size(wb,1) ~= n_storeys || size(wb,2) ~= n_bays || size(wul,1) ~= n_storeys|| size(wul,2) ~= n_cols ||...
            size(wlu,1) ~= n_storeys|| size(wlu,2) ~= n_cols|| size(Pl,1) ~= n_storeys || size(Pl,2) ~=n_cols ||...
            length(Rl) ~= n_cols || size(Ab,1) ~= n_storeys || size(Ab,2) ~= n_bays || size(Ac,1) ~= n_storeys ||...
            size(Ac,2) ~= n_cols || size(B,1) ~= n_storeys || size(B,2) ~= n_cols || size(Pu,1) ~= n_storeys ||...
            size(Pu,2) ~= n_cols)
        fprintf('\nError 100 - Invalid input dimensions')
        P_vec = -100;
        P_opt = -100;
        err_flag = -100;
        return
    end
    % ii. Confirm that prodS > 0 at Pl
    [prodS_zero,~] = prodS(E0,Ic,Hk,Pl,RL,RR,Rl,vc,wul,wlu,ones(n_storeys,n_cols),ones(n_storeys,n_cols),B,1,bool_inel,Ac,Fy,...
                        bool_shear,shear_assumption,kappa,v0,1,bool_axial,Ab,Lb);
    if prodS_zero < 0
        fprintf('\nError 133 - Input frame is unstable at minimum loads')
        P_opt = -133;
        P_vec = -133;
        err_flag = -133;
        return
    end
    
    
    %% C. Run global search
    fprintf('\nMessage - Beginning global search algorithm')
    prodS_val = prodS_zero;    % Tracks the determinant of the current optimal solution
    if MIN_OR_MAX == "MIN"
        opt_val = realmax;    % Tracks the current optimal value of sum(P)
    else
        opt_val = -realmax;    
    end
    opt_soln = [];            % Tracks the current optimal solution
    algo = 'interior-point';  % fmincon algorithm 'interior-point', 'sqp', or 'active-set'
    fminopt = optimoptions('fmincon', 'Algorithm', algo, 'Display', 'iter', ...
                        'MaxFunctionEvaluations', 100, 'MaxIterations', 10, ...
                        'FunctionTolerance', opt_tol*1000, 'ConstraintTolerance', 0*pi/2, 'StepTolerance', 0.05);
    hybrid_function = [];% {@fmincon,fminopt}; %[] or {@fmincon,fminopt}
    options = optimoptions('particleswarm','SwarmSize',swarm_size,'HybridFcn',hybrid_function... 
        ,'Display','iter','MaxIterations',max_iters,'MaxTime',3600); 
    
    % i. Full particle swarm algorithm performed on all loads
    n_vars = n_storeys*n_cols; % Let X be a list of pseudo-proportional load factors (0 corresponding to Pl and 1 corresponding to Pu)
    try
    [X, fval, exitflag, output] = particleswarm(...
        @(X) obj_func_decomp(X,[],Pl,Pu,E0,Ic,Rl,RL,RR,vc,wul,wlu,Hk,MIN_OR_MAX,K_tol,max_newton,B,bool_inel,Ac,Fy,n_multistart,...
                        bool_shear,shear_assumption,kappa,v0,bool_axial,Ab,Lb)...
        , n_vars, zeros(1,n_vars), ones(1,n_vars), options);
    catch
        fprintf('\nError 143 - Caught a fatal error during the global optimization')
    end
    opt_val = fval;
    [opt_soln,prodS_val,~,err_flag] = get_crit_loads_decomp_Brent(...
        X,[],Pl,Pu,E0,Ic,Rl,RL,RR,vc,wul,wlu,Hk,B,bool_inel,Ac,Fy,...
                    bool_shear,shear_assumption,kappa,v0,bool_axial,Ab,Lb);
    num_evals = num_evals + output.funccount;
    fprintf('\nMessage - Optimization completed with exit flag %.0d, with %.0d function evaluations and %.0d iterations', exitflag,output.funccount, output.iterations)
    fprintf('\nMessage - Current optimal solution has sumP = %.5e', opt_val)

    % ii. Only perform particle swarm on non-zero loads in the current optimal solution
    for repeat_PSO = 1:n_PSO-1
        non_zero_indices = find(abs(opt_soln-Pl)>opt_tol); % Returns of a vector of 1D indices of non-zero (or not at Pl) values in the optimal solution
        n_vars = length(non_zero_indices);                 % Reduced number of variables that will be optimized (note: sphere is dimension n-1)
        mask = abs(opt_soln-Pl)>opt_tol;                   % A matrix of booleans which will be multiplied by the transformed input variables
        if n_vars == 0
            fprintf('\nError 119 - Current optimal solution has no loads (all zero)')
            err_flag = -119;
        else
            try
            [X_mask, fval, exitflag, output] = particleswarm(...
                @(X_mask) obj_func_decomp(...
                X_mask,mask,Pl,Pu,E0,Ic,Rl,RL,RR,vc,wul,wlu,Hk,MIN_OR_MAX,K_tol,max_newton,B,bool_inel,Ac,Fy,n_multistart,...
                        bool_shear,shear_assumption,kappa,v0,bool_axial,Ab,Lb)...
                , n_vars, zeros(1,n_vars), ones(1,n_vars), options);
            catch
                fprintf('\nError 143 - Caught a fatal error during the global optimization')
            end
            if fval < opt_val
                opt_val = fval;
                [opt_soln,prodS_val,~,err_flag] = get_crit_loads_decomp_Brent(...
                    X_mask,mask,Pl,Pu,E0,Ic,Rl,RL,RR,vc,wul,wlu,Hk,B,bool_inel,Ac,Fy,...
                    bool_shear,shear_assumption,kappa,v0,bool_axial,Ab,Lb);
            end
            num_evals = num_evals + output.funccount;
            fprintf('\nMessage - Optimization completed with exit flag %.0d, with %.0d function evaluations and %.0d iterations', exitflag,output.funccount, output.iterations)
            fprintf('\nMessage - Current optimal solution has sumP = %.5e', opt_val)
        end
    end
    
    %% D. Output results
    P_opt = opt_val;    % Optimal value of sum P
    P_vec = opt_soln;   % Optimal combination of P
%     if abs(prodS_val) > K_tol
%         fprintf('\nError 117 - Final optimal solution does not satisfy prodS = 0')
%         err_flag = -117;
%     end
    fprintf('\nMessage - Final solution has sumP = %.5e and prodS = %.5e',opt_val,prodS_val)
    fprintf('\nMessage - Final solution exited with error flag: %.0f',err_flag)
    fprintf('\nMessage - Number of objective function evaluations: %.0f',num_evals)
    fprintf('\nMessage - Final solution is:')
    for i=1:n_storeys
        for j=1:n_cols
            fprintf('\nP(%.0f,%.0f) = %.5e',i,j,P_vec(i,j))
        end
    end
    
    % Calculate Nu for the optimal pseudo-proportional loading scenario
    if isempty(mask)
        mask = ones(size(Pl,1),size(Pl,2));
    end
    indices = find(mask==1);
    X_matrix = zeros(n_storeys,n_cols);
    for index = 1:length(indices)
        X_matrix(indices(index)) = X_mask(index);
    end
    vc_buckle = zeros(n_storeys,n_cols,n_storeys);
    for storey = 1:n_storeys
        vc_buckle(:,:,storey) = vc;
    end
    Nu = Iterate_Nu_Multi(E0,Ic,Hk,RL,RR,Rl,vc_buckle,... 
        ones(n_storeys,n_cols),ones(n_storeys,n_cols),1,...
        bool_inel,Ac,Fy,...
        1,Pl,Pu,X_mask,mask,1,wul,wlu,...
        bool_shear,shear_assumption,kappa,v0); 
    N_opt = zeros(n_storeys,n_cols);
    for i=1:n_storeys
        for j = 1:n_cols
            for ss = i:n_storeys
                N_opt(i,j) = N_opt(i,j) + P_vec(ss,j);
            end
        end
    end
    fprintf('\nMessage - Rotational buckling analysis:')
    for i=1:n_storeys
        for j=1:n_cols
            fprintf('\nNu(%.0f,%.0f) = %.5e',i,j,Nu(i,j))
        end
    end
    for i=1:n_storeys
        for j=1:n_cols
            fprintf('\n(N/Nu)(%.0f,%.0f) = %.5f',i,j,N_opt(i,j)/Nu(i,j))
        end
    end
    %fclose(fileID)

end

%% Objective function evaluation (Section 7.6.4)
% X = [n_storey * n_cols - 1] vector of input variables in spherical coordinates
% Pl, Pu = [n_storey * n_cols] lower and upper bound matrices for the loads
% X_mask = [n_vars] vector of spherical coordinates (for n_vars+1 loads)
function f = obj_func_decomp(...
    X_mask,mask,Pl,Pu,E0,Ic,Rl,RL,RR,vc,wul,wlu,Hk,MIN_OR_MAX,K_tol,max_newton,B, bool_inel,Ac,Fy,n_multistart,...
                        bool_shear,shear_assumption,kappa,v0,bool_axial,Ab,Lb)
    
    [soln,prodS_val,num_evals,err_flag] = ...
        get_crit_loads_decomp_Brent(...
        X_mask,mask,Pl,Pu,E0,Ic,Rl,RL,RR,vc,wul,wlu,Hk,B,bool_inel,Ac,Fy,...
                    bool_shear,shear_assumption,kappa,v0,bool_axial,Ab,Lb);
    f = sum(sum(soln));
    if MIN_OR_MAX == "MAX"
        f = -f;
    end
    
    % Reject the case if no solution found
    if err_flag == -131 || err_flag == -135 || err_flag == -140 || ...
            err_flag == -151 || err_flag == -152 || err_flag == -153 || err_flag == -160 || err_flag == -161 % no solution found
        f = Inf;
        return
    end
    
    % Reject the case if the determinant of the solution isn't zero
%     if abs(prodS_val) > K_tol
%         err_flag = -132;
%         f = Inf;
%         return
%         fprintf('\nError 132 - prodS = %.3e is not close enough to zero',prodS_val)
%     end
    
    % Reject the case if the bounds are not satisfied
    N = zeros(size(soln,1),size(soln,2));
    for i = 1:size(soln,1)
        for ii = i:size(soln,1)
            N(i,:) = N(i,:) + soln(ii,:);
        end
    end
    if sum(sum(soln>Pu)) > 0 || sum(sum(soln<Pl)) > 0
        if sum(sum(soln<Pl)) > 0
            fprintf('\nError 122a - Root P for which prodS = 0 does not fall within Pl bounds - rejecting case')
        end
        if sum(sum(soln>Pu)) > 0
            fprintf('\nError 122b - Root P for which prodS = 0 does not fall within Pu bounds - rejecting case')
        end
        err_flag = -122;
        f = Inf;
        return
    end
    
    %fprintf('\nObjective function evaluation completed with %.0d prodS evaluations and error flag %.0d', num_evals,err_flag)
end

%% Get the critical loads via Brent's method (Section 7.6.4)
function [soln,prodS_val,num_evals,err_flag] = get_crit_loads_decomp_Brent(...
    X,mask,Pl,Pu,E0,Ic,Rl,RL,RR,vc,wul,wlu,Hk,B,bool_inel,Ac,Fy,...
                        bool_shear,shear_assumption,kappa,v0,bool_axial,Ab,Lb)

    num_evals = 0;
    fprintf('\nBeginning root find algorithm for an instance of X')
    
    % 1. Check if X = zeros (will result in a no-search)
    if sum(sum(X)) == 0
        err_flag = -135;
        soln = -135;
        prodS_val = -135;
        fprintf('\nError 135 - Input X vector is all zeros - rejecting case')
        return;
    end

    % 2. Process mask = [] will let it search using all variables
    n_storeys = size(Pl,1);
    n_cols = size(Pl,2);
    if isempty(mask)
        mask = ones(size(Pl,1),size(Pl,2));
    end
    indices = find(mask==1);
    X_matrix = zeros(n_storeys,n_cols);
    for index = 1:length(indices)
        X_matrix(indices(index)) = X(index);
    end
    
    % 3. Get Nu for all columns given the current proportional loading scheme
    vc_buckle = zeros(n_storeys,n_cols,n_storeys);
    for storey = 1:n_storeys
        vc_buckle(:,:,storey) = vc;
    end
    Nu = Iterate_Nu_Multi(E0,Ic,Hk,RL,RR,Rl,vc_buckle,... 
        ones(n_storeys,n_cols),ones(n_storeys,n_cols),1,...
        bool_inel,Ac,Fy,...
        1,Pl,Pu,X,mask,1,wul,wlu,...
        bool_shear,shear_assumption,kappa,v0); % bool_decomp = 1 for decomposition method?? may need to validate
    for storey=1:n_storeys
        for col = 1:n_cols
            if Nu(storey,col)==-140
                err_flag = -140;
                soln = -140;
                prodS_val = -140;
                fprintf('\nError 140 - Nu(%.0d,%.0d) did not converge.',storey,col)
                return;
            end
            if Nu(storey,col)<0
                err_flag = -154;
                soln = -154;
                prodS_val = -154;
                fprintf('\nError 154 - Nu(i,j) returned an error or is not positive.')
                return;
            end
        end
    end
    
    % 4. Get the lowest value of lambda that will cause any Nu(i,j) to be reached in the frame for each storey
    storey_lu = realmax*ones(n_storeys,1);
    for storey = 1:n_storeys
        for col=1:n_cols
            sumPl_above = 0;
            sumPuPlX_above = 0;
            for ss = storey:n_storeys
                sumPl_above = sumPl_above + Pl(ss,col);
                sumPuPlX_above = sumPuPlX_above + (Pu(ss,col)-Pl(ss,col))*X_matrix(ss,col);
            end
            if sumPuPlX_above ~= 0
                lu = (Nu(storey,col) - sumPl_above)/sumPuPlX_above;
                if lu < storey_lu(storey)
                    storey_lu(storey) = lu;
                end
            end
        end
        fprintf('\nStorey %.0d, lambda_u  = %.5f',storey,storey_lu(storey))
        if storey_lu(storey) == realmax
            if sum(X_matrix(storey,:)) == 0
                fprintf('\nStorey %.0d is not loaded in current X case',storey)
            end
        end
    end
    min_lu = min(storey_lu);
    
    % 5. For each storey, find the value of lambda corresponding to Ki = 0 if it exists
    err_flag = 0;
    num_evals = 0;
    lr_min = min_lu; % Tracks the lowest lambda root (initial value is the lowest value of lu, which is an upper bound)
    for storey=1:n_storeys
        
        % 6. Ensure that at Ki > 0 at lambda = 0 (P = Pl). Otherwise the frame is already unstable at Pl.
        P = Pl;
        [~,Ki] = prodS(E0,Ic,Hk,P,RL,RR,Rl,vc,wul,wlu,...
            ones(n_storeys,n_cols),ones(n_storeys,n_cols),B,1,bool_inel,Ac,Fy,...
                        bool_shear,shear_assumption,kappa,v0,1,bool_axial,Ab,Lb);
        num_evals = num_evals + 1;
        if Ki(storey) < 0
            err_flag = -152;
            soln = -152;
            prodS_val = -152;
            fprintf('\nError 152 - Frame is unstable (Ki < 0) at Pl')
            return;
        end
        lambda_pos = 0;
        
        % 7. Rotational buckling load cannot possibly occur in the storey if its storey_lu is greater than the global min_lu
        % If this is the case, then check the value of Ki at lr_min.
        % If it is +ve then no root exists between lambda = 0 and lambda = lr_min.
        % If it is -ve then a root exists between lambda = 0 and lambda = lr_min.
        neg_found = 0;
        lambda_neg = realmax;
        if storey_lu(storey) > lr_min
            
            % 8. Check if K < 0 at lr_min - if so then root-find
            P = Pl;
            for index=1:length(indices)
                P(indices(index)) = Pl(indices(index)) + (Pu(indices(index))-Pl(indices(index)))*X(index)*lr_min;
            end
            [~,Ki] = prodS(E0,Ic,Hk,P,RL,RR,Rl,vc,wul,wlu,...
                ones(n_storeys,n_cols),ones(n_storeys,n_cols),B,1,bool_inel,Ac,Fy,...
                        bool_shear,shear_assumption,kappa,v0,1,bool_axial,Ab,Lb);
            num_evals = num_evals + 1;
            % If K > 0 or somewhere in the frame has yielded (isnan) then this storey does not govern stability. Skip to the next.
            if Ki(storey) > 0 || isnan(Ki(storey)) 
                fprintf('\nStorey %.0d found not to govern - skipping',storey)
                continue;
            else
                % Skip the following block since a point with K < 0 has already been found, but root find for K = 0
                neg_found = 1; 
                lambda_neg = lr_min;
            end
        end
        
        % 9. Check that at lambda_u, the value of Ki is negative. If not, decrease lambda by epsilon until a negative value is encountered
        % Theoretically, at lambda_u, Ki = -Inf, and increases with decreasing lambda
        eps = 1e-12;
        lu = storey_lu(storey);
        while lu > 0 && neg_found == 0
            P = Pl;
            for index=1:length(indices)
                P(indices(index)) = Pl(indices(index)) + (Pu(indices(index))-Pl(indices(index)))*X(index)*lu;
            end
            [~,Ki] = prodS(E0,Ic,Hk,P,RL,RR,Rl,vc,wul,wlu,...
                ones(n_storeys,n_cols),ones(n_storeys,n_cols),B,1,bool_inel,Ac,Fy,...
                        bool_shear,shear_assumption,kappa,v0,1,bool_axial,Ab,Lb);
            num_evals = num_evals + 1;
            if Ki(storey) < 0
                neg_found = 1;
                lambda_neg = lu;
                fprintf('\nStorey %.0d negative Ki found at lambda = %.5f',storey,lu)
                break;
            end
            lu = storey_lu(storey)-eps;
            eps = eps*2;
        end
        if neg_found == 0 % Point with negative Ki not found
            err_flag = -151;
            fprintf('\nWarning 151 - Could not find a point with Ki < 0 for root finding algorithm to start')
            
            % 10. Confirm that Nu has been reached in a column. If so, then assuming that Nu was calculated correctly
            % Ki = 0 should occur within storey_lu(storey) - eps or slightly above storey_lu(storey)
            % (it can be very hard to find this discontinuity because in some cases it occurs over a range of about 1E-11...)            
            lu = storey_lu(storey);
            P = Pl;
            for index=1:length(indices)
                P(indices(index)) = Pl(indices(index)) + (Pu(indices(index))-Pl(indices(index)))*X(index)*lu;
            end
            N = zeros(n_storeys,n_cols);
            for ss = 1:n_storeys
                for jj = 1:n_cols
                    for sss = ss:n_storeys
                        N(ss,jj) = N(ss,jj) + P(sss,jj);
                    end
                end
            end
            if min(min(abs(N-Nu))) < 0.1 % tolerance of N = Nu in Newtons
                % Since N = Nu and no root was found by decreasing lambda, assume that K = 0 occurs at the point N = Nu
                lr = storey_lu(storey);
                if lr < lr_min
                    lr_min = lr;
                end
                fprintf('\nStorey %.0d - N = Nu but no root found; assuming K = 0 at l = lu',storey)
                continue;
            else
                soln = -151;
                prodS_val = -151;
                fprintf('\nError 151 - Could not find a point with Ki < 0 for root finding algorithm to start and no column has N = Nu')
                return;
            end
        end
        
        % 11. Root find for lambda_r satisfying Ki = 0
        % By the intermediate value theorem, if a negative Ki and positive Ki are found, then the interval between these two points contains the root
        % Also, only a single root for Ki = 0 exists since we have decomposed into single storeys (Ki is monotonic in lambda)
        [lr,err_flag,ne] = Brents_Method(@(lambda) S_storey_prop(lambda,storey,Pl,Pu,X,mask,...
                    E0,Ic,Hk,RL,RR,Rl,vc,wul,wlu,ones(n_storeys,n_cols),ones(n_storeys,n_cols),...
                    B,bool_inel,Ac,Fy,...
                        bool_shear,shear_assumption,kappa,v0,1,bool_axial,Ab,Lb),...
                lambda_pos,lambda_neg,1e-12);
        num_evals = num_evals + ne;
        if err_flag == -160
            err_flag = -160;
            soln = -160;
            prodS_val = -160;
            fprintf('\nError 160 - Input to Brent''s Method did not contain points with opposite signs')
            return;
        elseif err_flag == -161
            err_flag = -161;
            soln = -161;
            prodS_val = -161;
            fprintf('\nError 161 - Brent''s method did not converge; function must have been discontinuous between input points')
            return;
        else
            fprintf('\nStorey %.0d Brent''s method succeeded with Ki = 0 at lambda = %.5f',storey,lr)
        end
        
        % Track the lowest found root of lambda - this will find the lowest lambda that satisfies Ki = 0 (i.e. the first buckling mode)
        if lr < lr_min
            lr_min = lr;
        end
    end
    
    % 12. Process results 
    if lr_min == realmax % No solution found
        err_flag = -153;
        soln = -153;
        prodS_val = -153;
        fprintf('\nError 153 - No solution found for getting critical loads via Brent''s method')
        return;
    else % Return solution - soln is the minimum loading configuration with prodS = 0 with the given proportional load factors
        soln = Pl;
        for index=1:length(indices)
            soln(indices(index)) = Pl(indices(index)) + (Pu(indices(index))-Pl(indices(index)))*X(index)*lr_min;
        end
        [prodS_val,~] = prodS(E0,Ic,Hk,soln,RL,RR,Rl,vc,wul,wlu,...
                ones(n_storeys,n_cols),ones(n_storeys,n_cols),B,1,bool_inel,Ac,Fy,...
                        bool_shear,shear_assumption,kappa,v0,1,bool_axial,Ab,Lb);
        num_evals = num_evals + 1;
        err_flag = 0;
        fprintf('\nRoot found with Brent''s method, prodS = %.3e, l = %.3f, P = %.3e',prodS_val,lr_min,sum(sum(soln)))
    end
    
end

%% Evaluate the decomposed lateral stiffness of the current storey based on the proportional load factor
function K_i = S_storey_prop(lambda,storey,Pl,Pu,X,mask,...
                    E0,Ic,Hk,RL,RR,Rl,vc,wul,wlu,su,sl,B,bool_inel,Ac,Fy,...
                        bool_shear,shear_assumption,kappa,v0,bool_buckle,bool_axial,Ab,Lb)
                
    % mask = [] will let it search using all variables
    n_storeys = size(Pl,1);
    n_cols = size(Pl,2);
    if isempty(mask)
        mask = ones(size(Pl,1),size(Pl,2));
    end
    indices = find(mask==1);
    X_matrix = zeros(n_storeys,n_cols);
    for index = 1:length(indices)
        X_matrix(indices(index)) = X(index);
    end

    P = Pl;
    for index=1:length(indices)
        P(indices(index)) = Pl(indices(index)) + (Pu(indices(index))-Pl(indices(index)))*X(index)*lambda;
    end
    [~,Ki] = prodS(E0,Ic,Hk,P,RL,RR,Rl,vc,wul,wlu,su,sl,B,1,bool_inel,Ac,Fy,...
                        bool_shear,shear_assumption,kappa,v0,bool_buckle,bool_axial,Ab,Lb);
    K_i = Ki(storey);
    
end

%% Use the decomposition method to calculate the product sum of storey lateral stiffnesses in the frame
function [prodK, Ki] = prodS(E0,Ic,Hk,P,RL,RR,Rl,vc,wul,wlu,su,sl,B,incl_cols,bool_inel,Ac,Fy,...
                        bool_shear,shear_assumption,kappa,v0,bool_buckle,bool_axial,Ab,Lb)

    % Bug fix - If any elements in P are very close to zero (within 10) there will be an error in calculating Ki - snap these values to zero 
    P(abs(P)<10.0) = 0;

    % Get end fixities and axial load coefficients of all columns
    [ru_m, rl_m, phi_m, omega_m] = Get_End_Fixities(E0,Ic,Hk,P,RL,RR,Rl,vc,wul,wlu,su,sl,incl_cols,bool_inel,Ac,Fy,0,...
                        bool_shear,shear_assumption,kappa,v0,bool_buckle);
    n_storeys = size(P,1);
    n_cols = size(P,2);
    
    % Get product sum of storey stiffnesses
    prodK = 1;
    Ki = zeros(n_storeys,1);
    for i = 1:n_storeys
        S_storey = 0;
        S_storey_sum_only = 0;
        
        for j = 1:n_cols
            
            % Calculate E 
            N = 0;
            for ss = i:1:n_storeys
                N = N + P(ss,j);
            end
            if bool_inel == 1
                E = Get_E_IE(N,Ac(i,j),E0,Fy);
                if bool_buckle == 1
                    e = 0.5*(E0/E-1);
                else
                    e = 0;
                end
                G = E0/(2*(1+v0) + 3*e);
            else
                E = E0;
                G = E0/(2*(1+v0));
            end
            EI = E*Ic(i,j);
            kAG = G*kappa*Ac(i,j);
            if bool_shear == 1
                eta = EI/Hk(i)^2/kAG;
            else %%bool_shear == 0
                eta = 0;
            end
            
            % Calculate beta
            ru = ru_m(i,j);
            rl = rl_m(i,j);
            phi = phi_m(i,j);
            omega = omega_m(i,j);
            wphi = phi*omega;
            beta_0 = (rl+ru+rl*ru)/(4-rl*ru);
            zeta_0 = 12*eta*beta_0;
            if phi > 0
                a1 = 3*(rl*(1-ru) + ru*(1-rl));
                a2 = 9*rl*ru - (1-rl)*(1-ru)*wphi^2;
                a3 = 18*rl*ru + a1*wphi^2;
                beta = 1/12*omega*phi^3 * ...
                        (a1*wphi*cos(wphi) + a2*sin(wphi)) /...
                        (18*rl*ru - a3*cos(wphi) + (a1-a2)*wphi*sin(wphi));
                zeta = (omega^2-1)*( 18*rl*ru*(1-cos(wphi)) + a1*wphi*sin(wphi) )/...
                        (18*rl*ru - a3*cos(wphi) + (a1-a2)*wphi*sin(wphi));
            else %% phi <= 0 -- no accounting for tension yet
                beta = beta_0;
                zeta = zeta_0;
            end
            
            % Calculate S
            Si = 12*EI/Hk(i)^3*beta/(1+zeta); % No Taylor series
            S_storey = S_storey + Si + B(i,j);
            S_storey_sum_only = S_storey_sum_only + Si + B(i,j);
            if bool_axial == 1 && j<n_cols % if approximating the effect of axial beam deformations
                kb = E0*Ab(i,j)/Lb(i,j);
                S_storey = 1/(1/S_storey + 1/kb);
            end
        end
        
        % Bug fix for axial deformations - cap Seq at sumS
        if bool_axial == 1 && S_storey > S_storey_sum_only
            S_storey = S_storey_sum_only;
        end

        %S_storey = S_storey + sum(B(i,:));  % moved inside the loop

        % Exact results can be obtained by only considering the bracing in the current storey - hence the following code is obsolete
        % This can be justified in that we don't consider the lateral stiffness of connecting columns either
        %if storey<n_storeys
        %    S_storey = S_storey - chi(storey)*sum(B(storey+1,:)); % Substract equivalent bracing stiffness from storey above
        %end
        %if storey>1
        %    S_storey = S_storey + 1/(chi(storey-1))*sum(B(storey-1,:)); % Add equivalent bracing stiffness from storey below
        %end

        % Product sum of storey lateral stiffness
        prodK = prodK * S_storey;
        Ki(i) = S_storey;
    end

end