% Calculates the column stiffnesses Si and frame stiffness sumS, given:
%   P = (n+1) vector of total axial load, N
%   Ll,Lm,Lu = (n+1) vectors of column segment lengths, m
%   ml,mm,mu = (n+1) vectors of column segment degradation factors
%   E0 = ambient modulus, N/m2
%   Ic = (n+1) vector of column inertia, m4
%   rl,ru = (n+1) vectors of column end fixity factors
%   taul, tauu = (n+1) vectors of column modulus modification factors
%   D0 = scalar or (n+1) vector of out-of-plumbness
%   Kb = scalar lateral brace stiffness

% Updated January 24, 2019 - Removed D0/D
function [phi,pl,pm,pu,psi,Si,sumS,D] = ...
    Get_S_3S_190124(P,Ll,Lm,Lu,ml,mm,mu,E0,Ic,rl,ru,taul,tauu,D0,Kb)
    
    % Verify input
    n = length(P)-1;
    if (length(Ic) ~= n+1 || length(Ll)~=n+1|| length(Lm)~=n+1|| length(Lu)~=n+1 || ...
            length(ml)~=n+1|| length(mm)~=n+1|| length(mu)~=n+1 ||...
            length(E0)~=1 || length(rl)~=n+1|| length(ru)~=n+1 ||...
            length(taul)~=n+1|| length(tauu)~=n+1)
        sumS = -100;
        Si = -100;
        phi = -100;pl = -100;pm = -100;pu = -100;
        psi = -100;
        return;
    end
    if (length(D0) == 1)
        D0 = ones(1,n+1)*D0;
    end 
    
    % Pre-processing
    L = Ll+Lm+Lu;
    %D = 1;
    
    % Primary load factors
    phi = zeros(1,n+1);
    pl = zeros(1,n+1); pm = zeros(1,n+1); pu = zeros(1,n+1);
    for i=1:n+1
        phi(i) = (P(i)/(E0*Ic(i)))^0.5 * L(i);
        pl(i) = phi(i) / ml(i)^0.5;
        pm(i) = phi(i) / mm(i)^0.5;
        pu(i) = phi(i) / mu(i)^0.5;
    end
    
    % Intermediate load variables
    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=pm.^2.*Sl.*Sm.*Su-pl.*pm.*Cl.*Cm.*Su-pm.*pu.*Sl.*Cm.*Cu-pl.*pu.*Cl.*Sm.*Cu;
    a2=pl.^2.*pm.*pu.*Sl.*Cm.*Cu+pl.*pm.^2.*pu.*Cl.*Sm.*Cu+pl.*pm.*pu.^2.*Cl.*Cm.*Su-pl.^2.*pu.^2.*Sl.*Sm.*Su;
    a3=pl.^2.*pm.*Sl.*Cm.*Su+pl.*pm.^2.*Cl.*Sm.*Su+pl.^2.*pu.*Sl.*Sm.*Cu-pl.*pm.*pu.*Cl.*Cm.*Cu;
    a4=pl.*pu.^2.*Cl.*Sm.*Su+pm.*pu.^2.*Sl.*Cm.*Su+pm.^2.*pu.*Sl.*Sm.*Cu-pl.*pm.*pu.*Cl.*Cm.*Cu;
    a5=a3+a4+2*pl.*pm.*pu;
     
    % Psi and lateral stiffness
    psi = zeros(1,n+1);
    Si = zeros(1,n+1);
    N = 0;
    for i=1:n+1
        psi(i) = (9*ru(i)*rl(i)*a5(i) -...
                    3*phi(i)^2*(rl(i)*(1-ru(i))*tauu(i)+ru(i)*(1-rl(i))*taul(i))*a1(i) ) /...
            ( phi(i)^4*taul(i)*tauu(i)*(1-rl(i))*(1-ru(i))*a1(i) + ...
                9*ru(i)*rl(i)*a2(i) - ...
                3*phi(i)^2*(rl(i)*(1-ru(i))*tauu(i)*a3(i)+ru(i)*(1-rl(i))*taul(i)*a4(i)) );
        Si(i)= phi(i)^2*E0*Ic(i)/L(i)^3 * ...
               ( 1/(psi(i)-1) ...%- D0(i)/D );
               );
        N = N + P(i)*D0(i)/L(i);
    end
    
    % Frame lateral stiffness
    sumS = sum(Si)+Kb;
    D = N/sumS;
    
end