|
3217
|
1 SUBROUTINE ZRATI(ZR, ZI, FNU, N, CYR, CYI, TOL) |
|
|
2 C***BEGIN PROLOGUE ZRATI |
|
|
3 C***REFER TO ZBESI,ZBESK,ZBESH |
|
|
4 C |
|
|
5 C ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD |
|
|
6 C RECURRENCE. THE STARTING INDEX IS DETERMINED BY FORWARD |
|
|
7 C RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B, |
|
|
8 C MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973, |
|
|
9 C BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER, |
|
|
10 C BY D. J. SOOKNE. |
|
|
11 C |
|
|
12 C***ROUTINES CALLED XZABS,ZDIV |
|
|
13 C***END PROLOGUE ZRATI |
|
|
14 C 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, XZABS |
|
|
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 = XZABS(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 = XZABS(P2R,P2I) |
|
|
46 AP1 = XZABS(P1R,P1I) |
|
|
47 C----------------------------------------------------------------------- |
|
|
48 C THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU |
|
|
49 C GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT |
|
|
50 C P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR |
|
|
51 C PREMATURELY. |
|
|
52 C----------------------------------------------------------------------- |
|
|
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 = XZABS(P2R,P2I) |
|
|
74 IF (AP1.LE.TEST) GO TO 10 |
|
|
75 IF (ITIME.EQ.2) GO TO 20 |
|
|
76 AK = XZABS(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 = XZABS(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 |
