strvh0 Function

  FUNCTION strvh0(xvalue) RESULT(fn_val)
    !! CALGO 757 Struve function \(\mathbf{H}_0(x)\).
    !
    ! DESCRIPTION:
    !     This function calculates the value of the Struve function
    !     of order 0, denoted H0(x), for the argument XVALUE, defined
    !
    !         STRVHO(x) = (2/pi) integral{0 to pi/2} sin(x cos(t)) dt
    !
    !     H0 also satisfies the second-order equation
    !
    !                 x*D(Df) + Df + x*f = 2x/pi
    !
    !     The code uses Chebyshev expansions whose coefficients are given to 20D.
    !
    !   ERROR RETURNS:
    !     As the asymptotic expansion of H0 involves the Bessel function
    !     of the second kind Y0, there is a problem for large x, since
    !     we cannot accurately calculate the value of Y0. An error message
    !     is printed and STRVH0 returns the value 0.0.
    !
    !   MACHINE-DEPENDENT CONSTANTS:
    !     NTERM1 - The no. of terms to be used in the array ARRH0. The
    !              recommended value is such that
    !                      ABS(ARRH0(NTERM1)) < EPS/100.
    !     NTERM2 - The no. of terms to be used in the array ARRH0A. The
    !              recommended value is such that
    !                      ABS(ARRH0A(NTERM2)) < EPS/100.
    !     NTERM3 - The no. of terms to be used in the array AY0ASP. The
    !              recommended value is such that
    !                      ABS(AY0ASP(NTERM3)) < EPS/100.
    !     NTERM4 - The no. of terms to be used in the array AY0ASQ. The
    !              recommended value is such that
    !                      ABS(AY0ASQ(NTERM4)) < EPS/100.
    !     XLOW - The value for which H0(x) = 2*x/pi to machine precision, if
    !            abs(x) < XLOW. The recommended value is
    !                      XLOW = 3 * SQRT(EPSNEG)
    !     XHIGH - The value above which we are unable to calculate Y0 with
    !             any reasonable accuracy. An error message is printed and
    !              STRVH0 returns the value 0.0. The recommended value is
    !                      XHIGH = 1/EPS.
    !
    !     For values of EPS and EPSNEG refer to the file MACHCON.TXT.

    REAL(wp), INTENT(IN) :: xvalue
    REAL(wp)             :: fn_val

    INTEGER  :: indsgn, nterm1, nterm2, nterm3, nterm4
    REAL(wp) :: h0as, t, x, xhigh, xlow, xmp4, xsq, y0p, y0q, y0val

    REAL(wp), PARAMETER :: zero = 0.0_wp, half = 0.5_wp, one = 1.0_wp, &
                           eight = 8.0_wp, eleven = 11.0_wp, twenty = 20.0_wp, &
                           onehun = 100.0_wp, sixtp5 = 60.5_wp, &
                           two62 = 262.0_wp, thr2p5 = 302.5_wp, &
                           piby4 = 0.78539816339744830962_wp, &
                           rt2bpi = 0.79788456080286535588_wp, &
                           twobpi = 0.63661977236758134308_wp

    REAL(wp), PARAMETER :: arrh0(0:19) = [ &
      0.28696487399013225740_wp, -0.25405332681618352305_wp, &
      0.20774026739323894439_wp, -0.20364029560386585140_wp, &
      0.12888469086866186016_wp, -0.4825632815622261202e-1_wp, &
      0.1168629347569001242e-1_wp, -0.198118135642418416e-2_wp, &
      0.24899138512421286e-3_wp, -0.2418827913785950e-4_wp, &
      0.187437547993431e-5_wp, -0.11873346074362e-6_wp, &
      0.626984943346e-8_wp, -0.28045546793e-9_wp, &
      0.1076941205e-10_wp, -0.35904793e-12_wp, &
      0.1049447e-13_wp, -0.27119e-15_wp, &
      0.624e-17_wp, -0.13e-18_wp &
    ]
    REAL(wp), PARAMETER :: arrh0a(0:20) = [ &
      1.99291885751992305515_wp, -0.384232668701456887e-2_wp, &
     -0.32871993712353050e-3_wp, -0.2941181203703409e-4_wp, &
     -0.267315351987066e-5_wp, -0.24681031075013e-6_wp, &
     -0.2295014861143e-7_wp, -0.215682231833e-8_wp, &
     -0.20303506483e-9_wp, -0.1934575509e-10_wp, &
     -0.182773144e-11_wp, -0.17768424e-12_wp, &
     -0.1643296e-13_wp, -0.171569e-14_wp, &
     -0.13368e-15_wp, -0.2077e-16_wp, &
      0.2e-19_wp, -0.55e-18_wp, &
      0.10e-18_wp, -0.4e-19_wp, &
      0.1e-19_wp &
    ]
    REAL(wp), PARAMETER :: ay0asp(0:12) = [ &
      1.99944639402398271568_wp, -0.28650778647031958e-3_wp, &
     -0.1005072797437620e-4_wp, -0.35835941002463e-6_wp, &
     -0.1287965120531e-7_wp, -0.46609486636e-9_wp, &
     -0.1693769454e-10_wp, -0.61852269e-12_wp, &
     -0.2261841e-13_wp, -0.83268e-15_wp, &
     -0.3042e-16_wp, -0.115e-17_wp, &
     -0.4e-19_wp &
    ]
    REAL(wp), PARAMETER :: ay0asq(0:13) = [ &
      1.99542681386828604092_wp, -0.236013192867514472e-2_wp, &
     -0.7601538908502966e-4_wp, -0.256108871456343e-5_wp, &
     -0.8750292185106e-7_wp, -0.304304212159e-8_wp, &
     -0.10621428314e-9_wp, -0.377371479e-11_wp, &
     -0.13213687e-12_wp, -0.488621e-14_wp, &
     -0.15809e-15_wp, -0.762e-17_wp, &
     -0.3e-19_wp, -0.3e-19_wp &
    ]

    ! Start computation
    x = xvalue
    indsgn = 1
    IF (x < zero) THEN
      x = -x
      indsgn = -1
    END IF

    ! Compute the machine-dependent constants.
    h0as = eps_wp
    xhigh = one / (2*eps_wp)

    ! Error test
    IF (ABS(xvalue) > xhigh) THEN
      fn_val = zero
      RETURN
    END IF

    ! continue with machine constants
    t = h0as / onehun
    IF (x <= eleven) THEN
      DO  nterm1 = 19, 0, -1
        IF (ABS(arrh0(nterm1)) > t) EXIT
      END DO
      y0p = SQRT(h0as)
      xlow = y0p + y0p + y0p
    ELSE
      DO  nterm2 = 20, 0, -1
        IF (ABS(arrh0a(nterm2)) > t) EXIT
      END DO
      DO  nterm3 = 12, 0, -1
        IF (ABS(ay0asp(nterm3)) > t) EXIT
      END DO
      DO  nterm4 = 13, 0, -1
        IF (ABS(ay0asq(nterm4)) > t) EXIT
      END DO
    END IF

    ! Code for abs(x) <= 11
    IF (x <= eleven) THEN
      IF (x < xlow) THEN
        fn_val = twobpi * x
      ELSE
        t = ((x*x)/sixtp5-half) - half
        fn_val = twobpi * x * cheval(nterm1,arrh0,t)
      END IF
    ELSE
      ! Code for abs(x) > 11
      xsq = x * x
      t = (two62-xsq) / (twenty+xsq)
      y0p = cheval(nterm3,ay0asp,t)
      y0q = cheval(nterm4,ay0asq,t) / (eight*x)
      xmp4 = x - piby4
      y0val = y0p * SIN(xmp4) - y0q * COS(xmp4)
      y0val = y0val * rt2bpi / SQRT(x)
      t = (thr2p5-xsq) / (sixtp5+xsq)
      h0as = twobpi * cheval(nterm2,arrh0a,t) / x
      fn_val = y0val + h0as
    END IF
    IF (indsgn == -1) fn_val = -fn_val
    RETURN
  END FUNCTION strvh0