function.f90

任意函数表达式求值模块 模块还提供了数学上常用的函数

         & 0.17760601686906713536d0, -0.07380429510868722524d0 /)
    real(fdp), dimension(0:111), parameter :: b = (/ &
         &-0.00000000002958527319d0, -0.00000000004799424787d0,&
         & 0.00000000070761365903d0, 0.00000000786916310954d0,&
         & 0.00000004406321746330d0, 0.00000010474207054310d0, &
         &-0.00000070359581447993d0, -0.00001068225180236166d0, &
         &-0.00007636033201484949d0, -0.00034318294473105478d0, &
         &-0.00050019381740765567d0, 0.00898034680285547482d0,&
         & 0.14783488851798565262d0, 0.01218096458136948754d0,  &
         &-0.00000000003790394288d0, -0.00000000083012324083d0, &
         &-0.00000000416344143065d0, -0.00000000159750317942d0,&
         & 0.00000011897479804282d0, 0.00000089462034665663d0,&
         & 0.00000279277855080310d0, -0.00000710267818373816d0, &
         &-0.00014158631248747538d0, -0.00089512296835766028d0, &
         &-0.00286419212921734238d0, 0.00448696141199596815d0,&
         & 0.16247276853497028694d0, 0.16806523426268420517d0,&
         & 0.00000000049508820408d0, 0.00000000240256475451d0, &
         &-0.00000000341287024760d0, -0.00000007926555593644d0, &
         &-0.00000027510687602104d0, 0.00000082850014912691d0,&
         & 0.00001202838962100801d0, 0.00004446440998573947d0, &
         &-0.00006345799532264893d0, -0.00153636349657759957d0, &
         &-0.00784582749467635199d0, -0.01091113270659651752d0,&
         & 0.15855078058220179592d0, 0.33112075987570150227d0, &
         &-0.00000000121381163034d0, -0.00000000185252445637d0,&
         & 0.00000002782181715752d0, 0.00000009983716245339d0, &
         &-0.00000051382328105461d0, -0.00000413045112612025d0,&
         & 0.00000069219807212556d0, 0.00011346881914173059d0,&
         & 0.00044857142622311310d0, -0.00075862241142011713d0, &
         &-0.01330170710405210993d0, -0.04337207915110066583d0,&
         & 0.10708263075506099964d0, 0.46937172516119821029d0,&
         & 0.00000000123954313286d0, -0.00000000131506473470d0, &
         &-0.00000003208158247443d0, 0.00000006936208476862d0,&
         & 0.00000096326321605564d0, -0.00000186813677726279d0, &
         &-0.00003296596487362508d0, -0.00001027532398838717d0,&
         & 0.00088296805857821394d0, 0.00288502916411875494d0, &
         &-0.00981551522729451135d0, -0.08175499367255420548d0, &
         &-0.01970902432425660780d0, 0.51958016003260732459d0, &
         &-0.00000000050293922567d0, 0.00000000259058191150d0,&
         & 0.00000000380069015768d0, -0.00000012579287605830d0,&
         & 0.00000029666593806065d0, 0.00000532199410651341d0, &
         &-0.00001466013674557952d0, -0.00021196179800800308d0,&
         & 0.00014889075157975250d0, 0.00598670133725335804d0,&
         & 0.00917434640154761191d0, -0.08592436550224288681d0, &
         &-0.19700894914535846940d0, 0.41202451505149651586d0,&
         & 0.00000000002689416111d0, -0.00000000060225798842d0,&
         & 0.00000000716632287642d0, -0.00000003038834206093d0, &
         &-0.00000049829955736697d0, 0.00000364678661079850d0,&
         & 0.00002599276311771812d0, -0.00016374249229498868d0, &
         &-0.00112385259950125883d0, 0.00342431039748960465d0,&
         & 0.03016663018029814938d0, -0.02432512430353342927d0, &
         &-0.31797371916576712163d0, 0.14418001571733803493d0,&
         & 0.00000000000194632210d0, -0.00000000013451327719d0,&
         & 0.00000000336581427068d0, 0.00000002163920016772d0, &
         &-0.00000050868502638387d0, -0.00000126562372408989d0,&
         & 0.00003554054033687039d0, 0.00007490048474688881d0, &
         &-0.00142380042008861330d0, -0.00355314274445060456d0,&
         & 0.03042001885994619102d0, 0.07365497405664882878d0, &
         &-0.26882214651298246524d0, -0.16578636480856581212d0 /)
    real(fdp), dimension(0:125), parameter :: c = (/&
         & 0.00000000014133441141d0, -0.00000000116239182574d0,&
         & 0.00000000560813675548d0, 0.00000004361111919906d0, &
         &-0.00000019980743680592d0, -0.00000623828267379045d0,&
         & 0.00002659779111113435d0, 0.00034429768949740887d0, &
         &-0.00128058451746064875d0, -0.01165906901727504569d0,&
         & 0.03800023545011044967d0, 0.13079665132249175616d0, &
         &-0.30099732306965462304d0, -0.19470500862950453349d0,&
         & 0.00000000002985070388d0, -0.00000000043935370028d0, &
         &-0.00000000270094988368d0, 0.00000005270139367180d0,&
         & 0.00000035596645009533d0, -0.00000559393055101660d0, &
         &-0.00002415386269153308d0, 0.00034989509210648774d0,&
         & 0.00098204687587121096d0, -0.01241793885175385487d0, &
         &-0.01398519840974389960d0, 0.16758045653582919005d0,&
         & 0.02375823895638961834d0, -0.33948059288191103829d0,&
         & 0.00000000003440543098d0, -0.00000000001492311224d0, &
         &-0.00000000550307983816d0, 0.00000000288861169685d0,&
         & 0.00000065975576340131d0, -0.00000064253433004282d0, &
         &-0.00005095354686815779d0, 0.00006355454993858760d0,&
         & 0.00231750171436134067d0, -0.00344140664586595266d0, &
         &-0.04796153689875432703d0, 0.06553727330877767204d0,&
         & 0.27409127395927545297d0, -0.17324243491898233567d0,&
         & 0.00000000001869953108d0, 0.00000000035678205287d0, &
         &-0.00000000323369443017d0, -0.00000004932837601739d0,&
         & 0.00000040636932259569d0, 0.00000455294558480472d0, &
         &-0.00003376596464782138d0, -0.00025758322585307069d0,&
         & 0.00167416927479443727d0, 0.00735300543429220468d0, &
         &-0.03904554680817849396d0, -0.07593187710651205994d0,&
         & 0.25912851048611625179d0, 0.11731328614820863082d0, &
         &-0.00000000001168499781d0, 0.00000000040491570620d0,&
         & 0.00000000176042185153d0, -0.00000005796810963324d0, &
         &-0.00000017664264701199d0, 0.00000565232278687159d0,&
         & 0.00001060939366068237d0, -0.00034382320897779425d0, &
         &-0.00028786871040915476d0, 0.01103696602501040460d0,&
         & 0.00105742483935856071d0, -0.13664188676516064192d0,&
         & 0.02616867939853747003d0, 0.27020510536578747599d0, &
         &-0.00000000003097140650d0, 0.00000000010275535029d0,&
         & 0.00000000507546062007d0, -0.00000001697584561546d0, &
         &-0.00000058222618994516d0, 0.00000192483899143098d0,&
         & 0.00004338881805420200d0, -0.00013714074807064326d0, &
         &-0.00184722933402281351d0, 0.00517360719982082990d0,&
         & 0.03611657815299529568d0, -0.07491163418624572802d0, &
         &-0.20317989938720766824d0, 0.17121062620272384486d0, &
         &-0.00000000002289344046d0, -0.00000000027762757821d0,&
         & 0.00000000391330214230d0, 0.00000003673304032817d0, &
         &-0.00000047258696909223d0, -0.00000323109113206792d0,&
         & 0.00003749284946160478d0, 0.00017039067180362504d0, &
         &-0.00172637004639815403d0, -0.00454152531086626043d0,&
         & 0.03717220530377418124d0, 0.04489395902785574534d0, &
         &-0.23370422835726857839d0, -0.06753037249787639679d0,&
         & 0.00000000000473299038d0, -0.00000000040489249115d0, &
         &-0.00000000056747261733d0, 0.00000005656135976496d0,&
         & 0.00000003564150059329d0, -0.00000534821220288168d0,&
         & 0.00000007581385461348d0, 0.00031190256463318867d0, &
         &-0.00014634564376714538d0, -0.00958213551550479850d0,&
         & 0.00624681472076521329d0, 0.11513529702572439703d0, &
         &-0.05794254714300082167d0, -0.22523211169118786537d0,&
         & 0.00000000002732444065d0, -0.00000000017726475495d0, &
         &-0.00000000437562701999d0, 0.00000002681536885560d0,&
         & 0.00000048799100503315d0, -0.00000276490187540644d0, &
         &-0.00003513357925824063d0, 0.00017693155138571443d0,&
         & 0.00144533233498513085d0, -0.00599106092630544975d0, &
         &-0.02759437856689911694d0, 0.07945362316083459396d0,&
         & 0.15383825653750118008d0, -0.17121430684466928734d0 /)
    real(fdp), dimension(0:51), parameter :: d = (/&
         & 1.059601355592185731d-14, -2.71150591218550377d-13,&
         & 8.6514809056201638d-12, -4.6264028554286627d-10,&
         & 5.0815403835647104d-8, -1.76722552048141208d-5,&
         & 0.16286750396763997378d0, 2.949651820598278873d-13, &
         &-8.818215611676125741d-12, 3.571119876162253451d-10, &
         &-2.631924120993717060d-8, 4.709502795656698909d-6, &
         &-5.208333333333283282d-3, 7.18344107717531977d-15, &
         &-2.51623725588410308d-13, 8.6017784918920604d-12, &
         &-4.6256876614290359d-10, 5.0815343220437937d-8, &
         &-1.76722551764941970d-5, 0.16286750396763433767d0,&
         & 2.2327570859680094777d-13, -8.464594853517051292d-12,&
         & 3.563766464349055183d-10, -2.631843986737892965d-8,&
         & 4.709502342288659410d-6, -5.2083333332278466225d-3,&
         & 5.15413392842889366d-15, -2.27740238380640162d-13,&
         & 8.4827767197609014d-12, -4.6224753682737618d-10,&
         & 5.0814848128929134d-8, -1.76722547638767480d-5,&
         & 0.16286750396748926663d0, 1.7316195320192170887d-13, &
         &-7.971122772293919646d-12, 3.544039469911895749d-10, &
         &-2.631443902081701081d-8, 4.709498228695400603d-6, &
         &-5.2083333315143653610d-3, 3.84653681453798517d-15, &
         &-2.04464520778789011d-13, 8.3089298605177838d-12, &
         &-4.6155016158412096d-10, 5.0813263696466650d-8, &
         &-1.76722528311426167d-5, 0.16286750396650065930d0,&
         & 1.3797879972460878797d-13, -7.448089381011684812d-12,&
         & 3.512733797106959780d-10, -2.630500895563592722d-8,&
         & 4.709483934775839193d-6, -5.2083333227940760113d-3 /)

    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)
       y = ((((((a(8) * t + a(9)) * t + a(10)) * t + a(11)) * t +&
            & a(12)) * t + a(13)) * t + a(14)) * t + a(15) + y *&
            & log(x)
    else if (x  <  4.5d0) then
       t = x - 4 / x
       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)
    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 * sin(x + theta)
    end if

    f = y
  end function f_besy0
  


  recursive function f_besy1(x) result(f)
    !
    ! Bessel function 2nd kind, order 1
    !

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

    real(fdp) :: t, v, y, theta
    integer :: k
    real(fdp), dimension(0:15), parameter :: a = (/ &
         &-0.00000000000009464990d0, 0.00000000002141432796d0, &
         &-0.00000000359779443470d0, 0.00000043173541397185d0, &
         &-0.00003453883313621946d0, 0.00165786399054057257d0, &
         &-0.03978873577297383381d0, 0.31830988618379067154d0, &
         &-0.00000000005571931463d0, 0.00000000893098810365d0, &
         &-0.00000099267351748541d0, 0.00007164268731546281d0, &
         &-0.00295530533604595764d0, 0.05434868816050718549d0, &
         &-0.19605709064623884320d0, -0.63661977236758134367d0 /)
    real(fdp), dimension(0:111), parameter :: b = (/&
         & 0.00000000002553685555d0, 0.00000000007384080449d0, &
         &-0.00000000029605446927d0, -0.00000000520817391060d0, &
         &-0.00000003676554001376d0, -0.00000015326525336270d0, &
         &-0.00000013967947147386d0, 0.00000388730181348638d0,&
         & 0.00004052240780599446d0, 0.00025133005806772436d0,&
         & 0.00109380629039004777d0, 0.00298157276348996763d0,&
         & 0.00047188983238582062d0, -0.19643830264444486066d0,&
         & 0.00000000001056916126d0, 0.00000000062160122556d0,&
         & 0.00000000385741197273d0, 0.00000000749146940600d0, &
         &-0.00000006306874425799d0, -0.00000068698336577990d0, &
         &-0.00000329785491054157d0, -0.00000557893716372629d0,&
         & 0.00004634569030169549d0, 0.00049086790937752985d0,&
         & 0.00256488729380268318d0, 0.00822571374211646414d0,&
         & 0.01094512145046992425d0, -0.19159562910398869317d0, &
         &-0.00000000040211480479d0, -0.00000000229721816531d0,&
         & 0.00000000062931855684d0, 0.00000006509879184236d0,&
         & 0.00000031589903600713d0, -0.00000003093448289001d0, &
         &-0.00000845532004628483d0, -0.00004930719004112743d0, &
         &-0.00009865914153735679d0, 0.00047059803263198826d0,&
         & 0.00472271367216891688d0, 0.01915518502418162917d0,&
         & 0.03722391361384158740d0, -0.16933230788966878110d0,&
         & 0.00000000110728466147d0, 0.00000000205608441891d0, &
         &-0.00000002536199750468d0, -0.00000011751665325851d0,&
         & 0.00000034045343850770d0, 0.00000436083390507752d0,&
         & 0.00000834725302600388d0, -0.00007154971448642437d0, &
         &-0.00051956165332380432d0, -0.00100811886858456098d0,&
         & 0.00436854257164588225d0, 0.03425644258611852849d0,&
         & 0.09073971172014934990d0, -0.10791597118914788657d0, &
         &-0.00000000117749768422d0, 0.00000000168403859797d0,&
         & 0.00000003471030171126d0, -0.00000006227816248757d0, &
         &-0.00000118405395385178d0, 0.00000014478514453927d0,&
         & 0.00003630744090017806d0, 0.00010481065165516454d0, &
         &-0.00052483993430622472d0, -0.00407014291708523893d0, &
         &-0.00574311514714191566d0, 0.03534768318503458091d0,&
         & 0.16539360869993020163d0, 0.01986240425870529979d0,&
         & 0.00000000040080505520d0, -0.00000000336044652853d0, &
         &-0.00000000092387130705d0, 0.00000017840823064933d0, &
         &-0.00000022643988836534d0, -0.00000808896250262948d0, &
         &-0.00000169403516892975d0, 0.00026617475073041589d0,&
         & 0.00072562705080786576d0, -0.00399208199997748377d0, &
         &-0.02399816920965700841d0, -0.00925574744557542440d0,&
         & 0.20082814022820666990d0, 0.21040564991987672576d0,&
         & 0.00000000006692975871d0, 0.00000000016722027050d0, &
         &-0.00000001433097981251d0, 0.00000005455077276306d0,&
         & 0.00000104954046999028d0, -0.00000339572219821681d0, &
         &-0.00005541577161092338d0, 0.00004344113487321887d0,&
         & 0.00184798594297387512d0, 0.00301041508468150450d0, &
         &-0.02789739141907399514d0, -0.09421628519129373841d0,&