k = 14 * (k - 8)
       y = ((((((((((((c(k) * t + c(k + 1)) * t + c(k + 2)) * t +&
            & c(k + 3)) * t + c(k + 4)) * t + c(k + 5)) * t + c(k&
            & + 6)) * t + c(k + 7)) * t + c(k + 8)) * t + c(k +&
            & 9)) * t + c(k + 10)) * t + c(k + 11)) * t + c(k +&
            & 12)) * t + c(k + 13)
    else
       v = 24 / w
       t = v * v
       k = 13 * (int(t))
       y = ((((((d(k) * t + d(k + 1)) * t + d(k + 2)) * t + d(k +&
            & 3)) * t + d(k + 4)) * t + d(k + 5)) * t + d(k + 6))&
            & * sqrt(v)
       theta = (((((d(k + 7) * t + d(k + 8)) * t + d(k + 9)) * t &
            &+ d(k + 10)) * t + d(k + 11)) * t + d(k + 12)) * v -&
            & pi4
       y = y * cos(w + theta)
    end if

    f = y
  end function f_besj0



  recursive function f_besj1(x) result(f)
    !
    ! Bessel function 1st kind, order 1
    !
    real(fdp), intent(in) :: x
    real(fdp) :: f
    
    integer :: k
    real(fdp) :: t, w, v, y, theta
    real(fdp), dimension(0:7), parameter :: a = (/ &
         &-0.00000000000014810349d0, 0.00000000003363594618d0, &
         &-0.00000000565140051697d0, 0.00000067816840144764d0, &
         &-0.00005425347222188379d0, 0.00260416666666662438d0, &
         &-0.06249999999999999799d0, 0.49999999999999999998d0 /)
    real(fdp), dimension(0:64), parameter :: b = (/&
         & 0.00000000000243721316d0, -0.00000000009400554763d0,&
         & 0.00000000306053389980d0, -0.00000008287270492518d0,&
         & 0.00000183020515991344d0, -0.00003219783841164382d0,&
         & 0.00043795830161515318d0, -0.00442952351530868999d0,&
         & 0.03157908273375945955d0, -0.14682160488052520107d0,&
         & 0.39309619054093640008d0, -0.47952808215101070280d0,&
         & 0.14148999344027125140d0, 0.00000000000182119257d0, &
         &-0.00000000006862117678d0, 0.00000000217327908360d0, &
         &-0.00000005693592917820d0, 0.00000120771046483277d0, &
         &-0.00002020151799736374d0, 0.00025745933218048448d0, &
         &-0.00238514907946126334d0, 0.01499220060892984289d0, &
         &-0.05707238494868888345d0, 0.10375225210588234727d0, &
         &-0.02721551202427354117d0, -0.06420643306727498985d0,&
         & 0.000000000001352611196d0, -0.000000000049706947875d0,&
         & 0.000000001527944986332d0, -0.000000038602878823401d0,&
         & 0.000000782618036237845d0, -0.000012349994748451100d0,&
         & 0.000145508295194426686d0, -0.001203649737425854162d0,&
         & 0.006299092495799005109d0, -0.016449840761170764763d0,&
         & 0.002106328565019748701d0, 0.058527410006860734650d0, &
         &-0.031896615709705053191d0, 0.000000000000997982124d0, &
         &-0.000000000035702556073d0, 0.000000001062332772617d0, &
         &-0.000000025779624221725d0, 0.000000496382962683556d0, &
         &-0.000007310776625173004d0, 0.000078028107569541842d0, &
         &-0.000550624088538081113d0, 0.002081442840335570371d0, &
         &-0.000771292652260286633d0, -0.019541271866742634199d0,&
         & 0.033361194224480445382d0, 0.017516628654559387164d0,&
         & 0.000000000000731050661d0, -0.000000000025404499912d0,&
         & 0.000000000729360079088d0, -0.000000016915375004937d0,&
         & 0.000000306748319652546d0, -0.000004151324014331739d0,&
         & 0.000038793392054271497d0, -0.000211180556924525773d0,&
         & 0.000274577195102593786d0, 0.003378676555289966782d0, &
         &-0.013842821799754920148d0, -0.002041834048574905921d0,&
         & 0.032167266073736023299d0 /)
    real(fdp), dimension(0:69), parameter :: c = (/ &
         &-0.00000000001185964494d0, 0.00000000039110295657d0,&
         & 0.00000000180385519493d0, -0.00000005575391345723d0, &
         &-0.00000018635897017174d0, 0.00000542738239401869d0,&
         & 0.00001181490114244279d0, -0.00033000319398521070d0, &
         &-0.00037717832892725053d0, 0.01070685852970608288d0,&
         & 0.00356629346707622489d0, -0.13524776185998074716d0,&
         & 0.00980725611657523952d0, 0.27312196367405374425d0, &
         &-0.00000000003029591097d0, 0.00000000009259293559d0,&
         & 0.00000000496321971223d0, -0.00000001518137078639d0, &
         &-0.00000057045127595547d0, 0.00000171237271302072d0,&
         & 0.00004271400348035384d0, -0.00012152454198713258d0, &
         &-0.00184155714921474963d0, 0.00462994691003219055d0,&
         & 0.03671737063840232452d0, -0.06863857568599167175d0, &
         &-0.21090395092505707655d0, 0.16126443075752985095d0, &
         &-0.00000000002197602080d0, -0.00000000027659100729d0,&
         & 0.00000000374295124827d0, 0.00000003684765777023d0, &
         &-0.00000045072801091574d0, -0.00000327941630669276d0,&
         & 0.00003571371554516300d0, 0.00017664005411843533d0, &
         &-0.00165119297594774104d0, -0.00485925381792986774d0,&
         & 0.03593306985381680131d0, 0.04997877588191962563d0, &
         &-0.22913866929783936544d0, -0.07885001422733148814d0,&
         & 0.00000000000516292316d0, -0.00000000039445956763d0, &
         &-0.00000000066220021263d0, 0.00000005511286218639d0,&
         & 0.00000005012579400780d0, -0.00000522111059203425d0, &
         &-0.00000134311394455105d0, 0.00030612891890766805d0, &
         &-0.00007103391195326182d0, -0.00949316714311443491d0,&
         & 0.00455036998246516948d0, 0.11540391585989614784d0, &
         &-0.04779493761902840455d0, -0.22837862066532347460d0,&
         & 0.00000000002697817493d0, -0.00000000016633326949d0, &
         &-0.00000000433134860350d0, 0.00000002508404686362d0,&
         & 0.00000048528284780984d0, -0.00000258267851112118d0, &
         &-0.00003521049080466759d0, 0.00016566324273339952d0,&
         & 0.00146474737522491617d0, -0.00565140892697147306d0, &
         &-0.02833882055679300400d0, 0.07580744376982855057d0,&
         & 0.16012275906960187978d0, -0.16548380461475971845d0 /)
    real(fdp), dimension(0:51), parameter :: d = (/ &
         &-1.272346002224188092d-14, 3.370464692346669075d-13, &
         &-1.144940314335484869d-11, 6.863141561083429745d-10, &
         &-9.491933932960924159d-8, 5.301676561445687562d-5,&
         & 0.1628675039676399740d0, -3.652982212914147794d-13,&
         & 1.151126750560028914d-11, -5.165585095674343486d-10,&
         & 4.657991250060549892d-8, -1.186794704692706504d-5,&
         & 1.562499999999994026d-2, -8.713069680903981555d-15,&
         & 3.140780373478474935d-13, -1.139089186076256597d-11,&
         & 6.862299023338785566d-10, -9.491926788274594674d-8,&
         & 5.301676558106268323d-5, 0.1628675039676466220d0, &
         &-2.792555727162752006d-13, 1.108650207651756807d-11, &
         &-5.156745588549830981d-10, 4.657894859077370979d-8, &
         &-1.186794650130550256d-5, 1.562499999987299901d-2, &
         &-6.304859171204770696d-15, 2.857249044208791652d-13, &
         &-1.124956921556753188d-11, 6.858482894906716661d-10, &
         &-9.491867953516898460d-8, 5.301676509057781574d-5,&
         & 0.1628675039678191167d0, -2.185193490132496053d-13,&
         & 1.048820673697426074d-11, -5.132819367467680132d-10,&
         & 4.657409437372994220d-8, -1.186794150862988921d-5,&
         & 1.562499999779270706d-2, -4.740417209792009850d-15,&
         & 2.578715253644144182d-13, -1.104148898414138857d-11,&
         & 6.850134201626289183d-10, -9.491678234174919640d-8,&
         & 5.301676277588728159d-5, 0.1628675039690033136d0, &
         &-1.755122057493842290d-13, 9.848723331445182397d-12, &
         &-5.094535425482245697d-10, 4.656255982268609304d-8, &
         &-1.186792402114394891d-5, 1.562499998712198636d-2 /)

    w = abs(x)

    if (w  <  1) then
       t = w * w
       y = (((((((a(0) * t + a(1)) * t + a(2)) * t + a(3)) * t +&
            & a(4)) * t + a(5)) * t + a(6)) * t + a(7)) * w
    else if (w  <  8.5d0) then
       t = w * w * 0.0625d0
       k = int(t)
       t = t - (k + 0.5d0)
       k = k * 13
       y = ((((((((((((b(k) * t + b(k + 1)) * t + b(k + 2)) * t +&
            & b(k + 3)) * t + b(k + 4)) * t + b(k + 5)) * t + b(k&
            & + 6)) * t + b(k + 7)) * t + b(k + 8)) * t + b(k +&
            & 9)) * t + b(k + 10)) * t + b(k + 11)) * t + b(k +&
            & 12)) * w
    else if (w  <  12.5d0) then
       k = int(w)
       t = w - (k + 0.5d0)
       k = 14 * (k - 8)
       y = ((((((((((((c(k) * t + c(k + 1)) * t + c(k + 2)) * t +&
            & c(k + 3)) * t + c(k + 4)) * t + c(k + 5)) * t + c(k&
            & + 6)) * t + c(k + 7)) * t + c(k + 8)) * t + c(k +&
            & 9)) * t + c(k + 10)) * t + c(k + 11)) * t + c(k +&
            & 12)) * t + c(k + 13)
    else
       v = 24 / w
       t = v * v
       k = 13 * (int(t))
       y = ((((((d(k) * t + d(k + 1)) * t + d(k + 2)) * t + d(k +&
            & 3)) * t + d(k + 4)) * t + d(k + 5)) * t + d(k + 6))&
            & * sqrt(v)
       theta = (((((d(k + 7) * t + d(k + 8)) * t + d(k + 9)) * t &
            &+ d(k + 10)) * t + d(k + 11)) * t + d(k + 12)) * v -&
            & pi4
       y = y * sin(w + theta)
    end if

    if (x  <  0) y = -y

    f = y
  end function f_besj1



  recursive function f_besjn(n, x) result(f)
    !
    ! Bessel function 1st kind, order n
    !

    real(fdp), intent(in) :: n, x
    real(fdp) :: f

    integer, parameter :: iexp=maxexponent(x)/2
    real(fdp), parameter :: iacc = 50.0d0
    integer :: j, jsum, m, nn
    real(fdp) :: ax, bj, bjm, bjp, summ, tox
    
    nn = nint(n)

    if (nn == 0) then
       f = f_besj0(x)
    else if (nn == 1) then
       f = f_besj1(x)
    else
       ax=abs(x)
       if (ax*ax <= 8.0d0*tiny(x)) then
          f = 0.0d0
       else if (ax > n) then
          tox = 2.0d0/ax
          bjm = f_besj0(ax)
          bj = f_besj1(ax)
          do j = 1, nn-1
             bjp = j*tox*bj-bjm
             bjm = bj
             bj = bjp
          end do
          f = bj
       else
          tox = 2.0d0/ax
          m = 2*((nn+int(sqrt(iacc*n)))/2)
          f = 0.0
          jsum = 0
          summ = 0.0
          bjp = 0.0
          bj = 1.0
          do j = m, 1, -1
             bjm = j*tox*bj-bjp
             bjp = bj
             bj = bjm
             if (exponent(bj) > iexp) then
                bj = scale(bj, -iexp)
                bjp = scale(bjp, -iexp)
                f = scale(f, -iexp)
                summ = scale(summ, -iexp)
             end if
             if (jsum /= 0) summ = summ+bj
             jsum = 1-jsum
             if (j == nn) f = bjp
          end do
          summ = 2.0d0*summ-bj
          f = f/summ
       end if
       f = merge(-f, f, x < 0.0d0 .and. mod(nn,2) == 1)
    end if
  end function f_besjn
  


  recursive function f_besy0(x) result(f)
    !
    ! Bessel function 2nd kind, order 0
    !

    real(fdp), intent(in) :: x
    real(fdp) :: f
    
    real(fdp) :: w, t, v, y, theta
    integer :: k
    real(fdp), dimension(0:15), parameter :: a = (/ &
         &-0.00000000000151249795d0, 0.00000000029979612902d0, &
         &-0.00000004317352912436d0, 0.00000431735413787068d0, &
         &-0.00027631066508933090d0, 0.00994718394324338940d0, &
         &-0.15915494309189533339d0, 0.63661977236758134306d0,&
         & 0.00000000000409490035d0, -0.00000000076925095943d0,&
         & 0.00000010358472550303d0, -0.00000949500519343105d0,&
         & 0.00053860266685948738d0, -0.01607396802593822992d0,&