|
3217
|
1 SUBROUTINE ZMLRI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL) |
|
|
2 C***BEGIN PROLOGUE ZMLRI |
|
|
3 C***REFER TO ZBESI,ZBESK |
|
|
4 C |
|
|
5 C ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE |
|
|
6 C MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES. |
|
|
7 C |
|
|
8 C***ROUTINES CALLED DGAMLN,D1MACH,XZABS,XZEXP,XZLOG,ZMLT |
|
|
9 C***END PROLOGUE ZMLRI |
|
|
10 C COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z |
|
|
11 DOUBLE PRECISION ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI, |
|
|
12 * CNORMR, CONEI, CONER, FKAP, FKK, FLAM, FNF, FNU, PTI, PTR, P1I, |
|
|
13 * P1R, P2I, P2R, RAZ, RHO, RHO2, RZI, RZR, SCLE, STI, STR, SUMI, |
|
|
14 * SUMR, TFNF, TOL, TST, YI, YR, ZEROI, ZEROR, ZI, ZR, DGAMLN, |
|
|
15 * D1MACH, XZABS |
|
|
16 INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N, NZ |
|
|
17 DIMENSION YR(N), YI(N) |
|
|
18 DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 / |
|
|
19 SCLE = D1MACH(1)/TOL |
|
|
20 NZ=0 |
|
|
21 AZ = XZABS(ZR,ZI) |
|
|
22 IAZ = INT(SNGL(AZ)) |
|
|
23 IFNU = INT(SNGL(FNU)) |
|
|
24 INU = IFNU + N - 1 |
|
|
25 AT = DBLE(FLOAT(IAZ)) + 1.0D0 |
|
|
26 RAZ = 1.0D0/AZ |
|
|
27 STR = ZR*RAZ |
|
|
28 STI = -ZI*RAZ |
|
|
29 CKR = STR*AT*RAZ |
|
|
30 CKI = STI*AT*RAZ |
|
|
31 RZR = (STR+STR)*RAZ |
|
|
32 RZI = (STI+STI)*RAZ |
|
|
33 P1R = ZEROR |
|
|
34 P1I = ZEROI |
|
|
35 P2R = CONER |
|
|
36 P2I = CONEI |
|
|
37 ACK = (AT+1.0D0)*RAZ |
|
|
38 RHO = ACK + DSQRT(ACK*ACK-1.0D0) |
|
|
39 RHO2 = RHO*RHO |
|
|
40 TST = (RHO2+RHO2)/((RHO2-1.0D0)*(RHO-1.0D0)) |
|
|
41 TST = TST/TOL |
|
|
42 C----------------------------------------------------------------------- |
|
|
43 C COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES |
|
|
44 C----------------------------------------------------------------------- |
|
|
45 AK = AT |
|
|
46 DO 10 I=1,80 |
|
|
47 PTR = P2R |
|
|
48 PTI = P2I |
|
|
49 P2R = P1R - (CKR*PTR-CKI*PTI) |
|
|
50 P2I = P1I - (CKI*PTR+CKR*PTI) |
|
|
51 P1R = PTR |
|
|
52 P1I = PTI |
|
|
53 CKR = CKR + RZR |
|
|
54 CKI = CKI + RZI |
|
|
55 AP = XZABS(P2R,P2I) |
|
|
56 IF (AP.GT.TST*AK*AK) GO TO 20 |
|
|
57 AK = AK + 1.0D0 |
|
|
58 10 CONTINUE |
|
|
59 GO TO 110 |
|
|
60 20 CONTINUE |
|
|
61 I = I + 1 |
|
|
62 K = 0 |
|
|
63 IF (INU.LT.IAZ) GO TO 40 |
|
|
64 C----------------------------------------------------------------------- |
|
|
65 C COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS |
|
|
66 C----------------------------------------------------------------------- |
|
|
67 P1R = ZEROR |
|
|
68 P1I = ZEROI |
|
|
69 P2R = CONER |
|
|
70 P2I = CONEI |
|
|
71 AT = DBLE(FLOAT(INU)) + 1.0D0 |
|
|
72 STR = ZR*RAZ |
|
|
73 STI = -ZI*RAZ |
|
|
74 CKR = STR*AT*RAZ |
|
|
75 CKI = STI*AT*RAZ |
|
|
76 ACK = AT*RAZ |
|
|
77 TST = DSQRT(ACK/TOL) |
|
|
78 ITIME = 1 |
|
|
79 DO 30 K=1,80 |
|
|
80 PTR = P2R |
|
|
81 PTI = P2I |
|
|
82 P2R = P1R - (CKR*PTR-CKI*PTI) |
|
|
83 P2I = P1I - (CKR*PTI+CKI*PTR) |
|
|
84 P1R = PTR |
|
|
85 P1I = PTI |
|
|
86 CKR = CKR + RZR |
|
|
87 CKI = CKI + RZI |
|
|
88 AP = XZABS(P2R,P2I) |
|
|
89 IF (AP.LT.TST) GO TO 30 |
|
|
90 IF (ITIME.EQ.2) GO TO 40 |
|
|
91 ACK = XZABS(CKR,CKI) |
|
|
92 FLAM = ACK + DSQRT(ACK*ACK-1.0D0) |
|
|
93 FKAP = AP/XZABS(P1R,P1I) |
|
|
94 RHO = DMIN1(FLAM,FKAP) |
|
|
95 TST = TST*DSQRT(RHO/(RHO*RHO-1.0D0)) |
|
|
96 ITIME = 2 |
|
|
97 30 CONTINUE |
|
|
98 GO TO 110 |
|
|
99 40 CONTINUE |
|
|
100 C----------------------------------------------------------------------- |
|
|
101 C BACKWARD RECURRENCE AND SUM NORMALIZING RELATION |
|
|
102 C----------------------------------------------------------------------- |
|
|
103 K = K + 1 |
|
|
104 KK = MAX0(I+IAZ,K+INU) |
|
|
105 FKK = DBLE(FLOAT(KK)) |
|
|
106 P1R = ZEROR |
|
|
107 P1I = ZEROI |
|
|
108 C----------------------------------------------------------------------- |
|
|
109 C SCALE P2 AND SUM BY SCLE |
|
|
110 C----------------------------------------------------------------------- |
|
|
111 P2R = SCLE |
|
|
112 P2I = ZEROI |
|
|
113 FNF = FNU - DBLE(FLOAT(IFNU)) |
|
|
114 TFNF = FNF + FNF |
|
|
115 BK = DGAMLN(FKK+TFNF+1.0D0,IDUM) - DGAMLN(FKK+1.0D0,IDUM) - |
|
|
116 * DGAMLN(TFNF+1.0D0,IDUM) |
|
|
117 BK = DEXP(BK) |
|
|
118 SUMR = ZEROR |
|
|
119 SUMI = ZEROI |
|
|
120 KM = KK - INU |
|
|
121 DO 50 I=1,KM |
|
|
122 PTR = P2R |
|
|
123 PTI = P2I |
|
|
124 P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI) |
|
|
125 P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI) |
|
|
126 P1R = PTR |
|
|
127 P1I = PTI |
|
|
128 AK = 1.0D0 - TFNF/(FKK+TFNF) |
|
|
129 ACK = BK*AK |
|
|
130 SUMR = SUMR + (ACK+BK)*P1R |
|
|
131 SUMI = SUMI + (ACK+BK)*P1I |
|
|
132 BK = ACK |
|
|
133 FKK = FKK - 1.0D0 |
|
|
134 50 CONTINUE |
|
|
135 YR(N) = P2R |
|
|
136 YI(N) = P2I |
|
|
137 IF (N.EQ.1) GO TO 70 |
|
|
138 DO 60 I=2,N |
|
|
139 PTR = P2R |
|
|
140 PTI = P2I |
|
|
141 P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI) |
|
|
142 P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI) |
|
|
143 P1R = PTR |
|
|
144 P1I = PTI |
|
|
145 AK = 1.0D0 - TFNF/(FKK+TFNF) |
|
|
146 ACK = BK*AK |
|
|
147 SUMR = SUMR + (ACK+BK)*P1R |
|
|
148 SUMI = SUMI + (ACK+BK)*P1I |
|
|
149 BK = ACK |
|
|
150 FKK = FKK - 1.0D0 |
|
|
151 M = N - I + 1 |
|
|
152 YR(M) = P2R |
|
|
153 YI(M) = P2I |
|
|
154 60 CONTINUE |
|
|
155 70 CONTINUE |
|
|
156 IF (IFNU.LE.0) GO TO 90 |
|
|
157 DO 80 I=1,IFNU |
|
|
158 PTR = P2R |
|
|
159 PTI = P2I |
|
|
160 P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI) |
|
|
161 P2I = P1I + (FKK+FNF)*(RZR*PTI+RZI*PTR) |
|
|
162 P1R = PTR |
|
|
163 P1I = PTI |
|
|
164 AK = 1.0D0 - TFNF/(FKK+TFNF) |
|
|
165 ACK = BK*AK |
|
|
166 SUMR = SUMR + (ACK+BK)*P1R |
|
|
167 SUMI = SUMI + (ACK+BK)*P1I |
|
|
168 BK = ACK |
|
|
169 FKK = FKK - 1.0D0 |
|
|
170 80 CONTINUE |
|
|
171 90 CONTINUE |
|
|
172 PTR = ZR |
|
|
173 PTI = ZI |
|
|
174 IF (KODE.EQ.2) PTR = ZEROR |
|
|
175 CALL XZLOG(RZR, RZI, STR, STI, IDUM) |
|
|
176 P1R = -FNF*STR + PTR |
|
|
177 P1I = -FNF*STI + PTI |
|
|
178 AP = DGAMLN(1.0D0+FNF,IDUM) |
|
|
179 PTR = P1R - AP |
|
|
180 PTI = P1I |
|
|
181 C----------------------------------------------------------------------- |
|
|
182 C THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW |
|
|
183 C IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES |
|
|
184 C----------------------------------------------------------------------- |
|
|
185 P2R = P2R + SUMR |
|
|
186 P2I = P2I + SUMI |
|
|
187 AP = XZABS(P2R,P2I) |
|
|
188 P1R = 1.0D0/AP |
|
|
189 CALL XZEXP(PTR, PTI, STR, STI) |
|
|
190 CKR = STR*P1R |
|
|
191 CKI = STI*P1R |
|
|
192 PTR = P2R*P1R |
|
|
193 PTI = -P2I*P1R |
|
|
194 CALL ZMLT(CKR, CKI, PTR, PTI, CNORMR, CNORMI) |
|
|
195 DO 100 I=1,N |
|
|
196 STR = YR(I)*CNORMR - YI(I)*CNORMI |
|
|
197 YI(I) = YR(I)*CNORMI + YI(I)*CNORMR |
|
|
198 YR(I) = STR |
|
|
199 100 CONTINUE |
|
|
200 RETURN |
|
|
201 110 CONTINUE |
|
|
202 NZ=-2 |
|
|
203 RETURN |
|
|
204 END |
