validate_quadrature Subroutine

   subroutine validate_quadrature(n_dirs, s_x, s_y, s_z, w)
      integer, intent(in) :: n_dirs
      real(rk), intent(in) :: s_x(:), s_y(:), s_z(:), w(:)

      integer :: m
      real(rk) :: w_sum, length
      real(rk) :: mxx, myy, mzz, mxy, mxz, myz
      real(rk), parameter :: tol_quad = 1.0e-5_rk

      w_sum = zero
      mxx = zero
      myy = zero
      mzz = zero
      mxy = zero
      mxz = zero
      myz = zero

      do m = 1, n_dirs
         length = sqrt(s_x(m)**2 + s_y(m)**2 + s_z(m)**2)
         if (abs(length - 1.0_rk) > tol_quad) then
            call fatal_error('dom', 'DOM quadrature direction vector does not have unit length')
         end if
         if (w(m) <= zero) then
            call fatal_error('dom', 'DOM quadrature has non-positive weights')
         end if

         w_sum = w_sum + w(m)
         mxx = mxx + w(m) * s_x(m)**2
         myy = myy + w(m) * s_y(m)**2
         mzz = mzz + w(m) * s_z(m)**2
         mxy = mxy + w(m) * s_x(m) * s_y(m)
         mxz = mxz + w(m) * s_x(m) * s_z(m)
         myz = myz + w(m) * s_y(m) * s_z(m)
      end do

      if (abs(w_sum - 4.0_rk*pi) > tol_quad) then
         call fatal_error('dom', 'DOM quadrature weights do not sum to 4*pi')
      end if

      if (abs(mxx - 4.0_rk*pi/3.0_rk) > tol_quad .or. &
          abs(myy - 4.0_rk*pi/3.0_rk) > tol_quad .or. &
          abs(mzz - 4.0_rk*pi/3.0_rk) > tol_quad) then
         call fatal_error('dom', 'DOM quadrature second moments are not isotropic')
      end if

      if (abs(mxy) > tol_quad .or. abs(mxz) > tol_quad .or. abs(myz) > tol_quad) then
         call fatal_error('dom', 'DOM quadrature cross moments are not zero')
      end if
   end subroutine validate_quadrature