% 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
function beta = Get_beta(phi,ru,rl)
    % Verify input
    n = length(phi)-1;
    if (length(phi) ~= n+1 || length(ru)~=n+1 || length(rl)~=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);
    for i=1:n+1
        if (phi(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))*phi(i)^2;
            a3(i) = 18*rl(i)*ru(i) + a1(i)*phi(i)^2;
            beta(i)= 1/12*phi(i)^3 * ...
                    (a1(i)*phi(i)*cos(phi(i)) + a2(i)*sin(phi(i))) /...
                    (18*rl(i)*ru(i) - a3(i)*cos(phi(i)) + (a1(i)-a2(i))*phi(i)*sin(phi(i)));
        else
            beta(i) = 1/gamma(i);
        end
    end
end