% Calculates chi 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 chi = Get_chi(phi,ru,rl)
    % Verify input
    n = length(phi)-1;
    if (length(phi) ~= n+1 || length(ru)~=n+1 || length(rl)~=n+1)
        chi = -100;
        return;
    end
    
    chi = zeros(1,n+1);
    a1 = zeros(1,n+1);
    a2 = zeros(1,n+1);
    a3 = zeros(1,n+1);
    for i=1:n+1
        if (phi > 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;
        chi(i)= 3*phi(i)^3*pi*sin(phi(i))*(ru(i)-rl(i)) / (pi^2-phi(i)^2) /...
                (18*rl(i)*ru(i) - a3(i)*cos(phi(i)) + (a1(i)-a2(i))*phi(i)*sin(phi(i)));
    
        end
    end
end