function.f90

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

         & 1.3365917359358069908d-14, -1.2932643065888544835d-13,&
         & 1.7450199447905602915d-12, 1.0419051209056979788d-11,&
         & 4.5804788198059832600d-10, 1.7442405450073548966d-8,&
         & 1.0059461453281292278d-6, 1.0729837434500161228d-4,&
         & 5.1503226940658446941d-2, 5.3771611477352308649d-14, &
         &-1.1396193006413731702d-12, 1.2858641335221653409d-11, &
         &-5.9802086004570057703d-11, 7.3666894305929510222d-10,&
         & 1.6731837150730356448d-8, 1.0070831435812128922d-6,&
         & 1.0729733111203704813d-4, 5.1503227360726294675d-2,&
         & 3.7819492084858931093d-14, -4.8600496888588034879d-13,&
         & 1.6898350504817224909d-12, 4.5884624327524255865d-11,&
         & 1.2521615963377513729d-10, 1.8959658437754727957d-8,&
         & 1.0020716710561353622d-6, 1.0730371198569275590d-4,&
         & 5.1503223833002307750d-2 /)
    real(fdp) :: w, t, y
    integer :: k

    w = abs(x)

    if (w  <  8.5d0) then
       t = w * w * 0.0625d0
       k = 13 * (int(t))
       y = (((((((((((a(k) * t + a(k + 1)) * t + a(k + 2)) * t + a(k &
            &+ 3)) * t + a(k + 4)) * t + a(k + 5)) * t + a(k + 6)) *&
            & t + a(k + 7)) * t + a(k + 8)) * t + a(k + 9)) * t + a(k&
            & + 10)) * t + a(k + 11)) * t + a(k + 12)
    else if (w  <  12.5d0) then
       k = int(w)
       t = w - k
       k = 14 * (k - 8)
       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
       t = 60 / w
       k = 9 * (int(t))
       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)) * sqrt(t) * exp(w)
    end if

    f = y
  end function f_besi0



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

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

    real(fdp), dimension(0:59), parameter :: a = (/&
         & 1.2787464404046789181d-10, 3.5705860060088241077d-9,&
         & 9.9611537619347335040d-8, 2.2395070088633043177d-6,&
         & 4.0312466928887462346d-5, 5.6437387840203722356d-4,&
         & 5.9259259312934746096d-3, 4.4444444443499008870d-2,&
         & 2.2222222222232042719d-1, 6.6666666666666139867d-1,&
         & 1.0000000000000001106d0, 4.9999999999999999962d-1,&
         & 1.7281952384448634449d-10, 3.0647204559976390130d-9,&
         & 1.0237662138842827028d-7, 2.2299494417341498163d-6,&
         & 4.0335364374929326943d-5, 5.6433440269141349899d-4,&
         & 5.9259754885893798654d-3, 4.4444399410880397870d-2,&
         & 2.2222225112835026730d-1, 6.6666665422146063244d-1,&
         & 1.0000000032274936821d0, 4.9999999961866867205d-1,&
         & 2.3216048939948030996d-10, 1.7443372702334489579d-9,&
         & 1.1596478963485415499d-7, 2.1446755518623035147d-6,&
         & 4.0697440347437076195d-5, 5.6324394900433192204d-4,&
         & 5.9283484996093060678d-3, 4.4440673899150997921d-2,&
         & 2.2222638016852657860d-1, 6.6666358151576732094d-1,&
         & 1.0000013834029985337d0, 4.9999971643129650249d-1,&
         & 3.1013758938255172562d-10, -8.4813676145611694984d-10,&
         & 1.5544980187411802596d-7, 1.7811109378708045726d-6,&
         & 4.2945322199060856985d-5, 5.5344850176852353639d-4,&
         & 5.9590327716950614802d-3, 4.4371611097707060659d-2,&
         & 2.2233578241986401111d-1, 6.6654747300463315310d-1,&
         & 1.0000756505206705927d0, 4.9997803664415994554d-1,&
         & 4.1214758313965020365d-10, -5.3613317735347429440d-9,&
         & 2.4661360807517345161d-7, 6.7144593918926723203d-7,&
         & 5.1988027944945587571d-5, 5.0165568586065803067d-4,&
         & 6.1717530047005289953d-3, 4.3745229577317251404d-2,&
         & 2.2363147971477747996d-1, 6.6475469131117660240d-1,&
         & 1.0015686689447547657d0, 4.9941120439785391891d-1 /)
    real(fdp), dimension(0:69), parameter :: b = (/&
         & 6.6324787943143095845d-8, 4.5125928898466638619d-7,&
         & 6.7937793134877246623d-6, 7.4580507871505926302d-5,&
         & 7.6866382927334005919d-4, 7.1185174803491859307d-3,&
         & 5.8721838073486424416d-2, 4.2473949281714196041d-1,&
         & 2.6396965606282079123d0, 13.710008536637016903d0,&
         & 57.158647688180932003d0, 179.46182892089389037d0,&
         & 377.57997362398478619d0, 399.87313678256009819d0,&
         & 1.7652713206027939711d-7, 1.1988179244834708057d-6,&
         & 1.8037851545747139231d-5, 1.9775785516370314656d-4,&
         & 2.0354870702829387283d-3, 1.8822164191032253600d-2,&
         & 1.5500485219010424263d-1, 1.1190100010560573210d0,&
         & 6.9391565185406617552d0, 35.948170579648649345d0,&
         & 149.41909525103032616d0, 467.42979492780642582d0,&
         & 979.04227423171290408d0, 1030.9147225169564443d0,&
         & 4.7022299276154507603d-7, 3.1878571710170115972d-6,&
         & 4.7940153875711448496d-5, 5.2496623508411440227d-4,&
         & 5.3968661134780824779d-3, 4.9837081920693776234d-2,&
         & 4.0979593830387765545d-1, 2.9533186922862948404d0,&
         & 18.278176130722516369d0, 94.476497150189121070d0,&
         & 391.66075612645333624d0, 1221.4182034643210345d0,&
         & 2548.6177980961291004d0, 2670.9883037012546541d0,&
         & 1.2535083724002034147d-6, 8.4845871420655708250d-6,&
         & 1.2753227372734042108d-4, 1.3950105363562648921d-3,&
         & 1.4325473993765291906d-2, 1.3212452778932829125d-1,&
         & 1.0849287786885151432d0, 7.8068089156260172673d0,&
         & 48.232254570679165833d0, 248.80659424902394371d0,&
         & 1029.0736929484210803d0, 3200.5629438795801652d0,&
         & 6656.7749162019607914d0, 6948.8586598121632302d0,&
         & 3.3439394490599745013d-6, 2.2600596902211837757d-5,&
         & 3.3955927589987356838d-4, 3.7105306061050972474d-3,&
         & 3.8065263634919156421d-2, 3.5068223415665236079d-1,&
         & 2.8760027832105027316d0, 20.665999500843274339d0,&
         & 127.47939148516390205d0, 656.43636874254000885d0,&
         & 2709.5242837932479920d0, 8407.1174233600734871d0,&
         & 17437.146284159740233d0, 18141.348781638831600d0 /)
    real(fdp), dimension(0:44), parameter :: c = (/ &
         &-2.8849790431465382128d-15, -3.5125350943844774657d-14, &
         &-7.4850867013707419750d-13, -1.8383904048277485153d-11, &
         &-5.7303556446977223342d-10, -2.4449502737311496525d-8, &
         &-1.6765373351766929724d-6, -3.2189516835265773471d-4,&
         & 5.1503226936425277377d-2, -5.8674306822281631119d-15, &
         &-9.4884898451194085565d-15, -8.5033865136600364340d-13, &
         &-1.8142997866945285736d-11, -5.7340238386338193949d-10, &
         &-2.4449138101742183665d-8, -1.6765375646678855842d-6, &
         &-3.2189516826945356325d-4, 5.1503226936412017608d-2, &
         &-1.4723362506764340882d-14, 1.3945147385179042899d-13, &
         &-1.9618041857586930923d-12, -1.3343606394065121821d-11, &
         &-5.8649674606973244159d-10, -2.4426060539669553778d-8, &
         &-1.6765631828366988006d-6, -3.2189515191449587253d-4,&
         & 5.1503226931820146445d-2, -5.8203519372580372987d-14,&
         & 1.2266326995309845825d-12, -1.3921625844526453237d-11,&
         & 6.2228025878281625469d-11, -8.8636681342142794023d-10, &
         &-2.3661241616744818608d-8, -1.6777870960740520557d-6, &
         &-3.2189402882677074318d-4, 5.1503226479551959376d-2, &
         &-4.5801527369223291722d-14, 6.7998819697143727209d-13, &
         &-4.1624857909290468421d-12, -3.2849009406112440998d-11, &
         &-3.2478275690431118270d-10, -2.5739209934053714983d-8, &
         &-1.6730566573215739195d-6, -3.2190010909008684076d-4,&
         & 5.1503229866932077150d-2 /)

    real(fdp) :: w, t, y
    integer :: k

    w = abs(x)

    if (w  <  8.5d0) then
       t = w * w * 0.0625d0
       k = 12 * (int(t))
       y = (((((((((((a(k) * t + a(k + 1)) * t + a(k + 2)) * t + a(k &
            &+ 3)) * t + a(k + 4)) * t + a(k + 5)) * t + a(k + 6)) *&
            & t + a(k + 7)) * t + a(k + 8)) * t + a(k + 9)) * t + a(k&
            & + 10)) * t + a(k + 11)) * w
    else if (w  <  12.5d0) then
       k = int(w)
       t = w - k
       k = 14 * (k - 8)
       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
       t = 60 / w
       k = 9 * (int(t))
       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)) * sqrt(t) * exp(w)
    end if

    if (x  <  0) y = -y

    f = y
  end function f_besi1



  recursive function f_besin(n, x) result(f)
    !
    ! Modified bessel function 1st kind, order n
    !

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

    integer, parameter :: iacc=40, iexp=maxexponent(x)/2
    integer :: j, m, nn
    real(fdp) :: bi, bim, bip, tox

    nn = nint(n)

    if (nn == 0) then
       f = f_besi0(x)
    else if (nn == 1) then
       f = f_besi1(x)
    else
       f=0.0

       if (x*x <= 8.0d0*tiny(x)) return

       tox=2.0d0/abs(x)
       bip=0.0d0
       bi=1.0d0
       m=2*((nn+int(sqrt(real(iacc*nn)))))

       do j=m,1,-1
          bim=bip+j*tox*bi
          bip=bi
          bi=bim

          if (exponent(bi) > iexp) then
             f=scale(f,-iexp)
             bi=scale(bi,-iexp)
             bip=scale(bip,-iexp)
          end if

          if (j == nn) f=bip
       end do

       f=f*f_besi0(x)/bi

       if (x < 0.0 .and. mod(nn,2) == 1) f=-f
    end if
  end function f_besin



  recursive function f_besk0(x) result(f)
    !
    ! Modified bessel function 2nd kind, order 0
    !

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

    real(fdp), dimension(0:15), parameter :: a = (/&
         & 2.4307270476772195953d-12, 4.7091666363785304370d-10,&
         & 6.7816861334344265568d-8, 6.7816840204737508252d-6,&
         & 4.3402777777915334676d-4, 1.5624999999999872796d-2,&
         & 2.5000000000000000448d-1, 9.9999999999999999997d-1,&
         & 6.5878327432224993071d-12, 1.2083308769932888218d-9,&
         & 1.6271062073716412046d-7, 1.4914719278555277887d-5,&
         & 8.4603509071212245667d-4, 2.5248929932162333910d-2,&
         & 2.7898287891460312491d-1, 1.1593151565841244874d-1 /)
    real(fdp), dimension(0:111), parameter :: b = (/ &
         &-4.6430702971053162197d-13, 1.0377936059563728230d-11, &
         &-1.0298475936392057807d-10, 5.3632747492333959219d-10, &
         &-2.1674628861036068105d-10, -2.3316071545820437669d-8,&
         & 2.2557819578691704059d-7, -9.2325694638587080009d-7, &
         &-3.3569097781613661759d-6, 8.7355061305812582974d-5, &
         &-6.8021202111645760475d-4, 2.7434654781323362319d-4,&
         & 1.0031787169953909561d-1, 4.2102443824070833334d-1,&
         & 4.1447451117883103686d-12, -3.4026589638576604315d-11,&
         & 9.3398790624638977468d-12, 1.5184181750799852630d-9, &
         &-1.1364911665083029464d-8, 2.0619457602095915719d-8,&
         & 3.0431018037572243630d-7, -2.9749736264474555510d-6,&
         & 8.0143661611467038568d-6, 8.0937525149549218398d-5, &
         &-1.0356346549612699886d-3,