RYBESL Subroutine

public subroutine RYBESL(X, ALPHA, NB, BY, NCALC)

Explanation of machine-dependent constants

beta = Radix for the floating-point system p = Number of significant base-beta digits in the significand of a floating-point number minexp = Smallest representable power of beta maxexp = Smallest power of beta that overflows EPS = beta (-p) DEL = Machine number below which sin(x)/x = 1; approximately SQRT(EPS). XMIN = Smallest acceptable argument for RBESY; approximately max(2*betaminexp,2/XINF), rounded up XINF = Largest positive machine number; approximately beta**maxexp THRESH = Lower bound for use of the asymptotic form; approximately AINT(-LOG10(EPS/2.0))+1.0 XLARGE = Upper bound on X; approximately 1/DEL, because the sine and cosine functions have lost about half of their precision at that point.

 Approximate values for some important machines are:

                    beta    p     minexp      maxexp      EPS

CRAY-1 (S.P.) 2 48 -8193 8191 3.55E-15 Cyber 180/185 under NOS (S.P.) 2 48 -975 1070 3.55E-15 IEEE (IBM/XT, SUN, etc.) (S.P.) 2 24 -126 128 5.96E-8 IEEE (IBM/XT, SUN, etc.) (D.P.) 2 53 -1022 1024 1.11D-16 IBM 3033 (D.P.) 16 14 -65 63 1.39D-17 VAX (S.P.) 2 24 -128 127 5.96E-8 VAX D-Format (D.P.) 2 56 -128 127 1.39D-17 VAX G-Format (D.P.) 2 53 -1024 1023 1.11D-16

                     DEL      XMIN      XINF     THRESH  XLARGE

CRAY-1 (S.P.) 5.0E-8 3.67E-2466 5.45E+2465 15.0E0 2.0E7 Cyber 180/855 under NOS (S.P.) 5.0E-8 6.28E-294 1.26E+322 15.0E0 2.0E7 IEEE (IBM/XT, SUN, etc.) (S.P.) 1.0E-4 2.36E-38 3.40E+38 8.0E0 1.0E4 IEEE (IBM/XT, SUN, etc.) (D.P.) 1.0D-8 4.46D-308 1.79D+308 16.0D0 1.0D8 IBM 3033 (D.P.) 1.0D-8 2.77D-76 7.23D+75 17.0D0 1.0D8 VAX (S.P.) 1.0E-4 1.18E-38 1.70E+38 8.0E0 1.0E4 VAX D-Format (D.P.) 1.0D-9 1.18D-38 1.70D+38 17.0D0 1.0D9 VAX G-Format (D.P.) 1.0D-8 2.23D-308 8.98D+307 16.0D0 1.0D8

Error returns

In case of an error, NCALC .NE. NB, and not all Y's are calculated to the desired accuracy.

NCALC .LT. -1: An argument is out of range. For example, NB .LE. 0, IZE is not 1 or 2, or IZE=1 and ABS(X) .GE. XMAX. In this case, BY(1) = 0.0, the remainder of the BY-vector is not calculated, and NCALC is set to MIN0(NB,0)-2 so that NCALC .NE. NB. NCALC = -1: Y(ALPHA,X) .GE. XINF. The requested function values are set to 0.0. 1 .LT. NCALC .LT. NB: Not all requested function values could be calculated accurately. BY(I) contains correct function values for I .LE. NCALC, and and the remaining NB-NCALC array elements contain 0.0.

Intrinsic functions required are:

 DBLE, EXP, INT, MAX, MIN, REAL, SQRT

Acknowledgement

This program draws heavily on Temme's Algol program for Y(a,x) and Y(a+1,x) and on Campbell's programs for Y_nu(x). Temme's scheme is used for x < THRESH, and Campbell's scheme is used in the asymptotic region. Segments of code from both sources have been translated into Fortran 77, merged, and heavily modified. Modifications include parameterization of machine dependencies, use of a new approximation for ln(gamma(x)), and built-in protection against over/underflow.

References: "Bessel functions J_nu(x) and Y_nu(x) of real order and real argument," Campbell, J. B., Comp. Phy. Comm. 18, 1979, pp. 133-142.

         "On the numerical evaluation of the ordinary
          Bessel function of the second kind," Temme,
          N. M., J. Comput. Phys. 21, 1976, pp. 343-350.

Latest modification: March 19, 1990

Modified by: W. J. Cody Applied Mathematics Division Argonne National Laboratory Argonne, IL 60439

Arguments