module rys use precision, only: dp use constants, only: pi use rys_lut, only: mxrys, maux, nauxs, maprys, xasymp, rts_hermit, wts_hermit, rtsaux, wtsaux implicit none real(dp), parameter :: eps_eig = 1.0e-14_dp real(dp), parameter :: pio4 = pi/4.0_dp type rys_root_t integer :: nroots = -1 real(dp) :: X real(dp) :: U(mxrys), W(mxrys) contains procedure, public :: evaluate, evaluate_testing end type rys_root_t private public rys_root_t contains subroutine evaluate(this) class(rys_root_t), intent(inout) :: this select case (this%nroots) case (1); call rys_rt1(this%x, this%U(1), this%W(1)) case (2); call rys_rt2(this%x, this%U(1:2), this%W(1:2)) case (3); call rys_rt3(this%x, this%U(1:3), this%W(1:3)) case (4); call rys_rt4(this%x, this%U(1:4), this%W(1:4)) case (5); call rys_rt5(this%x, this%U(1:5), this%W(1:5)) case (6:mxrys); call rys_general(this%x, this%U(1:this%nroots), this%W(1:this%nroots), this%nroots) case default error stop "internal error in rys" end select end subroutine evaluate subroutine evaluate_testing(this) class(rys_root_t), intent(inout) :: this this%U = 0.0_dp this%W = 0.0_dp call rys_general(this%x, this%U(1:this%nroots), this%W(1:this%nroots), this%nroots) end subroutine evaluate_testing subroutine rys_rt1(x, r, w) real(KIND=dp), intent(IN) :: & x real(KIND=dp), intent(OUT) :: & r real(KIND=dp), intent(OUT) :: & w real(KIND=dp) :: & f1, y, recx, e if (x <= 3.0e-07_dp) then r = (2.5e+00_dp-x)/(7.5e+00_dp-x) w = 1.0e+00_dp-x/3.0e+00_dp elseif (x <= 1.0e+00_dp) then f1 = ((((((((-8.36313918003957e-08_dp*x+ & 1.21222603512827e-06_dp)*x- & 1.15662609053481e-05_dp)*x+ & 9.25197374512647e-05_dp)*x- & 6.40994113129432e-04_dp)*x+ & 3.78787044215009e-03_dp)*x- & 1.85185172458485e-02_dp)*x+ & 7.14285713298222e-02_dp)*x- & 1.99999999997023e-01_dp)*x+ & 3.33333333333318e-01_dp w = 2*x*f1+exp(-x) r = f1/w elseif (x <= 3.0e+00_dp) then y = x-2.0e+00_dp f1 = ((((((((((-1.61702782425558e-10_dp*y+ & 1.96215250865776e-09_dp)*y- & 2.14234468198419e-08_dp)*y+ & 2.17216556336318e-07_dp)*y- & 1.98850171329371e-06_dp)*y+ & 1.62429321438911e-05_dp)*y- & 1.16740298039895e-04_dp)*y+ & 7.24888732052332e-04_dp)*y- & 3.79490003707156e-03_dp)*y+ & 1.61723488664661e-02_dp)*y- & 5.29428148329736e-02_dp)*y+ & 1.15702180856167e-01_dp w = 2*x*f1+exp(-x) r = f1/w elseif (x <= 5.0e+00_dp) then y = x-4.0e+00_dp f1 = ((((((((((-2.62453564772299e-11_dp*y+ & 3.24031041623823e-10_dp)*y- & 3.614965656163e-09_dp)*y+ & 3.760256799971e-08_dp)*y- & 3.553558319675e-07_dp)*y+ & 3.022556449731e-06_dp)*y- & 2.290098979647e-05_dp)*y+ & 1.526537461148e-04_dp)*y- & 8.81947375894379e-04_dp)*y+ & 4.33207949514611e-03_dp)*y- & 1.75257821619926e-02_dp)*y+ & 5.28406320615584e-02_dp w = 2*x*f1+exp(-x) r = f1/w elseif (x <= 10.0e+00_dp) then e = exp(-x) recx = 1.0e+00_dp/x w = ((((((4.6897511375022e-01_dp*recx- & 6.9955602298985e-01_dp)*recx+ & 5.3689283271887e-01_dp)*recx- & 3.2883030418398e-01_dp)*recx+ & 2.4645596956002e-01_dp)*recx- & 4.9984072848436e-01_dp)*recx- & 3.1501078774085e-06_dp)*e+sqrt(PIo4*recx) f1 = recx*(w-e)/2 r = f1/w elseif (x <= 15.0e+00_dp) then e = exp(-x) recx = 1.0e+00_dp/x w = (((-1.8784686463512e-01_dp*recx+ & 2.2991849164985e-01_dp)*recx- & 4.9893752514047e-01_dp)*recx- & 2.1916512131607e-05_dp)*e+sqrt(PIo4*recx) f1 = recx*(w-e)/2 r = f1/w elseif (x <= 33.0e+00_dp) then e = exp(-x) recx = 1.0e+00_dp/x w = ((1.9623264149430e-01_dp*recx- & 4.9695241464490e-01_dp)*recx- & 6.0156581186481e-05_dp)*e+sqrt(PIo4*recx) f1 = recx*(w-e)/2 r = f1/w else recx = 1.0e+00_dp/x w = sqrt(PIo4*recx) r = 0.5e+00_dp*recx end if r = r/(1.0e00_dp-r) end subroutine rys_rt1 !----------------------------------------------------------------- subroutine rys_rt2(x, r, w) real(KIND=dp), intent(IN) :: & x real(KIND=dp), intent(OUT) :: & r(2), w(2) real(KIND=dp), parameter :: & r12 = 2.75255128608411e-01_dp, & r22 = 2.72474487139158e+00_dp, & w22 = 9.17517095361369e-02_dp real(KIND=dp) :: & f1, y, recx, e, r1, r2, w1, w2 if (x <= 3.0e-07_dp) then r1 = 1.30693606237085e-01_dp-2.90430236082028e-02_dp*x r2 = 2.86930639376291e+00_dp-6.37623643058102e-01_dp*x w1 = 6.52145154862545e-01_dp-1.22713621927067e-01_dp*x w2 = 3.47854845137453e-01_dp-2.10619711404725e-01_dp*x elseif (x <= 1.0e+00_dp) then f1 = ((((((((-8.36313918003957e-08_dp*x+ & 1.21222603512827e-06_dp)*x- & 1.15662609053481e-05_dp)*x+ & 9.25197374512647e-05_dp)*x- & 6.40994113129432e-04_dp)*x+ & 3.78787044215009e-03_dp)*x- & 1.85185172458485e-02_dp)*x+ & 7.14285713298222e-02_dp)*x- & 1.99999999997023e-01_dp)*x+ & 3.33333333333318e-01_dp w1 = 2*x*f1+exp(-x) r1 = (((((((-2.35234358048491e-09_dp*x+ & 2.49173650389842e-08_dp)*x- & 4.558315364581e-08_dp)*x- & 2.447252174587e-06_dp)*x+ & 4.743292959463e-05_dp)*x- & 5.33184749432408e-04_dp)*x+ & 4.44654947116579e-03_dp)*x- & 2.90430236084697e-02_dp)*x+ & 1.30693606237085e-01_dp r2 = (((((((-2.47404902329170e-08_dp*x+ & 2.36809910635906e-07_dp)*x+ & 1.835367736310e-06_dp)*x- & 2.066168802076e-05_dp)*x- & 1.345693393936e-04_dp)*x- & 5.88154362858038e-05_dp)*x+ & 5.32735082098139e-02_dp)*x- & 6.37623643056745e-01_dp)*x+ & 2.86930639376289e+00_dp w2 = ((f1-w1)*r1+f1)*(1.0e+00_dp+r2)/(r2-r1) w1 = w1-w2 elseif (x <= 3.0e+00_dp) then y = x-2.0e+00_dp f1 = ((((((((((-1.61702782425558e-10_dp*y+ & 1.96215250865776e-09_dp)*y- & 2.14234468198419e-08_dp)*y+ & 2.17216556336318e-07_dp)*y- & 1.98850171329371e-06_dp)*y+ & 1.62429321438911e-05_dp)*y- & 1.16740298039895e-04_dp)*y+ & 7.24888732052332e-04_dp)*y- & 3.79490003707156e-03_dp)*y+ & 1.61723488664661e-02_dp)*y- & 5.29428148329736e-02_dp)*y+ & 1.15702180856167e-01_dp w1 = 2*x*f1+exp(-x) r1 = (((((((((-6.36859636616415e-12_dp*y+ & 8.47417064776270e-11_dp)*y- & 5.152207846962e-10_dp)*y- & 3.846389873308e-10_dp)*y+ & 8.472253388380e-08_dp)*y- & 1.85306035634293e-06_dp)*y+ & 2.47191693238413e-05_dp)*y- & 2.49018321709815e-04_dp)*y+ & 2.19173220020161e-03_dp)*y- & 1.63329339286794e-02_dp)*y+ & 8.68085688285261e-02_dp r2 = (((((((((1.45331350488343e-10_dp*y+ & 2.07111465297976e-09_dp)*y- & 1.878920917404e-08_dp)*y- & 1.725838516261e-07_dp)*y+ & 2.247389642339e-06_dp)*y+ & 9.76783813082564e-06_dp)*y- & 1.93160765581969e-04_dp)*y- & 1.58064140671893e-03_dp)*y+ & 4.85928174507904e-02_dp)*y- & 4.30761584997596e-01_dp)*y+ & 1.80400974537950e+00_dp w2 = ((f1-w1)*r1+f1)*(1.0e+00_dp+r2)/(r2-r1) w1 = w1-w2 elseif (x <= 5.0e+00_dp) then y = x-4.0e+00_dp f1 = ((((((((((-2.62453564772299e-11_dp*y+ & 3.24031041623823e-10_dp)*y- & 3.614965656163e-09_dp)*y+ & 3.760256799971e-08_dp)*y- & 3.553558319675e-07_dp)*y+ & 3.022556449731e-06_dp)*y- & 2.290098979647e-05_dp)*y+ & 1.526537461148e-04_dp)*y- & 8.81947375894379e-04_dp)*y+ & 4.33207949514611e-03_dp)*y- & 1.75257821619926e-02_dp)*y+ & 5.28406320615584e-02_dp w1 = 2*x*f1+exp(-x) r1 = ((((((((-4.11560117487296e-12_dp*y+ & 7.10910223886747e-11_dp)*y- & 1.73508862390291e-09_dp)*y+ & 5.93066856324744e-08_dp)*y- & 9.76085576741771e-07_dp)*y+ & 1.08484384385679e-05_dp)*y- & 1.12608004981982e-04_dp)*y+ & 1.16210907653515e-03_dp)*y- & 9.89572595720351e-03_dp)*y+ & 6.12589701086408e-02_dp r2 = (((((((((-1.80555625241001e-10_dp*y+ & 5.44072475994123e-10_dp)*y+ & 1.603498045240e-08_dp)*y- & 1.497986283037e-07_dp)*y- & 7.017002532106e-07_dp)*y+ & 1.85882653064034e-05_dp)*y- & 2.04685420150802e-05_dp)*y- & 2.49327728643089e-03_dp)*y+ & 3.56550690684281e-02_dp)*y- & 2.60417417692375e-01_dp)*y+ & 1.12155283108289e+00_dp w2 = ((f1-w1)*r1+f1)*(1.0e+00_dp+r2)/(r2-r1) w1 = w1-w2 elseif (x <= 10.0e+00_dp) then e = exp(-x) recx = 1.0e+00_dp/x w1 = ((((((4.6897511375022e-01_dp*recx- & 6.9955602298985e-01_dp)*recx+ & 5.3689283271887e-01_dp)*recx- & 3.2883030418398e-01_dp)*recx+ & 2.4645596956002e-01_dp)*recx- & 4.9984072848436e-01_dp)*recx- & 3.1501078774085e-06_dp)*e+sqrt(PIo4*recx) f1 = recx*(w1-e)/2 y = x-7.5e+00_dp r1 = (((((((((((((-1.43632730148572e-16_dp*y+ & 2.38198922570405e-16_dp)*y+ & 1.358319618800e-14_dp)*y- & 7.064522786879e-14_dp)*y- & 7.719300212748e-13_dp)*y+ & 7.802544789997e-12_dp)*y+ & 6.628721099436e-11_dp)*y- & 1.775564159743e-09_dp)*y+ & 1.713828823990e-08_dp)*y- & 1.497500187053e-07_dp)*y+ & 2.283485114279e-06_dp)*y- & 3.76953869614706e-05_dp)*y+ & 4.74791204651451e-04_dp)*y- & 4.60448960876139e-03_dp)*y+ & 3.72458587837249e-02_dp r2 = ((((((((((((2.48791622798900e-14_dp*y- & 1.36113510175724e-13_dp)*y- & 2.224334349799e-12_dp)*y+ & 4.190559455515e-11_dp)*y- & 2.222722579924e-10_dp)*y- & 2.624183464275e-09_dp)*y+ & 6.128153450169e-08_dp)*y- & 4.383376014528e-07_dp)*y- & 2.49952200232910e-06_dp)*y+ & 1.03236647888320e-04_dp)*y- & 1.44614664924989e-03_dp)*y+ & 1.35094294917224e-02_dp)*y- & 9.53478510453887e-02_dp)*y+ & 5.44765245686790e-01_dp w2 = ((f1-w1)*r1+f1)*(1.0e+00_dp+r2)/(r2-r1) w1 = w1-w2 elseif (x <= 15.0e+00_dp) then e = exp(-x) recx = 1.0e+00_dp/x w1 = (((-1.8784686463512e-01_dp*recx+ & 2.2991849164985e-01_dp)*recx- & 4.9893752514047e-01_dp)*recx- & 2.1916512131607e-05_dp)*e+ & sqrt(PIo4*recx) f1 = recx*(w1-e)/2 r1 = ((((-1.01041157064226e-05_dp*x+ & 1.19483054115173e-03_dp)*x- & 6.73760231824074e-02_dp)*x+ & 1.25705571069895e+00_dp)*x+ & (((-8.57609422987199e+03_dp*recx+ & 5.91005939591842e+03_dp)*recx- & 1.70807677109425e+03_dp)*recx+ & 2.64536689959503e+02_dp)*recx- & 2.38570496490846e+01_dp)*e+r12/(x-r12) r2 = (((3.39024225137123e-04_dp*x- & 9.34976436343509e-02_dp)*x- & 4.22216483306320e+00_dp)*x+ & (((-2.08457050986847e+03_dp*recx- & 1.04999071905664e+03_dp)*recx+ & 3.39891508992661e+02_dp)*recx- & 1.56184800325063e+02_dp)*recx+ & 8.00839033297501e+00_dp)*e+ & r22/(x-r22) w2 = ((f1-w1)*r1+f1)* & (1.0e+00_dp+r2)/(r2-r1) w1 = w1-w2 elseif (x <= 33.0e+00_dp) then e = exp(-x) recx = 1.0e+00_dp/x w1 = ((1.9623264149430e-01_dp*recx- & 4.9695241464490e-01_dp)*recx- & 6.0156581186481e-05_dp)*e+sqrt(PIo4*recx) f1 = recx*(w1-e)/2 r1 = ((((-1.14906395546354e-06_dp*x+ & 1.76003409708332e-04_dp)*x- & 1.71984023644904e-02_dp)*x- & 1.37292644149838e-01_dp)*x+ & (-4.75742064274859e+01_dp*recx+ & 9.21005186542857e+00_dp)*recx- & 2.31080873898939e-02_dp)*e+r12/(x-r12) r2 = (((3.64921633404158e-04_dp*x- & 9.71850973831558e-02_dp)*x- & 4.02886174850252e+00_dp)*x+ & (-1.35831002139173e+02_dp*recx- & 8.66891724287962e+01_dp)*recx+ & 2.98011277766958e+00_dp)*e+r22/(x-r22) w2 = ((f1-w1)*r1+f1)* & (1.0e+00_dp+r2)/(r2-r1) w1 = w1-w2 elseif (x <= 40.0e+00_dp) then w1 = sqrt(PIo4/x) e = exp(-x) r1 = (-8.78947307498880e-01_dp*x+ & 1.09243702330261e+01_dp)*e+r12/(x-r12) r2 = (-9.28903924275977e+00_dp*x+ & 8.10642367843811e+01_dp)*e+r22/(x-r22) w2 = (4.46857389308400e+00_dp*x- & 7.79250653461045e+01_dp)*e+w22*w1 w1 = w1-w2 else r1 = r12/(x-r12) r2 = r22/(x-r22) w1 = sqrt(PIo4/x) w2 = w22*w1 w1 = w1-w2 end if r = (/r1, r2/) w = (/w1, w2/) end subroutine rys_rt2 !----------------------------------------------------------------- subroutine rys_rt3(x, r, w) real(KIND=dp), intent(IN) :: & x real(KIND=dp), intent(OUT) :: & r(3), w(3) real(KIND=dp), parameter :: & r13 = 1.90163509193487e-01_dp, & r23 = 1.78449274854325e+00_dp, & w23 = 1.77231492083829e-01_dp, & r33 = 5.52534374226326e+00_dp, & w33 = 5.11156880411248e-03_dp real(KIND=dp) :: & f1, f2, y, recx, e, & a1, a2, t1, t2, t3, r1, r2, r3, w1, w2, w3 if (x <= 3.0e-07_dp) then r1 = 6.03769246832797e-02_dp- & 9.28875764357368e-03_dp*x r2 = 7.76823355931043e-01_dp- & 1.19511285527878e-01_dp*x r3 = 6.66279971938567e+00_dp- & 1.02504611068957e+00_dp*x w1 = 4.67913934572691e-01_dp- & 5.64876917232519e-02_dp*x w2 = 3.60761573048137e-01_dp- & 1.49077186455208e-01_dp*x w3 = 1.71324492379169e-01_dp- & 1.27768455150979e-01_dp*x elseif (x <= 1.0e+00_dp) then r1 = ((((((-5.10186691538870e-10_dp*x+ & 2.40134415703450e-08_dp)*x- & 5.01081057744427e-07_dp)*x+ & 7.58291285499256e-06_dp)*x- & 9.55085533670919e-05_dp)*x+ & 1.02893039315878e-03_dp)*x- & 9.28875764374337e-03_dp)*x+ & 6.03769246832810e-02_dp r2 = ((((((-1.29646524960555e-08_dp*x+ & 7.74602292865683e-08_dp)*x+ & 1.56022811158727e-06_dp)*x- & 1.58051990661661e-05_dp)*x- & 3.30447806384059e-04_dp)*x+ & 9.74266885190267e-03_dp)*x- & 1.19511285526388e-01_dp)*x+ & 7.76823355931033e-01_dp r3 = ((((((-9.28536484109606e-09_dp*x- & 3.02786290067014e-07_dp)*x- & 2.50734477064200e-06_dp)*x- & 7.32728109752881e-06_dp)*x+ & 2.44217481700129e-04_dp)*x+ & 4.94758452357327e-02_dp)*x- & 1.02504611065774e+00_dp)*x+ & 6.66279971938553e+00_dp f2 = ((((((((-7.60911486098850e-08_dp*x+ & 1.09552870123182e-06_dp)*x- & 1.03463270693454e-05_dp)*x+ & 8.16324851790106e-05_dp)*x- & 5.55526624875562e-04_dp)*x+ & 3.20512054753924e-03_dp)*x- & 1.51515139838540e-02_dp)*x+ & 5.55555554649585e-02_dp)*x- & 1.42857142854412e-01_dp)*x+ & 1.99999999999986e-01_dp e = exp(-x) f1 = (2*x*f2+e)/3.0e+00_dp w1 = 2*x*f1+e t1 = r1/(r1+1.0e+00_dp) t2 = r2/(r2+1.0e+00_dp) t3 = r3/(r3+1.0e+00_dp) a2 = f2-t1*f1 a1 = f1-t1*w1 w3 = (a2-t2*a1)/((t3-t2)*(t3-t1)) w2 = (t3*a1-a2)/((t3-t2)*(t2-t1)) w1 = w1-w2-w3 elseif (x <= 3.0e+00_dp) then y = x-2.0e+00_dp r1 = ((((((((1.44687969563318e-12_dp*y+ & 4.85300143926755e-12_dp)*y- & 6.55098264095516e-10_dp)*y+ & 1.56592951656828e-08_dp)*y- & 2.60122498274734e-07_dp)*y+ & 3.86118485517386e-06_dp)*y- & 5.13430986707889e-05_dp)*y+ & 6.03194524398109e-04_dp)*y- & 6.11219349825090e-03_dp)*y+ & 4.52578254679079e-02_dp r2 = (((((((6.95964248788138e-10_dp*y- & 5.35281831445517e-09_dp)*y- & 6.745205954533e-08_dp)*y+ & 1.502366784525e-06_dp)*y+ & 9.923326947376e-07_dp)*y- & 3.89147469249594e-04_dp)*y+ & 7.51549330892401e-03_dp)*y- & 8.48778120363400e-02_dp)*y+ & 5.73928229597613e-01_dp r3 = ((((((((-2.81496588401439e-10_dp*y+ & 3.61058041895031e-09_dp)*y+ & 4.53631789436255e-08_dp)*y- & 1.40971837780847e-07_dp)*y- & 6.05865557561067e-06_dp)*y- & 5.15964042227127e-05_dp)*y+ & 3.34761560498171e-05_dp)*y+ & 5.04871005319119e-02_dp)*y- & 8.24708946991557e-01_dp)*y+ & 4.81234667357205e+00_dp f2 = ((((((((((-1.48044231072140e-10_dp*y+ & 1.78157031325097e-09_dp)*y- & 1.92514145088973e-08_dp)*y+ & 1.92804632038796e-07_dp)*y- & 1.73806555021045e-06_dp)*y+ & 1.39195169625425e-05_dp)*y- & 9.74574633246452e-05_dp)*y+ & 5.83701488646511e-04_dp)*y- & 2.89955494844975e-03_dp)*y+ & 1.13847001113810e-02_dp)*y- & 3.23446977320647e-02_dp)*y+ & 5.29428148329709e-02_dp e = exp(-x) f1 = (2*x*f2+e)/3.0e+00_dp w1 = 2*x*f1+e t1 = r1/(r1+1.0e+00_dp) t2 = r2/(r2+1.0e+00_dp) t3 = r3/(r3+1.0e+00_dp) a2 = f2-t1*f1 a1 = f1-t1*w1 w3 = (a2-t2*a1)/((t3-t2)*(t3-t1)) w2 = (t3*a1-a2)/((t3-t2)*(t2-t1)) w1 = w1-w2-w3 elseif (x <= 5.0e+00_dp) then y = x-4.0e+00_dp r1 = (((((((1.44265709189601e-11_dp*y- & 4.66622033006074e-10_dp)*y+ & 7.649155832025e-09_dp)*y- & 1.229940017368e-07_dp)*y+ & 2.026002142457e-06_dp)*y- & 2.87048671521677e-05_dp)*y+ & 3.70326938096287e-04_dp)*y- & 4.21006346373634e-03_dp)*y+ & 3.50898470729044e-02_dp r2 = ((((((((-2.65526039155651e-11_dp*y+ & 1.97549041402552e-10_dp)*y+ & 2.15971131403034e-09_dp)*y- & 7.95045680685193e-08_dp)*y+ & 5.15021914287057e-07_dp)*y+ & 1.11788717230514e-05_dp)*y- & 3.33739312603632e-04_dp)*y+ & 5.30601428208358e-03_dp)*y- & 5.93483267268959e-02_dp)*y+ & 4.31180523260239e-01_dp r3 = ((((((((-3.92833750584041e-10_dp*y- & 4.16423229782280e-09_dp)*y+ & 4.42413039572867e-08_dp)*y+ & 6.40574545989551e-07_dp)*y- & 3.05512456576552e-06_dp)*y- & 1.05296443527943e-04_dp)*y- & 6.14120969315617e-04_dp)*y+ & 4.89665802767005e-02_dp)*y- & 6.24498381002855e-01_dp)*y+ & 3.36412312243724e+00_dp f2 = ((((((((((-2.36788772599074e-11_dp*y+ & 2.89147476459092e-10_dp)*y- & 3.18111322308846e-09_dp)*y+ & 3.25336816562485e-08_dp)*y- & 3.00873821471489e-07_dp)*y+ & 2.48749160874431e-06_dp)*y- & 1.81353179793672e-05_dp)*y+ & 1.14504948737066e-04_dp)*y- & 6.10614987696677e-04_dp)*y+ & 2.64584212770942e-03_dp)*y- & 8.66415899015349e-03_dp)*y+ & 1.75257821619922e-02_dp e = exp(-x) f1 = ((x+x)*f2+e)/3.0e+00_dp w1 = (x+x)*f1+e t1 = r1/(r1+1.0e+00_dp) t2 = r2/(r2+1.0e+00_dp) t3 = r3/(r3+1.0e+00_dp) a2 = f2-t1*f1 a1 = f1-t1*w1 w3 = (a2-t2*a1)/((t3-t2)*(t3-t1)) w2 = (t3*a1-a2)/((t3-t2)*(t2-t1)) w1 = w1-w2-w3 elseif (x <= 10.0e+00_dp) then e = exp(-x) recx = 1.0e+00_dp/x w1 = ((((((4.6897511375022e-01_dp*recx- & 6.9955602298985e-01_dp)*recx+ & 5.3689283271887e-01_dp)*recx- & 3.2883030418398e-01_dp)*recx+ & 2.4645596956002e-01_dp)*recx- & 4.9984072848436e-01_dp)*recx- & 3.1501078774085e-06_dp)*e+sqrt(PIo4*recx) f1 = recx*(w1-e)/2 f2 = recx*(f1+f1+f1-e)/2 y = x-7.5e+00_dp r1 = (((((((((((5.74429401360115e-16_dp*y+ & 7.11884203790984e-16_dp)*y- & 6.736701449826e-14_dp)*y- & 6.264613873998e-13_dp)*y+ & 1.315418927040e-11_dp)*y- & 4.23879635610964e-11_dp)*y+ & 1.39032379769474e-09_dp)*y- & 4.65449552856856e-08_dp)*y+ & 7.34609900170759e-07_dp)*y- & 1.08656008854077e-05_dp)*y+ & 1.77930381549953e-04_dp)*y- & 2.39864911618015e-03_dp)*y+ & 2.39112249488821e-02_dp r2 = (((((((((((1.13464096209120e-14_dp*y+ & 6.99375313934242e-15_dp)*y- & 8.595618132088e-13_dp)*y- & 5.293620408757e-12_dp)*y- & 2.492175211635e-11_dp)*y+ & 2.73681574882729e-09_dp)*y- & 1.06656985608482e-08_dp)*y- & 4.40252529648056e-07_dp)*y+ & 9.68100917793911e-06_dp)*y- & 1.68211091755327e-04_dp)*y+ & 2.69443611274173e-03_dp)*y- & 3.23845035189063e-02_dp)*y+ & 2.75969447451882e-01_dp r3 = ((((((((((((6.66339416996191e-15_dp*y+ & 1.84955640200794e-13_dp)*y- & 1.985141104444e-12_dp)*y- & 2.309293727603e-11_dp)*y+ & 3.917984522103e-10_dp)*y+ & 1.663165279876e-09_dp)*y- & 6.205591993923e-08_dp)*y+ & 8.769581622041e-09_dp)*y+ & 8.97224398620038e-06_dp)*y- & 3.14232666170796e-05_dp)*y- & 1.83917335649633e-03_dp)*y+ & 3.51246831672571e-02_dp)*y- & 3.22335051270860e-01_dp)*y+ & 1.73582831755430e+00_dp t1 = r1/(r1+1.0e+00_dp) t2 = r2/(r2+1.0e+00_dp) t3 = r3/(r3+1.0e+00_dp) a2 = f2-t1*f1 a1 = f1-t1*w1 w3 = (a2-t2*a1)/((t3-t2)*(t3-t1)) w2 = (t3*a1-a2)/((t3-t2)*(t2-t1)) w1 = w1-w2-w3 elseif (x <= 15.0e+00_dp) then e = exp(-x) recx = 1.0e+00_dp/x w1 = (((-1.8784686463512e-01_dp*recx+ & 2.2991849164985e-01_dp)*recx- & 4.9893752514047e-01_dp)*recx- & 2.1916512131607e-05_dp)*e+sqrt(PIo4*recx) f1 = recx*(w1-e)/2 f2 = recx*(f1+f1+f1-e)/2 y = x-12.5e+00_dp r1 = (((((((((((4.42133001283090e-16_dp*y- & 2.77189767070441e-15_dp)*y- & 4.084026087887e-14_dp)*y+ & 5.379885121517e-13_dp)*y+ & 1.882093066702e-12_dp)*y- & 8.67286219861085e-11_dp)*y+ & 7.11372337079797e-10_dp)*y- & 3.55578027040563e-09_dp)*y+ & 1.29454702851936e-07_dp)*y- & 4.14222202791434e-06_dp)*y+ & 8.04427643593792e-05_dp)*y- & 1.18587782909876e-03_dp)*y+ & 1.53435577063174e-02_dp r2 = (((((((((((6.85146742119357e-15_dp*y- & 1.08257654410279e-14_dp)*y- & 8.579165965128e-13_dp)*y+ & 6.642452485783e-12_dp)*y+ & 4.798806828724e-11_dp)*y- & 1.13413908163831e-09_dp)*y+ & 7.08558457182751e-09_dp)*y- & 5.59678576054633e-08_dp)*y+ & 2.51020389884249e-06_dp)*y- & 6.63678914608681e-05_dp)*y+ & 1.11888323089714e-03_dp)*y- & 1.45361636398178e-02_dp)*y+ & 1.65077877454402e-01_dp r3 = ((((((((((((3.20622388697743e-15_dp*y- & 2.73458804864628e-14_dp)*y- & 3.157134329361e-13_dp)*y+ & 8.654129268056e-12_dp)*y- & 5.625235879301e-11_dp)*y- & 7.718080513708e-10_dp)*y+ & 2.064664199164e-08_dp)*y- & 1.567725007761e-07_dp)*y- & 1.57938204115055e-06_dp)*y+ & 6.27436306915967e-05_dp)*y- & 1.01308723606946e-03_dp)*y+ & 1.13901881430697e-02_dp)*y- & 1.01449652899450e-01_dp)*y+ & 7.77203937334739e-01_dp t1 = r1/(r1+1.0e+00_dp) t2 = r2/(r2+1.0e+00_dp) t3 = r3/(r3+1.0e+00_dp) a2 = f2-t1*f1 a1 = f1-t1*w1 w3 = (a2-t2*a1)/((t3-t2)*(t3-t1)) w2 = (t3*a1-a2)/((t3-t2)*(t2-t1)) w1 = w1-w2-w3 elseif (x <= 20.0e+00_dp) then e = exp(-x) recx = 1.0e+00_dp/x w1 = ((1.9623264149430e-01_dp*recx- & 4.9695241464490e-01_dp)*recx- & 6.0156581186481e-05_dp)*e+sqrt(PIo4*recx) f1 = recx*(w1-e)/2 f2 = recx*(f1+f1+f1-e)/2 r1 = ((((((-2.43270989903742e-06_dp*x+ & 3.57901398988359e-04_dp)*x- & 2.34112415981143e-02_dp)*x+ & 7.81425144913975e-01_dp)*x- & 1.73209218219175e+01_dp)*x+ & 2.43517435690398e+02_dp)*x+ & (-1.97611541576986e+04_dp*recx+ & 9.82441363463929e+03_dp)*recx- & 2.07970687843258e+03_dp)*e+r13/(x-r13) r2 = (((((-2.62627010965435e-04_dp*x+ & 3.49187925428138e-02_dp)*x- & 3.09337618731880e+00_dp)*x+ & 1.07037141010778e+02_dp)*x- & 2.36659637247087e+03_dp)*x+ & ((-2.91669113681020e+06_dp*recx+ & 1.41129505262758e+06_dp)*recx- & 2.91532335433779e+05_dp)*recx+ & 3.35202872835409e+04_dp)*e+r23/(x-r23) r3 = (((((9.31856404738601e-05_dp*x- & 2.87029400759565e-02_dp)*x- & 7.83503697918455e-01_dp)*x- & 1.84338896480695e+01_dp)*x+ & 4.04996712650414e+02_dp)*x+ & (-1.89829509315154e+05_dp*recx+ & 5.11498390849158e+04_dp)*recx- & 6.88145821789955e+03_dp)*e+r33/(x-r33) t1 = r1/(r1+1.0e+00_dp) t2 = r2/(r2+1.0e+00_dp) t3 = r3/(r3+1.0e+00_dp) a2 = f2-t1*f1 a1 = f1-t1*w1 w3 = (a2-t2*a1)/((t3-t2)*(t3-t1)) w2 = (t3*a1-a2)/((t3-t2)*(t2-t1)) w1 = w1-w2-w3 elseif (x <= 33.0e+00_dp) then e = exp(-x) recx = 1.0e+00_dp/x w1 = ((1.9623264149430e-01_dp*recx- & 4.9695241464490e-01_dp)*recx- & 6.0156581186481e-05_dp)*e+sqrt(PIo4*recx) f1 = recx*(w1-e)/2 f2 = recx*(f1+f1+f1-e)/2 r1 = ((((-4.97561537069643e-04_dp*x- & 5.00929599665316e-02_dp)*x+ & 1.31099142238996e+00_dp)*x- & 1.88336409225481e+01_dp)*x- & 6.60344754467191e+02_dp*recx+ & 1.64931462413877e+02_dp)*e+r13/(x-r13) r2 = ((((-4.48218898474906e-03_dp*x- & 5.17373211334924e-01_dp)*x+ & 1.13691058739678e+01_dp)*x- & 1.65426392885291e+02_dp)*x- & 6.30909125686731e+03_dp*recx+ & 1.52231757709236e+03_dp)*e+r23/(x-r23) r3 = ((((-1.38368602394293e-02_dp*x- & 1.77293428863008e+00_dp)*x+ & 1.73639054044562e+01_dp)*x- & 3.57615122086961e+02_dp)*x- & 1.45734701095912e+04_dp*recx+ & 2.69831813951849e+03_dp)*e+r33/(x-r33) t1 = r1/(r1+1.0e+00_dp) t2 = r2/(r2+1.0e+00_dp) t3 = r3/(r3+1.0e+00_dp) a2 = f2-t1*f1 a1 = f1-t1*w1 w3 = (a2-t2*a1)/((t3-t2)*(t3-t1)) w2 = (t3*a1-a2)/((t3-t2)*(t2-t1)) w1 = w1-w2-w3 elseif (x <= 47.0e+00_dp) then w1 = sqrt(PIo4/x) e = exp(-x) r1 = ((-7.39058467995275e+00_dp*x+ & 3.21318352526305e+02_dp)*x- & 3.99433696473658e+03_dp)*e+r13/(x-r13) r2 = ((-7.38726243906513e+01_dp*x+ & 3.13569966333873e+03_dp)*x- & 3.86862867311321e+04_dp)*e+r23/(x-r23) r3 = ((-2.63750565461336e+02_dp*x+ & 1.04412168692352e+04_dp)*x- & 1.28094577915394e+05_dp)*e+r33/(x-r33) w3 = (((1.52258947224714e-01_dp*x- & 8.30661900042651e+00_dp)*x+ & 1.92977367967984e+02_dp)*x- & 1.67787926005344e+03_dp)*e+w33*w1 w2 = ((6.15072615497811e+01_dp*x- & 2.91980647450269e+03_dp)*x+ & 3.80794303087338e+04_dp)*e+w23*w1 w1 = w1-w2-w3 else r1 = r13/(x-r13) r2 = r23/(x-r23) r3 = r33/(x-r33) w1 = sqrt(PIo4/x) w2 = w23*w1 w3 = w33*w1 w1 = w1-w2-w3 end if r = (/r1, r2, r3/) w = (/w1, w2, w3/) end subroutine rys_rt3 !----------------------------------------------------------------- subroutine rys_rt4(x, r, w) real(KIND=dp), intent(IN) :: & x real(KIND=dp), intent(OUT) :: & r(4), w(4) real(KIND=dp), parameter :: & r14 = 1.45303521503316e-01_dp, & r24 = 1.33909728812636e+00_dp, & w24 = 2.34479815323517e-01_dp, & r34 = 3.92696350135829e+00_dp, & w34 = 1.92704402415764e-02_dp, & r44 = 8.58863568901199e+00_dp, & w44 = 2.25229076750736e-04_dp real(KIND=dp) :: & y, recx, e, r1, r2, r3, r4, w1, w2, w3, w4 if (x <= 3.0e-07_dp) then r1 = 3.48198973061471e-02_dp-4.09645850660395e-03_dp*x r2 = 3.81567185080042e-01_dp-4.48902570656719e-02_dp*x r3 = 1.73730726945891e+00_dp-2.04389090547327e-01_dp*x r4 = 1.18463056481549e+01_dp-1.39368301742312e+00_dp*x w1 = 3.62683783378362e-01_dp-3.13844305713928e-02_dp*x w2 = 3.13706645877886e-01_dp-8.98046242557724e-02_dp*x w3 = 2.22381034453372e-01_dp-1.29314370958973e-01_dp*x w4 = 1.01228536290376e-01_dp-8.28299075414321e-02_dp*x elseif (x <= 1.0e+00_dp) then r1 = ((((((-1.95309614628539e-10_dp*x+ & 5.19765728707592e-09_dp)*x- & 1.01756452250573e-07_dp)*x+ & 1.72365935872131e-06_dp)*x- & 2.61203523522184e-05_dp)*x+ & 3.52921308769880e-04_dp)*x- & 4.09645850658433e-03_dp)*x+ & 3.48198973061469e-02_dp r2 = (((((-1.89554881382342e-08_dp*x+ & 3.07583114342365e-07_dp)*x+ & 1.270981734393e-06_dp)*x- & 1.417298563884e-04_dp)*x+ & 3.226979163176e-03_dp)*x- & 4.48902570678178e-02_dp)*x+ & 3.81567185080039e-01_dp r3 = ((((((1.77280535300416e-09_dp*x+ & 3.36524958870615e-08_dp)*x- & 2.58341529013893e-07_dp)*x- & 1.13644895662320e-05_dp)*x- & 7.91549618884063e-05_dp)*x+ & 1.03825827346828e-02_dp)*x- & 2.04389090525137e-01_dp)*x+ & 1.73730726945889e+00_dp r4 = (((((-5.61188882415248e-08_dp*x- & 2.49480733072460e-07_dp)*x+ & 3.428685057114e-06_dp)*x+ & 1.679007454539e-04_dp)*x+ & 4.722855585715e-02_dp)*x- & 1.39368301737828e+00_dp)*x+ & 1.18463056481543e+01_dp w1 = ((((((-1.14649303201279e-08_dp*x+ & 1.88015570196787e-07_dp)*x- & 2.33305875372323e-06_dp)*x+ & 2.68880044371597e-05_dp)*x- & 2.94268428977387e-04_dp)*x+ & 3.06548909776613e-03_dp)*x- & 3.13844305680096e-02_dp)*x+ & 3.62683783378335e-01_dp w2 = ((((((((-4.11720483772634e-09_dp*x+ & 6.54963481852134e-08_dp)*x- & 7.20045285129626e-07_dp)*x+ & 6.93779646721723e-06_dp)*x- & 6.05367572016373e-05_dp)*x+ & 4.74241566251899e-04_dp)*x- & 3.26956188125316e-03_dp)*x+ & 1.91883866626681e-02_dp)*x- & 8.98046242565811e-02_dp)*x+ & 3.13706645877886e-01_dp w3 = ((((((((-3.41688436990215e-08_dp*x+ & 5.07238960340773e-07_dp)*x- & 5.01675628408220e-06_dp)*x+ & 4.20363420922845e-05_dp)*x- & 3.08040221166823e-04_dp)*x+ & 1.94431864731239e-03_dp)*x- & 1.02477820460278e-02_dp)*x+ & 4.28670143840073e-02_dp)*x- & 1.29314370962569e-01_dp)*x+ & 2.22381034453369e-01_dp w4 = (((((((((4.99660550769508e-09_dp*x- & 7.94585963310120e-08_dp)*x+ & 8.359072409485e-07_dp)*x- & 7.422369210610e-06_dp)*x+ & 5.763374308160e-05_dp)*x- & 3.86645606718233e-04_dp)*x+ & 2.18417516259781e-03_dp)*x- & 9.99791027771119e-03_dp)*x+ & 3.48791097377370e-02_dp)*x- & 8.28299075413889e-02_dp)*x+ & 1.01228536290376e-01_dp elseif (x <= 5.0e+00_dp) then y = x-3.0e+00_dp r1 = (((((((((-1.48570633747284e-15_dp*y- & 1.33273068108777e-13_dp)*y+ & 4.068543696670e-12_dp)*y- & 9.163164161821e-11_dp)*y+ & 2.046819017845e-09_dp)*y- & 4.03076426299031e-08_dp)*y+ & 7.29407420660149e-07_dp)*y- & 1.23118059980833e-05_dp)*y+ & 1.88796581246938e-04_dp)*y- & 2.53262912046853e-03_dp)*y+ & 2.51198234505021e-02_dp r2 = (((((((((1.35830583483312e-13_dp*y- & 2.29772605964836e-12_dp)*y- & 3.821500128045e-12_dp)*y+ & 6.844424214735e-10_dp)*y- & 1.048063352259e-08_dp)*y+ & 1.50083186233363e-08_dp)*y+ & 3.48848942324454e-06_dp)*y- & 1.08694174399193e-04_dp)*y+ & 2.08048885251999e-03_dp)*y- & 2.91205805373793e-02_dp)*y+ & 2.72276489515713e-01_dp r3 = (((((((((5.02799392850289e-13_dp*y+ & 1.07461812944084e-11_dp)*y- & 1.482277886411e-10_dp)*y- & 2.153585661215e-09_dp)*y+ & 3.654087802817e-08_dp)*y+ & 5.15929575830120e-07_dp)*y- & 9.52388379435709e-06_dp)*y- & 2.16552440036426e-04_dp)*y+ & 9.03551469568320e-03_dp)*y- & 1.45505469175613e-01_dp)*y+ & 1.21449092319186e+00_dp r4 = (((((((((-1.08510370291979e-12_dp*y+ & 6.41492397277798e-11_dp)*y+ & 7.542387436125e-10_dp)*y- & 2.213111836647e-09_dp)*y- & 1.448228963549e-07_dp)*y- & 1.95670833237101e-06_dp)*y- & 1.07481314670844e-05_dp)*y+ & 1.49335941252765e-04_dp)*y+ & 4.87791531990593e-02_dp)*y- & 1.10559909038653e+00_dp)*y+ & 8.09502028611780e+00_dp w1 = ((((((((((-4.65801912689961e-14_dp*y+ & 7.58669507106800e-13_dp)*y- & 1.186387548048e-11_dp)*y+ & 1.862334710665e-10_dp)*y- & 2.799399389539e-09_dp)*y+ & 4.148972684255e-08_dp)*y- & 5.933568079600e-07_dp)*y+ & 8.168349266115e-06_dp)*y- & 1.08989176177409e-04_dp)*y+ & 1.41357961729531e-03_dp)*y- & 1.87588361833659e-02_dp)*y+ & 2.89898651436026e-01_dp w2 = ((((((((((((-1.46345073267549e-14_dp*y+ & 2.25644205432182e-13_dp)*y- & 3.116258693847e-12_dp)*y+ & 4.321908756610e-11_dp)*y- & 5.673270062669e-10_dp)*y+ & 7.006295962960e-09_dp)*y- & 8.120186517000e-08_dp)*y+ & 8.775294645770e-07_dp)*y- & 8.77829235749024e-06_dp)*y+ & 8.04372147732379e-05_dp)*y- & 6.64149238804153e-04_dp)*y+ & 4.81181506827225e-03_dp)*y- & 2.88982669486183e-02_dp)*y+ & 1.56247249979288e-01_dp w3 = (((((((((((((9.06812118895365e-15_dp*y- & 1.40541322766087e-13_dp)*y+ & 1.919270015269e-12_dp)*y- & 2.605135739010e-11_dp)*y+ & 3.299685839012e-10_dp)*y- & 3.86354139348735e-09_dp)*y+ & 4.16265847927498e-08_dp)*y- & 4.09462835471470e-07_dp)*y+ & 3.64018881086111e-06_dp)*y- & 2.88665153269386e-05_dp)*y+ & 2.00515819789028e-04_dp)*y- & 1.18791896897934e-03_dp)*y+ & 5.75223633388589e-03_dp)*y- & 2.09400418772687e-02_dp)*y+ & 4.85368861938873e-02_dp w4 = ((((((((((((((-9.74835552342257e-16_dp*y+ & 1.57857099317175e-14_dp)*y- & 2.249993780112e-13_dp)*y+ & 3.173422008953e-12_dp)*y- & 4.161159459680e-11_dp)*y+ & 5.021343560166e-10_dp)*y- & 5.545047534808e-09_dp)*y+ & 5.554146993491e-08_dp)*y- & 4.99048696190133e-07_dp)*y+ & 3.96650392371311e-06_dp)*y- & 2.73816413291214e-05_dp)*y+ & 1.60106988333186e-04_dp)*y- & 7.64560567879592e-04_dp)*y+ & 2.81330044426892e-03_dp)*y- & 7.16227030134947e-03_dp)*y+ & 9.66077262223353e-03_dp elseif (x <= 10.0e+00_dp) then y = x-7.5e+00_dp r1 = (((((((((4.64217329776215e-15_dp*y- & 6.27892383644164e-15_dp)*y+ & 3.462236347446e-13_dp)*y- & 2.927229355350e-11_dp)*y+ & 5.090355371676e-10_dp)*y- & 9.97272656345253e-09_dp)*y+ & 2.37835295639281e-07_dp)*y- & 4.60301761310921e-06_dp)*y+ & 8.42824204233222e-05_dp)*y- & 1.37983082233081e-03_dp)*y+ & 1.66630865869375e-02_dp r2 = (((((((((2.93981127919047e-14_dp*y+ & 8.47635639065744e-13_dp)*y- & 1.446314544774e-11_dp)*y- & 6.149155555753e-12_dp)*y+ & 8.484275604612e-10_dp)*y- & 6.10898827887652e-08_dp)*y+ & 2.39156093611106e-06_dp)*y- & 5.35837089462592e-05_dp)*y+ & 1.00967602595557e-03_dp)*y- & 1.57769317127372e-02_dp)*y+ & 1.74853819464285e-01_dp r3 = ((((((((((2.93523563363000e-14_dp*y- & 6.40041776667020e-14_dp)*y- & 2.695740446312e-12_dp)*y+ & 1.027082960169e-10_dp)*y- & 5.822038656780e-10_dp)*y- & 3.159991002539e-08_dp)*y+ & 4.327249251331e-07_dp)*y+ & 4.856768455119e-06_dp)*y- & 2.54617989427762e-04_dp)*y+ & 5.54843378106589e-03_dp)*y- & 7.95013029486684e-02_dp)*y+ & 7.20206142703162e-01_dp r4 = (((((((((((-1.62212382394553e-14_dp*y+ & 7.68943641360593e-13_dp)*y+ & 5.764015756615e-12_dp)*y- & 1.380635298784e-10_dp)*y- & 1.476849808675e-09_dp)*y+ & 1.84347052385605e-08_dp)*y+ & 3.34382940759405e-07_dp)*y- & 1.39428366421645e-06_dp)*y- & 7.50249313713996e-05_dp)*y- & 6.26495899187507e-04_dp)*y+ & 4.69716410901162e-02_dp)*y- & 6.66871297428209e-01_dp)*y+ & 4.11207530217806e+00_dp w1 = ((((((((((-1.65995045235997e-15_dp*y+ & 6.91838935879598e-14_dp)*y- & 9.131223418888e-13_dp)*y+ & 1.403341829454e-11_dp)*y- & 3.672235069444e-10_dp)*y+ & 6.366962546990e-09_dp)*y- & 1.039220021671e-07_dp)*y+ & 1.959098751715e-06_dp)*y- & 3.33474893152939e-05_dp)*y+ & 5.72164211151013e-04_dp)*y- & 1.05583210553392e-02_dp)*y+ & 2.26696066029591e-01_dp w2 = ((((((((((((-3.57248951192047e-16_dp*y+ & 6.25708409149331e-15_dp)*y- & 9.657033089714e-14_dp)*y+ & 1.507864898748e-12_dp)*y- & 2.332522256110e-11_dp)*y+ & 3.428545616603e-10_dp)*y- & 4.698730937661e-09_dp)*y+ & 6.219977635130e-08_dp)*y- & 7.83008889613661e-07_dp)*y+ & 9.08621687041567e-06_dp)*y- & 9.86368311253873e-05_dp)*y+ & 9.69632496710088e-04_dp)*y- & 8.14594214284187e-03_dp)*y+ & 8.50218447733457e-02_dp w3 = (((((((((((((1.64742458534277e-16_dp*y- & 2.68512265928410e-15_dp)*y+ & 3.788890667676e-14_dp)*y- & 5.508918529823e-13_dp)*y+ & 7.555896810069e-12_dp)*y- & 9.69039768312637e-11_dp)*y+ & 1.16034263529672e-09_dp)*y- & 1.28771698573873e-08_dp)*y+ & 1.31949431805798e-07_dp)*y- & 1.23673915616005e-06_dp)*y+ & 1.04189803544936e-05_dp)*y- & 7.79566003744742e-05_dp)*y+ & 5.03162624754434e-04_dp)*y- & 2.55138844587555e-03_dp)*y+ & 1.13250730954014e-02_dp w4 = ((((((((((((((-1.55714130075679e-17_dp*y+ & 2.57193722698891e-16_dp)*y- & 3.626606654097e-15_dp)*y+ & 5.234734676175e-14_dp)*y- & 7.067105402134e-13_dp)*y+ & 8.793512664890e-12_dp)*y- & 1.006088923498e-10_dp)*y+ & 1.050565098393e-09_dp)*y- & 9.91517881772662e-09_dp)*y+ & 8.35835975882941e-08_dp)*y- & 6.19785782240693e-07_dp)*y+ & 3.95841149373135e-06_dp)*y- & 2.11366761402403e-05_dp)*y+ & 9.00474771229507e-05_dp)*y- & 2.78777909813289e-04_dp)*y+ & 5.26543779837487e-04_dp elseif (x <= 15.0e+00_dp) then y = x-12.5e+00_dp r1 = (((((((((((4.94869622744119e-17_dp*y+ & 8.03568805739160e-16_dp)*y- & 5.599125915431e-15_dp)*y- & 1.378685560217e-13_dp)*y+ & 7.006511663249e-13_dp)*y+ & 1.30391406991118e-11_dp)*y+ & 8.06987313467541e-11_dp)*y- & 5.20644072732933e-09_dp)*y+ & 7.72794187755457e-08_dp)*y- & 1.61512612564194e-06_dp)*y+ & 4.15083811185831e-05_dp)*y- & 7.87855975560199e-04_dp)*y+ & 1.14189319050009e-02_dp r2 = (((((((((((4.89224285522336e-16_dp*y+ & 1.06390248099712e-14_dp)*y- & 5.446260182933e-14_dp)*y- & 1.613630106295e-12_dp)*y+ & 3.910179118937e-12_dp)*y+ & 1.90712434258806e-10_dp)*y+ & 8.78470199094761e-10_dp)*y- & 5.97332993206797e-08_dp)*y+ & 9.25750831481589e-07_dp)*y- & 2.02362185197088e-05_dp)*y+ & 4.92341968336776e-04_dp)*y- & 8.68438439874703e-03_dp)*y+ & 1.15825965127958e-01_dp r3 = ((((((((((6.12419396208408e-14_dp*y+ & 1.12328861406073e-13_dp)*y- & 9.051094103059e-12_dp)*y- & 4.781797525341e-11_dp)*y+ & 1.660828868694e-09_dp)*y+ & 4.499058798868e-10_dp)*y- & 2.519549641933e-07_dp)*y+ & 4.977444040180e-06_dp)*y- & 1.25858350034589e-04_dp)*y+ & 2.70279176970044e-03_dp)*y- & 3.99327850801083e-02_dp)*y+ & 4.33467200855434e-01_dp r4 = (((((((((((4.63414725924048e-14_dp*y- & 4.72757262693062e-14_dp)*y- & 1.001926833832e-11_dp)*y+ & 6.074107718414e-11_dp)*y+ & 1.576976911942e-09_dp)*y- & 2.01186401974027e-08_dp)*y- & 1.84530195217118e-07_dp)*y+ & 5.02333087806827e-06_dp)*y+ & 9.66961790843006e-06_dp)*y- & 1.58522208889528e-03_dp)*y+ & 2.80539673938339e-02_dp)*y- & 2.78953904330072e-01_dp)*y+ & 1.82835655238235e+00_dp w4 = (((((((((((((2.90401781000996e-18_dp*y- & 4.63389683098251e-17_dp)*y+ & 6.274018198326e-16_dp)*y- & 8.936002188168e-15_dp)*y+ & 1.194719074934e-13_dp)*y- & 1.45501321259466e-12_dp)*y+ & 1.64090830181013e-11_dp)*y- & 1.71987745310181e-10_dp)*y+ & 1.63738403295718e-09_dp)*y- & 1.39237504892842e-08_dp)*y+ & 1.06527318142151e-07_dp)*y- & 7.27634957230524e-07_dp)*y+ & 4.12159381310339e-06_dp)*y- & 1.74648169719173e-05_dp)*y+ & 8.50290130067818e-05_dp w3 = ((((((((((((-4.19569145459480e-17_dp*y+ & 5.94344180261644e-16_dp)*y- & 1.148797566469e-14_dp)*y+ & 1.881303962576e-13_dp)*y- & 2.413554618391e-12_dp)*y+ & 3.372127423047e-11_dp)*y- & 4.933988617784e-10_dp)*y+ & 6.116545396281e-09_dp)*y- & 6.69965691739299e-08_dp)*y+ & 7.52380085447161e-07_dp)*y- & 8.08708393262321e-06_dp)*y+ & 6.88603417296672e-05_dp)*y- & 4.67067112993427e-04_dp)*y+ & 5.42313365864597e-03_dp w2 = ((((((((((-6.22272689880615e-15_dp*y+ & 1.04126809657554e-13_dp)*y- & 6.842418230913e-13_dp)*y+ & 1.576841731919e-11_dp)*y- & 4.203948834175e-10_dp)*y+ & 6.287255934781e-09_dp)*y- & 8.307159819228e-08_dp)*y+ & 1.356478091922e-06_dp)*y- & 2.08065576105639e-05_dp)*y+ & 2.52396730332340e-04_dp)*y- & 2.94484050194539e-03_dp)*y+ & 6.01396183129168e-02_dp recx = 1.0e+00_dp/x w1 = (((-1.8784686463512e-01_dp*recx+ & 2.2991849164985e-01_dp)*recx- & 4.9893752514047e-01_dp)*recx- & 2.1916512131607e-05_dp)*exp(-x)+ & sqrt(PIo4*recx)-w4-w3-w2 elseif (x <= 20.0e+00_dp) then w1 = sqrt(PIo4/x) y = x-17.5e+00_dp r1 = (((((((((((4.36701759531398e-17_dp*y- & 1.12860600219889e-16_dp)*y- & 6.149849164164e-15_dp)*y+ & 5.820231579541e-14_dp)*y+ & 4.396602872143e-13_dp)*y- & 1.24330365320172e-11_dp)*y+ & 6.71083474044549e-11_dp)*y+ & 2.43865205376067e-10_dp)*y+ & 1.67559587099969e-08_dp)*y- & 9.32738632357572e-07_dp)*y+ & 2.39030487004977e-05_dp)*y- & 4.68648206591515e-04_dp)*y+ & 8.34977776583956e-03_dp r2 = (((((((((((4.98913142288158e-16_dp*y- & 2.60732537093612e-16_dp)*y- & 7.775156445127e-14_dp)*y+ & 5.766105220086e-13_dp)*y+ & 6.432696729600e-12_dp)*y- & 1.39571683725792e-10_dp)*y+ & 5.95451479522191e-10_dp)*y+ & 2.42471442836205e-09_dp)*y+ & 2.47485710143120e-07_dp)*y- & 1.14710398652091e-05_dp)*y+ & 2.71252453754519e-04_dp)*y- & 4.96812745851408e-03_dp)*y+ & 8.26020602026780e-02_dp r3 = (((((((((((1.91498302509009e-15_dp*y+ & 1.48840394311115e-14_dp)*y- & 4.316925145767e-13_dp)*y+ & 1.186495793471e-12_dp)*y+ & 4.615806713055e-11_dp)*y- & 5.54336148667141e-10_dp)*y+ & 3.48789978951367e-10_dp)*y- & 2.79188977451042e-09_dp)*y+ & 2.09563208958551e-06_dp)*y- & 6.76512715080324e-05_dp)*y+ & 1.32129867629062e-03_dp)*y- & 2.05062147771513e-02_dp)*y+ & 2.88068671894324e-01_dp r4 = (((((((((((-5.43697691672942e-15_dp*y- & 1.12483395714468e-13_dp)*y+ & 2.826607936174e-12_dp)*y- & 1.266734493280e-11_dp)*y- & 4.258722866437e-10_dp)*y+ & 9.45486578503261e-09_dp)*y- & 5.86635622821309e-08_dp)*y- & 1.28835028104639e-06_dp)*y+ & 4.41413815691885e-05_dp)*y- & 7.61738385590776e-04_dp)*y+ & 9.66090902985550e-03_dp)*y- & 1.01410568057649e-01_dp)*y+ & 9.54714798156712e-01_dp w4 = ((((((((((((-7.56882223582704e-19_dp*y+ & 7.53541779268175e-18_dp)*y- & 1.157318032236e-16_dp)*y+ & 2.411195002314e-15_dp)*y- & 3.601794386996e-14_dp)*y+ & 4.082150659615e-13_dp)*y- & 4.289542980767e-12_dp)*y+ & 5.086829642731e-11_dp)*y- & 6.35435561050807e-10_dp)*y+ & 6.82309323251123e-09_dp)*y- & 5.63374555753167e-08_dp)*y+ & 3.57005361100431e-07_dp)*y- & 2.40050045173721e-06_dp)*y+ & 4.94171300536397e-05_dp w3 = (((((((((((-5.54451040921657e-17_dp*y+ & 2.68748367250999e-16_dp)*y+ & 1.349020069254e-14_dp)*y- & 2.507452792892e-13_dp)*y+ & 1.944339743818e-12_dp)*y- & 1.29816917658823e-11_dp)*y+ & 3.49977768819641e-10_dp)*y- & 8.67270669346398e-09_dp)*y+ & 1.31381116840118e-07_dp)*y- & 1.36790720600822e-06_dp)*y+ & 1.19210697673160e-05_dp)*y- & 1.42181943986587e-04_dp)*y+ & 4.12615396191829e-03_dp w2 = (((((((((((-1.86506057729700e-16_dp*y+ & 1.16661114435809e-15_dp)*y+ & 2.563712856363e-14_dp)*y- & 4.498350984631e-13_dp)*y+ & 1.765194089338e-12_dp)*y+ & 9.04483676345625e-12_dp)*y+ & 4.98930345609785e-10_dp)*y- & 2.11964170928181e-08_dp)*y+ & 3.98295476005614e-07_dp)*y- & 5.49390160829409e-06_dp)*y+ & 7.74065155353262e-05_dp)*y- & 1.48201933009105e-03_dp)*y+ & 4.97836392625268e-02_dp w1 = ((1.9623264149430e-01_dp/x- & 4.9695241464490e-01_dp)/x- & 6.0156581186481e-05_dp)*exp(-x)+ & w1-w2-w3-w4 elseif (x <= 25.0e+00_dp) then w1 = sqrt(PIo4/x) e = exp(-x) recx = 1.0e+00_dp/x r1 = ((((((-4.45711399441838e-05_dp*x+ & 1.27267770241379e-03_dp)*x- & 2.36954961381262e-01_dp)*x+ & 1.54330657903756e+01_dp)*x- & 5.22799159267808e+02_dp)*x+ & 1.05951216669313e+04_dp)*x+ & (-2.51177235556236e+06_dp*recx+ & 8.72975373557709e+05_dp)*recx- & 1.29194382386499e+05_dp)*e+r14/(x-r14) r2 = (((((-7.85617372254488e-02_dp*x+ & 6.35653573484868e+00_dp)*x- & 3.38296938763990e+02_dp)*x+ & 1.25120495802096e+04_dp)*x- & 3.16847570511637e+05_dp)*x+ & ((-1.02427466127427e+09_dp*recx+ & 3.70104713293016e+08_dp)*recx- & 5.87119005093822e+07_dp)*recx+ & 5.38614211391604e+06_dp)*e+r24/(x-r24) r3 = (((((-2.37900485051067e-01_dp*x+ & 1.84122184400896e+01_dp)*x- & 1.00200731304146e+03_dp)*x+ & 3.75151841595736e+04_dp)*x- & 9.50626663390130e+05_dp)*x+ & ((-2.88139014651985e+09_dp*recx+ & 1.06625915044526e+09_dp)*recx- & 1.72465289687396e+08_dp)*recx+ & 1.60419390230055e+07_dp)*e+r34/(x-r34) r4 = ((((((-6.00691586407385e-04_dp*x- & 3.64479545338439e-01_dp)*x+ & 1.57496131755179e+01_dp)*x- & 6.54944248734901e+02_dp)*x+ & 1.70830039597097e+04_dp)*x- & 2.90517939780207e+05_dp)*x+ & (+3.49059698304732e+07_dp*recx- & 1.64944522586065e+07_dp)*recx+ & 2.96817940164703e+06_dp)*e+r44/(x-r44) w4 = (((((((2.33766206773151e-07_dp*x- & 3.81542906607063e-05_dp)*x+ & 3.51416601267000e-03_dp)*x- & 1.66538571864728e-01_dp)*x+ & 4.80006136831847e+00_dp)*x- & 8.73165934223603e+01_dp)*x+ & 9.77683627474638e+02_dp)*x+ & 1.66000945117640e+04_dp*recx- & 6.14479071209961e+03_dp)*e+w44*w1 w3 = ((((((2.36392855180768e-04_dp*x- & 9.16785337967013e-03_dp)*x+ & 4.62186525041313e-01_dp)*x- & 1.96943786006540e+01_dp)*x+ & 4.99169195295559e+02_dp)*x- & 6.21419845845090e+03_dp)*x+ & ((+5.21445053212414e+07_dp*recx- & 1.34113464389309e+07_dp)*recx+ & 1.13673298305631e+06_dp)*recx- & 2.81501182042707e+03_dp)*e+w34*w1 w2 = ((((((7.29841848989391e-04_dp*x- & 3.53899555749875e-02_dp)*x+ & 2.07797425718513e+00_dp)*x- & 1.00464709786287e+02_dp)*x+ & 3.15206108877819e+03_dp)*x- & 6.27054715090012e+04_dp)*x+ & (+1.54721246264919e+07_dp*recx- & 5.26074391316381e+06_dp)*recx+ & 7.67135400969617e+05_dp)*e+w24*w1 w1 = ((1.9623264149430e-01_dp*recx- & 4.9695241464490e-01_dp)*recx- & 6.0156581186481e-05_dp)*e+ & w1-w2-w3-w4 elseif (x <= 35.0e+00_dp) then w1 = sqrt(PIo4/x) e = exp(-x) recx = 1.0e+00_dp/x r1 = ((((((-4.45711399441838e-05_dp*x+ & 1.27267770241379e-03_dp)*x- & 2.36954961381262e-01_dp)*x+ & 1.54330657903756e+01_dp)*x- & 5.22799159267808e+02_dp)*x+ & 1.05951216669313e+04_dp)*x+ & (-2.51177235556236e+06_dp*recx+ & 8.72975373557709e+05_dp)*recx- & 1.29194382386499e+05_dp)*e+r14/(x-r14) r2 = (((((-7.85617372254488e-02_dp*x+ & 6.35653573484868e+00_dp)*x- & 3.38296938763990e+02_dp)*x+ & 1.25120495802096e+04_dp)*x- & 3.16847570511637e+05_dp)*x+ & ((-1.02427466127427e+09_dp*recx+ & 3.70104713293016e+08_dp)*recx- & 5.87119005093822e+07_dp)*recx+ & 5.38614211391604e+06_dp)*e+r24/(x-r24) r3 = (((((-2.37900485051067e-01_dp*x+ & 1.84122184400896e+01_dp)*x- & 1.00200731304146e+03_dp)*x+ & 3.75151841595736e+04_dp)*x- & 9.50626663390130e+05_dp)*x+ & ((-2.88139014651985e+09_dp*recx+ & 1.06625915044526e+09_dp)*recx- & 1.72465289687396e+08_dp)*recx+ & 1.60419390230055e+07_dp)*e+r34/(x-r34) r4 = ((((((-6.00691586407385e-04_dp*x- & 3.64479545338439e-01_dp)*x+ & 1.57496131755179e+01_dp)*x- & 6.54944248734901e+02_dp)*x+ & 1.70830039597097e+04_dp)*x- & 2.90517939780207e+05_dp)*x+ & (+3.49059698304732e+07_dp*recx- & 1.64944522586065e+07_dp)*recx+ & 2.96817940164703e+06_dp)*e+r44/(x-r44) w4 = ((((((5.74245945342286e-06_dp*x- & 7.58735928102351e-05_dp)*x+ & 2.35072857922892e-04_dp)*x- & 3.78812134013125e-03_dp)*x+ & 3.09871652785805e-01_dp)*x- & 7.11108633061306e+00_dp)*x+ & 5.55297573149528e+01_dp)*e+w44*w1 w3 = ((((((2.36392855180768e-04_dp*x- & 9.16785337967013e-03_dp)*x+ & 4.62186525041313e-01_dp)*x- & 1.96943786006540e+01_dp)*x+ & 4.99169195295559e+02_dp)*x- & 6.21419845845090e+03_dp)*x+ & ((+5.21445053212414e+07_dp*recx- & 1.34113464389309e+07_dp)*recx+ & 1.13673298305631e+06_dp)*recx- & 2.81501182042707e+03_dp)*e+w34*w1 w2 = ((((((7.29841848989391e-04_dp*x- & 3.53899555749875e-02_dp)*x+ & 2.07797425718513e+00_dp)*x- & 1.00464709786287e+02_dp)*x+ & 3.15206108877819e+03_dp)*x- & 6.27054715090012e+04_dp)*x+ & (+1.54721246264919e+07_dp*recx- & 5.26074391316381e+06_dp)*recx+ & 7.67135400969617e+05_dp)*e+w24*w1 w1 = ((1.9623264149430e-01_dp*recx- & 4.9695241464490e-01_dp)*recx- & 6.0156581186481e-05_dp)*e+ & w1-w2-w3-w4 elseif (x <= 53.0e+00_dp) then w1 = sqrt(PIo4/x) e = exp(-x)*(x*x)**2 r4 = ((-2.19135070169653e-03_dp*x- & 1.19108256987623e-01_dp)*x- & 7.50238795695573e-01_dp)*e+r44/(x-r44) r3 = ((-9.65842534508637e-04_dp*x- & 4.49822013469279e-02_dp)*x+ & 6.08784033347757e-01_dp)*e+r34/(x-r34) r2 = ((-3.62569791162153e-04_dp*x- & 9.09231717268466e-03_dp)*x+ & 1.84336760556262e-01_dp)*e+r24/(x-r24) r1 = ((-4.07557525914600e-05_dp*x- & 6.88846864931685e-04_dp)*x+ & 1.74725309199384e-02_dp)*e+r14/(x-r14) w4 = ((5.76631982000990e-06_dp*x- & 7.89187283804890e-05_dp)*x+ & 3.28297971853126e-04_dp)*e+w44*w1 w3 = ((2.08294969857230e-04_dp*x- & 3.77489954837361e-03_dp)*x+ & 2.09857151617436e-02_dp)*e+w34*w1 w2 = ((6.16374517326469e-04_dp*x- & 1.26711744680092e-02_dp)*x+ & 8.14504890732155e-02_dp)*e+w24*w1 w1 = w1-w2-w3-w4 else r1 = r14/(x-r14) r2 = r24/(x-r24) r3 = r34/(x-r34) r4 = r44/(x-r44) w1 = sqrt(PIo4/x) w2 = w24*w1 w3 = w34*w1 w4 = w44*w1 w1 = w1-w2-w3-w4 end if r = (/r1, r2, r3, r4/) w = (/w1, w2, w3, w4/) end subroutine rys_rt4 !----------------------------------------------------------------- subroutine rys_rt5(x, r, w) real(KIND=dp), intent(IN) :: & x real(KIND=dp), intent(OUT) :: & r(5), w(5) real(KIND=dp), parameter :: & r15 = 1.17581320211778e-01_dp, & r25 = 1.07456201243690e+00_dp, & w25 = 2.70967405960535e-01_dp, & r35 = 3.08593744371754e+00_dp, & w35 = 3.82231610015404e-02_dp, & r45 = 6.41472973366203e+00_dp, & w45 = 1.51614186862443e-03_dp, & r55 = 1.18071894899717e+01_dp, & w55 = 8.62130526143657e-06_dp real(KIND=dp) :: y, e, r1, r2, r3, r4, r5, w1, w2, w3, w4, w5 if (x <= 3.0e-07_dp) then r1 = 2.26659266316985e-02_dp-2.15865967920897e-03_dp*x r2 = 2.31271692140903e-01_dp-2.20258754389745e-02_dp*x r3 = 8.57346024118836e-01_dp-8.16520023025515e-02_dp*x r4 = 2.97353038120346e+00_dp-2.83193369647137e-01_dp*x r5 = 1.84151859759051e+01_dp-1.75382723579439e+00_dp*x w1 = 2.95524224714752e-01_dp-1.96867576909777e-02_dp*x w2 = 2.69266719309995e-01_dp-5.61737590184721e-02_dp*x w3 = 2.19086362515981e-01_dp-9.71152726793658e-02_dp*x w4 = 1.49451349150580e-01_dp-1.02979262193565e-01_dp*x w5 = 6.66713443086877e-02_dp-5.73782817488315e-02_dp*x elseif (x <= 1.0e+00_dp) then r1 = ((((((-4.46679165328413e-11_dp*x+ & 1.21879111988031e-09_dp)*x- & 2.62975022612104e-08_dp)*x+ & 5.15106194905897e-07_dp)*x- & 9.27933625824749e-06_dp)*x+ & 1.51794097682482e-04_dp)*x- & 2.15865967920301e-03_dp)*x+ & 2.26659266316985e-02_dp r2 = ((((((1.93117331714174e-10_dp*x- & 4.57267589660699e-09_dp)*x+ & 2.48339908218932e-08_dp)*x+ & 1.50716729438474e-06_dp)*x- & 6.07268757707381e-05_dp)*x+ & 1.37506939145643e-03_dp)*x- & 2.20258754419939e-02_dp)*x+ & 2.31271692140905e-01_dp r3 = (((((4.84989776180094e-09_dp*x+ & 1.31538893944284e-07_dp)*x- & 2.766753852879e-06_dp)*x- & 7.651163510626e-05_dp)*x+ & 4.033058545972e-03_dp)*x- & 8.16520022916145e-02_dp)*x+ & 8.57346024118779e-01_dp r4 = ((((-2.48581772214623e-07_dp*x- & 4.34482635782585e-06_dp)*x- & 7.46018257987630e-07_dp)*x+ & 1.01210776517279e-02_dp)*x- & 2.83193369640005e-01_dp)*x+ & 2.97353038120345e+00_dp r5 = (((((-8.92432153868554e-09_dp*x+ & 1.77288899268988e-08_dp)*x+ & 3.040754680666e-06_dp)*x+ & 1.058229325071e-04_dp)*x+ & 4.596379534985e-02_dp)*x- & 1.75382723579114e+00_dp)*x+ & 1.84151859759049e+01_dp w1 = ((((((-2.03822632771791e-09_dp*x+ & 3.89110229133810e-08_dp)*x- & 5.84914787904823e-07_dp)*x+ & 8.30316168666696e-06_dp)*x- & 1.13218402310546e-04_dp)*x+ & 1.49128888586790e-03_dp)*x- & 1.96867576904816e-02_dp)*x+ & 2.95524224714749e-01_dp w2 = (((((((8.62848118397570e-09_dp*x- & 1.38975551148989e-07_dp)*x+ & 1.602894068228e-06_dp)*x- & 1.646364300836e-05_dp)*x+ & 1.538445806778e-04_dp)*x- & 1.28848868034502e-03_dp)*x+ & 9.38866933338584e-03_dp)*x- & 5.61737590178812e-02_dp)*x+ & 2.69266719309991e-01_dp w3 = ((((((((-9.41953204205665e-09_dp*x+ & 1.47452251067755e-07_dp)*x- & 1.57456991199322e-06_dp)*x+ & 1.45098401798393e-05_dp)*x- & 1.18858834181513e-04_dp)*x+ & 8.53697675984210e-04_dp)*x- & 5.22877807397165e-03_dp)*x+ & 2.60854524809786e-02_dp)*x- & 9.71152726809059e-02_dp)*x+ & 2.19086362515979e-01_dp w4 = ((((((((-3.84961617022042e-08_dp*x+ & 5.66595396544470e-07_dp)*x- & 5.52351805403748e-06_dp)*x+ & 4.53160377546073e-05_dp)*x- & 3.22542784865557e-04_dp)*x+ & 1.95682017370967e-03_dp)*x- & 9.77232537679229e-03_dp)*x+ & 3.79455945268632e-02_dp)*x- & 1.02979262192227e-01_dp)*x+ & 1.49451349150573e-01_dp w5 = (((((((((4.09594812521430e-09_dp*x- & 6.47097874264417e-08_dp)*x+ & 6.743541482689e-07_dp)*x- & 5.917993920224e-06_dp)*x+ & 4.531969237381e-05_dp)*x- & 2.99102856679638e-04_dp)*x+ & 1.65695765202643e-03_dp)*x- & 7.40671222520653e-03_dp)*x+ & 2.50889946832192e-02_dp)*x- & 5.73782817487958e-02_dp)*x+ & 6.66713443086877e-02_dp elseif (x <= 5.0e+00_dp) then y = x-3.0e+00_dp r1 = ((((((((-2.58163897135138e-14_dp*y+ & 8.14127461488273e-13_dp)*y- & 2.11414838976129e-11_dp)*y+ & 5.09822003260014e-10_dp)*y- & 1.16002134438663e-08_dp)*y+ & 2.46810694414540e-07_dp)*y- & 4.92556826124502e-06_dp)*y+ & 9.02580687971053e-05_dp)*y- & 1.45190025120726e-03_dp)*y+ & 1.73416786387475e-02_dp r2 = (((((((((1.04525287289788e-14_dp*y+ & 5.44611782010773e-14_dp)*y- & 4.831059411392e-12_dp)*y+ & 1.136643908832e-10_dp)*y- & 1.104373076913e-09_dp)*y- & 2.35346740649916e-08_dp)*y+ & 1.43772622028764e-06_dp)*y- & 4.23405023015273e-05_dp)*y+ & 9.12034574793379e-04_dp)*y- & 1.52479441718739e-02_dp)*y+ & 1.76055265928744e-01_dp r3 = (((((((((-6.89693150857911e-14_dp*y+ & 5.92064260918861e-13_dp)*y+ & 1.847170956043e-11_dp)*y- & 3.390752744265e-10_dp)*y- & 2.995532064116e-09_dp)*y+ & 1.57456141058535e-07_dp)*y- & 3.95859409711346e-07_dp)*y- & 9.58924580919747e-05_dp)*y+ & 3.23551502557785e-03_dp)*y- & 5.97587007636479e-02_dp)*y+ & 6.46432853383057e-01_dp r4 = ((((((((-3.61293809667763e-12_dp*y- & 2.70803518291085e-11_dp)*y+ & 8.83758848468769e-10_dp)*y+ & 1.59166632851267e-08_dp)*y- & 1.32581997983422e-07_dp)*y- & 7.60223407443995e-06_dp)*y- & 7.41019244900952e-05_dp)*y+ & 9.81432631743423e-03_dp)*y- & 2.23055570487771e-01_dp)*y+ & 2.21460798080643e+00_dp r5 = (((((((((7.12332088345321e-13_dp*y+ & 3.16578501501894e-12_dp)*y- & 8.776668218053e-11_dp)*y- & 2.342817613343e-09_dp)*y- & 3.496962018025e-08_dp)*y- & 3.03172870136802e-07_dp)*y+ & 1.50511293969805e-06_dp)*y+ & 1.37704919387696e-04_dp)*y+ & 4.70723869619745e-02_dp)*y- & 1.47486623003693e+00_dp)*y+ & 1.35704792175847e+01_dp w1 = (((((((((1.04348658616398e-13_dp*y- & 1.94147461891055e-12_dp)*y+ & 3.485512360993e-11_dp)*y- & 6.277497362235e-10_dp)*y+ & 1.100758247388e-08_dp)*y- & 1.88329804969573e-07_dp)*y+ & 3.12338120839468e-06_dp)*y- & 5.04404167403568e-05_dp)*y+ & 8.00338056610995e-04_dp)*y- & 1.30892406559521e-02_dp)*y+ & 2.47383140241103e-01_dp w2 = (((((((((((3.23496149760478e-14_dp*y- & 5.24314473469311e-13_dp)*y+ & 7.743219385056e-12_dp)*y- & 1.146022750992e-10_dp)*y+ & 1.615238462197e-09_dp)*y- & 2.15479017572233e-08_dp)*y+ & 2.70933462557631e-07_dp)*y- & 3.18750295288531e-06_dp)*y+ & 3.47425221210099e-05_dp)*y- & 3.45558237388223e-04_dp)*y+ & 3.05779768191621e-03_dp)*y- & 2.29118251223003e-02_dp)*y+ & 1.59834227924213e-01_dp w3 = ((((((((((((-3.42790561802876e-14_dp*y+ & 5.26475736681542e-13_dp)*y- & 7.184330797139e-12_dp)*y+ & 9.763932908544e-11_dp)*y- & 1.244014559219e-09_dp)*y+ & 1.472744068942e-08_dp)*y- & 1.611749975234e-07_dp)*y+ & 1.616487851917e-06_dp)*y- & 1.46852359124154e-05_dp)*y+ & 1.18900349101069e-04_dp)*y- & 8.37562373221756e-04_dp)*y+ & 4.93752683045845e-03_dp)*y- & 2.25514728915673e-02_dp)*y+ & 6.95211812453929e-02_dp w4 = (((((((((((((1.04072340345039e-14_dp*y- & 1.60808044529211e-13_dp)*y+ & 2.183534866798e-12_dp)*y- & 2.939403008391e-11_dp)*y+ & 3.679254029085e-10_dp)*y- & 4.23775673047899e-09_dp)*y+ & 4.46559231067006e-08_dp)*y- & 4.26488836563267e-07_dp)*y+ & 3.64721335274973e-06_dp)*y- & 2.74868382777722e-05_dp)*y+ & 1.78586118867488e-04_dp)*y- & 9.68428981886534e-04_dp)*y+ & 4.16002324339929e-03_dp)*y- & 1.28290192663141e-02_dp)*y+ & 2.22353727685016e-02_dp w5 = ((((((((((((((-8.16770412525963e-16_dp*y+ & 1.31376515047977e-14_dp)*y- & 1.856950818865e-13_dp)*y+ & 2.596836515749e-12_dp)*y- & 3.372639523006e-11_dp)*y+ & 4.025371849467e-10_dp)*y- & 4.389453269417e-09_dp)*y+ & 4.332753856271e-08_dp)*y- & 3.82673275931962e-07_dp)*y+ & 2.98006900751543e-06_dp)*y- & 2.00718990300052e-05_dp)*y+ & 1.13876001386361e-04_dp)*y- & 5.23627942443563e-04_dp)*y+ & 1.83524565118203e-03_dp)*y- & 4.37785737450783e-03_dp)*y+ & 5.36963805223095e-03_dp elseif (x <= 10.0e+00_dp) then y = x-7.5e+00_dp r1 = ((((((((-1.13825201010775e-14_dp*y+ & 1.89737681670375e-13_dp)*y- & 4.81561201185876e-12_dp)*y+ & 1.56666512163407e-10_dp)*y- & 3.73782213255083e-09_dp)*y+ & 9.15858355075147e-08_dp)*y- & 2.13775073585629e-06_dp)*y+ & 4.56547356365536e-05_dp)*y- & 8.68003909323740e-04_dp)*y+ & 1.22703754069176e-02_dp r2 = (((((((((-3.67160504428358e-15_dp*y+ & 1.27876280158297e-14_dp)*y- & 1.296476623788e-12_dp)*y+ & 1.477175434354e-11_dp)*y+ & 5.464102147892e-10_dp)*y- & 2.42538340602723e-08_dp)*y+ & 8.20460740637617e-07_dp)*y- & 2.20379304598661e-05_dp)*y+ & 4.90295372978785e-04_dp)*y- & 9.14294111576119e-03_dp)*y+ & 1.22590403403690e-01_dp r3 = (((((((((1.39017367502123e-14_dp*y- & 6.96391385426890e-13_dp)*y+ & 1.176946020731e-12_dp)*y+ & 1.725627235645e-10_dp)*y- & 3.686383856300e-09_dp)*y+ & 2.87495324207095e-08_dp)*y+ & 1.71307311000282e-06_dp)*y- & 7.94273603184629e-05_dp)*y+ & 2.00938064965897e-03_dp)*y- & 3.63329491677178e-02_dp)*y+ & 4.34393683888443e-01_dp r4 = ((((((((((-1.27815158195209e-14_dp*y+ & 1.99910415869821e-14_dp)*y+ & 3.753542914426e-12_dp)*y- & 2.708018219579e-11_dp)*y- & 1.190574776587e-09_dp)*y+ & 1.106696436509e-08_dp)*y+ & 3.954955671326e-07_dp)*y- & 4.398596059588e-06_dp)*y- & 2.01087998907735e-04_dp)*y+ & 7.89092425542937e-03_dp)*y- & 1.42056749162695e-01_dp)*y+ & 1.39964149420683e+00_dp r5 = ((((((((((-1.19442341030461e-13_dp*y- & 2.34074833275956e-12_dp)*y+ & 6.861649627426e-12_dp)*y+ & 6.082671496226e-10_dp)*y+ & 5.381160105420e-09_dp)*y- & 6.253297138700e-08_dp)*y- & 2.135966835050e-06_dp)*y- & 2.373394341886e-05_dp)*y+ & 2.88711171412814e-06_dp)*y+ & 4.85221195290753e-02_dp)*y- & 1.04346091985269e+00_dp)*y+ & 7.89901551676692e+00_dp w1 = (((((((((7.95526040108997e-15_dp*y- & 2.48593096128045e-13_dp)*y+ & 4.761246208720e-12_dp)*y- & 9.535763686605e-11_dp)*y+ & 2.225273630974e-09_dp)*y- & 4.49796778054865e-08_dp)*y+ & 9.17812870287386e-07_dp)*y- & 1.86764236490502e-05_dp)*y+ & 3.76807779068053e-04_dp)*y- & 8.10456360143408e-03_dp)*y+ & 2.01097936411496e-01_dp w2 = (((((((((((1.25678686624734e-15_dp*y- & 2.34266248891173e-14_dp)*y+ & 3.973252415832e-13_dp)*y- & 6.830539401049e-12_dp)*y+ & 1.140771033372e-10_dp)*y- & 1.82546185762009e-09_dp)*y+ & 2.77209637550134e-08_dp)*y- & 4.01726946190383e-07_dp)*y+ & 5.48227244014763e-06_dp)*y- & 6.95676245982121e-05_dp)*y+ & 8.05193921815776e-04_dp)*y- & 8.15528438784469e-03_dp)*y+ & 9.71769901268114e-02_dp w3 = ((((((((((((-8.20929494859896e-16_dp*y+ & 1.37356038393016e-14_dp)*y- & 2.022863065220e-13_dp)*y+ & 3.058055403795e-12_dp)*y- & 4.387890955243e-11_dp)*y+ & 5.923946274445e-10_dp)*y- & 7.503659964159e-09_dp)*y+ & 8.851599803902e-08_dp)*y- & 9.65561998415038e-07_dp)*y+ & 9.60884622778092e-06_dp)*y- & 8.56551787594404e-05_dp)*y+ & 6.66057194311179e-04_dp)*y- & 4.17753183902198e-03_dp)*y+ & 2.25443826852447e-02_dp w4 = ((((((((((((((-1.08764612488790e-17_dp*y+ & 1.85299909689937e-16_dp)*y- & 2.730195628655e-15_dp)*y+ & 4.127368817265e-14_dp)*y- & 5.881379088074e-13_dp)*y+ & 7.805245193391e-12_dp)*y- & 9.632707991704e-11_dp)*y+ & 1.099047050624e-09_dp)*y- & 1.15042731790748e-08_dp)*y+ & 1.09415155268932e-07_dp)*y- & 9.33687124875935e-07_dp)*y+ & 7.02338477986218e-06_dp)*y- & 4.53759748787756e-05_dp)*y+ & 2.41722511389146e-04_dp)*y- & 9.75935943447037e-04_dp)*y+ & 2.57520532789644e-03_dp w5 = (((((((((((((((7.28996979748849e-19_dp*y- & 1.26518146195173e-17_dp)*y+ & 1.886145834486e-16_dp)*y- & 2.876728287383e-15_dp)*y+ & 4.114588668138e-14_dp)*y- & 5.44436631413933e-13_dp)*y+ & 6.64976446790959e-12_dp)*y- & 7.44560069974940e-11_dp)*y+ & 7.57553198166848e-10_dp)*y- & 6.92956101109829e-09_dp)*y+ & 5.62222859033624e-08_dp)*y- & 3.97500114084351e-07_dp)*y+ & 2.39039126138140e-06_dp)*y- & 1.18023950002105e-05_dp)*y+ & 4.52254031046244e-05_dp)*y- & 1.21113782150370e-04_dp)*y+ & 1.75013126731224e-04_dp elseif (x <= 15.0e+00_dp) then y = x-12.5e+00_dp r1 = ((((((((((-4.16387977337393e-17_dp*y+ & 7.20872997373860e-16_dp)*y+ & 1.395993802064e-14_dp)*y+ & 3.660484641252e-14_dp)*y- & 4.154857548139e-12_dp)*y+ & 2.301379846544e-11_dp)*y- & 1.033307012866e-09_dp)*y+ & 3.997777641049e-08_dp)*y- & 9.35118186333939e-07_dp)*y+ & 2.38589932752937e-05_dp)*y- & 5.35185183652937e-04_dp)*y+ & 8.85218988709735e-03_dp r2 = ((((((((((-4.56279214732217e-16_dp*y+ & 6.24941647247927e-15_dp)*y+ & 1.737896339191e-13_dp)*y+ & 8.964205979517e-14_dp)*y- & 3.538906780633e-11_dp)*y+ & 9.561341254948e-11_dp)*y- & 9.772831891310e-09_dp)*y+ & 4.240340194620e-07_dp)*y- & 1.02384302866534e-05_dp)*y+ & 2.57987709704822e-04_dp)*y- & 5.54735977651677e-03_dp)*y+ & 8.68245143991948e-02_dp r3 = ((((((((((-2.52879337929239e-15_dp*y+ & 2.13925810087833e-14_dp)*y+ & 7.884307667104e-13_dp)*y- & 9.023398159510e-13_dp)*y- & 5.814101544957e-11_dp)*y- & 1.333480437968e-09_dp)*y- & 2.217064940373e-08_dp)*y+ & 1.643290788086e-06_dp)*y- & 4.39602147345028e-05_dp)*y+ & 1.08648982748911e-03_dp)*y- & 2.13014521653498e-02_dp)*y+ & 2.94150684465425e-01_dp r4 = ((((((((((-6.42391438038888e-15_dp*y+ & 5.37848223438815e-15_dp)*y+ & 8.960828117859e-13_dp)*y+ & 5.214153461337e-11_dp)*y- & 1.106601744067e-10_dp)*y- & 2.007890743962e-08_dp)*y+ & 1.543764346501e-07_dp)*y+ & 4.520749076914e-06_dp)*y- & 1.88893338587047e-04_dp)*y+ & 4.73264487389288e-03_dp)*y- & 7.91197893350253e-02_dp)*y+ & 8.60057928514554e-01_dp r5 = (((((((((((-2.24366166957225e-14_dp*y+ & 4.87224967526081e-14_dp)*y+ & 5.587369053655e-12_dp)*y- & 3.045253104617e-12_dp)*y- & 1.223983883080e-09_dp)*y- & 2.05603889396319e-09_dp)*y+ & 2.58604071603561e-07_dp)*y+ & 1.34240904266268e-06_dp)*y- & 5.72877569731162e-05_dp)*y- & 9.56275105032191e-04_dp)*y+ & 4.23367010370921e-02_dp)*y- & 5.76800927133412e-01_dp)*y+ & 3.87328263873381e+00_dp w1 = (((((((((8.98007931950169e-15_dp*y+ & 7.25673623859497e-14_dp)*y+ & 5.851494250405e-14_dp)*y- & 4.234204823846e-11_dp)*y+ & 3.911507312679e-10_dp)*y- & 9.65094802088511e-09_dp)*y+ & 3.42197444235714e-07_dp)*y- & 7.51821178144509e-06_dp)*y+ & 1.94218051498662e-04_dp)*y- & 5.38533819142287e-03_dp)*y+ & 1.68122596736809e-01_dp w2 = ((((((((((-1.05490525395105e-15_dp*y+ & 1.96855386549388e-14_dp)*y- & 5.500330153548e-13_dp)*y+ & 1.003849567976e-11_dp)*y- & 1.720997242621e-10_dp)*y+ & 3.533277061402e-09_dp)*y- & 6.389171736029e-08_dp)*y+ & 1.046236652393e-06_dp)*y- & 1.73148206795827e-05_dp)*y+ & 2.57820531617185e-04_dp)*y- & 3.46188265338350e-03_dp)*y+ & 7.03302497508176e-02_dp w3 = (((((((((((3.60020423754545e-16_dp*y- & 6.24245825017148e-15_dp)*y+ & 9.945311467434e-14_dp)*y- & 1.749051512721e-12_dp)*y+ & 2.768503957853e-11_dp)*y- & 4.08688551136506e-10_dp)*y+ & 6.04189063303610e-09_dp)*y- & 8.23540111024147e-08_dp)*y+ & 1.01503783870262e-06_dp)*y- & 1.20490761741576e-05_dp)*y+ & 1.26928442448148e-04_dp)*y- & 1.05539461930597e-03_dp)*y+ & 1.15543698537013e-02_dp w4 = (((((((((((((2.51163533058925e-18_dp*y- & 4.31723745510697e-17_dp)*y+ & 6.557620865832e-16_dp)*y- & 1.016528519495e-14_dp)*y+ & 1.491302084832e-13_dp)*y- & 2.06638666222265e-12_dp)*y+ & 2.67958697789258e-11_dp)*y- & 3.23322654638336e-10_dp)*y+ & 3.63722952167779e-09_dp)*y- & 3.75484943783021e-08_dp)*y+ & 3.49164261987184e-07_dp)*y- & 2.92658670674908e-06_dp)*y+ & 2.12937256719543e-05_dp)*y- & 1.19434130620929e-04_dp)*y+ & 6.45524336158384e-04_dp w5 = ((((((((((((((-1.29043630202811e-19_dp*y+ & 2.16234952241296e-18_dp)*y- & 3.107631557965e-17_dp)*y+ & 4.570804313173e-16_dp)*y- & 6.301348858104e-15_dp)*y+ & 8.031304476153e-14_dp)*y- & 9.446196472547e-13_dp)*y+ & 1.018245804339e-11_dp)*y- & 9.96995451348129e-11_dp)*y+ & 8.77489010276305e-10_dp)*y- & 6.84655877575364e-09_dp)*y+ & 4.64460857084983e-08_dp)*y- & 2.66924538268397e-07_dp)*y+ & 1.24621276265907e-06_dp)*y- & 4.30868944351523e-06_dp)*y+ & 9.94307982432868e-06_dp elseif (x <= 20.0e+00_dp) then y = x-17.5e+00_dp r1 = ((((((((((1.91875764545740e-16_dp*y+ & 7.8357401095707e-16_dp)*y- & 3.260875931644e-14_dp)*y- & 1.186752035569e-13_dp)*y+ & 4.275180095653e-12_dp)*y+ & 3.357056136731e-11_dp)*y- & 1.123776903884e-09_dp)*y+ & 1.231203269887e-08_dp)*y- & 3.99851421361031e-07_dp)*y+ & 1.45418822817771e-05_dp)*y- & 3.49912254976317e-04_dp)*y+ & 6.67768703938812e-03_dp r2 = ((((((((((2.02778478673555e-15_dp*y+ & 1.01640716785099e-14_dp)*y- & 3.385363492036e-13_dp)*y- & 1.615655871159e-12_dp)*y+ & 4.527419140333e-11_dp)*y+ & 3.853670706486e-10_dp)*y- & 1.184607130107e-08_dp)*y+ & 1.347873288827e-07_dp)*y- & 4.47788241748377e-06_dp)*y+ & 1.54942754358273e-04_dp)*y- & 3.55524254280266e-03_dp)*y+ & 6.44912219301603e-02_dp r3 = ((((((((((7.79850771456444e-15_dp*y+ & 6.00464406395001e-14_dp)*y- & 1.249779730869e-12_dp)*y- & 1.020720636353e-11_dp)*y+ & 1.814709816693e-10_dp)*y+ & 1.766397336977e-09_dp)*y- & 4.603559449010e-08_dp)*y+ & 5.863956443581e-07_dp)*y- & 2.03797212506691e-05_dp)*y+ & 6.31405161185185e-04_dp)*y- & 1.30102750145071e-02_dp)*y+ & 2.10244289044705e-01_dp r4 = (((((((((((-2.92397030777912e-15_dp*y+ & 1.94152129078465e-14_dp)*y+ & 4.859447665850e-13_dp)*y- & 3.217227223463e-12_dp)*y- & 7.484522135512e-11_dp)*y+ & 7.19101516047753e-10_dp)*y+ & 6.88409355245582e-09_dp)*y- & 1.44374545515769e-07_dp)*y+ & 2.74941013315834e-06_dp)*y- & 1.02790452049013e-04_dp)*y+ & 2.59924221372643e-03_dp)*y- & 4.35712368303551e-02_dp)*y+ & 5.62170709585029e-01_dp r5 = (((((((((((1.17976126840060e-14_dp*y+ & 1.24156229350669e-13_dp)*y- & 3.892741622280e-12_dp)*y- & 7.755793199043e-12_dp)*y+ & 9.492190032313e-10_dp)*y- & 4.98680128123353e-09_dp)*y- & 1.81502268782664e-07_dp)*y+ & 2.69463269394888e-06_dp)*y+ & 2.50032154421640e-05_dp)*y- & 1.33684303917681e-03_dp)*y+ & 2.29121951862538e-02_dp)*y- & 2.45653725061323e-01_dp)*y+ & 1.89999883453047e+00_dp w1 = ((((((((((1.74841995087592e-15_dp*y- & 6.95671892641256e-16_dp)*y- & 3.000659497257e-13_dp)*y+ & 2.021279817961e-13_dp)*y+ & 3.853596935400e-11_dp)*y+ & 1.461418533652e-10_dp)*y- & 1.014517563435e-08_dp)*y+ & 1.132736008979e-07_dp)*y- & 2.86605475073259e-06_dp)*y+ & 1.21958354908768e-04_dp)*y- & 3.86293751153466e-03_dp)*y+ & 1.45298342081522e-01_dp w2 = ((((((((((-1.11199320525573e-15_dp*y+ & 1.85007587796671e-15_dp)*y+ & 1.220613939709e-13_dp)*y+ & 1.275068098526e-12_dp)*y- & 5.341838883262e-11_dp)*y+ & 6.161037256669e-10_dp)*y- & 1.009147879750e-08_dp)*y+ & 2.907862965346e-07_dp)*y- & 6.12300038720919e-06_dp)*y+ & 1.00104454489518e-04_dp)*y- & 1.80677298502757e-03_dp)*y+ & 5.78009914536630e-02_dp w3 = ((((((((((-9.49816486853687e-16_dp*y+ & 6.67922080354234e-15_dp)*y+ & 2.606163540537e-15_dp)*y+ & 1.983799950150e-12_dp)*y- & 5.400548574357e-11_dp)*y+ & 6.638043374114e-10_dp)*y- & 8.799518866802e-09_dp)*y+ & 1.791418482685e-07_dp)*y- & 2.96075397351101e-06_dp)*y+ & 3.38028206156144e-05_dp)*y- & 3.58426847857878e-04_dp)*y+ & 8.39213709428516e-03_dp w4 = (((((((((((1.33829971060180e-17_dp*y- & 3.44841877844140e-16_dp)*y+ & 4.745009557656e-15_dp)*y- & 6.033814209875e-14_dp)*y+ & 1.049256040808e-12_dp)*y- & 1.70859789556117e-11_dp)*y+ & 2.15219425727959e-10_dp)*y- & 2.52746574206884e-09_dp)*y+ & 3.27761714422960e-08_dp)*y- & 3.90387662925193e-07_dp)*y+ & 3.46340204593870e-06_dp)*y- & 2.43236345136782e-05_dp)*y+ & 3.54846978585226e-04_dp w5 = (((((((((((((2.69412277020887e-20_dp*y- & 4.24837886165685e-19_dp)*y+ & 6.030500065438e-18_dp)*y- & 9.069722758289e-17_dp)*y+ & 1.246599177672e-15_dp)*y- & 1.56872999797549e-14_dp)*y+ & 1.87305099552692e-13_dp)*y- & 2.09498886675861e-12_dp)*y+ & 2.11630022068394e-11_dp)*y- & 1.92566242323525e-10_dp)*y+ & 1.62012436344069e-09_dp)*y- & 1.23621614171556e-08_dp)*y+ & 7.72165684563049e-08_dp)*y- & 3.59858901591047e-07_dp)*y+ & 2.43682618601000e-06_dp elseif (x <= 25.0e+00_dp) then y = x-22.5e+00_dp r1 = (((((((((-1.13927848238726e-15_dp*y+ & 7.39404133595713e-15_dp)*y+ & 1.445982921243e-13_dp)*y- & 2.676703245252e-12_dp)*y+ & 5.823521627177e-12_dp)*y+ & 2.17264723874381e-10_dp)*y+ & 3.56242145897468e-09_dp)*y- & 3.03763737404491e-07_dp)*y+ & 9.46859114120901e-06_dp)*y- & 2.30896753853196e-04_dp)*y+ & 5.24663913001114e-03_dp r2 = ((((((((((2.89872355524581e-16_dp*y- & 1.22296292045864e-14_dp)*y+ & 6.184065097200e-14_dp)*y+ & 1.649846591230e-12_dp)*y- & 2.729713905266e-11_dp)*y+ & 3.709913790650e-11_dp)*y+ & 2.216486288382e-09_dp)*y+ & 4.616160236414e-08_dp)*y- & 3.32380270861364e-06_dp)*y+ & 9.84635072633776e-05_dp)*y- & 2.30092118015697e-03_dp)*y+ & 5.00845183695073e-02_dp r3 = ((((((((((1.97068646590923e-15_dp*y- & 4.89419270626800e-14_dp)*y+ & 1.136466605916e-13_dp)*y+ & 7.546203883874e-12_dp)*y- & 9.635646767455e-11_dp)*y- & 8.295965491209e-11_dp)*y+ & 7.534109114453e-09_dp)*y+ & 2.699970652707e-07_dp)*y- & 1.42982334217081e-05_dp)*y+ & 3.78290946669264e-04_dp)*y- & 8.03133015084373e-03_dp)*y+ & 1.58689469640791e-01_dp r4 = ((((((((((1.33642069941389e-14_dp*y- & 1.55850612605745e-13_dp)*y- & 7.522712577474e-13_dp)*y+ & 3.209520801187e-11_dp)*y- & 2.075594313618e-10_dp)*y- & 2.070575894402e-09_dp)*y+ & 7.323046997451e-09_dp)*y+ & 1.851491550417e-06_dp)*y- & 6.37524802411383e-05_dp)*y+ & 1.36795464918785e-03_dp)*y- & 2.42051126993146e-02_dp)*y+ & 3.97847167557815e-01_dp r5 = ((((((((((-6.07053986130526e-14_dp*y+ & 1.04447493138843e-12_dp)*y- & 4.286617818951e-13_dp)*y- & 2.632066100073e-10_dp)*y+ & 4.804518986559e-09_dp)*y- & 1.835675889421e-08_dp)*y- & 1.068175391334e-06_dp)*y+ & 3.292234974141e-05_dp)*y- & 5.94805357558251e-04_dp)*y+ & 8.29382168612791e-03_dp)*y- & 9.93122509049447e-02_dp)*y+ & 1.09857804755042e+00_dp w1 = (((((((((-9.10338640266542e-15_dp*y+ & 1.00438927627833e-13_dp)*y+ & 7.817349237071e-13_dp)*y- & 2.547619474232e-11_dp)*y+ & 1.479321506529e-10_dp)*y+ & 1.52314028857627e-09_dp)*y+ & 9.20072040917242e-09_dp)*y- & 2.19427111221848e-06_dp)*y+ & 8.65797782880311e-05_dp)*y- & 2.82718629312875e-03_dp)*y+ & 1.28718310443295e-01_dp w2 = (((((((((5.52380927618760e-15_dp*y- & 6.43424400204124e-14_dp)*y- & 2.358734508092e-13_dp)*y+ & 8.261326648131e-12_dp)*y+ & 9.229645304956e-11_dp)*y- & 5.68108973828949e-09_dp)*y+ & 1.22477891136278e-07_dp)*y- & 2.11919643127927e-06_dp)*y+ & 4.23605032368922e-05_dp)*y- & 1.14423444576221e-03_dp)*y+ & 5.06607252890186e-02_dp w3 = (((((((((3.99457454087556e-15_dp*y- & 5.11826702824182e-14_dp)*y- & 4.157593182747e-14_dp)*y+ & 4.214670817758e-12_dp)*y+ & 6.705582751532e-11_dp)*y- & 3.36086411698418e-09_dp)*y+ & 6.07453633298986e-08_dp)*y- & 7.40736211041247e-07_dp)*y+ & 8.84176371665149e-06_dp)*y- & 1.72559275066834e-04_dp)*y+ & 7.16639814253567e-03_dp w4 = (((((((((((-2.14649508112234e-18_dp*y- & 2.45525846412281e-18_dp)*y+ & 6.126212599772e-16_dp)*y- & 8.526651626939e-15_dp)*y+ & 4.826636065733e-14_dp)*y- & 3.39554163649740e-13_dp)*y+ & 1.67070784862985e-11_dp)*y- & 4.42671979311163e-10_dp)*y+ & 6.77368055908400e-09_dp)*y- & 7.03520999708859e-08_dp)*y+ & 6.04993294708874e-07_dp)*y- & 7.80555094280483e-06_dp)*y+ & 2.85954806605017e-04_dp w5 = ((((((((((((-5.63938733073804e-21_dp*y+ & 6.92182516324628e-20_dp)*y- & 1.586937691507e-18_dp)*y+ & 3.357639744582e-17_dp)*y- & 4.810285046442e-16_dp)*y+ & 5.386312669975e-15_dp)*y- & 6.117895297439e-14_dp)*y+ & 8.441808227634e-13_dp)*y- & 1.18527596836592e-11_dp)*y+ & 1.36296870441445e-10_dp)*y- & 1.17842611094141e-09_dp)*y+ & 7.80430641995926e-09_dp)*y- & 5.97767417400540e-08_dp)*y+ & 1.65186146094969e-06_dp elseif (x <= 40.0e+00_dp) then w1 = sqrt(PIo4/x) e = exp(-x) r1 = ((((((((-1.73363958895356e-06_dp*x+ & 1.19921331441483e-04_dp)*x- & 1.59437614121125e-02_dp)*x+ & 1.13467897349442e+00_dp)*x- & 4.47216460864586e+01_dp)*x+ & 1.06251216612604e+03_dp)*x- & 1.52073917378512e+04_dp)*x+ & 1.20662887111273e+05_dp)*x- & 4.07186366852475e+05_dp)*e+r15/(x-r15) r2 = ((((((((-1.60102542621710e-05_dp*x+ & 1.10331262112395e-03_dp)*x- & 1.50043662589017e-01_dp)*x+ & 1.05563640866077e+01_dp)*x- & 4.10468817024806e+02_dp)*x+ & 9.62604416506819e+03_dp)*x- & 1.35888069838270e+05_dp)*x+ & 1.06107577038340e+06_dp)*x- & 3.51190792816119e+06_dp)*e+r25/(x-r25) r3 = ((((((((-4.48880032128422e-05_dp*x+ & 2.69025112122177e-03_dp)*x- & 4.01048115525954e-01_dp)*x+ & 2.78360021977405e+01_dp)*x- & 1.04891729356965e+03_dp)*x+ & 2.36985942687423e+04_dp)*x- & 3.19504627257548e+05_dp)*x+ & 2.34879693563358e+06_dp)*x- & 7.16341568174085e+06_dp)*e+r35/(x-r35) r4 = ((((((((-6.38526371092582e-05_dp*x- & 2.29263585792626e-03_dp)*x- & 7.65735935499627e-02_dp)*x+ & 9.12692349152792e+00_dp)*x- & 2.32077034386717e+02_dp)*x+ & 2.81839578728845e+02_dp)*x+ & 9.59529683876419e+04_dp)*x- & 1.77638956809518e+06_dp)*x+ & 1.02489759645410e+07_dp)*e+r45/(x-r45) r5 = ((((((((-3.59049364231569e-05_dp*x- & 2.25963977930044e-02_dp)*x+ & 1.12594870794668e+00_dp)*x- & 4.56752462103909e+01_dp)*x+ & 1.05804526830637e+03_dp)*x- & 1.16003199605875e+04_dp)*x- & 4.07297627297272e+04_dp)*x+ & 2.22215528319857e+06_dp)*x- & 1.61196455032613e+07_dp)*e+r55/(x-r55) w5 = (((((((((-4.61100906133970e-10_dp*x+ & 1.43069932644286e-07_dp)*x- & 1.63960915431080e-05_dp)*x+ & 1.15791154612838e-03_dp)*x- & 5.30573476742071e-02_dp)*x+ & 1.61156533367153e+00_dp)*x- & 3.23248143316007e+01_dp)*x+ & 4.12007318109157e+02_dp)*x- & 3.02260070158372e+03_dp)*x+ & 9.71575094154768e+03_dp)*e+w55*w1 w4 = (((((((((-2.40799435809950e-08_dp*x+ & 8.12621667601546e-06_dp)*x- & 9.04491430884113e-04_dp)*x+ & 6.37686375770059e-02_dp)*x- & 2.96135703135647e+00_dp)*x+ & 9.15142356996330e+01_dp)*x- & 1.86971865249111e+03_dp)*x+ & 2.42945528916947e+04_dp)*x- & 1.81852473229081e+05_dp)*x+ & 5.96854758661427e+05_dp)*e+w45*w1 w3 = ((((((((1.83574464457207e-05_dp*x- & 1.54837969489927e-03_dp)*x+ & 1.18520453711586e-01_dp)*x- & 6.69649981309161e+00_dp)*x+ & 2.44789386487321e+02_dp)*x- & 5.68832664556359e+03_dp)*x+ & 8.14507604229357e+04_dp)*x- & 6.55181056671474e+05_dp)*x+ & 2.26410896607237e+06_dp)*e+w35*w1 w2 = ((((((((2.77778345870650e-05_dp*x- & 2.22835017655890e-03_dp)*x+ & 1.61077633475573e-01_dp)*x- & 8.96743743396132e+00_dp)*x+ & 3.28062687293374e+02_dp)*x- & 7.65722701219557e+03_dp)*x+ & 1.10255055017664e+05_dp)*x- & 8.92528122219324e+05_dp)*x+ & 3.10638627744347e+06_dp)*e+w25*w1 w1 = w1-0.01962e+00_dp*e-w2-w3-w4-w5 elseif (x <= 59.0e+00_dp) then w1 = sqrt(PIo4/x) y = x**3 e = y*exp(-x) r1 = (((-2.43758528330205e-02_dp*x+ & 2.07301567989771e+00_dp)*x- & 6.45964225381113e+01_dp)*x+ & 7.14160088655470e+02_dp)*e+r15/(x-r15) r2 = (((-2.28861955413636e-01_dp*x+ & 1.93190784733691e+01_dp)*x- & 5.99774730340912e+02_dp)*x+ & 6.61844165304871e+03_dp)*e+r25/(x-r25) r3 = (((-6.95053039285586e-01_dp*x+ & 5.76874090316016e+01_dp)*x- & 1.77704143225520e+03_dp)*x+ & 1.95366082947811e+04_dp)*e+r35/(x-r35) r4 = (((-1.58072809087018e+00_dp*x+ & 1.27050801091948e+02_dp)*x- & 3.86687350914280e+03_dp)*x+ & 4.23024828121420e+04_dp)*e+r45/(x-r45) r5 = (((-3.33963830405396e+00_dp*x+ & 2.51830424600204e+02_dp)*x- & 7.57728527654961e+03_dp)*x+ & 8.21966816595690e+04_dp)*e+r55/(x-r55) e = y*e w5 = ((1.35482430510942e-08_dp*x- & 3.27722199212781e-07_dp)*x+ & 2.41522703684296e-06_dp)*e+w55*w1 w4 = ((1.23464092261605e-06_dp*x- & 3.55224564275590e-05_dp)*x+ & 3.03274662192286e-04_dp)*e+w45*w1 w3 = ((1.34547929260279e-05_dp*x- & 4.19389884772726e-04_dp)*x+ & 3.87706687610809e-03_dp)*e+w35*w1 w2 = ((2.09539509123135e-05_dp*x- & 6.87646614786982e-04_dp)*x+ & 6.68743788585688e-03_dp)*e+w25*w1 w1 = w1-w2-w3-w4-w5 else w1 = sqrt(PIo4/x) r1 = r15/(x-r15) r2 = r25/(x-r25) r3 = r35/(x-r35) r4 = r45/(x-r45) r5 = r55/(x-r55) w2 = w25*w1 w3 = w35*w1 w4 = w45*w1 w5 = w55*w1 w1 = w1-w2-w3-w4-w5 end if r = (/r1, r2, r3, r4, r5/) w = (/w1, w2, w3, w4, w5/) end subroutine rys_rt5 !----------------------------------------------------------------- subroutine rys_general(x, r, w, nroots) real(KIND=dp), intent(IN) :: x real(KIND=dp), intent(OUT) :: r(:), w(:) integer, intent(IN) :: nroots real(KIND=dp) :: beta(0:mxrys) real(KIND=dp) :: rgrid(maux), wgrid(maux) real(KIND=dp) :: scr1(maux), scr2(maux) integer :: naux, rmap naux = nauxs(nroots) rmap = maprys(nroots) if (x >= xasymp(nroots)) then r(:nroots) = rts_hermit(:nroots, nroots)/x w(:nroots) = wts_hermit(:nroots, nroots)/sqrt(x) else rgrid(:naux) = rtsaux(:naux, rmap) wgrid(:naux) = wtsaux(:naux, rmap)*exp(-x*rtsaux(:naux, rmap)) call discretized_stieltjes(nroots, naux, rgrid, wgrid, & r, beta, scr1, scr2) call golub_welsch(nroots, r, beta, eps_eig, w(1:nroots)) end if r(1:nroots) = r(1:nroots)/(1-r(1:nroots)) end subroutine rys_general !----------------------------------------------------------------- subroutine discretized_stieltjes(n, naux, r, w, alpha, beta, & p_old, p) integer, intent(IN) :: n, naux real(KIND=dp), intent(IN) :: r(*), w(*) real(KIND=dp), intent(OUT) :: alpha(*), beta(*) real(KIND=dp), intent(INOUT) :: p_old(*), p(*) #ifdef DEBUG real(KIND=dp), parameter :: & sqrt_tiny = sqrt(tiny(0.0e+00_dp)), & sqrt_huge = sqrt(huge(0.0e+00_dp)) #endif real(KIND=dp) :: tmp, pp_old, pp, xpp integer :: i, k #ifdef DEBUG if (n <= 0 .or. n > naux) then error stop "Inconsistent n or naux in discretized_stieltjes" end if #endif pp_old = sum(w(:naux)) pp = dot_product(w(:naux),r(:naux)) alpha(1) = pp/pp_old beta(1) = pp_old if (n == 1) return p_old(1:naux) = 0.0e+00_dp p(1:naux) = 1.0e+00_dp do k = 1, n-1 pp = 0 xpp = 0 do i = 1, naux tmp = p(i) p(i) = (r(i)-alpha(k))*p(i) - beta(k)*p_old(i) p_old(i) = tmp pp = pp + w(i)*p(i)*p(i) xpp = xpp + r(i)*w(i)*p(i)*p(i) end do #ifdef DEBUG if (abs(xpp) > sqrt_huge .or. abs(pp) < sqrt_tiny) then error stop "Numerical issues in discretized_stieltjes" end if #endif alpha(k+1) = xpp/pp beta(k+1) = pp/pp_old pp_old = pp end do end subroutine discretized_stieltjes !----------------------------------------------------------------- subroutine golub_welsch(n, alpha, beta, eps, weight) integer, parameter :: GW_MAXIT = 30 integer, intent(IN) :: n real(KIND=dp), intent(INOUT) :: alpha(*) real(KIND=dp), intent(INOUT) :: beta(0:*) real(KIND=dp), intent(OUT) :: weight(*) real(KIND=dp), intent(IN) :: eps integer :: i, j, l, m real(KIND=dp) :: b, c, f, g, p, r, s, mu0 #ifdef DEBUG if (any(beta(:n) < 0.0e+00_dp)) then error stop "Negative beta coef. in golub_welsch" end if #endif if (n == 1) then weight(1) = beta(0) return end if mu0 = beta(0) beta(1:n-1) = sqrt(beta(1:n-1)) beta(n) = 0 weight(1) = 1.0e+00_dp weight(2:n) = 0.0e+00_dp large: do l = 1, n j = 0 iter: do do m = l, n-1 if (abs(beta(m)) <= & eps*(abs(alpha(m))+abs(alpha(m+1)))) exit end do if (m == l) exit iter if (j == GW_MAXIT) then #ifdef DEBUG error stop "Golub-Welsch procedure can't converge" #endif return end if j = j+1 g = (alpha(l+1)-alpha(l))/(2.0e+00_dp*beta(l)) r = sqrt(g*g+1.0e+00_dp) g = alpha(m)-alpha(l)+beta(l)/(g+sign(r, g)) s = 1.0e+00_dp c = 1.0e+00_dp p = 0.0e+00_dp do i = m-1, l, -1 f = s*beta(i) b = c*beta(i) if (abs(f) < abs(g)) then s = f/g r = sqrt(s*s+1.0e+00_dp) beta(i+1) = g*r c = 1.0e+00_dp/r s = s*c else c = g/f r = sqrt(c*c+1.0e+00_dp) beta(i+1) = f*r s = 1.0e+00_dp/r c = c*s end if g = alpha(i+1)-p r = (alpha(i)-g)*s+2.0e+00_dp*c*b p = s*r alpha(i+1) = g+p g = c*r-b f = weight(i+1) weight(i+1) = s*weight(i) + c*f weight(i) = c*weight(i) - s*f end do alpha(l) = alpha(l)-p beta(l) = g beta(m) = 0.0e+00_dp end do iter end do large weight(1:n) = mu0*weight(1:n)*weight(1:n) end subroutine golub_welsch end module rys