% Calculates the critical sway load, Psw
%   E0 = scalar, modulus of steel
%   fy = scalar, yield stress of steel
%   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 = (n+1) vector of semi-bracing stiffness
%   bool_inel = scalar, "1" to consider tangent modulus, "0" otherwise
function Psw = Get_Psw_2001(E0,fy,Ac,Ic,Lc,ru0,rl0,Kb,lambda,bool_inel)
    % Verify input
    n = length(Ic)-1;
    if (length(Lc)~=n+1 || length(E0)~=1 || length(fy)~=1 || length(bool_inel)~=1||...
            length(ru0)~=n+1|| length(rl0)~=n+1 || length(Kb)~=n+1)
        Psw = -100;
        return;
    end
    
    % Initial guess
    ru = ru0;
    rl = rl0;
    K = Get_K(ru,rl);
    Ec = E0*ones(1,n+1);
    Pb = pi^2*Ec.*Ic./(K.*Lc).^2;
    Pmax = min(Pb./lambda);
    
    % Multistart
    n_trials = 50;
    lowest_Psw = nan;
    for trial = 1:n_trials

        P_val = Pmax*trial/n_trials;
        max_iter = 200;
        tol = 1e-3;
        step = -1;
        for num_iter = 1:max_iter

            % Get S
            P = P_val*lambda;
            if (bool_inel == 1) 
                Ec = Get_ET(P,Ac,E0,fy); %tangent modulus
                tau = Ec/E0;
                ru = ru0./((1-tau).*ru0 + tau);
                rl = rl0./((1-tau).*rl0 + tau);
            end
            if (find(Ec==0)>0)
                P_val = P_val*(1-0.05*num_iter/max_iter);
                continue;
            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 + sum(Kb);

            % Get dS
            P_step = (P_val + step)*lambda;
            if (bool_inel == 1)
                Ec = Get_ET(P_step,Ac,E0,fy); %tangent modulus
                tau = Ec/E0;
                ru = ru0./((1-tau).*ru0 + tau);
                rl = rl0./((1-tau).*rl0 + tau);
            end
            if (find(Ec==0)>0)
                P_val = P_val*(1-0.05*num_iter/max_iter);
                continue;
            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 + sum(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
        if num_iter >= max_iter
            Psw = -102;
        end
        
        % Update lowest Psw value found so far
        if Psw > 0 && isnan(lowest_Psw)
            lowest_Psw = Psw;
        end
        if Psw > 0 && Psw < lowest_Psw
            lowest_Psw = Psw;
        end
    end
    
    if isnan(lowest_Psw)
        Psw = -103;
    else
        Psw = lowest_Psw;
    end
 end

