box_get_grids Subroutine

  subroutine box_get_grids (aky, theta0, akx, ikx)
    use theta_grid, only: is_effectively_zero_shear, nperiod, shat
    use constants, only: twopi
    use theta_grid, only: nperiod
    use warning_helpers, only: is_zero, is_not_zero
    implicit none
    !> The \(k_x \rho\) grid
    real, dimension (:), intent (out) :: akx
    !> The \(k_y \rho\) grid
    real, dimension (:), intent (out) :: aky
    !> The \(\theta_0(k_x, k_y)\) grid
    real, dimension (:,:), intent (out) :: theta0
    !> Discrete kx wavenumber grid indices
    integer, dimension (:), intent (out) :: ikx

    real :: dkx, dky, ratio
    integer :: i, naky, ntheta0

    naky = size(aky)
    ntheta0 = size(akx)

    dky = twopi/ly

    if(is_effectively_zero_shear()) then   ! non-quantized b/c assumed to be periodic instead linked boundary conditions
       if (is_zero(x0)) then
          if (rtwist > 0) then
             ratio = rtwist
          else
             ratio = 1. / abs(rtwist)
          end if
          dkx = dky / ratio
       else
          if (x0 > 0.) then
             dkx = 1./x0
          else
             dkx = -x0
          end if
       end if
    else
       if (jtwist /= 0) then
          dkx = dky * twopi*(2*nperiod-1)*abs(shat)/real(jtwist)
       else
          dkx = dky
       end if
    endif

    do i = 1, naky
       aky(i) = real(i-1)*dky
    end do

    do i = 1, (ntheta0+1)/2
       ikx(i) = i-1
       akx(i) = real(i-1)*dkx
    end do

    do i = (ntheta0+1)/2+1, ntheta0
       ikx(i) = i-ntheta0-1
       akx(i) = real(i-ntheta0-1)*dkx
    end do

    if (is_not_zero(shat)) then
       ! Should the above use is_effectively_zero_shear?
       do i = 1, ntheta0
          theta0(i,1) = 0.0
          theta0(i,2:) = akx(i)/(aky(2:)*shat)
       end do
    else
       do i = 1, ntheta0
          theta0(i,1) = 0.0
          theta0(i,2:) = - akx(i)/aky(2:)   ! meaningless, so be careful
       end do
    end if

  end subroutine box_get_grids