% Calculates the column stiffnesses Si and frame stiffness sumS, given:
%   P = (n+1) vector of total axial load
%   Ec = (n+1) vector of column moduli
%   Ic = (n+1) vector of column inertia
%   Lc = (n+1) vector of column lengths
%   rl = (n+1) vector of lower end fixity factors
%   ru = (n+1) vector of upper end fixity factors
%   bool_shear = toggle for shear deformations (Timoshenko theory)
%   shear_assumption = 0 for shear acting on phi, 1 for shear acting on dy/dx

% February 15, 2019 - updated to only calculate zeta if considering shear
function [sumS,Si] = Get_S_Shear(P,Ec,Ic,Lc,rl,ru,kAG,bool_shear,shear_assumption)
    % Verify input
    n = length(Ec)-1;
    if (length(Ic) ~= n+1 || length(Lc)~=n+1 || length(rl)~=n+1 || length(ru) ~= n+1 ||...
        length(kAG)~= n+1 || length(P) ~=n+1)
        sumS = -100;
        Si = -100;
        return;
    end
    
    Si = zeros(1,n+1);
    w = zeros(1,n+1);
    zeta = zeros(1,n+1);
    phi = zeros(1,n+1);
    wphi = zeros(1,n+1);
    beta = zeros(1,n+1);
    gamma = Get_gamma(ru,rl);
    a1 = zeros(1,n+1);
    a2 = zeros(1,n+1);
    a3 = zeros(1,n+1);
    for i=1:n+1
        
        % Calculate omega and phi
        if (bool_shear == 1)
            if (shear_assumption == 0)
                w(i) = sqrt(1+P(i)/kAG(i));
            else
                w(i) = 1/sqrt(1-P(i)/kAG(i));
            end
        else
            w(i) = 1;
        end
        phi(i) = sqrt(P(i)/(Ec(i)*Ic(i)))*Lc(i);
        wphi(i) = w(i)*phi;
        
        % Calculate beta without shear
        if (wphi(i) > 0)
            a1(i) = 3*(rl(i)*(1-ru(i)) + ru(i)*(1-rl(i)));
            a2(i) = 9*rl(i)*ru(i) - (1-rl(i))*(1-ru(i))*wphi(i)^2;
            a3(i) = 18*rl(i)*ru(i) + a1(i)*wphi(i)^2;
            beta(i)= 1/12*w*phi(i)^3 * ...
                    (a1(i)*wphi(i)*cos(wphi(i)) + a2(i)*sin(wphi(i))) /...
                    (18*rl(i)*ru(i) - a3(i)*cos(wphi(i)) + (a1(i)-a2(i))*wphi(i)*sin(wphi(i)));
        else
            beta(i) = 1/gamma(i);
        end
        if (bool_shear == 1)
            zeta(i) = (w(i)^2-1)*(a1(i)*wphi(i)*sin(wphi(i)) + 18*rl(i)*ru(i)*(1-cos(wphi(i)))) /...
                (18*rl(i)*ru(i) - a3(i)*cos(wphi(i)) + (a1(i)-a2(i))*wphi(i)*sin(wphi(i)));
        else
            zeta(i) = 0;
        end
        
        Si(i)= 12*Ec(i)*Ic(i)/Lc(i)^3*beta(i)*1/(1+zeta(i));
    end
    sumS = sum(Si);
end