% Calculates the critical sway load, Psw
%   E0, fy = scalars, elastic modulus and yield stress
%   Ac = (n+1) vector of column area
%   Ic = (n+1) vector of column inertia
%   Lc = (n+1) vector of column lengths
%   ru = (n+1) vector of upper end fixity
%   rl = (n+1) vector of lower end fixity
%   Kb = scalar, total semi-bracing stiffness
%   P0 = scalar, initial guess load level
function Psw = Get_Psw(E0,fy0,Ac,Ic,Lc,ru,rl,Kb,lambda,P0)
    % Verify input
    n = length(Ac)-1;
    if (length(Ic) ~= n+1 || length(Lc)~=n+1 || length(Ac)~=n+1 ||...
            length(ru)~=n+1|| length(rl)~=n+1 || length(Kb)~=1)
        Psw = -100;
        return;
    end
    
    % Initial guess
    P_val = P0;
    
    max_iter = 100;
    tol = 1e-3;
    step = 1;
    for num_iter = 1:max_iter
        
        % Get S
        P = P_val*lambda;
        Ec = E0*ones(1,n+1);
        for i=1:n+1 % Inelastic modify
            Ec(i) = Get_E_IE(P(i),Ac(i),E0,fy0);
        end
        phi = Get_phi(P,Ec,Ic,Lc);
        beta = Get_beta(phi,ru,rl);
        [sumS,~] = Get_S_2001(P,Ec,Ic,Lc,beta);
        SS = sumS + Kb;
        
        % Get dS
        P_step = (P_val + step)*lambda;
        Ec = E0*ones(1,n+1);
        for i=1:n+1 % Inelastic modify
            Ec(i) = Get_E_IE(P_step(i),Ac(i),E0,fy0);
        end
        phi = Get_phi(P_step,Ec,Ic,Lc);
        beta = Get_beta(phi,ru,rl);
        [sumS_step,~] = Get_S_2001(P_step,Ec,Ic,Lc,beta);
        SS_step = sumS_step + Kb;
        dS = SS_step - SS;
        dSdP = dS/step;
        
        % Process and iterate
        P_next = P_val - SS/dSdP;
        if abs(P_next - P_val)<tol
            Psw = P_next;
            break;
        end
        if num_iter == max_iter
            Psw = -101;
            break;
        end
        P_val = P_next;
    end
end