% Calculates the critical sway load, Psw, at elevated temperatures
%   E0 = scalar, ambient modulus of elasticity (EC3)
%   fy = scalar, ambient yield stress (EC3)
%   Tc = {n+1) vector of column temperature
%   Ic = (n+1) vector of column inertia
%   Ac = (n+1) vector of column area
%   Lc = (n+1) vector of column lengths
%   ru0 = (n+1) vector of upper end fixity at ambient temperature
%   rl0 = (n+1) vector of lower end fixity at ambient temperature
%   Kb = (n+1) vector of total semi-bracing stiffness
%   lambda = (n+1) vector of proportional load factors
%   Tb = (n) vector of beam temperatures
%   Ab = (n) vector of beam axea
%   Lb = (n) vector of beam length
%   kappa = shear area factor
%   bool_shear = 1 to consider shear deformations
%   bool_axial = 1 to consider axial deformations
%   shear_assumption = 0 for Engesser (phi), 1 for Haringx (dydx)
function Psw = Get_Psw_shear_axial(E0,fy,Tc,Ic,Ac,Lc,ru0,rl0,Kb,lambda,Eb,Ab,Lb,kappa,bool_shear,bool_axial,shear_assumption)
    % Verify input
    n = length(Ic)-1;
    if (length(E0)~=1 || length(fy)~=1||length(Tc)~=n+1 ||length(Lc)~=n+1 || length(E0)~=1 || length(fy)~=1 || length(bool_shear)~=1|| length(bool_axial)~=1||...
            length(ru0)~=n+1|| length(rl0)~=n+1 || length(Kb)~=n+1 ||length(Eb) ~= n||length(Ab) ~= n||length(Lb) ~= n)
        Psw = -100;
        return;
    end
    
    kb = Eb.*Ab./Lb;
    nu = Get_nu(Tc);
    
    % 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;
            for i=1:n+1
                Ec(i) = Get_E(Tc(i),P(i),Ac(i),E0,fy); %tangent modulus
            end
            tau = Ec/E0;
            ru = ru0./((1-tau).*ru0 + tau);
            rl = rl0./((1-tau).*rl0 + tau);
            if (find(Ec==0)>0)
                P_val = P_val*(1-0.05*num_iter/max_iter);
                continue;
            end
            Gc = Ec./(2.*(1+nu)); % Assume tangent shear modulus via direct application of isotopic equation
            eta = Ec.*Ic./Lc.^2 ./ (kappa.*Ac.*Gc);
            phi = Get_phi(P,Ec,Ic,Lc);
            [beta_S,zeta] = Get_beta_shear(phi,ru,rl,eta,shear_assumption);
            beta_N = Get_beta(phi,ru,rl);
            if bool_shear == 0
                kc = 12.*Ec.*Ic./Lc.^3 .* beta_N;
            else
                kc = 12.*Ec.*Ic./Lc.^3 .* beta_S./(1+zeta);
            end
            if bool_axial == 0
                SS = sum(kc);
            else
                SS = Get_Seq_axial(kc,kb,Kb);
            end
            
            % Get dS
            P_step = (P_val + step)*lambda;
            for i=1:n+1
                Ec(i) = Get_E(Tc(i),P_step(i),Ac(i),E0,fy); %tangent modulus
            end
            tau = Ec/E0;
            ru = ru0./((1-tau).*ru0 + tau);
            rl = rl0./((1-tau).*rl0 + tau);
            if (find(Ec==0)>0)
                P_val = P_val*(1-0.05*num_iter/max_iter);
                continue;
            end
            Gc = Ec./(2.*(1+nu)); % Assume tangent shear modulus via direct application of isotopic equation
            phi = Get_phi(P_step,Ec,Ic,Lc);
            [beta_S,zeta] = Get_beta_shear(phi,ru,rl,eta,shear_assumption);
            beta_N = Get_beta(phi,ru,rl);
            if bool_shear == 0
                kc_step = 12.*Ec.*Ic./Lc.^3 .* beta_N;
            else
                kc_step = 12.*Ec.*Ic./Lc.^3 .* beta_S./(1+zeta);
            end
            if bool_axial == 0
                SS_step = sum(kc_step);
            else
                SS_step = Get_Seq_axial(kc_step,kb,Kb);
            end
            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;
        % Check for yield failure (the EC3 is elastic-perfectly plastic for
        % low temperatures) + assume fy = fy0,
        Py = min(Ac)*fy;
        P = Py*lambda-(1e-8);
        for i=1:n+1
            Ec(i) = Get_E(Tc(i),P(i),Ac(i),E0,fy); %tangent modulus
        end
        tau = Ec/E0;
        ru = ru0./((1-tau).*ru0 + tau);
        rl = rl0./((1-tau).*rl0 + tau);
        Gc = Ec./(2.*(1+nu)); % Assume tangent shear modulus via direct application of isotopic equation
        eta = Ec.*Ic./Lc.^2 ./ (kappa.*Ac.*Gc);
        phi = Get_phi(P,Ec,Ic,Lc);
        [beta_S,zeta] = Get_beta_shear(phi,ru,rl,eta,shear_assumption);
        beta_N = Get_beta(phi,ru,rl);
        if bool_shear == 0
            kc = 12.*Ec.*Ic./Lc.^3 .* beta_N;
        else
            kc = 12.*Ec.*Ic./Lc.^3 .* beta_S./(1+zeta);
        end
        if bool_axial == 0
            SS = sum(kc);
        else
            SS = Get_Seq_axial(kc,kb,Kb);
        end
        if SS > 0
            Psw = Py; % Yielding governs failure
        end
    else
        Psw = lowest_Psw;
    end
    
end

