% Calculates the stiffness coefficient, beta, for columns, given:
%   phi = (n+1) vector of the load factors
%   ru = (n+1) vector of upper end fixity factors
%   rl = (n+1) vector of lower end fixity factors
%   eta = (n+1) vector, shear flexibility coefficient
%   shear_assumption = 0 for phi, 1 for dy/dx
function [beta,zeta] = Get_beta_shear(phi,ru,rl,eta,shear_assumption)
    % Verify input
    n = length(phi)-1;
    if (length(phi) ~= n+1 || length(ru)~=n+1 || length(rl)~=n+1 || length(eta)~=n+1)
        beta = -100;
        return;
    end
    
    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);
    wphi = zeros(1,n+1);
    zeta = zeros(1,n+1);
    if shear_assumption == 0
        omega = (1+eta.*phi.^2).^0.5;
    else
        omega = (1-eta.*phi^2).^(-0.5);
    end
    for i=1:n+1
        if (phi(i) > 0)
            wphi(i) = omega(i)*phi(i);
            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*omega(i)*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)));
            zeta(i)= (omega(i)^2-1) * ...
                    ( 18*rl(i)*ru(i)*(1-cos(wphi(i))) + a1(i)*wphi(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);
            zeta(i) = 12*eta(i)*(rl(i)+ru(i)+rl(i)*ru(i))/(4-rl(i)*ru(i));
        end
    end
end