% Calculates the axial restraint force, H, and stiffness, ki, on a column, given:
%   Ll, Lm, Lu = (n+1) vector, lengths of column segments, m
%   Tl, Tm, Tu = (n+1) vector, temperatures of column segments, °C
%   ml, mm, mu = (n+1) vector, degradation factors of column segments
%   T0 = scalar, ambient temperature, °C
%   E0 = scalar, ambient modulus, N/m2
%   G = (n+1) vector, gravity loads on columns, N
%   Sb = (n) vector, beam vertical stiffnesses, N/m
%   Ac = (n+1) vector, column sectional areas, m^2

% Updated on Dec 18, 2018 - accounts for tangent modulus E(T)
function [ki, H] = GetAxialRestraint_181218(Ll,Lm,Lu,Tl,Tm,Tu,T0,E0,fy0,G,Sb,Ac)

    %Verify input
    n = length(Sb);
    if (length(T0)~=1 || length(E0) ~=1 || length(fy0) ~=1 || length(G) ~= (n+1) || ...
            length(Ac) ~= (n+1)|| ...
            length(Ll) ~= n+1 || length(Lm) ~= n+1 || length(Lu)~= n+1 ||...
            length(Tl) ~= n+1 || length(Tm) ~= n+1 || length(Tu)~= n+1)
        ki = -110;
        H = -110;
        return;
    end
    
    ki = zeros(1,n+1);
    Psi = zeros(1,n+1);
    %Eeq = zeros(1,n+1); % obsolete term
    H = zeros(1,n+1);
    P = zeros(1,n+1);
    L = Ll+Lm+Lu;
    for i=1:(n+1)
        
        % Get beam stiffness on columns
        if (i ~= 1)
            ki(i) = ki(i) + Sb(i-1);
        end
        if (i ~= n+1)
            ki(i) = ki(i) + Sb(i);
        end
        
        Psi(i) =(...
                    Ll(i)*(0.002*(Tl(i)^2-T0^2)+12*(Tl(i)-T0))+...
                    Lm(i)*(0.002*(Tm(i)^2-T0^2)+12*(Tm(i)-T0))+...
                    Lu(i)*(0.002*(Tu(i)^2-T0^2)+12*(Tu(i)-T0))...
                )*10^-6;
        %Eeq(i) = 1/(1/E0*(  Ll(i)/L(i)/ml(i)+... % obsolete term
        %                    Lm(i)/L(i)/mm(i)+...
        %                    Lu(i)/L(i)/mu(i)-...
        %                    1)...
        %                 );
        
        % Iterate to determine the axial load
        P(i) = G(i); % First iteration
        max_iters = 10;
        for n_iters = 1:max_iters
            psi_M = Lu(i) * Get_Strain(Tu(i), P(i), Ac(i), E0, fy0) +...
                    Lm(i) * Get_Strain(Tm(i), P(i), Ac(i), E0, fy0) +...
                    Ll(i) * Get_Strain(Tl(i), P(i), Ac(i), E0, fy0);
            P(i) = (G(i) + ki(i) * (Psi(i) + psi_M)) / (1 + ki(i) * L(i) / (E0 * Ac(i)));
        end
        
        H(i) = P(i) - G(i);
        %H(i) = (G(i) + ki(i)*Psi(i))/... % old formula
        %       ( 1+ki(i)*L(i)/Ac(i)/Eeq(i) )...
        %       - G(i);
        
        % Special case, Eeq = infinity % obsolete term
        %if (sum(Tl(i)+Tm(i)+Tu(i))==0) 
        %    H(i) = (G(i) + ki(i)*Psi(i))/...
        %       ( 1+0 )...
        %       - G(i);
        %end
        
        % Special case, Eeq = 0 % obsolete term
        %if (ml(i) == 0 || mm(i) == 0 || mu(i) == 0)
        %    H(i) = -G(i);
        %end
    end

end