anchor.diff:     cc(3) = em*sm3
   
   END SUBROUTINE BD_CrvExtractCrv
  +!------------------------------------------------------------------------------
  +!> This subroutine determines whether R is a valid rotation matrix and, if not,
  +!> replaces it with the closest orthongonal rotation matrix
  +SUBROUTINE BD_CheckRotMat(R, Rout, ErrStat, ErrMsg)
  +
  +   REAL(BDKi),       INTENT(IN   )  :: R(3,3)        !< Rotation Matrix input
  +   REAL(BDKi),       INTENT(  OUT)  :: Rout(3,3)     !< Rotation Matrix output
  +   INTEGER(IntKi),   INTENT(  OUT)  :: ErrStat       !< Error status of the operation
  +   CHARACTER(*),     INTENT(  OUT)  :: ErrMsg        !< Error message if ErrStat /= ErrID_None
  +
  +   !Local variables
  +   REAL(BDKi)                  :: tempmat(3,3)
  +   REAL(BDKi)                  :: S(3)         !mjs--these three are the SVD matrices (S is actually a vector)
  +   REAL(BDKi)                  :: U(3,3)
  +   REAL(BDKi)                  :: VT(3,3)
  +   INTEGER(IntKi)              :: lwork = 27   !mjs--from LAPACK: dgesvd doc page, lwork >= MAX(1,3*MIN(M,N) + MAX(M,N),5*MIN(M,N))
  +   REAL(BDKi), ALLOCATABLE     :: work(:)          ! where M x N is dimension of R, and lwork is the dimension of work
  +   REAL(BDKi)                  :: Rr(3,3)      !mjs--correccted rotation matrix
  +   LOGICAL                     :: ortho        !mjs--logical value indicating whether R is orthogonal
  +   INTEGER                     :: i            ! loop variable/case indicator
  +
  +   INTEGER(IntKi)              :: ErrStat2 ! Temporary Error status
  +   CHARACTER(ErrMsgLen)        :: ErrMsg2  ! Temporary Error message
  +   character(*), parameter     :: RoutineName = 'BD_CheckRotMat'
  +
  +   ! Initialize ErrStat
  +   ErrStat = ErrID_None
  +   ErrMsg  = ""
  +
  +   ! mjs--Start by determining if R is a valid rotation matrix using the properties:
  +      ! 1) the eigenvalues of an orthogonal matrix have complex modulus == 1
  +      ! 2) a valid rotation matrix must have determinant == +1
  +      ! i.e., the singular values == 1
  +
  +   ! mjs--If \f$ \underline{\underline{R}} \f$ is not a valid roatation tensor,
  +      ! compute \f$ \underline{\underline{Rr}} \f$, the nearest orthogonal tensor
  +      ! to \f$ \underline{\underline{R}} \f$
  +      ! this is done via computing SVD for \f$ \underline{\underline{R}} = USV^T \f$
  +      ! and setting \f$ \underline{\underline{Rr}} = UV^T \f$
  +      ! otherwise, assign \f$ \underline{\underline{Rr}}  = \underline{\underline{R}} \f$
  +
  +   allocate (work(lwork))
  +   tempmat = R !mjs--need this to handle inout nature of input for LAPACK_dgesvd
  +   call LAPACK_dgesvd('A', 'A', 3, 3, tempmat, S, U, VT, work, lwork, ErrStat2, ErrMsg2)
  +   CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
  +   if (ErrStat >= AbortErrLev) return
  +
  +   do i = 1, 3
  +      ortho = equalrealnos(S(i), 1.0_BDKi)
  +      if (.not. ortho) exit
  +   end do
  +
  +   ! mjs--if an invalid rotatation matrix is passed, set lowest error status, so as not to abort,
  +      ! but still inform user
  +   if (.not. ortho) then
  +      ErrStat = 1
  +      ErrMsg = 'Passed invalid rotation matrix--fixing via SVD'
  +      CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
  +      if (ErrStat >= AbortErrLev) return
  +      Rout = matmul(U, VT)
  +   else
  +      Rout = R
  +   end if
  +
  +END SUBROUTINE BD_CheckRotMat
   !------------------------------------------------------------------------------
   !> This subroutine generates n-point gauss-legendre quadrature points and weights
   !!
  
positive.code:    !! _Note:_ The above equation does not match what is in the March 2016 BD manual (it is the transpose).
     !! It does, however, match equation 5.17 in the March 2016 "BeamDyn User's Guide and Theory Manual"
  
     ! mjs--Start by determining if R is a valid rotation matrix using the properties:
        ! 1) the eigenvalues of an orthogonal matrix have complex modulus == 1
        ! 2) a valid rotation matrix must have determinant == +1
        ! i.e., the singular values == 1
  
     ! mjs--If \f$ \underline{\underline{R}} \f$ is not a valid roatation tensor,
        ! compute \f$ \underline{\underline{Rr}} \f$, the nearest orthogonal tensor
        ! to \f$ \underline{\underline{R}} \f$
        ! this is done via computing SVD for \f$ \underline{\underline{R}} = USV^T \f$
        ! and setting \f$ \underline{\underline{Rr}} = UV^T \f$
        ! otherwise, assign \f$ \underline{\underline{Rr}}  = \underline{\underline{R}} \f$,
  
  
     allocate (work(lwork))
     tempmat = R !mjs--need this to handle inout nature of input for LAPACK_gesvd
     call LAPACK_dgesvd('A', 'A', 3, 3, tempmat, S, U, VT, work, lwork, ErrStat2, ErrMsg2)
     CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
     if (ErrStat >= AbortErrLev) return
  
     do i = 1, 3
        ortho = equalrealnos(S(i), 1.0_BDKi)
        if (.not. ortho) exit
     end do
  
     ! mjs--if an invalid rotatation matrix is passed, set lowest error status, so as not to abort,
        ! but still inform user
     if (.not. ortho) then
        ErrStat = 1
        ErrMsg = 'Passed invalid rotation matrix--fixing via SVD'
        CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
        if (ErrStat >= AbortErrLev) return
        Rr = matmul(U, VT)
     else
        Rr = R
     end if
  
        !> Find the trace of the matrix
        !! \f$  T = t_{r0} \left( 3 c_0 - c_1^2 - c_2^2 - c_3^2 \right) \f$
     pivot = Rr(1,1) + Rr(2,2) + Rr(3,3)
  
        !> - Condition 1: \f$ \underline{\underline{R(3,3)}} > T \f$
        !! This implies that \f$ c_3^2 > c_0^2 \f$
     IF (Rr(3,3) .GT. pivot) THEN
        !pivot = Rr(3,3)
        !ipivot = 3
        sm0  = Rr(2,1) - Rr(1,2)                           !  4 c_0 c_3 t_{r0}
        sm1  = Rr(1,3) + Rr(3,1)                           !  4 c_1 c_3 t_{r0}
        sm2  = Rr(2,3) + Rr(3,2)                           !  4 c_2 c_3 t_{r0}
  
negative.code: 
        !> - Condition 3: \f$ \underline{\underline{R(1,1)}} > T \f$
        !! This implies that \f$ c_1^2 > c_0^2 \f$
     ELSEIF (i == 1) THEN
        sm0  = Rr(3,2) - Rr(2,3)                           !  4 c_0 c_1 t_{r0}
        sm1  = 1.0_BDKi + Rr(1,1) - Rr(2,2) - Rr(3,3)      !  4 c_1^2   t_{r0}
        sm2  = Rr(1,2) + Rr(2,1)                           !  4 c_1 c_2 t_{r0}
  