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