get_quadrature Subroutine

   subroutine get_quadrature(quad_type, n_dirs, s_x, s_y, s_z, w)
      character(len=*), intent(in) :: quad_type
      integer, intent(out) :: n_dirs
      real(rk), allocatable, intent(out) :: s_x(:), s_y(:), s_z(:), w(:)

      integer :: L, i, j, k, sx, sy, sz, idx
      real(rk), allocatable :: mu(:)
      real(rk) :: w_ijk

      select case (trim(adjustl(quad_type)))
      case ('s2', 'S2')
         n_dirs = 8
         L = 1
         allocate(mu(L))
         mu(1) = 0.5773502691896258_rk
      case ('s4', 'S4')
         n_dirs = 24
         L = 2
         allocate(mu(L))
         mu(1) = 0.3500212_rk
         mu(2) = 0.8688903_rk
      case ('s6', 'S6')
         n_dirs = 48
         L = 3
         allocate(mu(L))
         mu(1) = 0.2666355_rk
         mu(2) = 0.6815076_rk
         mu(3) = 0.9261808_rk
      case ('s8', 'S8')
         n_dirs = 80
         L = 4
         allocate(mu(L))
         mu(1) = 0.2182179_rk
         mu(2) = 0.5773503_rk
         mu(3) = 0.7867958_rk
         mu(4) = 0.9511897_rk
      case default
         call fatal_error('dom', 'Unsupported DOM quadrature type')
      end select

      allocate(s_x(n_dirs), s_y(n_dirs), s_z(n_dirs), w(n_dirs))
      idx = 0

      ! Generate combinations in the first octant where i+j+k=L+2, then reflect
      ! to all octants.
      do i = 1, L
         do j = 1, L
            do k = 1, L
               if (i + j + k /= L + 2) cycle

               w_ijk = quadrature_octant_weight(trim(adjustl(quad_type)), i, j, k)

               do sx = -1, 1, 2
                  do sy = -1, 1, 2
                     do sz = -1, 1, 2
                        idx = idx + 1
                        if (idx > n_dirs) then
                           call fatal_error('dom', 'DOM quadrature generated too many directions')
                        end if
                        s_x(idx) = real(sx, rk) * mu(i)
                        s_y(idx) = real(sy, rk) * mu(j)
                        s_z(idx) = real(sz, rk) * mu(k)
                        w(idx) = w_ijk * (pi / 2.0_rk)
                     end do
                  end do
               end do
            end do
         end do
      end do

      deallocate(mu)

      if (idx /= n_dirs) then
         call fatal_error('dom', 'DOM quadrature generated an unexpected number of directions')
      end if

      call validate_quadrature(n_dirs, s_x, s_y, s_z, w)
   end subroutine get_quadrature