GNU Octave: libcruft/amos/crati.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CRATI(Z, FNU, N, CY, TOL)
00002 C***BEGIN PROLOGUE  CRATI
00003 C***REFER TO  CBESI,CBESK,CBESH
00004 C
00005 C     CRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD
00006 C     RECURRENCE.  THE STARTING INDEX IS DETERMINED BY FORWARD
00007 C     RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B,
00008 C     MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973,
00009 C     BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER,
00010 C     BY D. J. SOOKNE.
00011 C
00012 C***ROUTINES CALLED  (NONE)
00013 C***END PROLOGUE  CRATI
00014       COMPLEX CDFNU, CONE, CY, CZERO, PT, P1, P2, RZ, T1, Z
00015       REAL AK, AMAGZ, AP1, AP2, ARG, AZ, DFNU, FDNU, FLAM, FNU, FNUP,
00016      * RAP1, RHO, TEST, TEST1, TOL
00017       INTEGER I, ID, IDNU, INU, ITIME, K, KK, MAGZ, N
00018       DIMENSION CY(N)
00019       DATA CZERO, CONE / (0.0E0,0.0E0), (1.0E0,0.0E0) /
00020       AZ = CABS(Z)
00021       INU = INT(FNU)
00022       IDNU = INU + N - 1
00023       FDNU = FLOAT(IDNU)
00024       MAGZ = INT(AZ)
00025       AMAGZ = FLOAT(MAGZ+1)
00026       FNUP = AMAX1(AMAGZ,FDNU)
00027       ID = IDNU - MAGZ - 1
00028       ITIME = 1
00029       K = 1
00030       RZ = (CONE+CONE)/Z
00031       T1 = CMPLX(FNUP,0.0E0)*RZ
00032       P2 = -T1
00033       P1 = CONE
00034       T1 = T1 + RZ
00035       IF (ID.GT.0) ID = 0
00036       AP2 = CABS(P2)
00037       AP1 = CABS(P1)
00038 C-----------------------------------------------------------------------
00039 C     THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNX
00040 C     GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT
00041 C     P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR
00042 C     PREMATURELY.
00043 C-----------------------------------------------------------------------
00044       ARG = (AP2+AP2)/(AP1*TOL)
00045       TEST1 = SQRT(ARG)
00046       TEST = TEST1
00047       RAP1 = 1.0E0/AP1
00048       P1 = P1*CMPLX(RAP1,0.0E0)
00049       P2 = P2*CMPLX(RAP1,0.0E0)
00050       AP2 = AP2*RAP1
00051    10 CONTINUE
00052       K = K + 1
00053       AP1 = AP2
00054       PT = P2
00055       P2 = P1 - T1*P2
00056       P1 = PT
00057       T1 = T1 + RZ
00058       AP2 = CABS(P2)
00059       IF (AP1.LE.TEST) GO TO 10
00060       IF (ITIME.EQ.2) GO TO 20
00061       AK = CABS(T1)*0.5E0
00062       FLAM = AK + SQRT(AK*AK-1.0E0)
00063       RHO = AMIN1(AP2/AP1,FLAM)
00064       TEST = TEST1*SQRT(RHO/(RHO*RHO-1.0E0))
00065       ITIME = 2
00066       GO TO 10
00067    20 CONTINUE
00068       KK = K + 1 - ID
00069       AK = FLOAT(KK)
00070       DFNU = FNU + FLOAT(N-1)
00071       CDFNU = CMPLX(DFNU,0.0E0)
00072       T1 = CMPLX(AK,0.0E0)
00073       P1 = CMPLX(1.0E0/AP2,0.0E0)
00074       P2 = CZERO
00075       DO 30 I=1,KK
00076         PT = P1
00077         P1 = RZ*(CDFNU+T1)*P1 + P2
00078         P2 = PT
00079         T1 = T1 - CONE
00080    30 CONTINUE
00081       IF (REAL(P1).NE.0.0E0 .OR. AIMAG(P1).NE.0.0E0) GO TO 40
00082       P1 = CMPLX(TOL,TOL)
00083    40 CONTINUE
00084       CY(N) = P2/P1
00085       IF (N.EQ.1) RETURN
00086       K = N - 1
00087       AK = FLOAT(K)
00088       T1 = CMPLX(AK,0.0E0)
00089       CDFNU = CMPLX(FNU,0.0E0)*RZ
00090       DO 60 I=2,N
00091         PT = CDFNU + T1*RZ + CY(K+1)
00092         IF (REAL(PT).NE.0.0E0 .OR. AIMAG(PT).NE.0.0E0) GO TO 50
00093         PT = CMPLX(TOL,TOL)
00094    50   CONTINUE
00095         CY(K) = CONE/PT
00096         T1 = T1 - CONE
00097         K = K - 1
00098    60 CONTINUE
00099       RETURN
00100       END