! Licensed under the ACM Software License Agreement
! Copyright © 1970–2012 Association for Computing Machinery (ACM)
! See ColSpecF LICENSE file for details.

module calgo_757
!* # CALGO 757
! Algorithm 757. 
!
! Procedures:
!
! - `strvh0`: Struve function \(\mathbf{H}_0(x)\)
! - `strvh1`: Struve function \(\mathbf{H}_1(x)\)
!
! Other CALGO 757 procedures **not** yet included in this module:
!
! - `abram0`, `abram1`, `abram2`: Abramowitz functions \(f_m(x) = 
!     \int_{0}^{\infty} t^m e^{-t^2 -\frac{x}{t}} \, dt, \, \text{for } m = 0, 1, 2\)
! - `airygi`: Modified Airy function \(\mathrm{Gi}(x)\)
! - `airyhi`: Modified Airy function \(\mathrm{Hi}(x)\)
! - `airint`: Integral of the Airy function \(\mathrm{Ai}(x)\) from \(0\) to \(x\)
! - `birint`: Integral of the Airy function \(\mathrm{Bi}(x)\) from \(0\) to \(x\)
! - `j0int`: Integral of the Bessel function \(J_0(x)\) from \(0\) to \(x\)
! - `y0int`: Integral of the Bessel function \(Y_0(x)\) from \(0\) to \(x\)
! - `i0int`: Integral of the modified Bessel function \(I_0(x)\) from \(0\) to \(x\)
! - `k0int`: Integral of the modified Bessel function \(K_0(x)\) from \(0\) to \(x\)
! - `debye1`, `debye2`, `debye3`, `debye4`: Debye functions \(D_n(x) = 
!     \frac{n}{x^n} \int_{0}^{x} \frac{t^n}{e^t-1} \, dt,
!     \, \text{for } n = 1, 2, 3, 4\)
! - `strvl0`: Modified Struve function \(\mathbf{L}_0(x)\)
! - `strvl1`: Modified Struve function \(\mathbf{L}_1(x)\)
! - `i0ml0`: Difference between the modified Bessel and Struve functions
!     \(I_0(x) - \mathbf{L}_0(x)\)
! - `i1ml1`: Difference between the modified Bessel and Struve functions
!     \(I_1(x) - \mathbf{L}_1(x)\)
! - `synch1`: Synchrotron radiation function
!     \(f_1(x) = x \int_{x}^{\infty} K_{5/3}(t) \, dt\)
! - `synch2`: Synchrotron radiation function
!     \(f_2(x) = x K_{2/3}(x)\)
! - `tran02`, `tran03`, ..., `tran09`: Transport integrals \(J_n(x) = 
!     \int_{0}^{x} \frac{t^n e^t}{(e^t-1)^2} \, dt,
!     \, \text{for } n = 2, 3, \ldots, 9\) 
! - `atnint`: Inverse-tangent integral \(\int_{0}^{x} \frac{\arctan t}{t} \, dt\)
! - `clausn`: Clausen's integral
!     \(-\int_{0}^{x} \ln{\lvert 2 \sin\frac{t}{2} \rvert} \, dt\)
! - `exp3`: \(\int_{0}^{x} e^{-t^3} \, dt\)
! - `goodst`: \(\int_{0}^{\infty} \frac{e^{-t^2}}{t + x} \, dt\)
! - `lobach`: Lobachevski's integral
!     \(-\int_{0}^{\infty} \log \lvert \cos t \rvert \, dt\)
! - `strom`: Stromgren's integral
!     \(\frac{15}{4 \pi^4} \int_{0}^{x} \frac{t^7 e^{2 t}}{(e^t - 1)^3} \, dt\)
!
! ## Author
! Allan J. MacLeod
!
! ## History
! - 1996-01-23 - Allan J. MacLeod
!     - Original code.
! - 2000-06-13 - Allan J. MacLeod
!     - F77 code distributed by ACM: <https://calgo.acm.org/757.zip>
! - 2001-10-10 - Alan Miller
!     - F90 code adaptation by Alan Miller:
!       <https://jblevins.org/mirror/amiller/toms757.zip>
! - 2025-07-27 - Rodrigo Castro (GitHub: rodpcastro)
!     - Retained only the subroutines `strvh0` and `strvh1`; additional subroutines
!       will be included as required.
!     - Replaced array constructor `(/.../)` by the less verbose `[...]`.
!     - Replaced `dp` (double precision) by `wp` (working precision).
!     - Replaced `EPSILON(0.0_dp)` by `eps_wp` (CSF working precision machine epsilon).
!     - Removed subroutine `errprn` that prints error messages, which are unnecessary
!       for `strvh0` and `strvh1`.
!
! ## References
! 1. Allan J. MacLeod. 1996. Algorithm 757: MISCFUN, a software package to compute
!    uncommon special functions. ACM Trans. Math. Softw. 22, 3 (Sept. 1996), 288–301.
!*   <https://doi.org/10.1145/232826.232846>

  use csf_kinds, only: wp
  use csf_numerror, only: eps_wp

  implicit none
  private
  public :: strvh0, strvh1

contains

  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


  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


  FUNCTION cheval(n, a, t) RESULT(fn_val)
    ! This function evaluates a Chebyshev series, using the Clenshaw method
    ! with Reinsch modification, as analysed in the paper by Oliver.
    !
    ! INPUT PARAMETERS
    !     N - INTEGER - The no. of terms in the sequence
    !     A - REAL(wp) ARRAY, dimension 0 to N - The coefficients of
    !         the Chebyshev series
    !     T - REAL(wp) - The value at which the series is to be evaluated
    !
    ! REFERENCES
    !      "An error analysis of the modified Clenshaw method for
    !       evaluating Chebyshev and Fourier series" J. Oliver,
    !       J.I.M.A., vol. 20, 1977, pp379-391

    INTEGER, INTENT(IN)  :: n
    REAL(wp), INTENT(IN) :: a(0:n)
    REAL(wp), INTENT(IN) :: t
    REAL(wp)             :: fn_val
    INTEGER  :: i
    REAL(wp) :: d1, d2, tt, u0, u1, u2
    REAL(wp), PARAMETER :: zero = 0.0_wp, half = 0.5_wp, test = 0.6_wp, two = 2.0_wp

    u1 = zero

    ! If ABS(T) < 0.6 use the standard Clenshaw method
    IF (ABS(t) < test) THEN
      u0 = zero
      tt = t + t
      DO  i = n, 0, -1
        u2 = u1
        u1 = u0
        u0 = tt * u1 + a(i) - u2
      END DO
      fn_val = (u0-u2) / two
    ELSE
      ! If ABS(T) >= 0.6 use the Reinsch modification
      d1 = zero
      
      ! T >= 0.6 code
      IF (t > zero) THEN
        tt = (t-half) - half
        tt = tt + tt
        DO  i = n, 0, -1
          d2 = d1
          u2 = u1
          d1 = tt * u2 + a(i) + d2
          u1 = d1 + u2
        END DO
        fn_val = (d1+d2) / two
      ELSE
        ! T <= -0.6 code
        tt = (t+half) + half
        tt = tt + tt
        DO  i = n, 0, -1
          d2 = d1
          u2 = u1
          d1 = tt * u2 + a(i) - d2
          u1 = d1 - u2
        END DO
        fn_val = (d1-d2) / two
      END IF
    END IF
    RETURN
  END FUNCTION cheval

end module calgo_757