RKBESL Subroutine

public subroutine RKBESL(X, ALPHA, NB, IZE, BK, NCALC)

Explanation of machine-dependent constants

beta = Radix for the floating-point system minexp = Smallest representable power of beta maxexp = Smallest power of beta that overflows EPS = The smallest positive floating-point number such that 1.0+EPS .GT. 1.0 XMAX = Upper limit on the magnitude of X when IZE=1; Solution to equation: W(X) * (1-1/8X+9/128X2) = betaminexp where W(X) = EXP(-X)SQRT(PI/2X) SQXMIN = Square root of betaminexp XINF = Largest positive machine number; approximately betamaxexp XMIN = Smallest positive machine number; approximately beta*minexp

 Approximate values for some important machines are:

                      beta       minexp      maxexp      EPS

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

                     SQXMIN       XINF        XMIN      XMAX

CRAY-1 (S.P.) 6.77E-1234 5.45E+2465 4.59E-2467 5674.858 Cyber 180/855 under NOS (S.P.) 1.77E-147 1.26E+322 3.14E-294 672.788 IEEE (IBM/XT, SUN, etc.) (S.P.) 1.08E-19 3.40E+38 1.18E-38 85.337 IEEE (IBM/XT, SUN, etc.) (D.P.) 1.49D-154 1.79D+308 2.23D-308 705.342 IBM 3033 (D.P.) 7.35D-40 7.23D+75 5.40D-79 177.852 VAX (S.P.) 5.42E-20 1.70E+38 2.94E-39 86.715 VAX D-Format (D.P.) 5.42D-20 1.70D+38 2.94D-39 86.715 VAX G-Format (D.P.) 7.46D-155 8.98D+307 5.57D-309 706.728

Error returns

In case of an error, NCALC .NE. NB, and not all K'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, the B-vector is not calculated, and NCALC is set to MIN0(NB,0)-2 so that NCALC .NE. NB. NCALC = -1: Either K(ALPHA,X) .GE. XINF or K(ALPHA+NB-1,X)/K(ALPHA+NB-2,X) .GE. XINF. In this case, the B-vector is not calculated. Note that again NCALC .NE. NB.

0 .LT. NCALC .LT. NB: Not all requested function values could be calculated accurately. BK(I) contains correct function values for I .LE. NCALC, and contains the ratios K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array.

Intrinsic functions required are:

 ABS, AINT, EXP, INT, LOG, MAX, MIN, SINH, SQRT

Acknowledgement

This program is based on a program written by J. B. Campbell (2) that computes values of the Bessel functions K of real argument and real order. Modifications include the addition of non-scaled functions, parameterization of machine dependencies, and the use of more accurate approximations for SINH and SIN.

References: "On Temme's Algorithm for the Modified Bessel Functions of the Third Kind," Campbell, J. B., TOMS 6(4), Dec. 1980, pp. 581-586.

         "A FORTRAN IV Subroutine for the Modified Bessel
          Functions of the Third Kind of Real Order and Real
          Argument," Campbell, J. B., Report NRC/ERB-925,
          National Research Council, Canada.

Latest modification: May 30, 1989

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

Arguments