anchor.diff:      alpha%theta_dot = vector(qdot(4:6))
       alpha%v_dot     = vector(qdot(7:9))
       alpha%omega_dot = vector(qdot(10:12))
  -    
  +
       if (ELASTIC) then
          ! set the elastic state variables in the body
  -!       alpha%qs        = vector(q(13:15))    ! IS_QS:IE_QE
  -!       alpha%qs_dot    = vector(qdot(13:15)) ! IS_QS:IE_QE
  +       !       alpha%qs        = vector(q(13:15))    ! IS_QS:IE_QE
  +       !       alpha%qs_dot    = vector(qdot(13:15)) ! IS_QS:IE_QE
       end if
   
     end subroutine set_body_state
  
positive.code:   end function get_angrate_from_cosines
  
    !-----------------------------------------------------------
    ! we use the 3-2-1 Euler angles, the S matrix does not depend
    ! on c3/s3, however we keep these as inputs in case we ever  
    ! want to change the Euler sequence.
    !-----------------------------------------------------------
    function get_angrate_dot(theta, dtheta) result(SMAT_DOT)
  
      real(dp)     :: theta(NUM_SPAT_DIM), dtheta(NUM_SPAT_DIM)
      type(matrix) :: SMAT_DOT
      real(dp)     :: c1, c2, c3, s1, s2, s3
  
      ! Compute the sin/cos of the Euler angles
      c1 = cos(theta(1))
      s1 = sin(theta(1))
  
  
negative.code:   !
    ! Input: theta as an array
    ! Output: CMAT of type MATRIX
    ! 
    ! Ref: Section 2.2 Eq. 18/19 Hughes
    !*******************************************************************!
    function get_rotation_from_angles_array(theta) result(CMAT)
  
  