anchor.diff:      type(matrix) :: SMAT_DOT
   
       real(dp)     :: theta(NUM_SPAT_DIM), dtheta(NUM_SPAT_DIM)
  -    
  +
       ! convert vec to array
       theta = array(thetain); dtheta=array(dthetain);   
   
  
positive.code: 
    end function get_angrate_dot_cosines
  
    ! ********************************************************
    ! routine that prints the state and properties of the body
    ! ********************************************************
    subroutine print_body(alpha)
  
      use dispmodule
  
negative.code:     ! Assuming the body axis to be at the centre of mass of the body
      ! Inertial or body frame?
      alpha%KE = 0.5_dp*(mass * alpha%v *  alpha%v &
           & + alpha%omega*alpha%J* alpha%omega) ! + coupling term
  
      ! potential energy due to position
      ! include strain energy later
  