1*DECK ZRATI
2      SUBROUTINE ZRATI (ZR, ZI, FNU, N, CYR, CYI, TOL)
3C***BEGIN PROLOGUE  ZRATI
4C***SUBSIDIARY
5C***PURPOSE  Subsidiary to ZBESH, ZBESI and ZBESK
6C***LIBRARY   SLATEC
7C***TYPE      ALL (CRATI-A, ZRATI-A)
8C***AUTHOR  Amos, D. E., (SNL)
9C***DESCRIPTION
10C
11C     ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD
12C     RECURRENCE.  THE STARTING INDEX IS DETERMINED BY FORWARD
13C     RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B,
14C     MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973,
15C     BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER,
16C     BY D. J. SOOKNE.
17C
18C***SEE ALSO  ZBESH, ZBESI, ZBESK
19C***ROUTINES CALLED  ZABS, ZDIV
20C***REVISION HISTORY  (YYMMDD)
21C   830501  DATE WRITTEN
22C   910415  Prologue converted to Version 4.0 format.  (BAB)
23C***END PROLOGUE  ZRATI
24      DOUBLE PRECISION AK, AMAGZ, AP1, AP2, ARG, AZ, CDFNUI, CDFNUR,
25     * CONEI, CONER, CYI, CYR, CZEROI, CZEROR, DFNU, FDNU, FLAM, FNU,
26     * FNUP, PTI, PTR, P1I, P1R, P2I, P2R, RAK, RAP1, RHO, RT2, RZI,
27     * RZR, TEST, TEST1, TOL, TTI, TTR, T1I, T1R, ZI, ZR, ZABS
28      INTEGER I, ID, IDNU, INU, ITIME, K, KK, MAGZ, N
29      DIMENSION CYR(N), CYI(N)
30      EXTERNAL ZABS
31      DATA CZEROR,CZEROI,CONER,CONEI,RT2/
32     1 0.0D0, 0.0D0, 1.0D0, 0.0D0, 1.41421356237309505D0 /
33C***FIRST EXECUTABLE STATEMENT  ZRATI
34      AZ = ZABS(ZR,ZI)
35      INU = FNU
36      IDNU = INU + N - 1
37      MAGZ = AZ
38      AMAGZ = MAGZ+1
39      FDNU = IDNU
40      FNUP = MAX(AMAGZ,FDNU)
41      ID = IDNU - MAGZ - 1
42      ITIME = 1
43      K = 1
44      PTR = 1.0D0/AZ
45      RZR = PTR*(ZR+ZR)*PTR
46      RZI = -PTR*(ZI+ZI)*PTR
47      T1R = RZR*FNUP
48      T1I = RZI*FNUP
49      P2R = -T1R
50      P2I = -T1I
51      P1R = CONER
52      P1I = CONEI
53      T1R = T1R + RZR
54      T1I = T1I + RZI
55      IF (ID.GT.0) ID = 0
56      AP2 = ZABS(P2R,P2I)
57      AP1 = ZABS(P1R,P1I)
58C-----------------------------------------------------------------------
59C     THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU
60C     GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT
61C     P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR
62C     PREMATURELY.
63C-----------------------------------------------------------------------
64      ARG = (AP2+AP2)/(AP1*TOL)
65      TEST1 = SQRT(ARG)
66      TEST = TEST1
67      RAP1 = 1.0D0/AP1
68      P1R = P1R*RAP1
69      P1I = P1I*RAP1
70      P2R = P2R*RAP1
71      P2I = P2I*RAP1
72      AP2 = AP2*RAP1
73   10 CONTINUE
74      K = K + 1
75      AP1 = AP2
76      PTR = P2R
77      PTI = P2I
78      P2R = P1R - (T1R*PTR-T1I*PTI)
79      P2I = P1I - (T1R*PTI+T1I*PTR)
80      P1R = PTR
81      P1I = PTI
82      T1R = T1R + RZR
83      T1I = T1I + RZI
84      AP2 = ZABS(P2R,P2I)
85      IF (AP1.LE.TEST) GO TO 10
86      IF (ITIME.EQ.2) GO TO 20
87      AK = ZABS(T1R,T1I)*0.5D0
88      FLAM = AK + SQRT(AK*AK-1.0D0)
89      RHO = MIN(AP2/AP1,FLAM)
90      TEST = TEST1*SQRT(RHO/(RHO*RHO-1.0D0))
91      ITIME = 2
92      GO TO 10
93   20 CONTINUE
94      KK = K + 1 - ID
95      AK = KK
96      T1R = AK
97      T1I = CZEROI
98      DFNU = FNU + (N-1)
99      P1R = 1.0D0/AP2
100      P1I = CZEROI
101      P2R = CZEROR
102      P2I = CZEROI
103      DO 30 I=1,KK
104        PTR = P1R
105        PTI = P1I
106        RAP1 = DFNU + T1R
107        TTR = RZR*RAP1
108        TTI = RZI*RAP1
109        P1R = (PTR*TTR-PTI*TTI) + P2R
110        P1I = (PTR*TTI+PTI*TTR) + P2I
111        P2R = PTR
112        P2I = PTI
113        T1R = T1R - CONER
114   30 CONTINUE
115      IF (P1R.NE.CZEROR .OR. P1I.NE.CZEROI) GO TO 40
116      P1R = TOL
117      P1I = TOL
118   40 CONTINUE
119      CALL ZDIV(P2R, P2I, P1R, P1I, CYR(N), CYI(N))
120      IF (N.EQ.1) RETURN
121      K = N - 1
122      AK = K
123      T1R = AK
124      T1I = CZEROI
125      CDFNUR = FNU*RZR
126      CDFNUI = FNU*RZI
127      DO 60 I=2,N
128        PTR = CDFNUR + (T1R*RZR-T1I*RZI) + CYR(K+1)
129        PTI = CDFNUI + (T1R*RZI+T1I*RZR) + CYI(K+1)
130        AK = ZABS(PTR,PTI)
131        IF (AK.NE.CZEROR) GO TO 50
132        PTR = TOL
133        PTI = TOL
134        AK = TOL*RT2
135   50   CONTINUE
136        RAK = CONER/AK
137        CYR(K) = RAK*PTR*RAK
138        CYI(K) = -RAK*PTI*RAK
139        T1R = T1R - CONER
140        K = K - 1
141   60 CONTINUE
142      RETURN
143      END
144