Mercurial > octave
changeset 3114:17579a02f0b3
[project @ 1997-11-28 20:04:59 by jwe]
| author | jwe |
|---|---|
| date | Fri, 28 Nov 1997 20:05:56 +0000 |
| parents | 61bb314b2c3d |
| children | 7e925ec34aeb |
| files | libcruft/ChangeLog libcruft/Makefile.in libcruft/specfun/ribesl.f libcruft/specfun/rjbesl.f libcruft/specfun/rkbesl.f libcruft/specfun/rybesl.f |
| diffstat | 6 files changed, 1811 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/ChangeLog Fri Nov 28 05:32:41 1997 +0000 +++ b/libcruft/ChangeLog Fri Nov 28 20:05:56 1997 +0000 @@ -1,3 +1,8 @@ +Fri Nov 28 14:05:12 1997 John W. Eaton <jwe@bevo.che.wisc.edu> + + * specfun: New subdirectory. + * Makefile.in (CRUFT_DIRS): Add specfun. + Wed Nov 26 01:49:47 1997 John W. Eaton <jwe@bevo.che.wisc.edu> * slatec-fn/d9gmit.f, slatec-fn/d9lgic.f, slatec-fn/d9lgit.f,
--- a/libcruft/Makefile.in Fri Nov 28 05:32:41 1997 +0000 +++ b/libcruft/Makefile.in Fri Nov 28 20:05:56 1997 +0000 @@ -25,7 +25,7 @@ # configure.in and run autoconf). CRUFT_DIRS = balgen blas dassl eispack fftpack lapack linpack \ - minpack misc odepack quadpack ranlib slatec-fn \ + minpack misc odepack quadpack ranlib slatec-fn specfun \ villad SUBDIRS = $(CRUFT_DIRS)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/specfun/ribesl.f Fri Nov 28 20:05:56 1997 +0000 @@ -0,0 +1,428 @@ + SUBROUTINE RIBESL(X,ALPHA,NB,IZE,B,NCALC) +C------------------------------------------------------------------- +C +C This routine calculates Bessel functions I SUB(N+ALPHA) (X) +C for non-negative argument X, and non-negative order N+ALPHA, +C with or without exponential scaling. +C +C +C Explanation of variables in the calling sequence +C +C X - Working precision non-negative real argument for which +C I's or exponentially scaled I's (I*EXP(-X)) +C are to be calculated. If I's are to be calculated, +C X must be less than EXPARG (see below). +C ALPHA - Working precision fractional part of order for which +C I's or exponentially scaled I's (I*EXP(-X)) are +C to be calculated. 0 .LE. ALPHA .LT. 1.0. +C NB - Integer number of functions to be calculated, NB .GT. 0. +C The first function calculated is of order ALPHA, and the +C last is of order (NB - 1 + ALPHA). +C IZE - Integer type. IZE = 1 if unscaled I's are to calculated, +C and 2 if exponentially scaled I's are to be calculated. +C B - Working precision output vector of length NB. If the routine +C terminates normally (NCALC=NB), the vector B contains the +C functions I(ALPHA,X) through I(NB-1+ALPHA,X), or the +C corresponding exponentially scaled functions. +C NCALC - Integer output variable indicating possible errors. +C Before using the vector B, the user should check that +C NCALC=NB, i.e., all orders have been calculated to +C the desired accuracy. See error returns below. +C +C +C******************************************************************* +C******************************************************************* +C +C Explanation of machine-dependent constants +C +C beta = Radix for the floating-point system +C minexp = Smallest representable power of beta +C maxexp = Smallest power of beta that overflows +C it = Number of bits in the mantissa of a working precision +C variable +C NSIG = Decimal significance desired. Should be set to +C INT(LOG10(2)*it+1). Setting NSIG lower will result +C in decreased accuracy while setting NSIG higher will +C increase CPU time without increasing accuracy. The +C truncation error is limited to a relative error of +C T=.5*10**(-NSIG). +C ENTEN = 10.0 ** K, where K is the largest integer such that +C ENTEN is machine-representable in working precision +C ENSIG = 10.0 ** NSIG +C RTNSIG = 10.0 ** (-K) for the smallest integer K such that +C K .GE. NSIG/4 +C ENMTEN = Smallest ABS(X) such that X/4 does not underflow +C XLARGE = Upper limit on the magnitude of X when IZE=2. Bear +C in mind that if ABS(X)=N, then at least N iterations +C of the backward recursion will be executed. The value +C of 10.0 ** 4 is used on every machine. +C EXPARG = Largest working precision argument that the library +C EXP routine can handle and upper limit on the +C magnitude of X when IZE=1; approximately +C LOG(beta**maxexp) +C +C +C Approximate values for some important machines are: +C +C beta minexp maxexp it +C +C CRAY-1 (S.P.) 2 -8193 8191 48 +C Cyber 180/855 +C under NOS (S.P.) 2 -975 1070 48 +C IEEE (IBM/XT, +C SUN, etc.) (S.P.) 2 -126 128 24 +C IEEE (IBM/XT, +C SUN, etc.) (D.P.) 2 -1022 1024 53 +C IBM 3033 (D.P.) 16 -65 63 14 +C VAX (S.P.) 2 -128 127 24 +C VAX D-Format (D.P.) 2 -128 127 56 +C VAX G-Format (D.P.) 2 -1024 1023 53 +C +C +C NSIG ENTEN ENSIG RTNSIG +C +C CRAY-1 (S.P.) 15 1.0E+2465 1.0E+15 1.0E-4 +C Cyber 180/855 +C under NOS (S.P.) 15 1.0E+322 1.0E+15 1.0E-4 +C IEEE (IBM/XT, +C SUN, etc.) (S.P.) 8 1.0E+38 1.0E+8 1.0E-2 +C IEEE (IBM/XT, +C SUN, etc.) (D.P.) 16 1.0D+308 1.0D+16 1.0D-4 +C IBM 3033 (D.P.) 5 1.0D+75 1.0D+5 1.0D-2 +C VAX (S.P.) 8 1.0E+38 1.0E+8 1.0E-2 +C VAX D-Format (D.P.) 17 1.0D+38 1.0D+17 1.0D-5 +C VAX G-Format (D.P.) 16 1.0D+307 1.0D+16 1.0D-4 +C +C +C ENMTEN XLARGE EXPARG +C +C CRAY-1 (S.P.) 1.84E-2466 1.0E+4 5677 +C Cyber 180/855 +C under NOS (S.P.) 1.25E-293 1.0E+4 741 +C IEEE (IBM/XT, +C SUN, etc.) (S.P.) 4.70E-38 1.0E+4 88 +C IEEE (IBM/XT, +C SUN, etc.) (D.P.) 8.90D-308 1.0D+4 709 +C IBM 3033 (D.P.) 2.16D-78 1.0D+4 174 +C VAX (S.P.) 1.17E-38 1.0E+4 88 +C VAX D-Format (D.P.) 1.17D-38 1.0D+4 88 +C VAX G-Format (D.P.) 2.22D-308 1.0D+4 709 +C +C******************************************************************* +C******************************************************************* +C +C Error returns +C +C In case of an error, NCALC .NE. NB, and not all I's are +C calculated to the desired accuracy. +C +C NCALC .LT. 0: An argument is out of range. For example, +C NB .LE. 0, IZE is not 1 or 2, or IZE=1 and ABS(X) .GE. EXPARG. +C In this case, the B-vector is not calculated, and NCALC is +C set to MIN0(NB,0)-1 so that NCALC .NE. NB. +C +C NB .GT. NCALC .GT. 0: Not all requested function values could +C be calculated accurately. This usually occurs because NB is +C much larger than ABS(X). In this case, B(N) is calculated +C to the desired accuracy for N .LE. NCALC, but precision +C is lost for NCALC .LT. N .LE. NB. If B(N) does not vanish +C for N .GT. NCALC (because it is too small to be represented), +C and B(N)/B(NCALC) = 10**(-K), then only the first NSIG-K +C significant figures of B(N) can be trusted. +C +C +C Intrinsic functions required are: +C +C DBLE, EXP, DGAMMA, GAMMA, INT, MAX, MIN, REAL, SQRT +C +C +C Acknowledgement +C +C This program is based on a program written by David J. +C Sookne (2) that computes values of the Bessel functions J or +C I of real argument and integer order. Modifications include +C the restriction of the computation to the I Bessel function +C of non-negative real argument, the extension of the computation +C to arbitrary positive order, the inclusion of optional +C exponential scaling, and the elimination of most underflow. +C An earlier version was published in (3). +C +C References: "A Note on Backward Recurrence Algorithms," Olver, +C F. W. J., and Sookne, D. J., Math. Comp. 26, 1972, +C pp 941-947. +C +C "Bessel Functions of Real Argument and Integer Order," +C Sookne, D. J., NBS Jour. of Res. B. 77B, 1973, pp +C 125-132. +C +C "ALGORITHM 597, Sequence of Modified Bessel Functions +C of the First Kind," Cody, W. J., Trans. Math. Soft., +C 1983, pp. 242-245. +C +C Latest modification: May 30, 1989 +C +C Modified by: W. J. Cody and L. Stoltz +C Applied Mathematics Division +C Argonne National Laboratory +C Argonne, IL 60439 +C +C------------------------------------------------------------------- + INTEGER IZE,K,L,MAGX,N,NB,NBMX,NCALC,NEND,NSIG,NSTART + DOUBLE PRECISION DGAMMA, + 1 ALPHA,B,CONST,CONV,EM,EMPAL,EMP2AL,EN,ENMTEN,ENSIG, + 2 ENTEN,EXPARG,FUNC,HALF,HALFX,ONE,P,PLAST,POLD,PSAVE,PSAVEL, + 3 RTNSIG,SUM,TEMPA,TEMPB,TEMPC,TEST,TOVER,TWO,X,XLARGE,ZERO + DIMENSION B(NB) +C------------------------------------------------------------------- +C Mathematical constants +C------------------------------------------------------------------- + DATA ONE,TWO,ZERO,HALF,CONST/1.0D0,2.0D0,0.0D0,0.5D0,1.585D0/ +C------------------------------------------------------------------- +C Machine-dependent parameters +C------------------------------------------------------------------- + DATA NSIG,XLARGE,EXPARG /16,1.0D4,709.0D0/ + DATA ENTEN,ENSIG,RTNSIG/1.0D308,1.0D16,1.0D-4/ + DATA ENMTEN/8.9D-308/ +C------------------------------------------------------------------- +C Statement functions for conversion +C------------------------------------------------------------------- + CONV(N) = DBLE(N) + FUNC(X) = DGAMMA(X) +C------------------------------------------------------------------- +C Check for X, NB, OR IZE out of range. +C------------------------------------------------------------------- + IF ((NB.GT.0) .AND. (X .GE. ZERO) .AND. + 1 (ALPHA .GE. ZERO) .AND. (ALPHA .LT. ONE) .AND. + 2 (((IZE .EQ. 1) .AND. (X .LE. EXPARG)) .OR. + 3 ((IZE .EQ. 2) .AND. (X .LE. XLARGE)))) THEN +C------------------------------------------------------------------- +C Use 2-term ascending series for small X +C------------------------------------------------------------------- + NCALC = NB + MAGX = INT(X) + IF (X .GE. RTNSIG) THEN +C------------------------------------------------------------------- +C Initialize the forward sweep, the P-sequence of Olver +C------------------------------------------------------------------- + NBMX = NB-MAGX + N = MAGX+1 + EN = CONV(N+N) + (ALPHA+ALPHA) + PLAST = ONE + P = EN / X +C------------------------------------------------------------------- +C Calculate general significance test +C------------------------------------------------------------------- + TEST = ENSIG + ENSIG + IF (2*MAGX .GT. 5*NSIG) THEN + TEST = SQRT(TEST*P) + ELSE + TEST = TEST / CONST**MAGX + END IF + IF (NBMX .GE. 3) THEN +C------------------------------------------------------------------- +C Calculate P-sequence until N = NB-1. Check for possible overflow. +C------------------------------------------------------------------- + TOVER = ENTEN / ENSIG + NSTART = MAGX+2 + NEND = NB - 1 + DO 100 K = NSTART, NEND + N = K + EN = EN + TWO + POLD = PLAST + PLAST = P + P = EN * PLAST/X + POLD + IF (P .GT. TOVER) THEN +C------------------------------------------------------------------- +C To avoid overflow, divide P-sequence by TOVER. Calculate +C P-sequence until ABS(P) .GT. 1. +C------------------------------------------------------------------- + TOVER = ENTEN + P = P / TOVER + PLAST = PLAST / TOVER + PSAVE = P + PSAVEL = PLAST + NSTART = N + 1 + 60 N = N + 1 + EN = EN + TWO + POLD = PLAST + PLAST = P + P = EN * PLAST/X + POLD + IF (P .LE. ONE) GO TO 60 + TEMPB = EN / X +C------------------------------------------------------------------- +C Calculate backward test, and find NCALC, the highest N +C such that the test is passed. +C------------------------------------------------------------------- + TEST = POLD*PLAST / ENSIG + TEST = TEST*(HALF-HALF/(TEMPB*TEMPB)) + P = PLAST * TOVER + N = N - 1 + EN = EN - TWO + NEND = MIN0(NB,N) + DO 80 L = NSTART, NEND + NCALC = L + POLD = PSAVEL + PSAVEL = PSAVE + PSAVE = EN * PSAVEL/X + POLD + IF (PSAVE*PSAVEL .GT. TEST) GO TO 90 + 80 CONTINUE + NCALC = NEND + 1 + 90 NCALC = NCALC - 1 + GO TO 120 + END IF + 100 CONTINUE + N = NEND + EN = CONV(N+N) + (ALPHA+ALPHA) +C------------------------------------------------------------------- +C Calculate special significance test for NBMX .GT. 2. +C------------------------------------------------------------------- + TEST = MAX(TEST,SQRT(PLAST*ENSIG)*SQRT(P+P)) + END IF +C------------------------------------------------------------------- +C Calculate P-sequence until significance test passed. +C------------------------------------------------------------------- + 110 N = N + 1 + EN = EN + TWO + POLD = PLAST + PLAST = P + P = EN * PLAST/X + POLD + IF (P .LT. TEST) GO TO 110 +C------------------------------------------------------------------- +C Initialize the backward recursion and the normalization sum. +C------------------------------------------------------------------- + 120 N = N + 1 + EN = EN + TWO + TEMPB = ZERO + TEMPA = ONE / P + EM = CONV(N) - ONE + EMPAL = EM + ALPHA + EMP2AL = (EM - ONE) + (ALPHA + ALPHA) + SUM = TEMPA * EMPAL * EMP2AL / EM + NEND = N - NB + IF (NEND .LT. 0) THEN +C------------------------------------------------------------------- +C N .LT. NB, so store B(N) and set higher orders to zero. +C------------------------------------------------------------------- + B(N) = TEMPA + NEND = -NEND + DO 130 L = 1, NEND + 130 B(N+L) = ZERO + ELSE + IF (NEND .GT. 0) THEN +C------------------------------------------------------------------- +C Recur backward via difference equation, calculating (but +C not storing) B(N), until N = NB. +C------------------------------------------------------------------- + DO 140 L = 1, NEND + N = N - 1 + EN = EN - TWO + TEMPC = TEMPB + TEMPB = TEMPA + TEMPA = (EN*TEMPB) / X + TEMPC + EM = EM - ONE + EMP2AL = EMP2AL - ONE + IF (N .EQ. 1) GO TO 150 + IF (N .EQ. 2) EMP2AL = ONE + EMPAL = EMPAL - ONE + SUM = (SUM + TEMPA*EMPAL) * EMP2AL / EM + 140 CONTINUE + END IF +C------------------------------------------------------------------- +C Store B(NB) +C------------------------------------------------------------------- + 150 B(N) = TEMPA + IF (NB .LE. 1) THEN + SUM = (SUM + SUM) + TEMPA + GO TO 230 + END IF +C------------------------------------------------------------------- +C Calculate and Store B(NB-1) +C------------------------------------------------------------------- + N = N - 1 + EN = EN - TWO + B(N) = (EN*TEMPA) / X + TEMPB + IF (N .EQ. 1) GO TO 220 + EM = EM - ONE + EMP2AL = EMP2AL - ONE + IF (N .EQ. 2) EMP2AL = ONE + EMPAL = EMPAL - ONE + SUM = (SUM + B(N)*EMPAL) * EMP2AL / EM + END IF + NEND = N - 2 + IF (NEND .GT. 0) THEN +C------------------------------------------------------------------- +C Calculate via difference equation and store B(N), until N = 2. +C------------------------------------------------------------------- + DO 200 L = 1, NEND + N = N - 1 + EN = EN - TWO + B(N) = (EN*B(N+1)) / X +B(N+2) + EM = EM - ONE + EMP2AL = EMP2AL - ONE + IF (N .EQ. 2) EMP2AL = ONE + EMPAL = EMPAL - ONE + SUM = (SUM + B(N)*EMPAL) * EMP2AL / EM + 200 CONTINUE + END IF +C------------------------------------------------------------------- +C Calculate B(1) +C------------------------------------------------------------------- + B(1) = TWO*EMPAL*B(2) / X + B(3) + 220 SUM = (SUM + SUM) + B(1) +C------------------------------------------------------------------- +C Normalize. Divide all B(N) by sum. +C------------------------------------------------------------------- + 230 IF (ALPHA .NE. ZERO) + 1 SUM = SUM * FUNC(ONE+ALPHA) * (X*HALF)**(-ALPHA) + IF (IZE .EQ. 1) SUM = SUM * EXP(-X) + TEMPA = ENMTEN + IF (SUM .GT. ONE) TEMPA = TEMPA * SUM + DO 260 N = 1, NB + IF (B(N) .LT. TEMPA) B(N) = ZERO + B(N) = B(N) / SUM + 260 CONTINUE + RETURN +C------------------------------------------------------------------- +C Two-term ascending series for small X. +C------------------------------------------------------------------- + ELSE + TEMPA = ONE + EMPAL = ONE + ALPHA + HALFX = ZERO + IF (X .GT. ENMTEN) HALFX = HALF * X + IF (ALPHA .NE. ZERO) TEMPA = HALFX**ALPHA /FUNC(EMPAL) + IF (IZE .EQ. 2) TEMPA = TEMPA * EXP(-X) + TEMPB = ZERO + IF ((X+ONE) .GT. ONE) TEMPB = HALFX * HALFX + B(1) = TEMPA + TEMPA*TEMPB / EMPAL + IF ((X .NE. ZERO) .AND. (B(1) .EQ. ZERO)) NCALC = 0 + IF (NB .GT. 1) THEN + IF (X .EQ. ZERO) THEN + DO 310 N = 2, NB + B(N) = ZERO + 310 CONTINUE + ELSE +C------------------------------------------------------------------- +C Calculate higher-order functions. +C------------------------------------------------------------------- + TEMPC = HALFX + TOVER = (ENMTEN + ENMTEN) / X + IF (TEMPB .NE. ZERO) TOVER = ENMTEN / TEMPB + DO 340 N = 2, NB + TEMPA = TEMPA / EMPAL + EMPAL = EMPAL + ONE + TEMPA = TEMPA * TEMPC + IF (TEMPA .LE. TOVER*EMPAL) TEMPA = ZERO + B(N) = TEMPA + TEMPA*TEMPB / EMPAL + IF ((B(N) .EQ. ZERO) .AND. (NCALC .GT. N)) + 1 NCALC = N-1 + 340 CONTINUE + END IF + END IF + END IF + ELSE + NCALC = MIN0(NB,0)-1 + END IF + RETURN +C---------- Last line of RIBESL ---------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/specfun/rjbesl.f Fri Nov 28 20:05:56 1997 +0000 @@ -0,0 +1,490 @@ + SUBROUTINE RJBESL(X, ALPHA, NB, B, NCALC) +C--------------------------------------------------------------------- +C This routine calculates Bessel functions J sub(N+ALPHA) (X) +C for non-negative argument X, and non-negative order N+ALPHA. +C +C +C Explanation of variables in the calling sequence. +C +C X - working precision non-negative real argument for which +C J's are to be calculated. +C ALPHA - working precision fractional part of order for which +C J's or exponentially scaled J'r (J*exp(X)) are +C to be calculated. 0 <= ALPHA < 1.0. +C NB - integer number of functions to be calculated, NB > 0. +C The first function calculated is of order ALPHA, and the +C last is of order (NB - 1 + ALPHA). +C B - working precision output vector of length NB. If RJBESL +C terminates normally (NCALC=NB), the vector B contains the +C functions J/ALPHA/(X) through J/NB-1+ALPHA/(X), or the +C corresponding exponentially scaled functions. +C NCALC - integer output variable indicating possible errors. +C Before using the vector B, the user should check that +C NCALC=NB, i.e., all orders have been calculated to +C the desired accuracy. See Error Returns below. +C +C +C******************************************************************* +C******************************************************************* +C +C Explanation of machine-dependent constants +C +C it = Number of bits in the mantissa of a working precision +C variable +C NSIG = Decimal significance desired. Should be set to +C INT(LOG10(2)*it+1). Setting NSIG lower will result +C in decreased accuracy while setting NSIG higher will +C increase CPU time without increasing accuracy. The +C truncation error is limited to a relative error of +C T=.5*10**(-NSIG). +C ENTEN = 10.0 ** K, where K is the largest integer such that +C ENTEN is machine-representable in working precision +C ENSIG = 10.0 ** NSIG +C RTNSIG = 10.0 ** (-K) for the smallest integer K such that +C K .GE. NSIG/4 +C ENMTEN = Smallest ABS(X) such that X/4 does not underflow +C XLARGE = Upper limit on the magnitude of X. If ABS(X)=N, +C then at least N iterations of the backward recursion +C will be executed. The value of 10.0 ** 4 is used on +C every machine. +C +C +C Approximate values for some important machines are: +C +C +C it NSIG ENTEN ENSIG +C +C CRAY-1 (S.P.) 48 15 1.0E+2465 1.0E+15 +C Cyber 180/855 +C under NOS (S.P.) 48 15 1.0E+322 1.0E+15 +C IEEE (IBM/XT, +C SUN, etc.) (S.P.) 24 8 1.0E+38 1.0E+8 +C IEEE (IBM/XT, +C SUN, etc.) (D.P.) 53 16 1.0D+308 1.0D+16 +C IBM 3033 (D.P.) 14 5 1.0D+75 1.0D+5 +C VAX (S.P.) 24 8 1.0E+38 1.0E+8 +C VAX D-Format (D.P.) 56 17 1.0D+38 1.0D+17 +C VAX G-Format (D.P.) 53 16 1.0D+307 1.0D+16 +C +C +C RTNSIG ENMTEN XLARGE +C +C CRAY-1 (S.P.) 1.0E-4 1.84E-2466 1.0E+4 +C Cyber 180/855 +C under NOS (S.P.) 1.0E-4 1.25E-293 1.0E+4 +C IEEE (IBM/XT, +C SUN, etc.) (S.P.) 1.0E-2 4.70E-38 1.0E+4 +C IEEE (IBM/XT, +C SUN, etc.) (D.P.) 1.0E-4 8.90D-308 1.0D+4 +C IBM 3033 (D.P.) 1.0E-2 2.16D-78 1.0D+4 +C VAX (S.P.) 1.0E-2 1.17E-38 1.0E+4 +C VAX D-Format (D.P.) 1.0E-5 1.17D-38 1.0D+4 +C VAX G-Format (D.P.) 1.0E-4 2.22D-308 1.0D+4 +C +C******************************************************************* +C******************************************************************* +C +C Error returns +C +C In case of an error, NCALC .NE. NB, and not all J's are +C calculated to the desired accuracy. +C +C NCALC .LT. 0: An argument is out of range. For example, +C NBES .LE. 0, ALPHA .LT. 0 or .GT. 1, or X is too large. +C In this case, B(1) is set to zero, the remainder of the +C B-vector is not calculated, and NCALC is set to +C MIN(NB,0)-1 so that NCALC .NE. NB. +C +C NB .GT. NCALC .GT. 0: Not all requested function values could +C be calculated accurately. This usually occurs because NB is +C much larger than ABS(X). In this case, B(N) is calculated +C to the desired accuracy for N .LE. NCALC, but precision +C is lost for NCALC .LT. N .LE. NB. If B(N) does not vanish +C for N .GT. NCALC (because it is too small to be represented), +C and B(N)/B(NCALC) = 10**(-K), then only the first NSIG-K +C significant figures of B(N) can be trusted. +C +C +C Intrinsic and other functions required are: +C +C ABS, AINT, COS, DBLE, GAMMA (or DGAMMA), INT, MAX, MIN, +C +C REAL, SIN, SQRT +C +C +C Acknowledgement +C +C This program is based on a program written by David J. Sookne +C (2) that computes values of the Bessel functions J or I of real +C argument and integer order. Modifications include the restriction +C of the computation to the J Bessel function of non-negative real +C argument, the extension of the computation to arbitrary positive +C order, and the elimination of most underflow. +C +C References: "A Note on Backward Recurrence Algorithms," Olver, +C F. W. J., and Sookne, D. J., Math. Comp. 26, 1972, +C pp 941-947. +C +C "Bessel Functions of Real Argument and Integer Order," +C Sookne, D. J., NBS Jour. of Res. B. 77B, 1973, pp +C 125-132. +C +C Latest modification: March 19, 1990 +C +C Author: W. J. Cody +C Applied Mathematics Division +C Argonne National Laboratory +C Argonne, IL 60439 +C +C--------------------------------------------------------------------- + INTEGER I,J,K,L,M,MAGX,N,NB,NBMX,NCALC,NEND,NSTART + DOUBLE PRECISION DGAMMA, + 1 ALPHA,ALPEM,ALP2EM,B,CAPP,CAPQ,CONV,EIGHTH,EM,EN,ENMTEN,ENSIG, + 2 ENTEN,FACT,FOUR,FUNC,GNU,HALF,HALFX,ONE,ONE30,P,PI2,PLAST, + 3 POLD,PSAVE,PSAVEL,RTNSIG,S,SUM,T,T1,TEMPA,TEMPB,TEMPC,TEST, + 4 THREE,THREE5,TOVER,TWO,TWOFIV,TWOPI1,TWOPI2,X,XC,XIN,XK,XLARGE, + 5 XM,VCOS,VSIN,Z,ZERO + DIMENSION B(NB), FACT(25) +C--------------------------------------------------------------------- +C Mathematical constants +C +C PI2 - 2 / PI +C TWOPI1 - first few significant digits of 2 * PI +C TWOPI2 - (2*PI - TWOPI) to working precision, i.e., +C TWOPI1 + TWOPI2 = 2 * PI to extra precision. +C--------------------------------------------------------------------- + DATA PI2, TWOPI1, TWOPI2 /0.636619772367581343075535D0,6.28125D0, + 1 1.935307179586476925286767D-3/ + DATA ZERO, EIGHTH, HALF, ONE /0.0D0,0.125D0,0.5D0,1.0D0/ + DATA TWO, THREE, FOUR, TWOFIV /2.0D0,3.0D0,4.0D0,25.0D0/ + DATA ONE30, THREE5 /130.0D0,35.0D0/ +C--------------------------------------------------------------------- +C Machine-dependent parameters +C--------------------------------------------------------------------- + DATA ENTEN, ENSIG, RTNSIG /1.0D38,1.0D17,1.0D-4/ + DATA ENMTEN, XLARGE /1.2D-37,1.0D4/ +C--------------------------------------------------------------------- +C Factorial(N) +C--------------------------------------------------------------------- + DATA FACT /1.0D0,1.0D0,2.0D0,6.0D0,24.0D0,1.2D2,7.2D2,5.04D3, + 1 4.032D4,3.6288D5,3.6288D6,3.99168D7,4.790016D8,6.2270208D9, + 2 8.71782912D10,1.307674368D12,2.0922789888D13,3.55687428096D14, + 3 6.402373705728D15,1.21645100408832D17,2.43290200817664D18, + 4 5.109094217170944D19,1.12400072777760768D21, + 5 2.585201673888497664D22,6.2044840173323943936D23/ +C--------------------------------------------------------------------- +C Statement functions for conversion and the gamma function. +C--------------------------------------------------------------------- + CONV(I) = DBLE(I) + FUNC(X) = DGAMMA(X) +C--------------------------------------------------------------------- +C Check for out of range arguments. +C--------------------------------------------------------------------- + MAGX = INT(X) + IF ((NB.GT.0) .AND. (X.GE.ZERO) .AND. (X.LE.XLARGE) + 1 .AND. (ALPHA.GE.ZERO) .AND. (ALPHA.LT.ONE)) + 2 THEN +C--------------------------------------------------------------------- +C Initialize result array to zero. +C--------------------------------------------------------------------- + NCALC = NB + DO 20 I=1,NB + B(I) = ZERO + 20 CONTINUE +C--------------------------------------------------------------------- +C Branch to use 2-term ascending series for small X and asymptotic +C form for large X when NB is not too large. +C--------------------------------------------------------------------- + IF (X.LT.RTNSIG) THEN +C--------------------------------------------------------------------- +C Two-term ascending series for small X. +C--------------------------------------------------------------------- + TEMPA = ONE + ALPEM = ONE + ALPHA + HALFX = ZERO + IF (X.GT.ENMTEN) HALFX = HALF*X + IF (ALPHA.NE.ZERO) + 1 TEMPA = HALFX**ALPHA/(ALPHA*FUNC(ALPHA)) + TEMPB = ZERO + IF ((X+ONE).GT.ONE) TEMPB = -HALFX*HALFX + B(1) = TEMPA + TEMPA*TEMPB/ALPEM + IF ((X.NE.ZERO) .AND. (B(1).EQ.ZERO)) NCALC = 0 + IF (NB .NE. 1) THEN + IF (X .LE. ZERO) THEN + DO 30 N=2,NB + B(N) = ZERO + 30 CONTINUE + ELSE +C--------------------------------------------------------------------- +C Calculate higher order functions. +C--------------------------------------------------------------------- + TEMPC = HALFX + TOVER = (ENMTEN+ENMTEN)/X + IF (TEMPB.NE.ZERO) TOVER = ENMTEN/TEMPB + DO 50 N=2,NB + TEMPA = TEMPA/ALPEM + ALPEM = ALPEM + ONE + TEMPA = TEMPA*TEMPC + IF (TEMPA.LE.TOVER*ALPEM) TEMPA = ZERO + B(N) = TEMPA + TEMPA*TEMPB/ALPEM + IF ((B(N).EQ.ZERO) .AND. (NCALC.GT.N)) + 1 NCALC = N-1 + 50 CONTINUE + END IF + END IF + ELSE IF ((X.GT.TWOFIV) .AND. (NB.LE.MAGX+1)) THEN +C--------------------------------------------------------------------- +C Asymptotic series for X .GT. 21.0. +C--------------------------------------------------------------------- + XC = SQRT(PI2/X) + XIN = (EIGHTH/X)**2 + M = 11 + IF (X.GE.THREE5) M = 8 + IF (X.GE.ONE30) M = 4 + XM = FOUR*CONV(M) +C--------------------------------------------------------------------- +C Argument reduction for SIN and COS routines. +C--------------------------------------------------------------------- + T = AINT(X/(TWOPI1+TWOPI2)+HALF) + Z = ((X-T*TWOPI1)-T*TWOPI2) - (ALPHA+HALF)/PI2 + VSIN = SIN(Z) + VCOS = COS(Z) + GNU = ALPHA + ALPHA + DO 80 I=1,2 + S = ((XM-ONE)-GNU)*((XM-ONE)+GNU)*XIN*HALF + T = (GNU-(XM-THREE))*(GNU+(XM-THREE)) + CAPP = S*T/FACT(2*M+1) + T1 = (GNU-(XM+ONE))*(GNU+(XM+ONE)) + CAPQ = S*T1/FACT(2*M+2) + XK = XM + K = M + M + T1 = T + DO 70 J=2,M + XK = XK - FOUR + S = ((XK-ONE)-GNU)*((XK-ONE)+GNU) + T = (GNU-(XK-THREE))*(GNU+(XK-THREE)) + CAPP = (CAPP+ONE/FACT(K-1))*S*T*XIN + CAPQ = (CAPQ+ONE/FACT(K))*S*T1*XIN + K = K - 2 + T1 = T + 70 CONTINUE + CAPP = CAPP + ONE + CAPQ = (CAPQ+ONE)*(GNU*GNU-ONE)*(EIGHTH/X) + B(I) = XC*(CAPP*VCOS-CAPQ*VSIN) + IF (NB.EQ.1) GO TO 300 + T = VSIN + VSIN = -VCOS + VCOS = T + GNU = GNU + TWO + 80 CONTINUE +C--------------------------------------------------------------------- +C If NB .GT. 2, compute J(X,ORDER+I) I = 2, NB-1 +C--------------------------------------------------------------------- + IF (NB .GT. 2) THEN + GNU = ALPHA + ALPHA + TWO + DO 90 J=3,NB + B(J) = GNU*B(J-1)/X - B(J-2) + GNU = GNU + TWO + 90 CONTINUE + END IF +C--------------------------------------------------------------------- +C Use recurrence to generate results. First initialize the +C calculation of P*S. +C--------------------------------------------------------------------- + ELSE + NBMX = NB - MAGX + N = MAGX + 1 + EN = CONV(N+N) + (ALPHA+ALPHA) + PLAST = ONE + P = EN/X +C--------------------------------------------------------------------- +C Calculate general significance test. +C--------------------------------------------------------------------- + TEST = ENSIG + ENSIG + IF (NBMX .GE. 3) THEN +C--------------------------------------------------------------------- +C Calculate P*S until N = NB-1. Check for possible overflow. +C--------------------------------------------------------------------- + TOVER = ENTEN/ENSIG + NSTART = MAGX + 2 + NEND = NB - 1 + EN = CONV(NSTART+NSTART) - TWO + (ALPHA+ALPHA) + DO 130 K=NSTART,NEND + N = K + EN = EN + TWO + POLD = PLAST + PLAST = P + P = EN*PLAST/X - POLD + IF (P.GT.TOVER) THEN +C--------------------------------------------------------------------- +C To avoid overflow, divide P*S by TOVER. Calculate P*S until +C ABS(P) .GT. 1. +C--------------------------------------------------------------------- + TOVER = ENTEN + P = P/TOVER + PLAST = PLAST/TOVER + PSAVE = P + PSAVEL = PLAST + NSTART = N + 1 + 100 N = N + 1 + EN = EN + TWO + POLD = PLAST + PLAST = P + P = EN*PLAST/X - POLD + IF (P.LE.ONE) GO TO 100 + TEMPB = EN/X +C--------------------------------------------------------------------- +C Calculate backward test and find NCALC, the highest N such that +C the test is passed. +C--------------------------------------------------------------------- + TEST = POLD*PLAST*(HALF-HALF/(TEMPB*TEMPB)) + TEST = TEST/ENSIG + P = PLAST*TOVER + N = N - 1 + EN = EN - TWO + NEND = MIN(NB,N) + DO 110 L=NSTART,NEND + POLD = PSAVEL + PSAVEL = PSAVE + PSAVE = EN*PSAVEL/X - POLD + IF (PSAVE*PSAVEL.GT.TEST) THEN + NCALC = L - 1 + GO TO 190 + END IF + 110 CONTINUE + NCALC = NEND + GO TO 190 + END IF + 130 CONTINUE + N = NEND + EN = CONV(N+N) + (ALPHA+ALPHA) +C--------------------------------------------------------------------- +C Calculate special significance test for NBMX .GT. 2. +C--------------------------------------------------------------------- + TEST = MAX(TEST,SQRT(PLAST*ENSIG)*SQRT(P+P)) + END IF +C--------------------------------------------------------------------- +C Calculate P*S until significance test passes. +C--------------------------------------------------------------------- + 140 N = N + 1 + EN = EN + TWO + POLD = PLAST + PLAST = P + P = EN*PLAST/X - POLD + IF (P.LT.TEST) GO TO 140 +C--------------------------------------------------------------------- +C Initialize the backward recursion and the normalization sum. +C--------------------------------------------------------------------- + 190 N = N + 1 + EN = EN + TWO + TEMPB = ZERO + TEMPA = ONE/P + M = 2*N - 4*(N/2) + SUM = ZERO + EM = CONV(N/2) + ALPEM = (EM-ONE) + ALPHA + ALP2EM = (EM+EM) + ALPHA + IF (M .NE. 0) SUM = TEMPA*ALPEM*ALP2EM/EM + NEND = N - NB + IF (NEND .GT. 0) THEN +C--------------------------------------------------------------------- +C Recur backward via difference equation, calculating (but not +C storing) B(N), until N = NB. +C--------------------------------------------------------------------- + DO 200 L=1,NEND + N = N - 1 + EN = EN - TWO + TEMPC = TEMPB + TEMPB = TEMPA + TEMPA = (EN*TEMPB)/X - TEMPC + M = 2 - M + IF (M .NE. 0) THEN + EM = EM - ONE + ALP2EM = (EM+EM) + ALPHA + IF (N.EQ.1) GO TO 210 + ALPEM = (EM-ONE) + ALPHA + IF (ALPEM.EQ.ZERO) ALPEM = ONE + SUM = (SUM+TEMPA*ALP2EM)*ALPEM/EM + END IF + 200 CONTINUE + END IF +C--------------------------------------------------------------------- +C Store B(NB). +C--------------------------------------------------------------------- + 210 B(N) = TEMPA + IF (NEND .GE. 0) THEN + IF (NB .LE. 1) THEN + ALP2EM = ALPHA + IF ((ALPHA+ONE).EQ.ONE) ALP2EM = ONE + SUM = SUM + B(1)*ALP2EM + GO TO 250 + ELSE +C--------------------------------------------------------------------- +C Calculate and store B(NB-1). +C--------------------------------------------------------------------- + N = N - 1 + EN = EN - TWO + B(N) = (EN*TEMPA)/X - TEMPB + IF (N.EQ.1) GO TO 240 + M = 2 - M + IF (M .NE. 0) THEN + EM = EM - ONE + ALP2EM = (EM+EM) + ALPHA + ALPEM = (EM-ONE) + ALPHA + IF (ALPEM.EQ.ZERO) ALPEM = ONE + SUM = (SUM+B(N)*ALP2EM)*ALPEM/EM + END IF + END IF + END IF + NEND = N - 2 + IF (NEND .NE. 0) THEN +C--------------------------------------------------------------------- +C Calculate via difference equation and store B(N), until N = 2. +C--------------------------------------------------------------------- + DO 230 L=1,NEND + N = N - 1 + EN = EN - TWO + B(N) = (EN*B(N+1))/X - B(N+2) + M = 2 - M + IF (M .NE. 0) THEN + EM = EM - ONE + ALP2EM = (EM+EM) + ALPHA + ALPEM = (EM-ONE) + ALPHA + IF (ALPEM.EQ.ZERO) ALPEM = ONE + SUM = (SUM+B(N)*ALP2EM)*ALPEM/EM + END IF + 230 CONTINUE + END IF +C--------------------------------------------------------------------- +C Calculate B(1). +C--------------------------------------------------------------------- + B(1) = TWO*(ALPHA+ONE)*B(2)/X - B(3) + 240 EM = EM - ONE + ALP2EM = (EM+EM) + ALPHA + IF (ALP2EM.EQ.ZERO) ALP2EM = ONE + SUM = SUM + B(1)*ALP2EM +C--------------------------------------------------------------------- +C Normalize. Divide all B(N) by sum. +C--------------------------------------------------------------------- + 250 IF ((ALPHA+ONE).NE.ONE) + 1 SUM = SUM*FUNC(ALPHA)*(X*HALF)**(-ALPHA) + TEMPA = ENMTEN + IF (SUM.GT.ONE) TEMPA = TEMPA*SUM + DO 260 N=1,NB + IF (ABS(B(N)).LT.TEMPA) B(N) = ZERO + B(N) = B(N)/SUM + 260 CONTINUE + END IF +C--------------------------------------------------------------------- +C Error return -- X, NB, or ALPHA is out of range. +C--------------------------------------------------------------------- + ELSE + B(1) = ZERO + NCALC = MIN(NB,0) - 1 + END IF +C--------------------------------------------------------------------- +C Exit +C--------------------------------------------------------------------- + 300 RETURN +C ---------- Last line of RJBESL ---------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/specfun/rkbesl.f Fri Nov 28 20:05:56 1997 +0000 @@ -0,0 +1,465 @@ + SUBROUTINE RKBESL(X,ALPHA,NB,IZE,BK,NCALC) +C------------------------------------------------------------------- +C +C This FORTRAN 77 routine calculates modified Bessel functions +C of the second kind, K SUB(N+ALPHA) (X), for non-negative +C argument X, and non-negative order N+ALPHA, with or without +C exponential scaling. +C +C Explanation of variables in the calling sequence +C +C Description of output values .. +C +C X - Working precision non-negative real argument for which +C K's or exponentially scaled K's (K*EXP(X)) +C are to be calculated. If K's are to be calculated, +C X must not be greater than XMAX (see below). +C ALPHA - Working precision fractional part of order for which +C K's or exponentially scaled K's (K*EXP(X)) are +C to be calculated. 0 .LE. ALPHA .LT. 1.0. +C NB - Integer number of functions to be calculated, NB .GT. 0. +C The first function calculated is of order ALPHA, and the +C last is of order (NB - 1 + ALPHA). +C IZE - Integer type. IZE = 1 if unscaled K's are to be calculated, +C and 2 if exponentially scaled K's are to be calculated. +C BK - Working precision output vector of length NB. If the +C routine terminates normally (NCALC=NB), the vector BK +C contains the functions K(ALPHA,X), ... , K(NB-1+ALPHA,X), +C or the corresponding exponentially scaled functions. +C If (0 .LT. NCALC .LT. NB), BK(I) contains correct function +C values for I .LE. NCALC, and contains the ratios +C K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array. +C NCALC - Integer output variable indicating possible errors. +C Before using the vector BK, the user should check that +C NCALC=NB, i.e., all orders have been calculated to +C the desired accuracy. See error returns below. +C +C +C******************************************************************* +C******************************************************************* +C +C Explanation of machine-dependent constants +C +C beta = Radix for the floating-point system +C minexp = Smallest representable power of beta +C maxexp = Smallest power of beta that overflows +C EPS = The smallest positive floating-point number such that +C 1.0+EPS .GT. 1.0 +C XMAX = Upper limit on the magnitude of X when IZE=1; Solution +C to equation: +C W(X) * (1-1/8X+9/128X**2) = beta**minexp +C where W(X) = EXP(-X)*SQRT(PI/2X) +C SQXMIN = Square root of beta**minexp +C XINF = Largest positive machine number; approximately +C beta**maxexp +C XMIN = Smallest positive machine number; approximately +C beta**minexp +C +C +C Approximate values for some important machines are: +C +C beta minexp maxexp EPS +C +C CRAY-1 (S.P.) 2 -8193 8191 7.11E-15 +C Cyber 180/185 +C under NOS (S.P.) 2 -975 1070 3.55E-15 +C IEEE (IBM/XT, +C SUN, etc.) (S.P.) 2 -126 128 1.19E-7 +C IEEE (IBM/XT, +C SUN, etc.) (D.P.) 2 -1022 1024 2.22D-16 +C IBM 3033 (D.P.) 16 -65 63 2.22D-16 +C VAX (S.P.) 2 -128 127 5.96E-8 +C VAX D-Format (D.P.) 2 -128 127 1.39D-17 +C VAX G-Format (D.P.) 2 -1024 1023 1.11D-16 +C +C +C SQXMIN XINF XMIN XMAX +C +C CRAY-1 (S.P.) 6.77E-1234 5.45E+2465 4.59E-2467 5674.858 +C Cyber 180/855 +C under NOS (S.P.) 1.77E-147 1.26E+322 3.14E-294 672.788 +C IEEE (IBM/XT, +C SUN, etc.) (S.P.) 1.08E-19 3.40E+38 1.18E-38 85.337 +C IEEE (IBM/XT, +C SUN, etc.) (D.P.) 1.49D-154 1.79D+308 2.23D-308 705.342 +C IBM 3033 (D.P.) 7.35D-40 7.23D+75 5.40D-79 177.852 +C VAX (S.P.) 5.42E-20 1.70E+38 2.94E-39 86.715 +C VAX D-Format (D.P.) 5.42D-20 1.70D+38 2.94D-39 86.715 +C VAX G-Format (D.P.) 7.46D-155 8.98D+307 5.57D-309 706.728 +C +C******************************************************************* +C******************************************************************* +C +C Error returns +C +C In case of an error, NCALC .NE. NB, and not all K's are +C calculated to the desired accuracy. +C +C NCALC .LT. -1: An argument is out of range. For example, +C NB .LE. 0, IZE is not 1 or 2, or IZE=1 and ABS(X) .GE. +C XMAX. In this case, the B-vector is not calculated, +C and NCALC is set to MIN0(NB,0)-2 so that NCALC .NE. NB. +C NCALC = -1: Either K(ALPHA,X) .GE. XINF or +C K(ALPHA+NB-1,X)/K(ALPHA+NB-2,X) .GE. XINF. In this case, +C the B-vector is not calculated. Note that again +C NCALC .NE. NB. +C +C 0 .LT. NCALC .LT. NB: Not all requested function values could +C be calculated accurately. BK(I) contains correct function +C values for I .LE. NCALC, and contains the ratios +C K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array. +C +C +C Intrinsic functions required are: +C +C ABS, AINT, EXP, INT, LOG, MAX, MIN, SINH, SQRT +C +C +C Acknowledgement +C +C This program is based on a program written by J. B. Campbell +C (2) that computes values of the Bessel functions K of real +C argument and real order. Modifications include the addition +C of non-scaled functions, parameterization of machine +C dependencies, and the use of more accurate approximations +C for SINH and SIN. +C +C References: "On Temme's Algorithm for the Modified Bessel +C Functions of the Third Kind," Campbell, J. B., +C TOMS 6(4), Dec. 1980, pp. 581-586. +C +C "A FORTRAN IV Subroutine for the Modified Bessel +C Functions of the Third Kind of Real Order and Real +C Argument," Campbell, J. B., Report NRC/ERB-925, +C National Research Council, Canada. +C +C Latest modification: May 30, 1989 +C +C Modified by: W. J. Cody and L. Stoltz +C Applied Mathematics Division +C Argonne National Laboratory +C Argonne, IL 60439 +C +C------------------------------------------------------------------- + INTEGER I,IEND,ITEMP,IZE,J,K,M,MPLUS1,NB,NCALC + DOUBLE PRECISION + 1 A,ALPHA,BLPHA,BK,BK1,BK2,C,D,DM,D1,D2,D3,ENU,EPS,ESTF,ESTM, + 2 EX,FOUR,F0,F1,F2,HALF,ONE,P,P0,Q,Q0,R,RATIO,S,SQXMIN,T,TINYX, + 3 TWO,TWONU,TWOX,T1,T2,WMINF,X,XINF,XMAX,XMIN,X2BY4,ZERO + DIMENSION BK(1),P(8),Q(7),R(5),S(4),T(6),ESTM(6),ESTF(7) +C--------------------------------------------------------------------- +C Mathematical constants +C A = LOG(2.D0) - Euler's constant +C D = SQRT(2.D0/PI) +C--------------------------------------------------------------------- + DATA HALF,ONE,TWO,ZERO/0.5D0,1.0D0,2.0D0,0.0D0/ + DATA FOUR,TINYX/4.0D0,1.0D-10/ + DATA A/ 0.11593151565841244881D0/,D/0.797884560802865364D0/ +C--------------------------------------------------------------------- +C Machine dependent parameters +C--------------------------------------------------------------------- + DATA EPS/2.22D-16/,SQXMIN/1.49D-154/,XINF/1.79D+308/ + DATA XMIN/2.23D-308/,XMAX/705.342D0/ +C--------------------------------------------------------------------- +C P, Q - Approximation for LOG(GAMMA(1+ALPHA))/ALPHA +C + Euler's constant +C Coefficients converted from hex to decimal and modified +C by W. J. Cody, 2/26/82 +C R, S - Approximation for (1-ALPHA*PI/SIN(ALPHA*PI))/(2.D0*ALPHA) +C T - Approximation for SINH(Y)/Y +C--------------------------------------------------------------------- + DATA P/ 0.805629875690432845D00, 0.204045500205365151D02, + 1 0.157705605106676174D03, 0.536671116469207504D03, + 2 0.900382759291288778D03, 0.730923886650660393D03, + 3 0.229299301509425145D03, 0.822467033424113231D00/ + DATA Q/ 0.294601986247850434D02, 0.277577868510221208D03, + 1 0.120670325591027438D04, 0.276291444159791519D04, + 2 0.344374050506564618D04, 0.221063190113378647D04, + 3 0.572267338359892221D03/ + DATA R/-0.48672575865218401848D+0, 0.13079485869097804016D+2, + 1 -0.10196490580880537526D+3, 0.34765409106507813131D+3, + 2 0.34958981245219347820D-3/ + DATA S/-0.25579105509976461286D+2, 0.21257260432226544008D+3, + 1 -0.61069018684944109624D+3, 0.42269668805777760407D+3/ + DATA T/ 0.16125990452916363814D-9, 0.25051878502858255354D-7, + 1 0.27557319615147964774D-5, 0.19841269840928373686D-3, + 2 0.83333333333334751799D-2, 0.16666666666666666446D+0/ + DATA ESTM/5.20583D1, 5.7607D0, 2.7782D0, 1.44303D1, 1.853004D2, + 1 9.3715D0/ + DATA ESTF/4.18341D1, 7.1075D0, 6.4306D0, 4.25110D1, 1.35633D0, + 1 8.45096D1, 2.0D1/ +C--------------------------------------------------------------------- + EX = X + ENU = ALPHA + NCALC = MIN(NB,0)-2 + IF ((NB .GT. 0) .AND. ((ENU .GE. ZERO) .AND. (ENU .LT. ONE)) + 1 .AND. ((IZE .GE. 1) .AND. (IZE .LE. 2)) .AND. + 2 ((IZE .NE. 1) .OR. (EX .LE. XMAX)) .AND. + 3 (EX .GT. ZERO)) THEN + K = 0 + IF (ENU .LT. SQXMIN) ENU = ZERO + IF (ENU .GT. HALF) THEN + K = 1 + ENU = ENU - ONE + END IF + TWONU = ENU+ENU + IEND = NB+K-1 + C = ENU*ENU + D3 = -C + IF (EX .LE. ONE) THEN +C--------------------------------------------------------------------- +C Calculation of P0 = GAMMA(1+ALPHA) * (2/X)**ALPHA +C Q0 = GAMMA(1-ALPHA) * (X/2)**ALPHA +C--------------------------------------------------------------------- + D1 = ZERO + D2 = P(1) + T1 = ONE + T2 = Q(1) + DO 10 I = 2,7,2 + D1 = C*D1+P(I) + D2 = C*D2+P(I+1) + T1 = C*T1+Q(I) + T2 = C*T2+Q(I+1) + 10 CONTINUE + D1 = ENU*D1 + T1 = ENU*T1 + F1 = LOG(EX) + F0 = A+ENU*(P(8)-ENU*(D1+D2)/(T1+T2))-F1 + Q0 = EXP(-ENU*(A-ENU*(P(8)+ENU*(D1-D2)/(T1-T2))-F1)) + F1 = ENU*F0 + P0 = EXP(F1) +C--------------------------------------------------------------------- +C Calculation of F0 = +C--------------------------------------------------------------------- + D1 = R(5) + T1 = ONE + DO 20 I = 1,4 + D1 = C*D1+R(I) + T1 = C*T1+S(I) + 20 CONTINUE + IF (ABS(F1) .LE. HALF) THEN + F1 = F1*F1 + D2 = ZERO + DO 30 I = 1,6 + D2 = F1*D2+T(I) + 30 CONTINUE + D2 = F0+F0*F1*D2 + ELSE + D2 = SINH(F1)/ENU + END IF + F0 = D2-ENU*D1/(T1*P0) + IF (EX .LE. TINYX) THEN +C-------------------------------------------------------------------- +C X.LE.1.0E-10 +C Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X) +C-------------------------------------------------------------------- + BK(1) = F0+EX*F0 + IF (IZE .EQ. 1) BK(1) = BK(1)-EX*BK(1) + RATIO = P0/F0 + C = EX*XINF + IF (K .NE. 0) THEN +C-------------------------------------------------------------------- +C Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X), +C ALPHA .GE. 1/2 +C-------------------------------------------------------------------- + NCALC = -1 + IF (BK(1) .GE. C/RATIO) GO TO 500 + BK(1) = RATIO*BK(1)/EX + TWONU = TWONU+TWO + RATIO = TWONU + END IF + NCALC = 1 + IF (NB .EQ. 1) GO TO 500 +C-------------------------------------------------------------------- +C Calculate K(ALPHA+L,X)/K(ALPHA+L-1,X), L = 1, 2, ... , NB-1 +C-------------------------------------------------------------------- + NCALC = -1 + DO 80 I = 2,NB + IF (RATIO .GE. C) GO TO 500 + BK(I) = RATIO/EX + TWONU = TWONU+TWO + RATIO = TWONU + 80 CONTINUE + NCALC = 1 + GO TO 420 + ELSE +C-------------------------------------------------------------------- +C 1.0E-10 .LT. X .LE. 1.0 +C-------------------------------------------------------------------- + C = ONE + X2BY4 = EX*EX/FOUR + P0 = HALF*P0 + Q0 = HALF*Q0 + D1 = -ONE + D2 = ZERO + BK1 = ZERO + BK2 = ZERO + F1 = F0 + F2 = P0 + 100 D1 = D1+TWO + D2 = D2+ONE + D3 = D1+D3 + C = X2BY4*C/D2 + F0 = (D2*F0+P0+Q0)/D3 + P0 = P0/(D2-ENU) + Q0 = Q0/(D2+ENU) + T1 = C*F0 + T2 = C*(P0-D2*F0) + BK1 = BK1+T1 + BK2 = BK2+T2 + IF ((ABS(T1/(F1+BK1)) .GT. EPS) .OR. + 1 (ABS(T2/(F2+BK2)) .GT. EPS)) GO TO 100 + BK1 = F1+BK1 + BK2 = TWO*(F2+BK2)/EX + IF (IZE .EQ. 2) THEN + D1 = EXP(EX) + BK1 = BK1*D1 + BK2 = BK2*D1 + END IF + WMINF = ESTF(1)*EX+ESTF(2) + END IF + ELSE IF (EPS*EX .GT. ONE) THEN +C-------------------------------------------------------------------- +C X .GT. ONE/EPS +C-------------------------------------------------------------------- + NCALC = NB + BK1 = ONE / (D*SQRT(EX)) + DO 110 I = 1, NB + BK(I) = BK1 + 110 CONTINUE + GO TO 500 + ELSE +C-------------------------------------------------------------------- +C X .GT. 1.0 +C-------------------------------------------------------------------- + TWOX = EX+EX + BLPHA = ZERO + RATIO = ZERO + IF (EX .LE. FOUR) THEN +C-------------------------------------------------------------------- +C Calculation of K(ALPHA+1,X)/K(ALPHA,X), 1.0 .LE. X .LE. 4.0 +C-------------------------------------------------------------------- + D2 = AINT(ESTM(1)/EX+ESTM(2)) + M = INT(D2) + D1 = D2+D2 + D2 = D2-HALF + D2 = D2*D2 + DO 120 I = 2,M + D1 = D1-TWO + D2 = D2-D1 + RATIO = (D3+D2)/(TWOX+D1-RATIO) + 120 CONTINUE +C-------------------------------------------------------------------- +C Calculation of I(|ALPHA|,X) and I(|ALPHA|+1,X) by backward +C recurrence and K(ALPHA,X) from the wronskian +C-------------------------------------------------------------------- + D2 = AINT(ESTM(3)*EX+ESTM(4)) + M = INT(D2) + C = ABS(ENU) + D3 = C+C + D1 = D3-ONE + F1 = XMIN + F0 = (TWO*(C+D2)/EX+HALF*EX/(C+D2+ONE))*XMIN + DO 130 I = 3,M + D2 = D2-ONE + F2 = (D3+D2+D2)*F0 + BLPHA = (ONE+D1/D2)*(F2+BLPHA) + F2 = F2/EX+F1 + F1 = F0 + F0 = F2 + 130 CONTINUE + F1 = (D3+TWO)*F0/EX+F1 + D1 = ZERO + T1 = ONE + DO 140 I = 1,7 + D1 = C*D1+P(I) + T1 = C*T1+Q(I) + 140 CONTINUE + P0 = EXP(C*(A+C*(P(8)-C*D1/T1)-LOG(EX)))/EX + F2 = (C+HALF-RATIO)*F1/EX + BK1 = P0+(D3*F0-F2+F0+BLPHA)/(F2+F1+F0)*P0 + IF (IZE .EQ. 1) BK1 = BK1*EXP(-EX) + WMINF = ESTF(3)*EX+ESTF(4) + ELSE +C-------------------------------------------------------------------- +C Calculation of K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X), by backward +C recurrence, for X .GT. 4.0 +C-------------------------------------------------------------------- + DM = AINT(ESTM(5)/EX+ESTM(6)) + M = INT(DM) + D2 = DM-HALF + D2 = D2*D2 + D1 = DM+DM + DO 160 I = 2,M + DM = DM-ONE + D1 = D1-TWO + D2 = D2-D1 + RATIO = (D3+D2)/(TWOX+D1-RATIO) + BLPHA = (RATIO+RATIO*BLPHA)/DM + 160 CONTINUE + BK1 = ONE/((D+D*BLPHA)*SQRT(EX)) + IF (IZE .EQ. 1) BK1 = BK1*EXP(-EX) + WMINF = ESTF(5)*(EX-ABS(EX-ESTF(7)))+ESTF(6) + END IF +C-------------------------------------------------------------------- +C Calculation of K(ALPHA+1,X) from K(ALPHA,X) and +C K(ALPHA+1,X)/K(ALPHA,X) +C-------------------------------------------------------------------- + BK2 = BK1+BK1*(ENU+HALF-RATIO)/EX + END IF +C-------------------------------------------------------------------- +C Calculation of 'NCALC', K(ALPHA+I,X), I = 0, 1, ... , NCALC-1, +C K(ALPHA+I,X)/K(ALPHA+I-1,X), I = NCALC, NCALC+1, ... , NB-1 +C-------------------------------------------------------------------- + NCALC = NB + BK(1) = BK1 + IF (IEND .EQ. 0) GO TO 500 + J = 2-K + IF (J .GT. 0) BK(J) = BK2 + IF (IEND .EQ. 1) GO TO 500 + M = MIN(INT(WMINF-ENU),IEND) + DO 190 I = 2,M + T1 = BK1 + BK1 = BK2 + TWONU = TWONU+TWO + IF (EX .LT. ONE) THEN + IF (BK1 .GE. (XINF/TWONU)*EX) GO TO 195 + GO TO 187 + ELSE + IF (BK1/EX .GE. XINF/TWONU) GO TO 195 + END IF + 187 CONTINUE + BK2 = TWONU/EX*BK1+T1 + ITEMP = I + J = J+1 + IF (J .GT. 0) BK(J) = BK2 + 190 CONTINUE + 195 M = ITEMP + IF (M .EQ. IEND) GO TO 500 + RATIO = BK2/BK1 + MPLUS1 = M+1 + NCALC = -1 + DO 410 I = MPLUS1,IEND + TWONU = TWONU+TWO + RATIO = TWONU/EX+ONE/RATIO + J = J+1 + IF (J .GT. 1) THEN + BK(J) = RATIO + ELSE + IF (BK2 .GE. XINF/RATIO) GO TO 500 + BK2 = RATIO*BK2 + END IF + 410 CONTINUE + NCALC = MAX(MPLUS1-K,1) + IF (NCALC .EQ. 1) BK(1) = BK2 + IF (NB .EQ. 1) GO TO 500 + 420 J = NCALC+1 + DO 430 I = J,NB + IF (BK(NCALC) .GE. XINF/BK(I)) GO TO 500 + BK(I) = BK(NCALC)*BK(I) + NCALC = I + 430 CONTINUE + END IF + 500 RETURN +C---------- Last line of RKBESL ---------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/specfun/rybesl.f Fri Nov 28 20:05:56 1997 +0000 @@ -0,0 +1,422 @@ + SUBROUTINE RYBESL(X,ALPHA,NB,BY,NCALC) +C---------------------------------------------------------------------- +C +C This routine calculates Bessel functions Y SUB(N+ALPHA) (X) +C for non-negative argument X, and non-negative order N+ALPHA. +C +C +C Explanation of variables in the calling sequence +C +C X - Working precision non-negative real argument for which +C Y's are to be calculated. +C ALPHA - Working precision fractional part of order for which +C Y's are to be calculated. 0 .LE. ALPHA .LT. 1.0. +C NB - Integer number of functions to be calculated, NB .GT. 0. +C The first function calculated is of order ALPHA, and the +C last is of order (NB - 1 + ALPHA). +C BY - Working precision output vector of length NB. If the +C routine terminates normally (NCALC=NB), the vector BY +C contains the functions Y(ALPHA,X), ... , Y(NB-1+ALPHA,X), +C If (0 .LT. NCALC .LT. NB), BY(I) contains correct function +C values for I .LE. NCALC, and contains the ratios +C Y(ALPHA+I-1,X)/Y(ALPHA+I-2,X) for the rest of the array. +C NCALC - Integer output variable indicating possible errors. +C Before using the vector BY, the user should check that +C NCALC=NB, i.e., all orders have been calculated to +C the desired accuracy. See error returns below. +C +C +C******************************************************************* +C******************************************************************* +C +C Explanation of machine-dependent constants +C +C beta = Radix for the floating-point system +C p = Number of significant base-beta digits in the +C significand of a floating-point number +C minexp = Smallest representable power of beta +C maxexp = Smallest power of beta that overflows +C EPS = beta ** (-p) +C DEL = Machine number below which sin(x)/x = 1; approximately +C SQRT(EPS). +C XMIN = Smallest acceptable argument for RBESY; approximately +C max(2*beta**minexp,2/XINF), rounded up +C XINF = Largest positive machine number; approximately +C beta**maxexp +C THRESH = Lower bound for use of the asymptotic form; approximately +C AINT(-LOG10(EPS/2.0))+1.0 +C XLARGE = Upper bound on X; approximately 1/DEL, because the sine +C and cosine functions have lost about half of their +C precision at that point. +C +C +C Approximate values for some important machines are: +C +C beta p minexp maxexp EPS +C +C CRAY-1 (S.P.) 2 48 -8193 8191 3.55E-15 +C Cyber 180/185 +C under NOS (S.P.) 2 48 -975 1070 3.55E-15 +C IEEE (IBM/XT, +C SUN, etc.) (S.P.) 2 24 -126 128 5.96E-8 +C IEEE (IBM/XT, +C SUN, etc.) (D.P.) 2 53 -1022 1024 1.11D-16 +C IBM 3033 (D.P.) 16 14 -65 63 1.39D-17 +C VAX (S.P.) 2 24 -128 127 5.96E-8 +C VAX D-Format (D.P.) 2 56 -128 127 1.39D-17 +C VAX G-Format (D.P.) 2 53 -1024 1023 1.11D-16 +C +C +C DEL XMIN XINF THRESH XLARGE +C +C CRAY-1 (S.P.) 5.0E-8 3.67E-2466 5.45E+2465 15.0E0 2.0E7 +C Cyber 180/855 +C under NOS (S.P.) 5.0E-8 6.28E-294 1.26E+322 15.0E0 2.0E7 +C IEEE (IBM/XT, +C SUN, etc.) (S.P.) 1.0E-4 2.36E-38 3.40E+38 8.0E0 1.0E4 +C IEEE (IBM/XT, +C SUN, etc.) (D.P.) 1.0D-8 4.46D-308 1.79D+308 16.0D0 1.0D8 +C IBM 3033 (D.P.) 1.0D-8 2.77D-76 7.23D+75 17.0D0 1.0D8 +C VAX (S.P.) 1.0E-4 1.18E-38 1.70E+38 8.0E0 1.0E4 +C VAX D-Format (D.P.) 1.0D-9 1.18D-38 1.70D+38 17.0D0 1.0D9 +C VAX G-Format (D.P.) 1.0D-8 2.23D-308 8.98D+307 16.0D0 1.0D8 +C +C******************************************************************* +C******************************************************************* +C +C Error returns +C +C In case of an error, NCALC .NE. NB, and not all Y's are +C calculated to the desired accuracy. +C +C NCALC .LT. -1: An argument is out of range. For example, +C NB .LE. 0, IZE is not 1 or 2, or IZE=1 and ABS(X) .GE. +C XMAX. In this case, BY(1) = 0.0, the remainder of the +C BY-vector is not calculated, and NCALC is set to +C MIN0(NB,0)-2 so that NCALC .NE. NB. +C NCALC = -1: Y(ALPHA,X) .GE. XINF. The requested function +C values are set to 0.0. +C 1 .LT. NCALC .LT. NB: Not all requested function values could +C be calculated accurately. BY(I) contains correct function +C values for I .LE. NCALC, and and the remaining NB-NCALC +C array elements contain 0.0. +C +C +C Intrinsic functions required are: +C +C DBLE, EXP, INT, MAX, MIN, REAL, SQRT +C +C +C Acknowledgement +C +C This program draws heavily on Temme's Algol program for Y(a,x) +C and Y(a+1,x) and on Campbell's programs for Y_nu(x). Temme's +C scheme is used for x < THRESH, and Campbell's scheme is used +C in the asymptotic region. Segments of code from both sources +C have been translated into Fortran 77, merged, and heavily modified. +C Modifications include parameterization of machine dependencies, +C use of a new approximation for ln(gamma(x)), and built-in +C protection against over/underflow. +C +C References: "Bessel functions J_nu(x) and Y_nu(x) of real +C order and real argument," Campbell, J. B., +C Comp. Phy. Comm. 18, 1979, pp. 133-142. +C +C "On the numerical evaluation of the ordinary +C Bessel function of the second kind," Temme, +C N. M., J. Comput. Phys. 21, 1976, pp. 343-350. +C +C Latest modification: March 19, 1990 +C +C Modified by: W. J. Cody +C Applied Mathematics Division +C Argonne National Laboratory +C Argonne, IL 60439 +C +C---------------------------------------------------------------------- + INTEGER I,K,NA,NB,NCALC + DOUBLE PRECISION + 1 ALFA,ALPHA,AYE,B,BY,C,CH,COSMU,D,DEL,DEN,DDIV,DIV,DMU,D1,D2, + 2 E,EIGHT,EN,ENU,EN1,EPS,EVEN,EX,F,FIVPI,G,GAMMA,H,HALF,ODD, + 3 ONBPI,ONE,ONE5,P,PA,PA1,PI,PIBY2,PIM5,Q,QA,QA1,Q0,R,S,SINMU, + 4 SQ2BPI,TEN9,TERM,THREE,THRESH,TWO,TWOBYX,X,XINF,XLARGE,XMIN, + 5 XNA,X2,YA,YA1,ZERO + DIMENSION BY(NB),CH(21) +C---------------------------------------------------------------------- +C Mathematical constants +C FIVPI = 5*PI +C PIM5 = 5*PI - 15 +C ONBPI = 1/PI +C PIBY2 = PI/2 +C SQ2BPI = SQUARE ROOT OF 2/PI +C---------------------------------------------------------------------- + DATA ZERO,HALF,ONE,TWO,THREE/0.0D0,0.5D0,1.0D0,2.0D0,3.0D0/ + DATA EIGHT,ONE5,TEN9/8.0D0,15.0D0,1.9D1/ + DATA FIVPI,PIBY2/1.5707963267948966192D1,1.5707963267948966192D0/ + DATA PI,SQ2BPI/3.1415926535897932385D0,7.9788456080286535588D-1/ + DATA PIM5,ONBPI/7.0796326794896619231D-1,3.1830988618379067154D-1/ +C---------------------------------------------------------------------- +C Machine-dependent constants +C---------------------------------------------------------------------- + DATA DEL,XMIN,XINF,EPS/1.0D-8,4.46D-308,1.79D308,1.11D-16/ + DATA THRESH,XLARGE/16.0D0,1.0D8/ +C---------------------------------------------------------------------- +C Coefficients for Chebyshev polynomial expansion of +C 1/gamma(1-x), abs(x) .le. .5 +C---------------------------------------------------------------------- + DATA CH/-0.67735241822398840964D-23,-0.61455180116049879894D-22, + 1 0.29017595056104745456D-20, 0.13639417919073099464D-18, + 2 0.23826220476859635824D-17,-0.90642907957550702534D-17, + 3 -0.14943667065169001769D-14,-0.33919078305362211264D-13, + 4 -0.17023776642512729175D-12, 0.91609750938768647911D-11, + 5 0.24230957900482704055D-09, 0.17451364971382984243D-08, + 6 -0.33126119768180852711D-07,-0.86592079961391259661D-06, + 7 -0.49717367041957398581D-05, 0.76309597585908126618D-04, + 8 0.12719271366545622927D-02, 0.17063050710955562222D-02, + 9 -0.76852840844786673690D-01,-0.28387654227602353814D+00, + A 0.92187029365045265648D+00/ +C---------------------------------------------------------------------- + EX = X + ENU = ALPHA + IF ((NB .GT. 0) .AND. (X .GE. XMIN) .AND. (EX .LT. XLARGE) + 1 .AND. (ENU .GE. ZERO) .AND. (ENU .LT. ONE)) THEN + XNA = AINT(ENU+HALF) + NA = INT(XNA) + IF (NA .EQ. 1) ENU = ENU - XNA + IF (ENU .EQ. -HALF) THEN + P = SQ2BPI/SQRT(EX) + YA = P * SIN(EX) + YA1 = -P * COS(EX) + ELSE IF (EX .LT. THREE) THEN +C---------------------------------------------------------------------- +C Use Temme's scheme for small X +C---------------------------------------------------------------------- + B = EX * HALF + D = -LOG(B) + F = ENU * D + E = B**(-ENU) + IF (ABS(ENU) .LT. DEL) THEN + C = ONBPI + ELSE + C = ENU / SIN(ENU*PI) + END IF +C---------------------------------------------------------------------- +C Computation of sinh(f)/f +C---------------------------------------------------------------------- + IF (ABS(F) .LT. ONE) THEN + X2 = F*F + EN = TEN9 + S = ONE + DO 80 I = 1, 9 + S = S*X2/EN/(EN-ONE)+ONE + EN = EN - TWO + 80 CONTINUE + ELSE + S = (E - ONE/E) * HALF / F + END IF +C---------------------------------------------------------------------- +C Computation of 1/gamma(1-a) using Chebyshev polynomials +C---------------------------------------------------------------------- + X2 = ENU*ENU*EIGHT + AYE = CH(1) + EVEN = ZERO + ALFA = CH(2) + ODD = ZERO + DO 40 I = 3, 19, 2 + EVEN = -(AYE+AYE+EVEN) + AYE = -EVEN*X2 - AYE + CH(I) + ODD = -(ALFA+ALFA+ODD) + ALFA = -ODD*X2 - ALFA + CH(I+1) + 40 CONTINUE + EVEN = (EVEN*HALF+AYE)*X2 - AYE + CH(21) + ODD = (ODD+ALFA)*TWO + GAMMA = ODD*ENU + EVEN +C---------------------------------------------------------------------- +C End of computation of 1/gamma(1-a) +C---------------------------------------------------------------------- + G = E * GAMMA + E = (E + ONE/E) * HALF + F = TWO*C*(ODD*E+EVEN*S*D) + E = ENU*ENU + P = G*C + Q = ONBPI / G + C = ENU*PIBY2 + IF (ABS(C) .LT. DEL) THEN + R = ONE + ELSE + R = SIN(C)/C + END IF + R = PI*C*R*R + C = ONE + D = - B*B + H = ZERO + YA = F + R*Q + YA1 = P + EN = ZERO + 100 EN = EN + ONE + IF (ABS(G/(ONE+ABS(YA))) + 1 + ABS(H/(ONE+ABS(YA1))) .GT. EPS) THEN + F = (F*EN+P+Q)/(EN*EN-E) + C = C * D/EN + P = P/(EN-ENU) + Q = Q/(EN+ENU) + G = C*(F+R*Q) + H = C*P - EN*G + YA = YA + G + YA1 = YA1+H + GO TO 100 + END IF + YA = -YA + YA1 = -YA1/B + ELSE IF (EX .LT. THRESH) THEN +C---------------------------------------------------------------------- +C Use Temme's scheme for moderate X +C---------------------------------------------------------------------- + C = (HALF-ENU)*(HALF+ENU) + B = EX + EX + E = (EX*ONBPI*COS(ENU*PI)/EPS) + E = E*E + P = ONE + Q = -EX + R = ONE + EX*EX + S = R + EN = TWO + 200 IF (R*EN*EN .LT. E) THEN + EN1 = EN+ONE + D = (EN-ONE+C/EN)/S + P = (EN+EN-P*D)/EN1 + Q = (-B+Q*D)/EN1 + S = P*P + Q*Q + R = R*S + EN = EN1 + GO TO 200 + END IF + F = P/S + P = F + G = -Q/S + Q = G + 220 EN = EN - ONE + IF (EN .GT. ZERO) THEN + R = EN1*(TWO-P)-TWO + S = B + EN1*Q + D = (EN-ONE+C/EN)/(R*R+S*S) + P = D*R + Q = D*S + E = F + ONE + F = P*E - G*Q + G = Q*E + P*G + EN1 = EN + GO TO 220 + END IF + F = ONE + F + D = F*F + G*G + PA = F/D + QA = -G/D + D = ENU + HALF -P + Q = Q + EX + PA1 = (PA*Q-QA*D)/EX + QA1 = (QA*Q+PA*D)/EX + B = EX - PIBY2*(ENU+HALF) + C = COS(B) + S = SIN(B) + D = SQ2BPI/SQRT(EX) + YA = D*(PA*S+QA*C) + YA1 = D*(QA1*S-PA1*C) + ELSE +C---------------------------------------------------------------------- +C Use Campbell's asymptotic scheme. +C---------------------------------------------------------------------- + NA = 0 + D1 = AINT(EX/FIVPI) + I = INT(D1) + DMU = ((EX-ONE5*D1)-D1*PIM5)-(ALPHA+HALF)*PIBY2 + IF (I-2*(I/2) .EQ. 0) THEN + COSMU = COS(DMU) + SINMU = SIN(DMU) + ELSE + COSMU = -COS(DMU) + SINMU = -SIN(DMU) + END IF + DDIV = EIGHT * EX + DMU = ALPHA + DEN = SQRT(EX) + DO 350 K = 1, 2 + P = COSMU + COSMU = SINMU + SINMU = -P + D1 = (TWO*DMU-ONE)*(TWO*DMU+ONE) + D2 = ZERO + DIV = DDIV + P = ZERO + Q = ZERO + Q0 = D1/DIV + TERM = Q0 + DO 310 I = 2, 20 + D2 = D2 + EIGHT + D1 = D1 - D2 + DIV = DIV + DDIV + TERM = -TERM*D1/DIV + P = P + TERM + D2 = D2 + EIGHT + D1 = D1 - D2 + DIV = DIV + DDIV + TERM = TERM*D1/DIV + Q = Q + TERM + IF (ABS(TERM) .LE. EPS) GO TO 320 + 310 CONTINUE + 320 P = P + ONE + Q = Q + Q0 + IF (K .EQ. 1) THEN + YA = SQ2BPI * (P*COSMU-Q*SINMU) / DEN + ELSE + YA1 = SQ2BPI * (P*COSMU-Q*SINMU) / DEN + END IF + DMU = DMU + ONE + 350 CONTINUE + END IF + IF (NA .EQ. 1) THEN + H = TWO*(ENU+ONE)/EX + IF (H .GT. ONE) THEN + IF (ABS(YA1) .GT. XINF/H) THEN + H = ZERO + YA = ZERO + END IF + END IF + H = H*YA1 - YA + YA = YA1 + YA1 = H + END IF +C---------------------------------------------------------------------- +C Now have first one or two Y's +C---------------------------------------------------------------------- + BY(1) = YA + BY(2) = YA1 + IF (YA1 .EQ. ZERO) THEN + NCALC = 1 + ELSE + AYE = ONE + ALPHA + TWOBYX = TWO/EX + NCALC = 2 + DO 400 I = 3, NB + IF (TWOBYX .LT. ONE) THEN + IF (ABS(BY(I-1))*TWOBYX .GE. XINF/AYE) + 1 GO TO 450 + ELSE + IF (ABS(BY(I-1)) .GE. XINF/AYE/TWOBYX ) + 1 GO TO 450 + END IF + BY(I) = TWOBYX*AYE*BY(I-1) - BY(I-2) + AYE = AYE + ONE + NCALC = NCALC + 1 + 400 CONTINUE + END IF + 450 DO 460 I = NCALC+1, NB + BY(I) = ZERO + 460 CONTINUE + ELSE + BY(1) = ZERO + NCALC = MIN(NB,0) - 1 + END IF + 900 RETURN +C---------- Last line of RYBESL ---------- + END