anchor.diff:  
       c3 = cos(theta(3))
       s3 = sin(theta(3))
  -    
  +
  +    SMAT_DOT = get_angrate_dot_cosines( c1, c2, c3, s1, s2, s3, dtheta)
  +
  +  end function get_angrate_dot_array
  +
  +
  +  !-----------------------------------------------------------
  +  ! Returns the time derivative of angular rate matrix
  +  !
  +  ! 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.
  +  !
  +  ! Input : The sines and cosines of euler angles
  +  ! Output: SMAT_DOT
  +  !
  +  !-----------------------------------------------------------
  +  function get_angrate_dot_cosines( c1, c2, c3, s1, s2, s3, dtheta) &
  +       &result(SMAT_DOT)
  +
  +    real(dp)     :: dtheta(NUM_SPAT_DIM) 
  +    real(dp)     :: c1, c2, c3, s1, s2, s3
  +    type(matrix) :: SMAT_DOT
  +
       SMAT_DOT= matrix( (/ 0.0_dp,  0.0_dp,  -c2*dtheta(2), &
            & 0.0_dp, -s1*dtheta(1), c1*c2*dtheta(1)-s1*s2*dtheta(2),&
  -         & 0.0_dp, -c1*dtheta(1), -s1*c2*dtheta(1)-c1*s2*dtheta(2) /))
  -
  -  end function get_angrate_dot
  -
  +         & 0.0_dp, -c1*dtheta(1), -s1*c2*dtheta(1)-c1*s2*dtheta(2)/))
  +
  +  end function get_angrate_dot_cosines
   
     ! ********************************************************
     ! routine that prints the state and properties of the body
  
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: 
      ! reaction torque
      alpha%gr = zeroV !alpha%J*alpha%omega_dot !? check
  
      ! translational + rotational KE = 0.5mV^2 + 0.5*wJw
      ! (*) between two vectors does a dot product (see utils.f90)
      ! Assuming the body axis to be at the centre of mass of the body
  