crati.f
111 lines
| 3.2 KiB
| text/x-fortran
|
FortranFixedLexer
| r1568 | *DECK CRATI | |||
| SUBROUTINE CRATI (Z, FNU, N, CY, TOL) | ||||
| C***BEGIN PROLOGUE CRATI | ||||
| C***SUBSIDIARY | ||||
| C***PURPOSE Subsidiary to CBESH, CBESI and CBESK | ||||
| C***LIBRARY SLATEC | ||||
| C***TYPE ALL (CRATI-A, ZRATI-A) | ||||
| C***AUTHOR Amos, D. E., (SNL) | ||||
| C***DESCRIPTION | ||||
| C | ||||
| C CRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD | ||||
| C RECURRENCE. THE STARTING INDEX IS DETERMINED BY FORWARD | ||||
| C RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B, | ||||
| C MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973, | ||||
| C BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER, | ||||
| C BY D. J. SOOKNE. | ||||
| C | ||||
| C***SEE ALSO CBESH, CBESI, CBESK | ||||
| C***ROUTINES CALLED (NONE) | ||||
| C***REVISION HISTORY (YYMMDD) | ||||
| C 830501 DATE WRITTEN | ||||
| C 910415 Prologue converted to Version 4.0 format. (BAB) | ||||
| C***END PROLOGUE CRATI | ||||
| COMPLEX CDFNU, CONE, CY, CZERO, PT, P1, P2, RZ, T1, Z | ||||
| REAL AK, AMAGZ, AP1, AP2, ARG, AZ, DFNU, FDNU, FLAM, FNU, FNUP, | ||||
| * RAP1, RHO, TEST, TEST1, TOL | ||||
| INTEGER I, ID, IDNU, INU, ITIME, K, KK, MAGZ, N | ||||
| DIMENSION CY(N) | ||||
| DATA CZERO, CONE / (0.0E0,0.0E0), (1.0E0,0.0E0) / | ||||
| C***FIRST EXECUTABLE STATEMENT CRATI | ||||
| AZ = ABS(Z) | ||||
| INU = FNU | ||||
| IDNU = INU + N - 1 | ||||
| FDNU = IDNU | ||||
| MAGZ = AZ | ||||
| AMAGZ = MAGZ+1 | ||||
| FNUP = MAX(AMAGZ,FDNU) | ||||
| ID = IDNU - MAGZ - 1 | ||||
| ITIME = 1 | ||||
| K = 1 | ||||
| RZ = (CONE+CONE)/Z | ||||
| T1 = CMPLX(FNUP,0.0E0)*RZ | ||||
| P2 = -T1 | ||||
| P1 = CONE | ||||
| T1 = T1 + RZ | ||||
| IF (ID.GT.0) ID = 0 | ||||
| AP2 = ABS(P2) | ||||
| AP1 = ABS(P1) | ||||
| C----------------------------------------------------------------------- | ||||
| C THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNX | ||||
| C GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT | ||||
| C P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR | ||||
| C PREMATURELY. | ||||
| C----------------------------------------------------------------------- | ||||
| ARG = (AP2+AP2)/(AP1*TOL) | ||||
| TEST1 = SQRT(ARG) | ||||
| TEST = TEST1 | ||||
| RAP1 = 1.0E0/AP1 | ||||
| P1 = P1*CMPLX(RAP1,0.0E0) | ||||
| P2 = P2*CMPLX(RAP1,0.0E0) | ||||
| AP2 = AP2*RAP1 | ||||
| 10 CONTINUE | ||||
| K = K + 1 | ||||
| AP1 = AP2 | ||||
| PT = P2 | ||||
| P2 = P1 - T1*P2 | ||||
| P1 = PT | ||||
| T1 = T1 + RZ | ||||
| AP2 = ABS(P2) | ||||
| IF (AP1.LE.TEST) GO TO 10 | ||||
| IF (ITIME.EQ.2) GO TO 20 | ||||
| AK = ABS(T1)*0.5E0 | ||||
| FLAM = AK + SQRT(AK*AK-1.0E0) | ||||
| RHO = MIN(AP2/AP1,FLAM) | ||||
| TEST = TEST1*SQRT(RHO/(RHO*RHO-1.0E0)) | ||||
| ITIME = 2 | ||||
| GO TO 10 | ||||
| 20 CONTINUE | ||||
| KK = K + 1 - ID | ||||
| AK = KK | ||||
| DFNU = FNU + (N-1) | ||||
| CDFNU = CMPLX(DFNU,0.0E0) | ||||
| T1 = CMPLX(AK,0.0E0) | ||||
| P1 = CMPLX(1.0E0/AP2,0.0E0) | ||||
| P2 = CZERO | ||||
| DO 30 I=1,KK | ||||
| PT = P1 | ||||
| P1 = RZ*(CDFNU+T1)*P1 + P2 | ||||
| P2 = PT | ||||
| T1 = T1 - CONE | ||||
| 30 CONTINUE | ||||
| IF (REAL(P1).NE.0.0E0 .OR. AIMAG(P1).NE.0.0E0) GO TO 40 | ||||
| P1 = CMPLX(TOL,TOL) | ||||
| 40 CONTINUE | ||||
| CY(N) = P2/P1 | ||||
| IF (N.EQ.1) RETURN | ||||
| K = N - 1 | ||||
| AK = K | ||||
| T1 = CMPLX(AK,0.0E0) | ||||
| CDFNU = CMPLX(FNU,0.0E0)*RZ | ||||
| DO 60 I=2,N | ||||
| PT = CDFNU + T1*RZ + CY(K+1) | ||||
| IF (REAL(PT).NE.0.0E0 .OR. AIMAG(PT).NE.0.0E0) GO TO 50 | ||||
| PT = CMPLX(TOL,TOL) | ||||
| 50 CONTINUE | ||||
| CY(K) = CONE/PT | ||||
| T1 = T1 - CONE | ||||
| K = K - 1 | ||||
| 60 CONTINUE | ||||
| RETURN | ||||
| END | ||||
