% Calculates the stiffness coefficient, beta, for columns, given:
%   phi = (n+1) vector of the load factors (IN TENSION!)
%   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_tens(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);

    for i=1:n+1
        if (phi(i) > 0)
            if shear_assumption == 0
                omega = (1-eta(i)*phi(i)^2)^0.5;
            else
                omega = (1+eta(i)*phi(i)^2)^(-0.5);
            end
            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)*(exp(2*wphi(i))+1) + a2(i)*(exp(2*wphi(i))-1)) /...
                    (36*rl(i)*ru(i)*exp(wphi(i)) - a3(i)*(exp(2*wphi(i))+1) + (a2(i)-a1(i))*wphi(i)*(exp(2*wphi(i))-1));
            zeta(i)= (omega(i)^2-1) * ...
                    (36*rl(i)*ru(i)*exp(wphi(i)) - 18*rl(i)*ru(i)*(exp(2*wphi(i))+1) - a1(i)*wphi(i)*(exp(2*wphi(i))-1)  ) /...
                    (36*rl(i)*ru(i)*exp(wphi(i)) - a3(i)*(exp(2*wphi(i))+1) + (a2(i)-a1(i))*wphi(i)*(exp(2*wphi(i))-1));
        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