anchor.diff:      double precision :: omegasq
   
       integer          :: i, j, k
  -    double precision :: theta2, theta3
  +    double precision :: theta2, theta3, omega(3)
   
       phi(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) = -phi(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3))    
  +
  +    omega = get_omega()
   
       ! To ensure that this process is fully general for any orientation of
       ! the rotation axis with respect to the stellar configuration, 
  
positive.code:     use meth_params_module, only: NVAR, URHO
      use prob_params_module, only: center
      use probdata_module, only: c, d_vector, rloc
  
      implicit none
      
  
negative.code:     ! The rotational potential is then given by:
      ! (1/2) * | omega x r |^2 = (1/2) |omega|^2 |r|^2 sin^2(theta)
  
      theta2 = acos( dot_product(omega,d_vector(:,2)) / (sqrt(sum(omega**2)) * sqrt(sum(d_vector(:,2)**2))) )
      theta3 = acos( dot_product(omega,d_vector(:,3)) / (sqrt(sum(omega**2)) * sqrt(sum(d_vector(:,3)**2))) )
  
      do k = rloc(3,2), rloc(3,2) + 1
         do j = rloc(2,2), rloc(2,2) + 1
  