%% Returns all of the end fixity factors and axial load coefficients of each column in a multi-storey frame given the following inputs:
% E0 = scalar, ambient modulus of elasticity
% Ic = [n_storeys × n_cols] moment of inertia of columns
% Hk = [n_storeys × 1] storey heights
% P = [n_storeys × n_cols] applied gravity loads
% RL = [n_storeys × n_bays] beam left end connection rotational stiffness
% RR = [n_storeys × n_bays] beam right end connection rotational stiffness
% Rl = [1 × n_cols] ground connection rotational stiffness
% vc = [n_storeys × n_cols] ratio of upper to lower connection rotations for each column
% wul = [n_storeys × n_cols] ratio of chord rotation to lower end rotation for each column
% wlu = [n_storeys × n_cols] ratio of chord rotation to upper end rotation for each column
% su = [n_storeys × n_cols] upper end splice fixity of columns (1 for continuous splice)
% sl = [n_storeys × n_cols] lower end splice fixity of columns (1 for continuous splice)
% incl_cols = toggle 1 to include column contribution to end fixity, 0 to exclude
% bool_inel = toggle 1 to consider tangent modulus adjustment, 0 to exclude
% bool_decomp = X OBSOLETE MAY 2, 2019 SHOULD ALWAYS BE OFF toggle 1 to decompose all other storeys, 0 to remain
% bool_shear = toggle for shear deformations (Timoshenko theory)
% shear_assumption = 0 for shear acting on phi, 1 for shear acting on dy/dx
% Ac = [n_storeys × n_cols] area of columns
% v0 = scalar, Poisson's ratio
% bool_buckle = boolean, e = 0 if false, e = Becque(2016) if true

function [ru, rl, phi, omega] = Get_End_Fixities(E0,Ic,Hk,P,RL,RR,Rl,vc,wul,wlu,su,sl,incl_cols,...
    bool_inel,Ac,Fy,bool_decomp,... % inelasticity
    bool_shear,shear_assumption,kappa,v0,bool_buckle) % shear

    % Setup the solution
    n_storeys = size(P,1);
    n_cols = size(P,2);
    ru = -109*ones(n_storeys,n_cols);
    rl = -109*ones(n_storeys,n_cols);
    phi = -109*ones(n_storeys,n_cols);
    
    if bool_shear == 1
        omega = -109*ones(n_storeys,n_cols);
    else
        omega = ones(n_storeys,n_cols);
    end
    
    % Cycle through each column of the frame
    for storey = 1:n_storeys
        for column = 1:n_cols

            % I. 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);
            if bool_shear == 1
                if shear_assumption == 0
                    omega(storey,column) = sqrt(1+N/kAG);
                else %shear_assumption == 1
                    omega(storey,column) = 1/sqrt(1-N/kAG);
                end
            end
            phi(storey,column) = sqrt(N/EI)*Hk(storey);

            % II. Calculate end fixities
            % a) Lower end fixity, rl - move up
            RCl_next = Rl(column); % Tracks the lower rotational restraint of the current local storey ss to the next storey
            for ss = 1:storey

                % Determine local storey properties
                N_ss = 0;
                for sss = ss:1:n_storeys
                    N_ss = N_ss + P(sss,column);
                end
                if bool_inel == 1
                    E = Get_E_IE(N_ss,Ac(ss,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_ss = E*Ic(ss,column);
                kAG_ss = G*kappa*Ac(ss,column);

                % Determine lower end fixity factor of local storey
                Rl_ss = RCl_next; % Update to the current storey
                rl_ss = 1/(1+3*EI_ss/(Rl_ss*Hk(ss)));

                % Calculate the effective rotational restraint on the next storey
                if incl_cols == 1
                    if bool_decomp == 0 % seems to work for Nu
                        RCl_next = R_prime(su(ss,column),sl(ss,column),N_ss,EI_ss,Hk(ss),1/vc(ss,column),wlu(ss,column),...
                            bool_shear,kAG_ss,shear_assumption); 
                    else %bool_decomp == 1, seems to work for decomposition NOPE - MAY 2, 2019
                        RCl_next = R_prime(su(ss,column),rl_ss,N_ss,EI_ss,Hk(ss),1/vc(ss,column),wlu(ss,column),...
                            bool_shear,kAG_ss,shear_assumption);
                    end
                else
                    RCl_next = 0;
                end
                
                % Add the upper beam connection of the current local storey to the RCl of the next storey
                Rbu_ss = 0;
                if column > 1 % Right end of beams
                    Rbu_ss = Rbu_ss + RR(ss,column-1);
                end
                if column < n_cols % Left end of beams
                    Rbu_ss = Rbu_ss + RL(ss,column);
                end
                RCl_next = RCl_next + Rbu_ss;
                if ss == storey
                    rl(storey,column) = 1/(1+3*EI_ss/(Rl_ss*Hk(ss)));
                end
            end

            % b) Upper end fixity, ru - move down
            RCu_next = 0; % Tracks the upper rotational restraint of the current local storey ss to the next storey
            for ss = n_storeys:-1:storey

                % Add the upper beam connection contribution of the current local storey
                Ru_ss = RCu_next; % Update to the current storey
                Rbu_ss = 0;
                if column > 1 % Right end of beams
                    Rbu_ss = Rbu_ss + RR(ss,column-1);
                end
                if column < n_cols % Left end of beams
                    Rbu_ss = Rbu_ss + RL(ss,column);
                end
                Ru_ss = Ru_ss + Rbu_ss;

                % Determine local storey properties
                N_ss = 0;
                for sss = ss:1:n_storeys
                    N_ss = N_ss + P(sss,column);
                end
                if bool_inel == 1
                    E = Get_E_IE(N_ss,Ac(ss,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_ss = E*Ic(ss,column);
                kAG_ss = G*kappa*Ac(ss,column);

                % Determine upper end fixity factor of local storey
                ru_ss = 1/(1+3*EI_ss/(Ru_ss*Hk(ss)));

                % Calculate the effective rotational restraint on the next storey
                if incl_cols == 1
                    if bool_decomp == 0 % seems to work for Nu
                        RCu_next = R_prime(sl(ss,column),su(ss,column),N_ss,EI_ss,Hk(ss),vc(ss,column),wul(ss,column),...
                            bool_shear,kAG_ss,shear_assumption);
                    else %bool_decomp == 1, seems to work for decomposition NOPE - MAY 2, 2019
                        RCu_next = R_prime(sl(ss,column),ru_ss,N_ss,EI_ss,Hk(ss),vc(ss,column),wul(ss,column),...
                            bool_shear,kAG_ss,shear_assumption);
                    end
                else
                    RCu_next = 0;
                end
                
                if ss == storey
                    ru(storey,column) = 1/(1+3*EI_ss/(Ru_ss*Hk(ss)));
                end
            end
        end
    end
end