% Calculates the rotational buckling limit, Pu, for a three-segment column
%   ru, rl = scalars, end fixity factors of column
%   tauu, taul = scalars, modulus modification factors of column
%   E0 = scalar, ambient modulus, N/m2
%   I = scalar, moment of inertia of column, m4
%   Lu, Lm, Ll = scalars, segment lengths in column, m
%   mu, mm, ml = scalars, degradation factors in segments of column

% Translated from VBA script Aug 17, 2018 ~TM
% Modiifed to account for tangent modulus Sep 18, 2018 ~TM
function [Pb] = Solve_Pb_3S(E0, fy0, I, A, Tu, Tm, Tl, con_l, con_u, Rl, Ru,...
                            rl0, ru0, Lu, Lm, Ll)

    % Calculate rotational resistance of top connection
    L = Lu + Lm + Ll;
    
    % Iteration settings
    maxIter = 30;
    numIter = 0;
    converge = "False";
    error_tol = 0.0001;
    
    % Iteration variables
    P_curr = 0;
    P_eps = 1;
    P_next = 0;
    e_curr = 0;
    de = 1;
    Pcr_min = 1E+100;
    found_a_Pcr = "False";
        
    % Initial values
    tl = 0.00001;
    
    % Check if it yields at zero load
    [~,~,~,~,~,~,~,~,yield] = Preprocess(Lu,Lm,Ll,Tu,Tm,Tl,0.0001,A,I,E0,fy0,con_l,con_u,Rl,Ru,rl0,ru0);
    if (yield == 1)
        Pb = 0;
        return
    end
    
    n_trials = 20;
    for k = 1:n_trials
        P_curr = 10 ^ (10*k/n_trials); % Initial guess
        numIter = 0;
        converge = "False";
        
        while (numIter < maxIter && converge == "False")
            
            [ru,rl,tauu,taul,mu,mm,ml,R_l,yield] = Preprocess(Lu,Lm,Ll,Tu,Tm,Tl,P_curr,A,I,E0,fy0,con_l,con_u,Rl,Ru,rl0,ru0);
            if (yield == 1) 
                P_curr = P_curr * 0.95;
                continue;
            end
            % Use Newton's method to zero out the error
            e_curr = GetError(P_curr, ru, rl, tauu, taul, E0, I, Lu, Lm, Ll, mu, mm, ml, L, R_l, tl);
            if (abs(e_curr) < error_tol) 
                % Found a critical load - is it the minimum one?
                converge = "True";
                found_a_Pcr = "True";
                if (P_curr > 0 && P_curr < Pcr_min) 
                    Pcr_min = P_curr; % Keep track of the lowest critical load
                end
                break;
            elseif e_curr == -100  %yl constraint not satisfied
                numIter = maxIter + 1;
                break;
            end
            
            % Compute derivative
            P_step = P_curr + P_eps;
            [ru,rl,tauu,taul,mu,mm,ml,R_l,yield] = Preprocess(Lu,Lm,Ll,Tu,Tm,Tl,P_step,A,I,E0,fy0,con_l,con_u,Rl,Ru,rl0,ru0);
            if (yield == 1) 
                P_curr = P_curr * 0.95;
                continue;
            end
            de = GetError(P_step, ru, rl, tauu, taul, E0, I, Lu, Lm, Ll, mu, mm, ml, L, R_l, tl)...
                     - e_curr;
            de_dP = de / P_eps;
            
            % Update P to the next iterate
            if de_dP == 0 
                numIter = maxIter + 1;
                break;
            end
            P_next = P_curr - e_curr / de_dP;
            if P_next < 0  % Avoid P < 0
                P_next = 0.00001;
            end
            if P_next > 1E+15 
                numIter = maxIter + 1;
                break;
            end
            P_curr = P_next;
            numIter = numIter + 1;
        end
        
    end
    
    if (found_a_Pcr == "True") 
        Pb = Pcr_min;
    else
        Pb = -101; %Error -101: No converged solutions found
    end 
    
end 

function [err] = GetError(P, ru, rl, tauu, taul, E0, I, Lu, Lm, Ll,...
    mu, mm, ml, L, R_l, tl)

        % Given values
        yl = 0; %Fix for elastic buckling!
        
        % Helper variables
        pu = L * (P / (E0 * mu * I)) ^ 0.5;
        pm = L * (P / (E0 * mm * I)) ^ 0.5;
        pl = L * (P / (E0 * ml * I)) ^ 0.5;
        Cu = cos(pu * Lu / L);
        Cm = cos(pm * Lm / L);
        Cl = cos(pl * Ll / L);
        Su = sin(pu * Lu / L);
        Sm = sin(pm * Lm / L);
        Sl = sin(pl * Ll / L);
        a1 = pl * pm * pu * (1 - Cl * Cm * Cu) + ...
             pu * pl ^ 2 * Sl * Sm * Cu + ...
             pm * pl ^ 2 * Sl * Cm * Su + ...
             pm ^ 2 * pl * Cl * Sm * Su;
        a2 = pl * pm * Cl * Cm * Su + ...
             pl * pu * Cl * Sm * Cu + ...
             pm * pu * Sl * Cm * Cu - ...
             pm ^ 2 * Sl * Sm * Su;
        a3 = pl * pm * pu + pm ^ 2 * Sl * Sm * Su - ...
             pl * pm * Cl * Cm * Su - ...
             pl * pu * Cl * Sm * Cu - ...
             pm * pu * Sl * Cm * Cu;
             
        % Shear force required to satisfy yl = 0 boundary condition
        S = -((R_l * a1 + L * P * a2) * tl + yl * P * pl * pm * pu) / (L * a3);
        
        % Based on P, calculate deflections, rotations and moments at stations
        Ml = R_l * tl;
        y1 = yl + Ml / P - Ml / P * Cl + (tl - S / P) / (pl / L) * Sl + S * Ll / P;
        t1 = Ml * pl / (P * L) * Sl + (tl - S / P) * Cl + S / P;
        M1 = Ml * Cl - (tl - S / P) * (P * L / pl) * Sl;
        y2 = y1 + M1 / P - M1 / P * Cm + (t1 - S / P) / (pm / L) * Sm + S * Lm / P;
        t2 = M1 * pm / (P * L) * Sm + (t1 - S / P) * Cm + S / P;
        M2 = M1 * Cm - (t1 - S / P) * (P * L / pm) * Sm;
        yu = y2 + M2 / P - M2 / P * Cu + (t2 - S / P) / (pu / L) * Su + S * Lu / P;
        tu = M2 * pu / (P * L) * Su + (t2 - S / P) * Cu + S / P;
        Mu = M2 * Cu - (t2 - S / P) * (P * L / pu) * Su;
        
        % Check if yu = 0 is satisfied - if not then I coded something wrong
        if (abs(yu) > 0.00000001) 
            err = -100; % Error -100: Coding error - y1 constraint not satisfied
        else
            R_u_crit = -Mu / tu;
            ru_crit = 1 / (1 + 3 * E0 * I / (tauu * L * R_u_crit));
            err = ru_crit - ru;
        end
end
%%
function [ru,rl,tauu,taul,mu,mm,ml,R_l,yield] = Preprocess(Lu,Lm,Ll,Tu,Tm,Tl,P,A,I,E0,fy0,con_l,con_u,Rl,Ru,rl0,ru0)
   % Begin preprocessing
    yield = 0;
    Eu = Get_E(Tu, P, A, E0, fy0);
    Em = Get_E(Tm, P, A, E0, fy0);
    El = Get_E(Tl, P, A, E0, fy0);
    mu = Eu / E0;
    mm = Em / E0;
    ml = El / E0;
    if (mu == 0 || mm == 0 || ml == 0)
        yield = 1;
        ru = -1; rl = -1; tauu = -1; taul = -1; R_l = -1;
        return
    end
    L = Lu+Lm+Ll;
    tauu = 1 / mu * (Lu / L) ^ 3 + 1 / mm * ((Lu + Lm) ^ 3 - Lu ^ 3) / L ^ 3 + 1 / ml * (1 - (Lu + Lm) ^ 3 / L ^ 3);
    taul = 1 / ml * (Ll / L) ^ 3 + 1 / mm * ((Ll + Lm) ^ 3 - Ll ^ 3) / L ^ 3 + 1 / mu * (1 - (Ll + Lm) ^ 3 / L ^ 3);
    if con_l == 1 
        if Rl == 0 
            rl = 0;
        else
            rl = 1 / (1 + 3 * E0 * I / L / Rl / taul);
        end
    else
        if rl0 == 0 
            rl = 0;
        else
            rl = 1 / (1 + (1 - rl0) / (taul * rl0));
        end
    end
    if con_u == 1 
        if Ru == 0 
            ru = 0;
        else
            ru = 1 / (1 + 3 * E0 * I / L / Ru / tauu);
        end 
    else
        if ru0 == 0 
            ru = 0;
        else
            ru = 1 / (1 + (1 - ru0) / (tauu * ru0));
        end 
    end
    if rl == 1 
        rl = 0.999999;
    end
    R_l = 3 * E0 * I * rl / (tauu * L * (1 - rl));
    % End preprocessing    
end