& 0.09950212445529771576d0, 0.37496975838746404080d0,&
         & 0.00000000004528583204d0, 0.00000000066861092459d0, &
         &-0.00000000924648168570d0, -0.00000008002058622489d0,&
         & 0.00000087578676979688d0, 0.00000628456150216647d0, &
         &-0.00004276239654834506d0, -0.00034630842915911896d0,&
         & 0.00089431202169021474d0, 0.01086380107342193819d0,&
         & 0.00145563493724746238d0, -0.14193383817937981478d0, &
         &-0.15148562296986494725d0, 0.35720232689066512346d0 /)
    real(fdp), dimension(0:125), parameter :: c= (/&
         & 0.00000000032016204980d0, -0.00000000177431944127d0,&
         & 0.00000001367520951871d0, -0.00000006171755212447d0, &
         &-0.00000043602779932091d0, 0.00000179827173442894d0,&
         & 0.00004990624965075363d0, -0.00018618453816525107d0, &
         &-0.00206578613620093339d0, 0.00640292258731756063d0,&
         & 0.04663627606907821596d0, -0.11400070635033154332d0, &
         &-0.26159330264498333979d0, 0.30099732306965462346d0, &
         &-0.00000000000059774708d0, -0.00000000038827648808d0,&
         & 0.00000000527269271362d0, 0.00000002971059814012d0, &
         &-0.00000052701406514027d0, -0.00000320369808906187d0,&
         & 0.00004475144442564492d0, 0.00016907703884518866d0, &
         &-0.00209937055264007564d0, -0.00491023437935628693d0,&
         & 0.04967175540701545140d0, 0.04195559522923170315d0, &
         &-0.33516091307165838036d0, -0.02375823895638961835d0,&
         & 0.00000000000283380103d0, -0.00000000044715450235d0,&
         & 0.00000000017695199610d0, 0.00000006053379840070d0, &
         &-0.00000002888550814408d0, -0.00000593780185020349d0,&
         & 0.00000514027455732102d0, 0.00035667482807472357d0, &
         &-0.00038132729962607730d0, -0.01158750857180657934d0,&
         & 0.01376562658346365929d0, 0.14388461069626297876d0, &
         &-0.13107454661755534290d0, -0.27409127395927545296d0,&
         & 0.00000000002683697058d0, -0.00000000024302161327d0, &
         &-0.00000000430151236236d0, 0.00000003557058903199d0,&
         & 0.00000049328952592926d0, -0.00000365732389065382d0, &
         &-0.00003642356546467712d0, 0.00023636175253326716d0,&
         & 0.00154549935517002108d0, -0.00837084637397210913d0, &
         &-0.02941202173717025197d0, 0.11713664042453548043d0,&
         & 0.15186375421302413107d0, -0.25912851048611625179d0,&
         & 0.00000000002914434604d0, 0.00000000015185121893d0, &
         &-0.00000000488084673389d0, -0.00000001936460341189d0,&
         & 0.00000057968735781303d0, 0.00000158978381365925d0, &
         &-0.00004521858314881102d0, -0.00007426575562367423d0,&
         & 0.00206293925392279862d0, 0.00143934355204571641d0, &
         &-0.04414786410004317489d0, -0.00317227451807568106d0,&
         & 0.27328377353032129600d0, -0.02616867939853747003d0,&
         & 0.00000000000644550655d0, 0.00000000040249587279d0, &
         &-0.00000000123789833341d0, -0.00000005582997578770d0,&
         & 0.00000016975984093144d0, 0.00000524003568623094d0, &
         &-0.00001539871212028104d0, -0.00030372172637669852d0,&
         & 0.00082284448843625174d0, 0.00923614667011392611d0, &
         &-0.02069442879928366384d0, -0.10834973445898588438d0,&
         & 0.14982326837249145873d0, 0.20317989938720766823d0, &
         &-0.00000000002115691215d0, 0.00000000029752035314d0,&
         & 0.00000000334739862269d0, -0.00000004304625868398d0, &
         &-0.00000036733494871207d0, 0.00000425328270524107d0,&
         & 0.00002584872967637478d0, -0.00026244994622929809d0, &
         &-0.00102234403086242668d0, 0.00863185023199066935d0,&
         & 0.01816610124346617164d0, -0.11151661591132254182d0, &
         &-0.08978791805571149950d0, 0.23370422835726857838d0, &
         &-0.00000000002957977086d0, -0.00000000006150420279d0,&
         & 0.00000000488089472193d0, 0.00000000624218182859d0, &
         &-0.00000056561995267846d0, -0.00000032077350100271d0,&
         & 0.00004278569848964830d0, -0.00000053069698280036d0, &
         &-0.00187141538785600230d0, 0.00073172821883575326d0,&
         & 0.03832854206202077374d0, -0.01874044416229564036d0, &
         &-0.23027059405144880640d0, 0.05794254714300082167d0, &
         &-0.00000000001203179184d0, -0.00000000035509887475d0,&
         & 0.00000000213620090322d0, 0.00000004813181550803d0, &
         &-0.00000026815627351126d0, -0.00000439191902440609d0,&
         & 0.00002211921535574542d0, 0.00024593505480524700d0, &
         &-0.00106158930833741898d0, -0.00722666167492552731d0,&
         & 0.02396424370522244155d0, 0.08278313570069734843d0, &
         &-0.15890724632166919294d0, -0.15383825653750118007d0 /)
    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 /)

    if (x  <  0.85d0) then
       t = x * x
       y = (((((((a(0) * t + a(1)) * t + a(2)) * t + a(3)) * t +&
            & a(4)) * t + a(5)) * t + a(6)) * t + a(7)) * x
       y = (((((((a(8) * t + a(9)) * t + a(10)) * t + a(11)) * t +&
            & a(12)) * t + a(13)) * t + a(14)) * t + a(15)) / x + y *&
            & log(x)
    else if (x  <  4.15d0) then
       v = 4 / x
       t = x - v
       k = int(t + 4)
       t = t - (k - 3.5d0)
       k = k * 14
       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)) * t +&
            & b(k + 13)) * v
    else if (x  <  12.5d0) then
       k = int(x)
       t = x - (k + 0.5d0)
       k = 14 * (k - 4)
       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 / x
       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(x + theta))
    end if

    f = y

  end function f_besy1



  recursive function f_besyn(n, x) result(f)
    !
    ! Bessel function 2nd kind, order n
    !

    real(fdp), intent(in) :: n, x
    real(fdp) :: f
    
    integer :: j, nn
    real(fdp) :: by, bym, byp, tox

    nn = nint(n)
    if (nn == 0) then
       f = f_besy0(x)
    else if (nn == 1) then
       f = f_besy1(x)
    else
       tox = 2.0d0/x
       by = f_besy1(x)
       bym = f_besy0(x)
       do j = 1, nn-1
          byp = j*tox*by-bym
          bym = by
          by = byp
       end do
       f = by
    end if
  end function f_besyn

  
  
  recursive function f_besi0(x) result(f)
    !
    ! Modified bessel function 1st kind, order 0
    !

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

    real(fdp), dimension(0:64), parameter :: a= (/&
         & 8.5246820682016865877d-11, 2.5966600546497407288d-9,&
         & 7.9689994568640180274d-8, 1.9906710409667748239d-6,&
         & 4.0312469446528002532d-5, 6.4499871606224265421d-4,&
         & 7.9012345761930579108d-3, 7.1111111109207045212d-2,&
         & 4.4444444444472490900d-1, 1.7777777777777532045d0,&
         & 4.0000000000000011182d0, 3.9999999999999999800d0,&
         & 1.0000000000000000001d0, 1.1520919130377195927d-10,&
         & 2.2287613013610985225d-9, 8.1903951930694585113d-8,&
         & 1.9821560631611544984d-6, 4.0335461940910133184d-5,&
         & 6.4495330974432203401d-4, 7.9013012611467520626d-3,&
         & 7.1111038160875566622d-2, 4.4444450319062699316d-1,&
         & 1.7777777439146450067d0, 4.0000000132337935071d0,&
         & 3.9999999968569015366d0, 1.0000000003426703174d0,&
         & 1.5476870780515238488d-10, 1.2685004214732975355d-9,&
         & 9.2776861851114223267d-8, 1.9063070109379044378d-6,&
         & 4.0698004389917945832d-5, 6.4370447244298070713d-4,&
         & 7.9044749458444976958d-3, 7.1105052411749363882d-2,&
         & 4.4445280640924755082d-1, 1.7777694934432109713d0,&
         & 4.0000055808824003386d0, 3.9999977081165740932d0,&
         & 1.0000004333949319118d0, 2.0675200625006793075d-10, &
         &-6.1689554705125681442d-10, 1.2436765915401571654d-7,&
         & 1.5830429403520613423d-6, 4.2947227560776583326d-5,&
         & 6.3249861665073441312d-4, 7.9454472840953930811d-3,&
         & 7.0994327785661860575d-2, 4.4467219586283000332d-1,&
         & 1.7774588182255374745d0, 4.0003038986252717972d0,&
         & 3.9998233869142057195d0, 1.0000472932961288324d0,&
         & 2.7475684794982708655d-10, -3.8991472076521332023d-9,&
         & 1.9730170483976049388d-7, 5.9651531561967674521d-7,&
         & 5.1992971474748995357d-5, 5.7327338675433770752d-4,&
         & 8.2293143836530412024d-3, 6.9990934858728039037d-2,&
         & 4.4726764292723985087d-1, 1.7726685170014087784d0,&
         & 4.0062907863712704432d0, 3.9952750700487845355d0,&
         & 1.0016354346654179322d0 /)
    real(fdp), dimension(0:69), parameter :: b = (/&
         & 6.7852367144945531383d-8, 4.6266061382821826854d-7,&
         & 6.9703135812354071774d-6, 7.6637663462953234134d-5,&
         & 7.9113515222612691636d-4, 7.3401204731103808981d-3,&
         & 6.0677114958668837046d-2, 4.3994941411651569622d-1,&
         & 2.7420017097661750609d0, 14.289661921740860534d0,&
         & 59.820609640320710779d0, 188.78998681199150629d0,&
         & 399.87313678256011180d0, 427.56411572180478514d0,&
         & 1.8042097874891098754d-7, 1.2277164312044637357d-6,&
         & 1.8484393221474274861d-5, 2.0293995900091309208d-4,&
         & 2.0918539850246207459d-3, 1.9375315654033949297d-2,&
         & 1.5985869016767185908d-1, 1.1565260527420641724d0,&
         & 7.1896341224206072113d0, 37.354773811947484532d0,&
         & 155.80993164266268457d0, 489.52113711585409180d0,&
         & 1030.9147225169564806d0, 1093.5883545113746958d0,&
         & 4.8017305613187493564d-7, 3.2613178439123800740d-6,&
         & 4.9073137508166159639d-5, 5.3806506676487583755d-4,&
         & 5.5387918291051866561d-3, 5.1223717488786549025d-2,&
         & 4.2190298621367914765d-1, 3.0463625987357355872d0,&
         & 18.895299447327733204d0, 97.915189029455461554d0,&
         & 407.13940115493494659d0, 1274.3088990480582632d0,&
         & 2670.9883037012547506d0, 2815.7166284662544712d0,&
         & 1.2789926338424623394d-6, 8.6718263067604918916d-6,&
         & 1.3041508821299929489d-4, 1.4282247373727478920d-3,&
         & 1.4684070635768789378d-2, 1.3561403190404185755d-1,&
         & 1.1152592585977393953d0, 8.0387088559465389038d0,&
         & 49.761318895895479206d0, 257.26842323135291380d0,&
         & 1066.8543146269566231d0, 3328.3874581009636362d0,&
         & 6948.8586598121634874d0, 7288.4893398212481055d0,&
         & 3.4093503681970328930d-6, 2.3079025203103376076d-5,&
         & 3.4691373283901830239d-4, 3.7949949772229085450d-3,&
         & 3.8974209677945602145d-2, 3.5949483804148783710d-1,&
         & 2.9522878893539528226d0, 21.246564609514287056d0,&
         & 131.28727387146173141d0, 677.38107093296675421d0,&
         & 2802.3724744545046518d0, 8718.5731420798254081d0,&
         & 18141.348781638832286d0, 18948.925349296308859d0 /)
    real(fdp), dimension(0:44), parameter :: c = (/&
         & 2.5568678676452702768d-15, 3.0393953792305924324d-14,&
         & 6.3343751991094840009d-13, 1.5041298011833009649d-11,&
         & 4.4569436918556541414d-10, 1.7463930514271679510d-8,&
         & 1.0059224011079852317d-6, 1.0729838945088577089d-4,&
         & 5.1503226936425277380d-2, 5.2527963991711562216d-15,&
         & 7.2021184814210056410d-15, 7.2561421229904797156d-13,&
         & 1.4823121466731042510d-11, 4.4602670450376245434d-10,&
         & 1.7463600061788679671d-8, 1.0059226091322347560d-6,&
         & 1.0729838937545111487d-4, 5.1503226936437300716d-2,&