eikonal Subroutine

  pure subroutine eikonal(theta, invR, trip, qval, bi, seik, dsdthet, dpsidrp)
    use constants, only: twopi
    implicit none
    real, dimension(-ntgrid:), intent(in) :: theta, invR, trip
    real, dimension(-ntgrid:), intent(out) :: seik, dsdthet
    real, intent(in out) :: qval
    real, intent(in) :: bi
    real, intent(out) :: dpsidrp

    dsdthet = -bi * invR**2 / trip
    call integrate(dsdthet, theta, seik, ntgrid)

    if (local_eq) then
       dpsidrp = -(seik(nth) - seik(-nth)) / (twopi * qval)
       seik = seik / dpsidrp
       dsdthet = dsdthet / dpsidrp
    else
       dpsidrp = 1.
       qval = -(seik(nth) - seik(-nth)) / (twopi)
    end if
  end subroutine eikonal