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

module calgo_715
!* # CALGO 715
! Algorithm 715. 
!
! Procedures:
!
! - `caljy0`: Bessel functions \(J_0(x)\) and \(Y_0(x)\)
! - `caljy1`: Bessel functions \(J_1(x)\) and \(Y_1(x)\)
!
! Other CALGO 715 procedures **not** yet included in this module:
!
! - `calcei`: Exponential integrals \(\mathrm{Ei}(x)\) and \(\mathrm{E}_1(x)\)
! - `gamma`: Gamma function \(\Gamma(x)\)
! - `algama`: Log Gamma function \(\ln\Gamma(x)\)
! - `psi`:  Digamma function \(\psi(x)\)
! - `calerf`: Error function \(\mathrm{erf}(x)\)
! - `daw`: Dawson's integral \(e^{-x^2} \int_{0}^{x} e^{t^2} dt\)
! - `anorm`: Normal distribution
!    \(\frac{1}{\sqrt{2\pi}} \int_{-\infty}^{x} e^{\frac{-t^2}{2}} dt\)
! - `rjbesl`: Bessel function \(J_{\nu}(x)\)
! - `rybesl`: Bessel function \(Y_{\nu}(x)\)
! - `calci0`: Bessel function \(I_0(x)\)
! - `calci1`: Bessel function \(I_1(x)\)
! - `ribesl`: Bessel function \(I_{\nu}(x)\)
! - `calck0`: Bessel function \(K_0(x)\)
! - `calck1`: Bessel function \(K_1(x)\)
! - `rkbesl`: Bessel function \(K_{\nu}(x)\)
!
! ## Author
! W. J. Cody
!
! ## History
! - 1992-03-15 - W. J. Cody
!     - Original code.
! - 2000-03-30 - W. J. Cody
!     - F77 code distributed by ACM: <https://calgo.acm.org/715.zip>
! - 2003-01-14 - Alan Miller
!     - F90 code adaptation by Alan Miller:
!       <https://jblevins.org/mirror/amiller/toms715.zip>
! - 2025-06-17 - Rodrigo Castro (GitHub: rodpcastro)
!     - Retained only subroutines `caljy0` and `caljy1`; additional subroutines will
!       be included as required.
!     - Replaced `dp` (double precision) by `wp` (working precision).
!     - Replaced array constructor `(/.../)` by the less verbose `[...]`.
!     - `xmax` is set such that `xmax` \( > \frac{2}{\pi \epsilon^2}\) , where
!       \(\epsilon\) is the machine epsilon.
!     - `xinf` is set to `huge(0.0_wp)`.
!
! ## References
! 1. W. J. Cody. 1993. Algorithm 715: SPECFUN–a portable FORTRAN package of special
!    function routines and test drivers. ACM Trans. Math. Softw. 19, 1 (March 1993),
!*   22–30. <https://doi.org/10.1145/151271.151273>

  use csf_kinds, only: wp

  implicit none
  private
  public :: caljy0, caljy1

contains

  SUBROUTINE caljy0(arg, result, jint)
    !* CALGO 715 Bessel functions \(J_0(x)\) and \(Y_0(x)\).
    !
    ! \(x \in \mathbb{R}\) for \(J_0(x)\)
    ! 
    ! \(\lbrace x \in \mathbb{R} \mid x \gt 0 \rbrace\) for \(Y_0(x)\)
    !
    ! To obtain:
    ! 
    ! - `y` = \(J_0(x)\), call `caljy0(x, y, 0)`
    ! - `y` = \(Y_0(x)\), call `caljy0(x, y, 1)`
    !*

    !-------------------------------------------------------------------
    ! This packet computes zero-order Bessel functions of the first and
    ! second kind (J0 and Y0), for real arguments X, where 0 < X <= XMAX
    ! for Y0, and |X| <= XMAX for J0.  It contains two function-type
    ! subprograms, BESJ0  and  BESY0, and one subroutine-type
    ! subprogram, CALJY0.  The calling statements for the primary
    ! entries are:
    !
    ! Y = BESJ0(X)
    !
    ! and
    !
    ! Y = BESY0(X),
    !
    ! where the entry points correspond to the functions J0(X) and Y0(X),
    ! respectively. The routine CALJY0 is intended for internal packet
    ! use only, all computations within the packet being concentrated in
    ! this one routine. The function subprograms invoke CALJY0 with
    ! the statement
    !
    ! CALL CALJY0(ARG, RESULT, JINT),
    !
    ! where the parameter usage is as follows:
    !
    ! Function     Parameters for CALJY0
    ! call         ARG                    RESULT    JINT
    !
    ! BESJ0(ARG)   |ARG| .LE. XMAX        J0(ARG)   0
    ! BESY0(ARG)   0 .LT. ARG .LE. XMAX   Y0(ARG)   1
    !
    ! The main computation uses unpublished minimax rational
    ! approximations for X .LE. 8.0, and an approximation from the
    ! book  Computer Approximations  by Hart, et. al., Wiley and Sons,
    ! New York, 1968, for arguments larger than 8.0  Part of this
    ! transportable packet is patterned after the machine-dependent
    ! FUNPACK program BESJ0(X), but cannot match that version for
    ! efficiency or accuracy. This version uses rational functions
    ! that are theoretically accurate to at least 18 significant decimal
    ! digits for X <= 8, and at least 18 decimal places for X > 8. The
    ! accuracy achieved depends on the arithmetic system, the compiler,
    ! the intrinsic functions, and proper selection of the machine-
    ! dependent constants.
    !
    !-------------------------------------------------------------------
    !
    ! The following machine-dependent constants must be declared in
    ! DATA statements. IEEE values are provided as a default.
    !
    ! XINF   = largest positive machine number
    ! XMAX   = largest acceptable argument. The functions AINT, SIN
    !          and COS must perform properly for ABS(X) .LE. XMAX.
    !          We recommend that XMAX be a small integer multiple of
    !          sqrt(1/eps), where eps is the smallest positive number
    !          such that  1+eps > 1.
    ! XSMALL = positive argument such that 1.0-(X/2)**2 = 1.0
    !          to machine precision for all ABS(X) .LE. XSMALL.
    !          We recommend that  XSMALL < sqrt(eps)/beta, where beta
    !          is the floating-point radix (usually 2 or 16).
    !
    ! Approximate values for some important machines are
    !
    !                       eps       XMAX      XSMALL    XINF
    !
    ! CDC 7600      (S.P.)  7.11E-15  1.34E+08  2.98E-08  1.26E+322
    ! CRAY-1        (S.P.)  7.11E-15  1.34E+08  2.98E-08  5.45E+2465
    ! IBM PC (8087) (S.P.)  5.96E-08  8.19E+03  1.22E-04  3.40E+38
    ! IBM PC (8087) (D.P.)  1.11D-16  2.68D+08  3.72D-09  1.79D+308
    ! IBM 195       (D.P.)  2.22D-16  6.87D+09  9.09D-13  7.23D+75
    ! UNIVAC 1108   (D.P.)  1.73D-18  4.30D+09  2.33D-10  8.98D+307
    ! VAX 11/780    (D.P.)  1.39D-17  1.07D+09  9.31D-10  1.70D+38
    !
    !-------------------------------------------------------------------
    !
    ! Error Returns
    !
    ! The program returns the value zero for X .GT. XMAX, and returns
    ! -XINF when BESLY0 is called with a negative or zero argument.
    !-------------------------------------------------------------------

    real(wp), INTENT(IN)  :: arg
    real(wp), INTENT(OUT) :: result
    INTEGER, INTENT(IN)   :: jint

    INTEGER  :: i
    real(wp) :: ax, down, prod, resj, r0, r1, up, w, wsq, xden, xnum, xy, z, zsq
    !-------------------------------------------------------------------
    ! Mathematical constants
    ! CONS = ln(.5) + Euler's gamma
    !-------------------------------------------------------------------
    real(wp), PARAMETER :: zero = 0.0_wp, one = 1.0_wp, three = 3.0_wp, &
                           four = 4.0_wp, eight = 8.0_wp, five5 = 5.5_wp, &
                           sixty4 = 64.0_wp, oneov8 = 0.125_wp, &
                           p17 = 0.1716_wp, two56 = 256.0_wp, &
                           cons = -1.1593151565841244881e-1_wp, &
                           pi2 = 6.3661977236758134308e-1_wp, &
                           twopi = 6.2831853071795864769_wp, &
                           twopi1 = 6.28125_wp, twopi2 = 1.9353071795864769253e-3_wp
    !-------------------------------------------------------------------
    ! Machine-dependent constants
    ! Default values:
    ! XMAX = 2.68e+8_wp, XSMALL = 3.72e-9_wp, XINF = 1.79e+308_wp
    !-------------------------------------------------------------------
    real(wp), PARAMETER :: XMAX = 1.0e+32_wp, XSMALL = 3.72e-9_wp, XINF = huge(0.0_wp)
    !-------------------------------------------------------------------
    ! Zeroes of Bessel functions
    !-------------------------------------------------------------------
    real(wp), PARAMETER :: XJ0 = 2.4048255576957727686_wp, &
                           XJ1 = 5.5200781102863106496_wp, &
                           XY0 = 8.9357696627916752158e-1_wp, &
                           XY1 = 3.9576784193148578684_wp, &
                           XY2 = 7.0860510603017726976_wp, &
                           XJ01 =  616.0_wp, XJ02 = -1.4244423042272313784e-3_wp, &
                           XJ11 = 1413.0_wp, XJ12 =  5.4686028631064959660e-4_wp, &
                           XY01 =  228.0_wp, XY02 =  2.9519662791675215849e-3_wp, &
                           XY11 = 1013.0_wp, XY12 =  6.4716931485786837568e-4_wp, &
                           XY21 = 1814.0_wp, XY22 =  1.1356030177269762362e-4_wp
    !-------------------------------------------------------------------
    ! Coefficients for rational approximation to ln(x/a)
    !--------------------------------------------------------------------
    real(wp), PARAMETER :: PLG(4) = [ &
      -2.4562334077563243311e+1_wp, 2.3642701335621505212e+2_wp, &
      -5.4989956895857911039e+2_wp, 3.5687548468071500413e+2_wp &
    ]
    real(wp), PARAMETER :: QLG(4) = [ &
      -3.5553900764052419184e+1_wp, 1.9400230218539473193e+2_wp, &
      -3.3442903192607538956e+2_wp, 1.7843774234035750207e+2_wp &
    ]
    !-------------------------------------------------------------------
    ! Coefficients for rational approximation of
    ! J0(X) / (X**2 - XJ0**2), XSMALL < |X| <= 4.0
    !--------------------------------------------------------------------
    real(wp), PARAMETER :: PJ0(7) = [ &
      6.6302997904833794242e+6_wp, -6.2140700423540120665e+8_wp, &
      2.7282507878605942706e+10_wp, -4.1298668500990866786e+11_wp, &
      -1.2117036164593528341e-1_wp, 1.0344222815443188943e+2_wp, &
      -3.6629814655107086448e+4_wp &
    ]
    real(wp), PARAMETER :: QJ0(5) = [ &
      4.5612696224219938200e+5_wp, 1.3985097372263433271e+8_wp, &
      2.6328198300859648632e+10_wp, 2.3883787996332290397e+12_wp, &
      9.3614022392337710626e+2_wp &
    ]
    !-------------------------------------------------------------------
    ! Coefficients for rational approximation of
    ! J0(X) / (X**2 - XJ1**2), 4.0 < |X| <= 8.0
    !-------------------------------------------------------------------
    real(wp), PARAMETER :: PJ1(8) = [ &
      4.4176707025325087628e+3_wp, 1.1725046279757103576e+4_wp, &
      1.0341910641583726701e+4_wp, -7.2879702464464618998e+3_wp, &
      -1.2254078161378989535e+4_wp, -1.8319397969392084011e+3_wp, &
      4.8591703355916499363e+1_wp, 7.4321196680624245801e+2_wp &
    ]
    real(wp), PARAMETER :: QJ1(7) = [ &
      3.3307310774649071172e+2_wp, -2.9458766545509337327e+3_wp, &
      1.8680990008359188352e+4_wp, -8.4055062591169562211e+4_wp, &
      2.4599102262586308984e+5_wp, -3.5783478026152301072e+5_wp, &
      -2.5258076240801555057e+1_wp &
    ]
    !-------------------------------------------------------------------
    ! Coefficients for rational approximation of
    ! (Y0(X) - 2 LN(X/XY0) J0(X)) / (X**2 - XY0**2), XSMALL < |X| <= 3.0
    !--------------------------------------------------------------------
    real(wp), PARAMETER :: PY0(6) = [ &
      1.0102532948020907590e+4_wp, -2.1287548474401797963e+6_wp, &
      2.0422274357376619816e+8_wp, -8.3716255451260504098e+9_wp, &
      1.0723538782003176831e+11_wp, -1.8402381979244993524e+1_wp &
    ]
    real(wp), PARAMETER :: QY0(5) = [ &
      6.6475986689240190091e+2_wp, 2.3889393209447253406e+5_wp, &
      5.5662956624278251596e+7_wp, 8.1617187777290363573e+9_wp, &
      5.8873865738997033405e+11_wp &
    ]
    !-------------------------------------------------------------------
    ! Coefficients for rational approximation of
    ! (Y0(X) - 2 LN(X/XY1) J0(X)) / (X**2 - XY1**2), 3.0 < |X| <= 5.5
    !--------------------------------------------------------------------
    real(wp), PARAMETER :: PY1(7) = [ &
      -1.4566865832663635920e+4_wp, 4.6905288611678631510e+6_wp, &
      -6.9590439394619619534e+8_wp, 4.3600098638603061642e+10_wp, &
      -5.5107435206722644429e+11_wp, -2.2213976967566192242e+13_wp, &
      1.7427031242901594547e+1_wp &
    ]
    real(wp), PARAMETER :: QY1(6) = [ &
      8.3030857612070288823e+2_wp, 4.0669982352539552018e+5_wp, &
      1.3960202770986831075e+8_wp, 3.4015103849971240096e+10_wp, &
      5.4266824419412347550e+12_wp, 4.3386146580707264428e+14_wp &
    ]
    !-------------------------------------------------------------------
    ! Coefficients for rational approximation of
    ! (Y0(X) - 2 LN(X/XY2) J0(X)) / (X**2 - XY2**2), 5.5 < |X| <= 8.0
    !--------------------------------------------------------------------
    real(wp), PARAMETER :: PY2(8) = [ &
      2.1363534169313901632e+4_wp, -1.0085539923498211426e+7_wp, &
      2.1958827170518100757e+9_wp, -1.9363051266772083678e+11_wp, &
      -1.2829912364088687306e+11_wp, 6.7016641869173237784e+14_wp, &
      -8.0728726905150210443e+15_wp, -1.7439661319197499338e+1_wp &
    ]
    real(wp), PARAMETER :: QY2(7) = [ &
      8.7903362168128450017e+2_wp, 5.3924739209768057030e+5_wp, &
      2.4727219475672302327e+8_wp, 8.6926121104209825246e+10_wp, &
      2.2598377924042897629e+13_wp, 3.9272425569640309819e+15_wp, &
      3.4563724628846457519e+17_wp &
    ]
    !-------------------------------------------------------------------
    ! Coefficients for Hart's approximation, |X| > 8.0
    !-------------------------------------------------------------------
    real(wp), PARAMETER :: P0(6) = [ &
      3.4806486443249270347e+3_wp, 2.1170523380864944322e+4_wp, &
      4.1345386639580765797e+4_wp, 2.2779090197304684302e+4_wp, &
      8.8961548424210455236e-1_wp, 1.5376201909008354296e+2_wp &
    ]
    real(wp), PARAMETER :: Q0(5) = [ &
      3.5028735138235608207e+3_wp, 2.1215350561880115730e+4_wp, &
      4.1370412495510416640e+4_wp, 2.2779090197304684318e+4_wp, &
      1.5711159858080893649e+2_wp &
    ]
    real(wp), PARAMETER :: P1(6) = [ &
      -2.2300261666214198472e+1_wp, -1.1183429920482737611e+2_wp, &
      -1.8591953644342993800e+2_wp, -8.9226600200800094098e+1_wp, &
      -8.8033303048680751817e-3_wp, -1.2441026745835638459_wp &
    ]
    real(wp), PARAMETER :: Q1(5) = [ &
      1.4887231232283756582e+3_wp, 7.2642780169211018836e+3_wp, &
      1.1951131543434613647e+4_wp, 5.7105024128512061905e+3_wp, &
      9.0593769594993125859e+1_wp &
    ]
    !-------------------------------------------------------------------
    ! Check for error conditions
    !-------------------------------------------------------------------
    ax = ABS(arg)
    IF (jint == 1 .AND. arg <= zero) THEN
      result = -xinf
      GO TO 80
    ELSE IF (ax > xmax) THEN
      result = zero
      GO TO 80
    END IF
    IF (ax <= eight) THEN
      IF (ax <= xsmall) THEN
        IF (jint == 0) THEN
          result = one
        ELSE
          result = pi2 * (LOG(ax)+cons)
        END IF
        GO TO 80
      END IF
    !-------------------------------------------------------------------
    ! Calculate J0 for appropriate interval, preserving
    ! accuracy near the zero of J0
    !-------------------------------------------------------------------
      zsq = ax * ax
      IF (ax <= four) THEN
        xnum = (pj0(5)*zsq + pj0(6)) * zsq + pj0(7)
        xden = zsq + qj0(5)
        DO  i = 1, 4
          xnum = xnum * zsq + pj0(i)
          xden = xden * zsq + qj0(i)
        END DO
        prod = ((ax-xj01/two56)-xj02) * (ax+xj0)
      ELSE
        wsq = one - zsq / sixty4
        xnum = pj1(7) * wsq + pj1(8)
        xden = wsq + qj1(7)
        DO  i = 1, 6
          xnum = xnum * wsq + pj1(i)
          xden = xden * wsq + qj1(i)
        END DO
        prod = (ax+xj1) * ((ax-xj11/two56)-xj12)
      END IF
      result = prod * xnum / xden
      IF (jint == 0) GO TO 80
    !-------------------------------------------------------------------
    ! Calculate Y0.  First find  RESJ = pi/2 ln(x/xn) J0(x),
    ! where xn is a zero of Y0
    !-------------------------------------------------------------------
      IF (ax <= three) THEN
        up = (ax-xy01/two56) - xy02
        xy = xy0
      ELSE IF (ax <= five5) THEN
        up = (ax-xy11/two56) - xy12
        xy = xy1
      ELSE
        up = (ax-xy21/two56) - xy22
        xy = xy2
      END IF
      down = ax + xy
      IF (ABS(up) < p17*down) THEN
        w = up / down
        wsq = w * w
        xnum = plg(1)
        xden = wsq + qlg(1)
        DO  i = 2, 4
          xnum = xnum * wsq + plg(i)
          xden = xden * wsq + qlg(i)
        END DO
        resj = pi2 * result * w * xnum / xden
      ELSE
        resj = pi2 * result * LOG(ax/xy)
      END IF
    !-------------------------------------------------------------------
    ! Now calculate Y0 for appropriate interval, preserving
    ! accuracy near the zero of Y0
    !-------------------------------------------------------------------
      IF (ax <= three) THEN
        xnum = py0(6) * zsq + py0(1)
        xden = zsq + qy0(1)
        DO  i = 2, 5
          xnum = xnum * zsq + py0(i)
          xden = xden * zsq + qy0(i)
        END DO
      ELSE IF (ax <= five5) THEN
        xnum = py1(7) * zsq + py1(1)
        xden = zsq + qy1(1)
        DO  i = 2, 6
          xnum = xnum * zsq + py1(i)
          xden = xden * zsq + qy1(i)
        END DO
      ELSE
        xnum = py2(8) * zsq + py2(1)
        xden = zsq + qy2(1)
        DO  i = 2, 7
          xnum = xnum * zsq + py2(i)
          xden = xden * zsq + qy2(i)
        END DO
      END IF
      result = resj + up * down * xnum / xden
    ELSE
    !-------------------------------------------------------------------
    ! Calculate J0 or Y0 for |ARG| > 8.0
    !-------------------------------------------------------------------
      z = eight / ax
      w = ax / twopi
      w = AINT(w) + oneov8
      w = (ax-w*twopi1) - w * twopi2
      zsq = z * z
      xnum = p0(5) * zsq + p0(6)
      xden = zsq + q0(5)
      up = p1(5) * zsq + p1(6)
      down = zsq + q1(5)
      DO  i = 1, 4
        xnum = xnum * zsq + p0(i)
        xden = xden * zsq + q0(i)
        up = up * zsq + p1(i)
        down = down * zsq + q1(i)
      END DO
      r0 = xnum / xden
      r1 = up / down
      IF (jint == 0) THEN
        result = SQRT(pi2/ax) * (r0*COS(w) - z*r1*SIN(w))
      ELSE
        result = SQRT(pi2/ax) * (r0*SIN(w) + z*r1*COS(w))
      END IF
    END IF
  80 RETURN
  END SUBROUTINE caljy0


  SUBROUTINE caljy1(arg, result, jint)
    !* CALGO 715 Bessel functions \(J_1(x)\) and \(Y_1(x)\).
    !
    ! \(x \in \mathbb{R}\) for \(J_1(x)\)
    ! 
    ! \(\lbrace x \in \mathbb{R} \mid x \gt 0 \rbrace\) for \(Y_1(x)\)
    !
    ! To obtain:
    ! 
    ! - `y` = \(J_1(x)\), call `caljy1(x, y, 0)`
    ! - `y` = \(Y_1(x)\), call `caljy1(x, y, 1)`
    !*

    !-------------------------------------------------------------------
    ! This packet computes first-order Bessel functions of the first and
    ! second kind (J1 and Y1), for real arguments X, where 0 < X <= XMAX
    ! for Y1, and |X| <= XMAX for J1. It contains two function-type
    ! subprograms, BESJ1 and BESY1, and one subroutine-type
    ! subprogram, CALJY1. The calling statements for the primary
    ! entries are:
    !
    ! Y = BESJ1(X)
    !
    !   and
    !
    ! Y = BESY1(X),
    !
    ! where the entry points correspond to the functions J1(X) and Y1(X),
    ! respectively. The routine CALJY1 is intended for internal packet
    ! use only, all computations within the packet being concentrated in
    ! this one routine. The function subprograms invoke CALJY1  with
    ! the statement
    !
    ! CALL CALJY1(ARG, RESULT, JINT),
    !
    ! where the parameter usage is as follows:
    !
    ! Function     Parameters for CALJY1
    ! call         ARG                    RESULT    JINT
    !
    ! BESJ1(ARG)   |ARG| .LE. XMAX        J1(ARG)   0
    ! BESY1(ARG)   0 .LT. ARG .LE. XMAX   Y1(ARG)   1
    !
    ! The main computation uses unpublished minimax rational
    ! approximations for X .LE. 8.0, and an approximation from the
    ! book Computer Approximations  by Hart, et. al., Wiley and Sons,
    ! New York, 1968, for arguments larger than 8.0 Part of this
    ! transportable packet is patterned after the machine-dependent
    ! FUNPACK program BESJ1(X), but cannot match that version for
    ! efficiency or accuracy. This version uses rational functions
    ! that are theoretically accurate to at least 18 significant decimal
    ! digits for X <= 8, and at least 18 decimal places for X > 8. The
    ! accuracy achieved depends on the arithmetic system, the compiler,
    ! the intrinsic functions, and proper selection of the machine-
    ! dependent constants.
    !
    !-------------------------------------------------------------------
    !
    ! The following machine-dependent constants must be declared in
    ! DATA statements.  IEEE values are provided as a default.
    !
    ! XINF   = largest positive machine number
    ! XMAX   = largest acceptable argument. The functions AINT, SIN
    !          and COS must perform properly for ABS(X) .LE. XMAX.
    !          We recommend that XMAX be a small integer multiple of
    !          sqrt(1/eps), where eps is the smallest positive number
    !          such that  1+eps > 1.
    ! XSMALL = positive argument such that 1.0-(1/2)(X/2)**2 = 1.0
    !          to machine precision for all ABS(X) .LE. XSMALL.
    !          We recommend that  XSMALL < sqrt(eps)/beta, where beta
    !          is the floating-point radix (usually 2 or 16).
    !
    ! Approximate values for some important machines are
    !
    !                       eps       XMAX      XSMALL    XINF
    !
    ! CDC 7600      (S.P.)  7.11E-15  1.34E+08  2.98E-08  1.26E+322
    ! CRAY-1        (S.P.)  7.11E-15  1.34E+08  2.98E-08  5.45E+2465
    ! IBM PC (8087) (S.P.)  5.96E-08  8.19E+03  1.22E-04  3.40E+38
    ! IBM PC (8087) (D.P.)  1.11D-16  2.68D+08  3.72D-09  1.79D+308
    ! IBM 195       (D.P.)  2.22D-16  6.87D+09  9.09D-13  7.23D+75
    ! UNIVAC 1108   (D.P.)  1.73D-18  4.30D+09  2.33D-10  8.98D+307
    ! VAX 11/780    (D.P.)  1.39D-17  1.07D+09  9.31D-10  1.70D+38
    !
    !--------------------------------------------------------------------
    !
    ! Error Returns
    !
    ! The program returns the value zero for  X .GT. XMAX, and returns
    ! -XINF when BESLY1 is called with a negative or zero argument.
    !--------------------------------------------------------------------

    real(wp), INTENT(IN)  :: arg
    real(wp), INTENT(OUT) :: result
    INTEGER, INTENT(IN)   :: jint

    INTEGER  :: i
    real(wp) :: ax, down, prod, resj, r0, r1, up, w, wsq, xden, xnum, xy, z, zsq
    !-------------------------------------------------------------------
    ! Mathematical constants
    !-------------------------------------------------------------------
    real(wp), PARAMETER :: eight = 8.0_wp, four = 4.0_wp, half = 0.5_wp, &
                           throv8 = 0.375_wp, p17 = 0.1716_wp, &
                           pi2 = 6.3661977236758134308e-1_wp, zero = 0.0_wp, &
                           twopi = 6.2831853071795864769_wp, twopi1 = 6.28125_wp, &
                           twopi2 = 1.9353071795864769253e-3_wp, two56 = 256.0_wp, &
                           rtpi2 = 7.9788456080286535588e-1_wp
    !-------------------------------------------------------------------
    ! Machine-dependent constants
    ! Default values:
    ! XMAX = 2.68e+8_wp, XSMALL = 3.72e-9_wp, XINF = 1.79e+308_wp
    !-------------------------------------------------------------------
    real(wp), PARAMETER :: XMAX = 1.0e+32_wp, XSMALL = 3.72e-9_wp, XINF = huge(0.0_wp)
    !-------------------------------------------------------------------
    ! Zeroes of Bessel functions
    !-------------------------------------------------------------------
    real(wp), PARAMETER :: & 
      XJ0 = 3.8317059702075123156_wp, XJ1 = 7.0155866698156187535_wp, &
      XY0 = 2.1971413260310170351_wp, XY1 = 5.4296810407941351328_wp, &
      XJ01 = 981.0_wp, XJ02 = -3.2527979248768438556e-4_wp, &
      XJ11 = 1796.0_wp, XJ12 = -3.8330184381246462950e-5_wp, &
      XY01 = 562.0_wp, XY02 = 1.8288260310170351490e-3_wp, &
      XY11 = 1390.0_wp, XY12 = -6.4592058648672279948e-6_wp
    !-------------------------------------------------------------------
    ! Coefficients for rational approximation to ln(x/a)
    !--------------------------------------------------------------------
    real(wp), PARAMETER :: PLG(4) = [ &
      -2.4562334077563243311e+1_wp, 2.3642701335621505212e+2_wp, &
      -5.4989956895857911039e+2_wp, 3.5687548468071500413e+2_wp &
    ]
    real(wp), PARAMETER :: QLG(4) = [ &
      -3.5553900764052419184e+1_wp, 1.9400230218539473193e+2_wp, &
      -3.3442903192607538956e+2_wp, 1.7843774234035750207e+2_wp &
    ]
    !-------------------------------------------------------------------
    ! Coefficients for rational approximation of
    ! J1(X) / (X * (X**2 - XJ0**2)), XSMALL < |X| <= 4.0
    !--------------------------------------------------------------------
    real(wp), PARAMETER :: PJ0(7) = [ &
      9.8062904098958257677e+5_wp, -1.1548696764841276794e+8_wp, &
      6.6781041261492395835e+9_wp, -1.4258509801366645672e+11_wp, &
      -4.4615792982775076130e+3_wp, 1.0650724020080236441e+1_wp, &
      -1.0767857011487300348e-2_wp &
    ]
    real(wp), PARAMETER :: QJ0(5) = [ &
      5.9117614494174794095e+5_wp, 2.0228375140097033958e+8_wp, &
      4.2091902282580133541e+10_wp, 4.1868604460820175290e+12_wp, &
      1.0742272239517380498e+3_wp &
    ]
    !-------------------------------------------------------------------
    ! Coefficients for rational approximation of
    ! J1(X) / (X * (X**2 - XJ1**2)), 4.0 < |X| <= 8.0
    !-------------------------------------------------------------------
    real(wp), PARAMETER :: PJ1(8) = [ &
      4.6179191852758252280_wp, -7.1329006872560947377e+3_wp, &
      4.5039658105749078904e+6_wp, -1.4437717718363239107e+9_wp, &
      2.3569285397217157313e+11_wp, -1.6324168293282543629e+13_wp, &
      1.1357022719979468624e+14_wp, 1.0051899717115285432e+15_wp &
    ]
    real(wp), PARAMETER :: QJ1(7) = [ &
      1.1267125065029138050e+6_wp, 6.4872502899596389593e+8_wp, &
      2.7622777286244082666e+11_wp, 8.4899346165481429307e+13_wp, &
      1.7128800897135812012e+16_wp, 1.7253905888447681194e+18_wp, &
      1.3886978985861357615e+3_wp &
    ]
    !-------------------------------------------------------------------
    ! Coefficients for rational approximation of
    ! (Y1(X) - 2 LN(X/XY0) J1(X)) / (X**2 - XY0**2), XSMALL < |X| <= 4.0
    !--------------------------------------------------------------------
    real(wp), PARAMETER :: PY0(7) = [ &
      2.2157953222280260820e+5_wp, -5.9157479997408395984e+7_wp, &
      7.2144548214502560419e+9_wp, -3.7595974497819597599e+11_wp, &
      5.4708611716525426053e+12_wp, 4.0535726612579544093e+13_wp, &
      -3.1714424660046133456e+2_wp &
    ]
    real(wp), PARAMETER :: QY0(6) = [ &
      8.2079908168393867438e+2_wp, 3.8136470753052572164e+5_wp, &
      1.2250435122182963220e+8_wp, 2.7800352738690585613e+10_wp, &
      4.1272286200406461981e+12_wp, 3.0737873921079286084e+14_wp &
    ]
    !--------------------------------------------------------------------
    ! Coefficients for rational approximation of
    ! (Y1(X) - 2 LN(X/XY1) J1(X)) / (X**2 - XY1**2), 4.0 < |X| <= 8.0
    !--------------------------------------------------------------------
    real(wp), PARAMETER :: PY1(9) = [ &
      1.9153806858264202986e+6_wp, -1.1957961912070617006e+9_wp, &
      3.7453673962438488783e+11_wp, -5.9530713129741981618e+13_wp, &
      4.0686275289804744814e+15_wp, -2.3638408497043134724e+16_wp, &
      -5.6808094574724204577e+18_wp, 1.1514276357909013326e+19_wp, &
      -1.2337180442012953128e+3_wp &
    ]
    real(wp), PARAMETER :: QY1(8) = [ &
      1.2855164849321609336e+3_wp, 1.0453748201934079734e+6_wp, &
      6.3550318087088919566e+8_wp, 3.0221766852960403645e+11_wp, &
      1.1187010065856971027e+14_wp, 3.0837179548112881950e+16_wp, &
      5.6968198822857178911e+18_wp, 5.3321844313316185697e+20_wp &
    ]
    !-------------------------------------------------------------------
    ! Coefficients for Hart's approximation, |X| > 8.0
    !-------------------------------------------------------------------
    real(wp), PARAMETER :: P0(6) = [ &
      -1.0982405543459346727e+5_wp, -1.5235293511811373833e+6_wp, &
      -6.6033732483649391093e+6_wp, -9.9422465050776411957e+6_wp, &
      -4.4357578167941278571e+6_wp, -1.6116166443246101165e+3_wp &
    ]
    real(wp), PARAMETER :: Q0(6) = [ &
      -1.0726385991103820119e+5_wp,-1.5118095066341608816e+6_wp, &
      -6.5853394797230870728e+6_wp,-9.9341243899345856590e+6_wp, &
      -4.4357578167941278568e+6_wp,-1.4550094401904961825e+3_wp &
    ]
    real(wp), PARAMETER :: P1(6) = [ &
      1.7063754290207680021e+3_wp, 1.8494262873223866797e+4_wp, &
      6.6178836581270835179e+4_wp, 8.5145160675335701966e+4_wp, &
      3.3220913409857223519e+4_wp, 3.5265133846636032186e+1_wp &
    ]
    real(wp), PARAMETER :: Q1(6) = [ &
      3.7890229745772202641e+4_wp, 4.0029443582266975117e+5_wp, &
      1.4194606696037208929e+6_wp, 1.8194580422439972989e+6_wp, &
      7.0871281941028743574e+5_wp, 8.6383677696049909675e+2_wp &
    ]
    !-------------------------------------------------------------------
    ! Check for error conditions
    !-------------------------------------------------------------------
    ax = ABS(arg)
    IF (jint == 1 .AND. (arg <= zero .OR. (arg < half .AND. ax*xinf < pi2))) THEN
      result = -xinf
      GO TO 80
    ELSE IF (ax > xmax) THEN
      result = zero
      GO TO 80
    END IF
    IF (ax > eight) THEN
      GO TO 60
    ELSE IF (ax <= xsmall) THEN
      IF (jint == 0) THEN
        result = arg * half
      ELSE
        result = -pi2 / ax
      END IF
      GO TO 80
    END IF
    !-------------------------------------------------------------------
    ! Calculate J1 for appropriate interval, preserving accuracy near
    ! the zero of J1
    !-------------------------------------------------------------------
    zsq = ax * ax
    IF (ax <= four) THEN
      xnum = (pj0(7)*zsq+pj0(6)) * zsq + pj0(5)
      xden = zsq + qj0(5)
      DO  i = 1, 4
        xnum = xnum * zsq + pj0(i)
        xden = xden * zsq + qj0(i)
      END DO
      prod = arg * ((ax-xj01/two56)-xj02) * (ax+xj0)
    ELSE
      xnum = pj1(1)
      xden = (zsq+qj1(7)) * zsq + qj1(1)
      DO  i = 2, 6
        xnum = xnum * zsq + pj1(i)
        xden = xden * zsq + qj1(i)
      END DO
      xnum = xnum * (ax-eight) * (ax+eight) + pj1(7)
      xnum = xnum * (ax-four) * (ax+four) + pj1(8)
      prod = arg * ((ax-xj11/two56)-xj12) * (ax+xj1)
    END IF
    result = prod * (xnum/xden)
    IF (jint == 0) GO TO 80
    !-------------------------------------------------------------------
    ! Calculate Y1.  First find  RESJ = pi/2 ln(x/xn) J1(x),
    ! where xn is a zero of Y1
    !-------------------------------------------------------------------
    IF (ax <= four) THEN
      up = (ax-xy01/two56) - xy02
      xy = xy0
    ELSE
      up = (ax-xy11/two56) - xy12
      xy = xy1
    END IF
    down = ax + xy
    IF (ABS(up) < p17*down) THEN
      w = up / down
      wsq = w * w
      xnum = plg(1)
      xden = wsq + qlg(1)
      DO  i = 2, 4
        xnum = xnum * wsq + plg(i)
        xden = xden * wsq + qlg(i)
      END DO
      resj = pi2 * result * w * xnum / xden
    ELSE
      resj = pi2 * result * LOG(ax/xy)
    END IF
    !-------------------------------------------------------------------
    ! Now calculate Y1 for appropriate interval, preserving
    ! accuracy near the zero of Y1
    !-------------------------------------------------------------------
    IF (ax <= four) THEN
      xnum = py0(7) * zsq + py0(1)
      xden = zsq + qy0(1)
      DO  i = 2, 6
        xnum = xnum * zsq + py0(i)
        xden = xden * zsq + qy0(i)
      END DO
    ELSE
      xnum = py1(9) * zsq + py1(1)
      xden = zsq + qy1(1)
      DO  i = 2, 8
        xnum = xnum * zsq + py1(i)
        xden = xden * zsq + qy1(i)
      END DO
    END IF
    result = resj + (up*down/ax) * xnum / xden
    GO TO 80
    !-------------------------------------------------------------------
    ! Calculate J1 or Y1 for |ARG| > 8.0
    !-------------------------------------------------------------------
  60 z = eight / ax
    w = AINT(ax/twopi) + throv8
    w = (ax-w*twopi1) - w * twopi2
    zsq = z * z
    xnum = p0(6)
    xden = zsq + q0(6)
    up = p1(6)
    down = zsq + q1(6)
    DO  i = 1, 5
      xnum = xnum * zsq + p0(i)
      xden = xden * zsq + q0(i)
      up = up * zsq + p1(i)
      down = down * zsq + q1(i)
    END DO
    r0 = xnum / xden
    r1 = up / down
    IF (jint == 0) THEN
      result = (rtpi2/SQRT(ax)) * (r0*COS(w)-z*r1*SIN(w))
    ELSE
      result = (rtpi2/SQRT(ax)) * (r0*SIN(w)+z*r1*COS(w))
    END IF
    IF ((jint == 0) .AND. (arg < zero)) result = -result
  80 RETURN
  END SUBROUTINE caljy1

end module calgo_715
