module rys_deriv
! Analytic X-derivatives of the closed-form Rys roots/weights
! rys_rt1..rys_rt5, via complex-step differentiation of type-substituted
! clones of the value routines in rys.F90. Gate 1 of the native
! nuclear-attraction Hessian. Generated mechanically (no hand transcription).
  use precision, only: dp
  use constants, only: pi
  implicit none
  real(dp), parameter :: pio4 = pi/4.0_dp
  private
  public :: rys_rt1_d, rys_rt2_d, rys_rt3_d, rys_rt4_d, rys_rt5_d
contains
  subroutine rys_rt1_d(x, dr, dw)
    ! Complex-step analytic derivatives du/dX, dw/dX of rys_rt1 (1 root).
    real(KIND=dp), intent(IN) :: x
    real(KIND=dp), intent(OUT) :: dr(1), dw(1)
    complex(KIND=dp) :: xc, rc, wc
    real(KIND=dp), parameter :: hcs = 1.0e-20_dp
    xc = cmplx(x, hcs, KIND=dp)
    call rys_rt1_csval(xc, rc, wc)
    dr(1) = aimag(rc)/hcs
    dw(1) = aimag(wc)/hcs
  end subroutine rys_rt1_d

  subroutine rys_rt1_csval(x, r, w)

    complex(KIND=dp), intent(IN) :: &
      x
    complex(KIND=dp), intent(OUT) :: &
      r
    complex(KIND=dp), intent(OUT) :: &
      w

    complex(KIND=dp) :: &
      f1, y, recx, e

    if (real(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 (real(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 (real(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 (real(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 (real(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 (real(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 (real(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_csval

  subroutine rys_rt2_d(x, dr, dw)
    ! Complex-step analytic derivatives du_t/dX, dw_t/dX of rys_rt2.
    ! dr(t)=du_t/dX, dw(t)=dw_t/dX, to machine precision (no cancellation).
    real(KIND=dp), intent(IN) :: x
    real(KIND=dp), intent(OUT) :: dr(2), dw(2)
    complex(KIND=dp) :: xc, rc(2), wc(2)
    real(KIND=dp), parameter :: hcs = 1.0e-20_dp
    xc = cmplx(x, hcs, KIND=dp)
    call rys_rt2_csval(xc, rc, wc)
    dr = aimag(rc)/hcs
    dw = aimag(wc)/hcs
  end subroutine rys_rt2_d

  subroutine rys_rt2_csval(x, r, w)

    complex(KIND=dp), intent(IN) :: &
      x
    complex(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

    complex(KIND=dp) :: &
      f1, y, recx, e, r1, r2, w1, w2

    if (real(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 (real(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 (real(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 (real(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 (real(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 (real(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 (real(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 (real(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_csval

  subroutine rys_rt3_d(x, dr, dw)
    ! Complex-step analytic derivatives du_t/dX, dw_t/dX of rys_rt3.
    ! dr(t)=du_t/dX, dw(t)=dw_t/dX, to machine precision (no cancellation).
    real(KIND=dp), intent(IN) :: x
    real(KIND=dp), intent(OUT) :: dr(3), dw(3)
    complex(KIND=dp) :: xc, rc(3), wc(3)
    real(KIND=dp), parameter :: hcs = 1.0e-20_dp
    xc = cmplx(x, hcs, KIND=dp)
    call rys_rt3_csval(xc, rc, wc)
    dr = aimag(rc)/hcs
    dw = aimag(wc)/hcs
  end subroutine rys_rt3_d

  subroutine rys_rt3_csval(x, r, w)

    complex(KIND=dp), intent(IN) :: &
      x
    complex(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

    complex(KIND=dp) :: &
      f1, f2, y, recx, e, &
      a1, a2, t1, t2, t3, r1, r2, r3, w1, w2, w3

    if (real(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 (real(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 (real(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 (real(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 (real(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 (real(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 (real(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 (real(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 (real(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_csval

  subroutine rys_rt4_d(x, dr, dw)
    ! Complex-step analytic derivatives du_t/dX, dw_t/dX of rys_rt4.
    ! dr(t)=du_t/dX, dw(t)=dw_t/dX, to machine precision (no cancellation).
    real(KIND=dp), intent(IN) :: x
    real(KIND=dp), intent(OUT) :: dr(4), dw(4)
    complex(KIND=dp) :: xc, rc(4), wc(4)
    real(KIND=dp), parameter :: hcs = 1.0e-20_dp
    xc = cmplx(x, hcs, KIND=dp)
    call rys_rt4_csval(xc, rc, wc)
    dr = aimag(rc)/hcs
    dw = aimag(wc)/hcs
  end subroutine rys_rt4_d

  subroutine rys_rt4_csval(x, r, w)

    complex(KIND=dp), intent(IN) :: &
      x
    complex(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

    complex(KIND=dp) :: &
      y, recx, e, r1, r2, r3, r4, w1, w2, w3, w4

    if (real(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 (real(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 (real(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 (real(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 (real(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 (real(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 (real(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 (real(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 (real(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_csval

  subroutine rys_rt5_d(x, dr, dw)
    ! Complex-step analytic derivatives du_t/dX, dw_t/dX of rys_rt5.
    ! dr(t)=du_t/dX, dw(t)=dw_t/dX, to machine precision (no cancellation).
    real(KIND=dp), intent(IN) :: x
    real(KIND=dp), intent(OUT) :: dr(5), dw(5)
    complex(KIND=dp) :: xc, rc(5), wc(5)
    real(KIND=dp), parameter :: hcs = 1.0e-20_dp
    xc = cmplx(x, hcs, KIND=dp)
    call rys_rt5_csval(xc, rc, wc)
    dr = aimag(rc)/hcs
    dw = aimag(wc)/hcs
  end subroutine rys_rt5_d

  subroutine rys_rt5_csval(x, r, w)

    complex(KIND=dp), intent(IN) :: &
      x
    complex(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

    complex(KIND=dp) :: y, e, r1, r2, r3, r4, r5, w1, w2, w3, w4, w5

    if (real(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 (real(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 (real(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 (real(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 (real(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 (real(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 (real(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 (real(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 (real(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_csval

end module rys_deriv
