anchor.diff:    use fundamental_constants_module, only: Gconst, c_light, parsec
     use probdata_module, only: gw_dist, axis_1, axis_2, axis_3
     use meth_params_module, only: rot_axis
  +  use prob_params_module, only: dim
   
     implicit none
   
  
positive.code: 
       h(:,:) = h(:,:) * TWO * Gconst / (c_light**4 * r)
  
       ! If rot_axis == 3, then h_+ = h_{11} = -h_{22} and h_x = h_{12} = h_{21}.
       ! Analogous statements hold along the other axes.
  
       ! We are adding here so that this calculation makes sense on multiple levels.
  
       if (dir .eq. axis_1) then
  
          h_plus_1  = h_plus_1  + h(axis_2,axis_2)
          h_cross_1 = h_cross_1 + h(axis_2,axis_3)
  
       else if (dir .eq. axis_2) then
  
          h_plus_2  = h_plus_2  + h(axis_3,axis_3)
          h_cross_2 = h_cross_2 + h(axis_3,axis_1)
  
       else if (dir .eq. axis_3) then
  
          h_plus_3  = h_plus_3  + h(axis_1,axis_1)
          h_cross_3 = h_cross_3 + h(axis_1,axis_2)
  
       endif
  
  
negative.code: 
             g = inertial_rotation(g, time)
  
             ! Absorb the factor of 2 outside the integral into the zone mass, for efficiency.
  
             dm = TWO * rho(i,j,k) * vol(i,j,k)
  
             if (dim .eq. 3) then
  
                do m = 1, 3
                   do l = 1, 3
                      dQtt(l,m) = dQtt(l,m) + dM * (vel(l) * vel(m) + pos(l) * g(m))
                   enddo
                enddo
  
             else
  
                ! For axisymmetric coordinates we need to be careful here.
                ! We want to calculate the quadrupole tensor in terms of
                ! Cartesian coordinates but our coordinates are cylindrical (R, z).
                ! What we can do is to first express the Cartesian coordinates
                ! as (x, y, z) = (R cos(phi), R sin(phi), z). Then we can integrate
                ! out the phi coordinate for each component. The off-diagonal components
                ! all then vanish automatically. The on-diagonal components xx and yy
                ! pick up a factor of cos**2(phi) which when integrated from (0, 2*pi)
                ! yields pi. Note that we're going to choose that the cylindrical z axis
                ! coincides with the Cartesian x-axis, which is our default choice.
  
                ! We also need to then divide by the volume by 2*pi since
                ! it has already been integrated out.
  
                dm = dm / (TWO * M_PI)
  
                dQtt(1,1) = dQtt(1,1) + dm * (TWO * M_PI) * (vel(2)**2 + pos(2) * g(2))
                dQtt(2,2) = dQtt(2,2) + dm * M_PI * (vel(1)**2 + pos(1) * g(1))
                dQtt(3,3) = dQtt(3,3) + dm * M_PI * (vel(1)**2 + pos(1) * g(1))
  
             endif
  
          enddo
       enddo
  