% Calculates the end rigidity contributions of the end connections, given:
%   E0 = scalar, ambient modulus
%   Ib = scalar, moment of inertia of beam (constant)
%   rN = scalar, near end fixity factor
%   rF = scalar, far end fixity factor
%   tauN = scalar, near end tau coefficient
%   tauF = scalar, far end tau coefficient
%   LN, LM, LF = scalar, near, middle and far end segment lengths
%   mN, mM, mF = scalar, near, middle, and far end modulus factors
%   v = scalar, ratio between end rotations

% January 11, 2019 - Fixed RR to take 1/v instead of v, and rN rF on
% denominator
function [RL, RR] = Get_R_3S(E0, Ib, rN, rF, tauN, tauF, LN,LM,LF,mN,mM,mF, v) 
    % Verify input
    if (length(E0) ~= 1 || length(Ib) ~= 1 || length(rN) ~= 1 || length(rF)~=1 ||...
            length(v) ~= 1 || length(tauN)~=1 || length(tauF)~=1 ||...
            length(LN) ~= 1 || length(LM) ~= 1 || length(LF)~=1 ||...
            length(mN) ~= 1 || length(mM) ~= 1|| length(mF)~= 1)
        R = -100;
        return;
    end
    L = LN+LM+LF;
    
    lb_NN = 1/L^3* ( LN^3*mM*mF + LM^3*mN*mF + LF^3*mN*mM + ... 
		  3*LN*LM^2*mN*mF + 3*LN*LF^2*mN*mM+3*LM*LN^2*mN*mF +...
		  3*LF*LM^2*mN*mM + 3*LF*LN^2*mN*mM+3*LM*LF^2*mM*mN +...
		  6*LN*LM*LF*mN*mM );
      
    lb_FF = 1/L^3* ( LN^3*mM*mF + LM^3*mN*mF + LF^3*mN*mM + ... 
		  3*LN*LM^2*mM*mF + 3*LN*LF^2*mM*mF+3*LM*LN^2*mM*mF+...
		  3*LF*LM^2*mN*mF + 3*LF*LN^2*mF*mM+3*LM*LF^2*mN*mF+...
		  6*LN*LM*LF*mF*mM );

    lb_NF = 1/L^3* ( LN^3*mM*mF + LM^3*mN*mF + LF^3*mN*mM+...
              3*LN*LM^2*mN*mF + 3*LN*LF^2*mN*mM+3*LM*LN^2*mM*mF+...
              3*LF*LM^2*mN*mF + 3*LF*LN^2*mF*mM+3*LM*LF^2*mN*mM+...
              6*LN*LM*LF*mN*mF );
          
    lb_FN = lb_NF;

    kN = 1/L^4*( ( LN^4*mM*mF + LM^4*mN*mF+LF^4*mN*mM) +...
              3*(LN*LM^3*mM*mF + LN*LF^3*mM*mF+LM*LF^3*mN*mF) +...
              4*(LM*LN^3*mM*mF + LF*LN^3*mM*mF + LF*LM^3*mN*mF) +...
              (LN*LM^3*mN*mF + LN*LF^3*mN*mM + LM*LF^3*mN*mM)+...
              6*(LN^2*LM^2*mM*mF+LN^2*LF^2*mM*mF + LM^2*LF^2*mN*mF)+...
              9*(LN*LM^2*LF*mM*mF + LN*LM*LF^2*mM*mF)+...
              12*(LN^2*LM*LF*mM*mF)+...
              3*(LN*LM^2*LF*mN*mF+LN*LM*LF^2*mN*mF) );

    kF = 1/L^4*( ( LN^4*mM*mF + LM^4*mN*mF+LF^4*mN*mM) +...
              3*(LM*LN^3*mN*mF + LF*LN^3*mN*mM+LF*LM^3*mN*mM) +...
              4*(LN*LM^3*mN*mF + LN*LF^3*mN*mM + LM*LF^3*mN*mM) +...
              (LM*LN^3*mM*mF + LF*LN^3*mM*mF + LF*LM^3*mN*mF)+...
              6*(LN^2*LM^2*mN*mF+LN^2*LF^2*mN*mM + LM^2*LF^2*mN*mM)+...
              9*(LN*LM^2*LF*mN*mM + LF*LM*LN^2*mM*mN)+...
              12*(LF^2*LM*LN*mM*mN)+...
              3*(LN*LM^2*LF*mN*mF+LF*LM*LN^2*mN*mF) );

    lb_A = mN*mM*mF*tauN*tauF;
    lb_B = 4*tauF *(kN - tauN*mN*mM*mF);
    lb_C = 4*tauN *(kF - tauF*mN*mM*mF);
    lb_D = 4*( tauF*kN + tauN*kF - tauN*tauF*mN*mM*mF  ) -...
               3/L^4*(LN^4*mM*mF/mN+LM^4*mN*mF/mM+LF^4*mN*mM/mF+...
               4*LN*LM^3*mF + 4*LN*LF^3*mM + 4*LM*LN^3*mF +...
               4*LM*LF^3*mN + 4*LF*LN^3*mM + 4*LF*LM^3*mN +...
               6*LN^2*LM^2*mF + 6*LN^2*LF^2*mM + 6*LM^2*LF^2*mN +...
               12*LN^2*LM*LF*mM + 12*LN*LM^2*LF*mM + 12*LN*LM*LF^2*mM );
           
    gN = tauF*mN*mM*mF*(1-rF) + lb_NN*rF;
    gF = tauN*mN*mM*mF*(1-rN) + lb_FF*rN; 
    
    RL = 6*E0*Ib*rN/L * (2*gN+v*lb_NF*rF) / (4*lb_A+rN*lb_B + rF*lb_C - rN*rF*lb_D);
    RR = 6*E0*Ib*rF/L * (2*gF+1/v*lb_FN*rN) / (4*lb_A+rN*lb_B + rF*lb_C - rF*rN*lb_D);
end