% Modified February 4, 2019 - fix v ~ TM

% Calculates the end rotational stiffness given:
%   Eb = (n) vector of beam moduli
%   Ib = (n) vector of beam inertia
%   Lb = (n) vector of beam lengths
%   rL = (n) vector of beam left end fixity factors (assume known)
%   rR = (n) vector of beam right end fixity factors (assume known)
%   v = scalar ratio of beam end rotations
function [RL,RR] = Get_R_190204(Eb,Ib,Lb,rL,rR,v)
    % Verify input
    n = length(Eb);
    if (length(Ib) ~= n || ...
            length(rL) ~= n ||...
            length(rR) ~= n || length(v) ~= n)
        RL = -100;
        RR = -100;
        return;
    end
    
    RL = zeros(1,n);
    RR = zeros(1,n);
    for j=1:n
        % Left end rigidity, Nm/rad
        RL(j) = 6*rL(j)/(4-rL(j)*rR(j))*Eb(j)*Ib(j)/Lb(j)*(2+v(j)*rR(j));
        
        % Right end rigidity, Nm/rad
        RR(j) = 6*rR(j)/(4-rL(j)*rR(j))*Eb(j)*Ib(j)/Lb(j)*(2+1/v(j)*rL(j));
    end
    
end