CALGO 757 Struve function .
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(in) | :: | xvalue |
FUNCTION strvh1(xvalue) RESULT(fn_val) !! CALGO 757 Struve function \(\mathbf{H}_1(x)\). ! ! DESCRIPTION: ! This function calculates the value of the Struve function ! of order 1, denoted H1(x), for the argument XVALUE, defined as ! 2 ! STRVH1(x) = (2x/pi) integral{0 to pi/2} sin( x cos(t))*sin t dt ! ! H1 also satisfies the second-order differential equation ! ! 2 2 2 2 ! x * D f + x * Df + (x - 1)f = 2x / pi ! ! The code uses Chebyshev expansions with the coefficients given to 20D. ! ! ERROR RETURNS: ! As the asymptotic expansion of H1 involves the Bessel function ! of the second kind Y1, there is a problem for large x, since ! we cannot accurately calculate the value of Y1. An error message ! is printed and STRVH1 returns the value 0.0. ! ! MACHINE-DEPENDENT CONSTANTS: ! NTERM1 - The no. of terms to be used in the array ARRH1. The ! recommended value is such that ! ABS(ARRH1(NTERM1)) < EPS/100. ! NTERM2 - The no. of terms to be used in the array ARRH1A. The ! recommended value is such that ! ABS(ARRH1A(NTERM2)) < EPS/100. ! NTERM3 - The no. of terms to be used in the array AY1ASP. The ! recommended value is such that ! ABS(AY1ASP(NTERM3)) < EPS/100. ! NTERM4 - The no. of terms to be used in the array AY1ASQ. The ! recommended value is such that ! ABS(AY1ASQ(NTERM4)) < EPS/100. ! XLOW1 - The value of x, below which H1(x) set to zero, if ! abs(x)<XLOW1. The recommended value is ! XLOW1 = 1.5 * SQRT(XMIN) ! XLOW2 - The value for which H1(x) = 2*x*x/pi to machine precision, if ! abs(x) < XLOW2. The recommended value is ! XLOW2 = SQRT(15*EPSNEG) ! XHIGH - The value above which we are unable to calculate Y1 with ! any reasonable accuracy. An error message is printed and ! STRVH1 returns the value 0.0. The recommended value is ! XHIGH = 1/EPS. ! ! For values of EPS, EPSNEG and XMIN refer to the file MACHCON.TXT. REAL(wp), INTENT(IN) :: xvalue REAL(wp) :: fn_val INTEGER :: nterm1, nterm2, nterm3, nterm4 REAL(wp) :: h1as, t, x, xhigh, xlow1, xlow2, xm3p4, xsq, y1p, y1q, y1val REAL(wp), PARAMETER :: zero = 0.0_wp, half = 0.5_wp, eight = 8.0_wp, & nine = 9.0_wp, fiften = 15.0_wp, twenty = 20.0_wp, & onehun = 100.0_wp, fortp5 = 40.5_wp, & one82 = 182.0_wp, tw02p5 = 202.5_wp, & rt2bpi = 0.79788456080286535588_wp, & thpby4 = 2.35619449019234492885_wp, & twobpi = 0.63661977236758134308_wp REAL(wp), PARAMETER :: arrh1(0:17) = [ & 0.17319061083675439319_wp, -0.12606917591352672005_wp, & 0.7908576160495357500e-1_wp, -0.3196493222321870820e-1_wp, & 0.808040581404918834e-2_wp, -0.136000820693074148e-2_wp, & 0.16227148619889471e-3_wp, -0.1442352451485929e-4_wp, & 0.99219525734072e-6_wp, -0.5441628049180e-7_wp, & 0.243631662563e-8_wp, -0.9077071338e-10_wp, & 0.285926585e-11_wp, -0.7716975e-13_wp, & 0.180489e-14_wp, -0.3694e-16_wp, & 0.67e-18_wp, -0.1e-19_wp & ] REAL(wp), PARAMETER :: arrh1a(0:21) = [ & 2.01083504951473379407_wp, 0.592218610036099903e-2_wp, & 0.55274322698414130e-3_wp, 0.5269873856311036e-4_wp, & 0.506374522140969e-5_wp, 0.49028736420678e-6_wp, & 0.4763540023525e-7_wp, 0.465258652283e-8_wp, & 0.45465166081e-9_wp, 0.4472462193e-10_wp, & 0.437308292e-11_wp, 0.43568368e-12_wp, & 0.4182190e-13_wp, 0.441044e-14_wp, & 0.36391e-15_wp, 0.5558e-16_wp, & -0.4e-19_wp, 0.163e-17_wp, & -0.34e-18_wp, 0.13e-18_wp, & -0.4e-19_wp, 0.1e-19_wp & ] REAL(wp), PARAMETER :: ay1asp(0:14) = [ & 2.00135240045889396402_wp, 0.71104241596461938e-3_wp, & 0.3665977028232449e-4_wp, 0.191301568657728e-5_wp, & 0.10046911389777e-6_wp, 0.530401742538e-8_wp, & 0.28100886176e-9_wp, 0.1493886051e-10_wp, & 0.79578420e-12_wp, 0.4252363e-13_wp, & 0.227195e-14_wp, 0.12216e-15_wp, & 0.650e-17_wp, 0.36e-18_wp, & 0.2e-19_wp & ] REAL(wp), PARAMETER :: ay1asq(0:15) = [ & 5.99065109477888189116_wp, -0.489593262336579635e-2_wp, & -0.23238321307070626e-3_wp, -0.1144734723857679e-4_wp, & -0.57169926189106e-6_wp, -0.2895516716917e-7_wp, & -0.147513345636e-8_wp, -0.7596537378e-10_wp, & -0.390658184e-11_wp, -0.20464654e-12_wp, & -0.1042636e-13_wp, -0.57702e-15_wp, & -0.2550e-16_wp, -0.210e-17_wp, & 0.2e-19_wp, -0.2e-19_wp & ] ! Start computation x = ABS(xvalue) ! Compute the machine-dependent constants. xhigh = (half+half) / (2*eps_wp) ! Error test IF (x > xhigh) THEN fn_val = zero RETURN END IF ! continue with machine constants h1as = eps_wp t = h1as / onehun IF (x <= nine) THEN DO nterm1 = 17, 0, -1 IF (ABS(arrh1(nterm1)) > t) EXIT END DO xlow1 = half * SQRT(TINY(zero)) xlow1 = xlow1 + xlow1 + xlow1 xlow2 = SQRT(fiften*h1as) ELSE DO nterm2 = 21, 0, -1 IF (ABS(arrh1a(nterm2)) > t) EXIT END DO DO nterm3 = 14, 0, -1 IF (ABS(ay1asp(nterm3)) > t) EXIT END DO DO nterm4 = 15, 0, -1 IF (ABS(ay1asq(nterm4)) > t) EXIT END DO END IF ! Code for abs(x) <= 9 IF (x <= nine) THEN IF (x < xlow1) THEN fn_val = zero ELSE xsq = x * x IF (x < xlow2) THEN fn_val = twobpi * xsq ELSE t = (xsq/fortp5-half) - half fn_val = twobpi * xsq * cheval(nterm1,arrh1,t) END IF END IF ELSE ! Code for abs(x) > 9 xsq = x * x t = (one82-xsq) / (twenty+xsq) y1p = cheval(nterm3,ay1asp,t) y1q = cheval(nterm4,ay1asq,t) / (eight*x) xm3p4 = x - thpby4 y1val = y1p * SIN(xm3p4) + y1q * COS(xm3p4) y1val = y1val * rt2bpi / SQRT(x) t = (tw02p5-xsq) / (fortp5+xsq) h1as = twobpi * cheval(nterm2,arrh1a,t) fn_val = y1val + h1as END IF RETURN END FUNCTION strvh1