J1 Function

    pure function J1(x) result(val)
      !> Bessel function J1(x)
      implicit none
      real(wp), intent(in) :: x
      real(wp)             :: val
      ! work variables
      real(wp) :: ax,xx,z,y
      real(wp), parameter :: r(6) = [ 72362614232.0_wp,-7895059235.0_wp, 242396853.1_wp, &
                                     -2972611.439_wp,   15704.48260_wp, -30.16036606_wp]
      real(wp), parameter :: s(6) = [144725228442.0_wp, 2300535178.0_wp, 18583304.74_wp, &
                                     99447.43394_wp,    376.9991397_wp,  1.0_wp]
      real(wp), parameter :: p(5) = [1.0_wp, 0.183105E-2_wp,     -0.3516396496E-4_wp, &
                                             0.2457520174E-5_wp, -0.240337019E-6_wp]
      real(wp), parameter :: q(5) = [ .04687499995_wp,-.2002690873E-3_wp,.8449199096E-5_wp,&
                                     -.88228987E-6_wp,.105787412E-6_wp]

      if (abs(x).lt.8.0_wp) then
         y   = x**2
         val = x*(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-2.356194491_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))))))*sign(1.0_wp,x)
      end if
      return
    end function J1