% Calculates the second-order rotational stiffness contribution of the near end of the member
% N = axial force; +ve for compression, -ve for tension !!
% EI = EI of the member
% L = length of the member
% rN, rF = near and far end connection fixities of the member (1 for a fixed connection, 0 for a pinned connection)
% vji = tF/tN, ratio of the end connection rotations
% wji = (yF-yN)/(L*tN), where yF and yN are the far and near end deflections of the member
% bool_shear = toggle for shear deformations (Timoshenko theory)
% shear_assumption = 0 for shear acting on phi, 1 for shear acting on dy/dx
function RP = R_prime(rN,rF,N,EI,L,vji,wji,...
    bool_shear,kAG,shear_assumption) % shear 

    if bool_shear == 1
        if shear_assumption == 0
            omega = sqrt(1+N/kAG);
        else %shear_assumption == 1
            omega = 1/sqrt(1-N/kAG);
        end
    else
        omega = 1;
    end
    
    wphi = omega* L*(abs(N)/EI)^0.5;
    if wphi == 0
        if bool_shear == 1
            n = EI/L^2/kAG;
        else
            n = 0;
        end
        RP = 6*EI/L * rN/(4-rN*rF+12*n*(rN+rF+rN*rF)) * (2+6*n*rF+vji*rF*(1-6*n) - (2+rF)*wji);
    elseif N > 0 % compression
        a1 = 3*(rN*(1-rF)+rF*(1-rN));
        a2 = 9*rN*rF - (1-rN)*(1-rF)*wphi^2;
        a3 = 18*rN*rF + a1*wphi^2;
        DEN = 18*rN*rF + (a1-a2)*wphi*sin(wphi) - a3*cos(wphi);
        zeta = (omega^2-1)*(18*rN*rF*(1-cos(wphi)) + a1*wphi*sin(wphi))/DEN;
        G1 = (omega^2-1)*(3*rF*sin(wphi))/( (1-rF)*wphi^2*sin(wphi) + 3*rF*(sin(wphi)-wphi*cos(wphi))  );
        G2 = (omega^2-1)*sin(wphi)/(sin(wphi)-wphi);
        NUM = (1-rF)*wphi^2*sin(wphi)*(1+G1-wji) + 3*rF*(sin(wphi)*(1+G1-vji*(1+G2))-wphi*cos(wphi)*(1+G1-wji)+wphi*(vji*(1+G2)-wji));
        RP = 3*EI/L * rN*wphi * NUM/DEN / (1+zeta);
    else % N < 0 % tension - NOTE - shear not implemented!
        if (bool_shear == 1)
            disp('Shear not implemented for tension!')
        end
        b1 = 3*(rN*(1-rF)+rF*(1-rN));
        b2 = 9*rN*rF + (1-rN)*(1-rF)*wphi^2;
        b3 = 18*rN*rF - b1*wphi^2;
        DEN = 36*rN*rF*exp(wphi) - b3*(exp(2*wphi)+1) + wphi*(b2-b1)*(exp(2*wphi)-1);
        NUM = (1-rF)*wphi^2*(exp(2*wphi)-1)*(1-wji) - 3*rF*((exp(2*wphi)-1)*(1-vji) - wphi*(1+exp(2*wphi))*(1-wji) + 2*wphi*exp(wphi)*(vji-wji));
        RP = 3*EI/L * rN*wphi * NUM/DEN;
    end

end