J0 Function

    pure function J0(x) result(val)
      !> Bessel function J0(x)
      implicit none
      real(wp), intent(in) :: x
      real(wp)             :: val
      ! work variables
      real(wp) :: ax,xx,z,y
      real(wp), parameter :: r(6) = [ 57568490574.0_wp,-13362590354.0_wp, 651619640.7_wp,&
                                     -11214424.18_wp,   77392.33017_wp,  -184.9052456_wp]
      real(wp), parameter :: s(6) = [ 57568490411.0_wp, 1029532985.0_wp, 9494680.718_wp, &
                                      59272.64853_wp,   267.8532712_wp,  1.0_wp]
      real(wp), parameter :: p(5) = [ 1.0_wp,-0.1098628627E-2_wp, 0.2734510407E-4_wp,&
                                             -0.2073370639E-5_wp, 0.2093887211E-6_wp]
      real(wp), parameter :: q(5) = [-.1562499995E-1_wp,.1430488765E-3_wp,-.6911147651E-5_wp,&
                                      .7621095161E-6_wp,-.934945152E-7_wp]

      if (abs(x).lt.8.0_wp) then
         y   = x**2
         val = (r(1)+y*(r(2)+y*(r(3)+y*(r(4)+y*(r(5)+y*r(6)))))) / &
               (s(1)+y*(s(2)+y*(s(3)+y*(s(4)+y*(s(5)+y*s(6))))))
      else
         ax  = abs(x)
         z   = 8.0_wp/ax
         y   = z**2
         xx  = ax-.785398164_wp
         val = sqrt(.636619772_wp/ax)*(cos(xx)*(p(1)+y*(p(2)+y*(p(3)+y*(p(4)+y*p(5))))) - &
               z*sin(xx)*(q(1)+y*(q(2)+y*(q(3)+y*(q(4)+y*q(5))))))
      end if
      return
    end function J0