|
3217
|
1 SUBROUTINE ZBKNU(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, |
|
|
2 * ALIM) |
|
|
3 C***BEGIN PROLOGUE ZBKNU |
|
|
4 C***REFER TO ZBESI,ZBESK,ZAIRY,ZBESH |
|
|
5 C |
|
|
6 C ZBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE. |
|
|
7 C |
|
|
8 C***ROUTINES CALLED DGAMLN,I1MACH,D1MACH,ZKSCL,ZSHCH,ZUCHK,XZABS,ZDIV, |
|
|
9 C XZEXP,XZLOG,ZMLT,XZSQRT |
|
|
10 C***END PROLOGUE ZBKNU |
|
|
11 C |
|
|
12 DOUBLE PRECISION AA, AK, ALIM, ASCLE, A1, A2, BB, BK, BRY, CAZ, |
|
|
13 * CBI, CBR, CC, CCHI, CCHR, CKI, CKR, COEFI, COEFR, CONEI, CONER, |
|
|
14 * CRSCR, CSCLR, CSHI, CSHR, CSI, CSR, CSRR, CSSR, CTWOR, |
|
|
15 * CZEROI, CZEROR, CZI, CZR, DNU, DNU2, DPI, ELIM, ETEST, FC, FHS, |
|
|
16 * FI, FK, FKS, FMUI, FMUR, FNU, FPI, FR, G1, G2, HPI, PI, PR, PTI, |
|
|
17 * PTR, P1I, P1R, P2I, P2M, P2R, QI, QR, RAK, RCAZ, RTHPI, RZI, |
|
|
18 * RZR, R1, S, SMUI, SMUR, SPI, STI, STR, S1I, S1R, S2I, S2R, TM, |
|
|
19 * TOL, TTH, T1, T2, YI, YR, ZI, ZR, DGAMLN, D1MACH, XZABS, ELM, |
|
|
20 * CELMR, ZDR, ZDI, AS, ALAS, HELIM, CYR, CYI |
|
|
21 INTEGER I, IFLAG, INU, K, KFLAG, KK, KMAX, KODE, KODED, N, NZ, |
|
|
22 * IDUM, I1MACH, J, IC, INUB, NW |
|
|
23 DIMENSION YR(N), YI(N), CC(8), CSSR(3), CSRR(3), BRY(3), CYR(2), |
|
|
24 * CYI(2) |
|
|
25 C COMPLEX Z,Y,A,B,RZ,SMU,FU,FMU,F,FLRZ,CZ,S1,S2,CSH,CCH |
|
|
26 C COMPLEX CK,P,Q,COEF,P1,P2,CBK,PT,CZERO,CONE,CTWO,ST,EZ,CS,DK |
|
|
27 C |
|
|
28 DATA KMAX / 30 / |
|
|
29 DATA CZEROR,CZEROI,CONER,CONEI,CTWOR,R1/ |
|
|
30 1 0.0D0 , 0.0D0 , 1.0D0 , 0.0D0 , 2.0D0 , 2.0D0 / |
|
|
31 DATA DPI, RTHPI, SPI ,HPI, FPI, TTH / |
|
|
32 1 3.14159265358979324D0, 1.25331413731550025D0, |
|
|
33 2 1.90985931710274403D0, 1.57079632679489662D0, |
|
|
34 3 1.89769999331517738D0, 6.66666666666666666D-01/ |
|
|
35 DATA CC(1), CC(2), CC(3), CC(4), CC(5), CC(6), CC(7), CC(8)/ |
|
|
36 1 5.77215664901532861D-01, -4.20026350340952355D-02, |
|
|
37 2 -4.21977345555443367D-02, 7.21894324666309954D-03, |
|
|
38 3 -2.15241674114950973D-04, -2.01348547807882387D-05, |
|
|
39 4 1.13302723198169588D-06, 6.11609510448141582D-09/ |
|
|
40 C |
|
|
41 CAZ = XZABS(ZR,ZI) |
|
|
42 CSCLR = 1.0D0/TOL |
|
|
43 CRSCR = TOL |
|
|
44 CSSR(1) = CSCLR |
|
|
45 CSSR(2) = 1.0D0 |
|
|
46 CSSR(3) = CRSCR |
|
|
47 CSRR(1) = CRSCR |
|
|
48 CSRR(2) = 1.0D0 |
|
|
49 CSRR(3) = CSCLR |
|
|
50 BRY(1) = 1.0D+3*D1MACH(1)/TOL |
|
|
51 BRY(2) = 1.0D0/BRY(1) |
|
|
52 BRY(3) = D1MACH(2) |
|
|
53 NZ = 0 |
|
|
54 IFLAG = 0 |
|
|
55 KODED = KODE |
|
|
56 RCAZ = 1.0D0/CAZ |
|
|
57 STR = ZR*RCAZ |
|
|
58 STI = -ZI*RCAZ |
|
|
59 RZR = (STR+STR)*RCAZ |
|
|
60 RZI = (STI+STI)*RCAZ |
|
|
61 INU = INT(SNGL(FNU+0.5D0)) |
|
|
62 DNU = FNU - DBLE(FLOAT(INU)) |
|
|
63 IF (DABS(DNU).EQ.0.5D0) GO TO 110 |
|
|
64 DNU2 = 0.0D0 |
|
|
65 IF (DABS(DNU).GT.TOL) DNU2 = DNU*DNU |
|
|
66 IF (CAZ.GT.R1) GO TO 110 |
|
|
67 C----------------------------------------------------------------------- |
|
|
68 C SERIES FOR CABS(Z).LE.R1 |
|
|
69 C----------------------------------------------------------------------- |
|
|
70 FC = 1.0D0 |
|
|
71 CALL XZLOG(RZR, RZI, SMUR, SMUI, IDUM) |
|
|
72 FMUR = SMUR*DNU |
|
|
73 FMUI = SMUI*DNU |
|
|
74 CALL ZSHCH(FMUR, FMUI, CSHR, CSHI, CCHR, CCHI) |
|
|
75 IF (DNU.EQ.0.0D0) GO TO 10 |
|
|
76 FC = DNU*DPI |
|
|
77 FC = FC/DSIN(FC) |
|
|
78 SMUR = CSHR/DNU |
|
|
79 SMUI = CSHI/DNU |
|
|
80 10 CONTINUE |
|
|
81 A2 = 1.0D0 + DNU |
|
|
82 C----------------------------------------------------------------------- |
|
|
83 C GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU) |
|
|
84 C----------------------------------------------------------------------- |
|
|
85 T2 = DEXP(-DGAMLN(A2,IDUM)) |
|
|
86 T1 = 1.0D0/(T2*FC) |
|
|
87 IF (DABS(DNU).GT.0.1D0) GO TO 40 |
|
|
88 C----------------------------------------------------------------------- |
|
|
89 C SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU) |
|
|
90 C----------------------------------------------------------------------- |
|
|
91 AK = 1.0D0 |
|
|
92 S = CC(1) |
|
|
93 DO 20 K=2,8 |
|
|
94 AK = AK*DNU2 |
|
|
95 TM = CC(K)*AK |
|
|
96 S = S + TM |
|
|
97 IF (DABS(TM).LT.TOL) GO TO 30 |
|
|
98 20 CONTINUE |
|
|
99 30 G1 = -S |
|
|
100 GO TO 50 |
|
|
101 40 CONTINUE |
|
|
102 G1 = (T1-T2)/(DNU+DNU) |
|
|
103 50 CONTINUE |
|
|
104 G2 = (T1+T2)*0.5D0 |
|
|
105 FR = FC*(CCHR*G1+SMUR*G2) |
|
|
106 FI = FC*(CCHI*G1+SMUI*G2) |
|
|
107 CALL XZEXP(FMUR, FMUI, STR, STI) |
|
|
108 PR = 0.5D0*STR/T2 |
|
|
109 PI = 0.5D0*STI/T2 |
|
|
110 CALL ZDIV(0.5D0, 0.0D0, STR, STI, PTR, PTI) |
|
|
111 QR = PTR/T1 |
|
|
112 QI = PTI/T1 |
|
|
113 S1R = FR |
|
|
114 S1I = FI |
|
|
115 S2R = PR |
|
|
116 S2I = PI |
|
|
117 AK = 1.0D0 |
|
|
118 A1 = 1.0D0 |
|
|
119 CKR = CONER |
|
|
120 CKI = CONEI |
|
|
121 BK = 1.0D0 - DNU2 |
|
|
122 IF (INU.GT.0 .OR. N.GT.1) GO TO 80 |
|
|
123 C----------------------------------------------------------------------- |
|
|
124 C GENERATE K(FNU,Z), 0.0D0 .LE. FNU .LT. 0.5D0 AND N=1 |
|
|
125 C----------------------------------------------------------------------- |
|
|
126 IF (CAZ.LT.TOL) GO TO 70 |
|
|
127 CALL ZMLT(ZR, ZI, ZR, ZI, CZR, CZI) |
|
|
128 CZR = 0.25D0*CZR |
|
|
129 CZI = 0.25D0*CZI |
|
|
130 T1 = 0.25D0*CAZ*CAZ |
|
|
131 60 CONTINUE |
|
|
132 FR = (FR*AK+PR+QR)/BK |
|
|
133 FI = (FI*AK+PI+QI)/BK |
|
|
134 STR = 1.0D0/(AK-DNU) |
|
|
135 PR = PR*STR |
|
|
136 PI = PI*STR |
|
|
137 STR = 1.0D0/(AK+DNU) |
|
|
138 QR = QR*STR |
|
|
139 QI = QI*STR |
|
|
140 STR = CKR*CZR - CKI*CZI |
|
|
141 RAK = 1.0D0/AK |
|
|
142 CKI = (CKR*CZI+CKI*CZR)*RAK |
|
|
143 CKR = STR*RAK |
|
|
144 S1R = CKR*FR - CKI*FI + S1R |
|
|
145 S1I = CKR*FI + CKI*FR + S1I |
|
|
146 A1 = A1*T1*RAK |
|
|
147 BK = BK + AK + AK + 1.0D0 |
|
|
148 AK = AK + 1.0D0 |
|
|
149 IF (A1.GT.TOL) GO TO 60 |
|
|
150 70 CONTINUE |
|
|
151 YR(1) = S1R |
|
|
152 YI(1) = S1I |
|
|
153 IF (KODED.EQ.1) RETURN |
|
|
154 CALL XZEXP(ZR, ZI, STR, STI) |
|
|
155 CALL ZMLT(S1R, S1I, STR, STI, YR(1), YI(1)) |
|
|
156 RETURN |
|
|
157 C----------------------------------------------------------------------- |
|
|
158 C GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE |
|
|
159 C----------------------------------------------------------------------- |
|
|
160 80 CONTINUE |
|
|
161 IF (CAZ.LT.TOL) GO TO 100 |
|
|
162 CALL ZMLT(ZR, ZI, ZR, ZI, CZR, CZI) |
|
|
163 CZR = 0.25D0*CZR |
|
|
164 CZI = 0.25D0*CZI |
|
|
165 T1 = 0.25D0*CAZ*CAZ |
|
|
166 90 CONTINUE |
|
|
167 FR = (FR*AK+PR+QR)/BK |
|
|
168 FI = (FI*AK+PI+QI)/BK |
|
|
169 STR = 1.0D0/(AK-DNU) |
|
|
170 PR = PR*STR |
|
|
171 PI = PI*STR |
|
|
172 STR = 1.0D0/(AK+DNU) |
|
|
173 QR = QR*STR |
|
|
174 QI = QI*STR |
|
|
175 STR = CKR*CZR - CKI*CZI |
|
|
176 RAK = 1.0D0/AK |
|
|
177 CKI = (CKR*CZI+CKI*CZR)*RAK |
|
|
178 CKR = STR*RAK |
|
|
179 S1R = CKR*FR - CKI*FI + S1R |
|
|
180 S1I = CKR*FI + CKI*FR + S1I |
|
|
181 STR = PR - FR*AK |
|
|
182 STI = PI - FI*AK |
|
|
183 S2R = CKR*STR - CKI*STI + S2R |
|
|
184 S2I = CKR*STI + CKI*STR + S2I |
|
|
185 A1 = A1*T1*RAK |
|
|
186 BK = BK + AK + AK + 1.0D0 |
|
|
187 AK = AK + 1.0D0 |
|
|
188 IF (A1.GT.TOL) GO TO 90 |
|
|
189 100 CONTINUE |
|
|
190 KFLAG = 2 |
|
|
191 A1 = FNU + 1.0D0 |
|
|
192 AK = A1*DABS(SMUR) |
|
|
193 IF (AK.GT.ALIM) KFLAG = 3 |
|
|
194 STR = CSSR(KFLAG) |
|
|
195 P2R = S2R*STR |
|
|
196 P2I = S2I*STR |
|
|
197 CALL ZMLT(P2R, P2I, RZR, RZI, S2R, S2I) |
|
|
198 S1R = S1R*STR |
|
|
199 S1I = S1I*STR |
|
|
200 IF (KODED.EQ.1) GO TO 210 |
|
|
201 CALL XZEXP(ZR, ZI, FR, FI) |
|
|
202 CALL ZMLT(S1R, S1I, FR, FI, S1R, S1I) |
|
|
203 CALL ZMLT(S2R, S2I, FR, FI, S2R, S2I) |
|
|
204 GO TO 210 |
|
|
205 C----------------------------------------------------------------------- |
|
|
206 C IFLAG=0 MEANS NO UNDERFLOW OCCURRED |
|
|
207 C IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH |
|
|
208 C KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD |
|
|
209 C RECURSION |
|
|
210 C----------------------------------------------------------------------- |
|
|
211 110 CONTINUE |
|
|
212 CALL XZSQRT(ZR, ZI, STR, STI) |
|
|
213 CALL ZDIV(RTHPI, CZEROI, STR, STI, COEFR, COEFI) |
|
|
214 KFLAG = 2 |
|
|
215 IF (KODED.EQ.2) GO TO 120 |
|
|
216 IF (ZR.GT.ALIM) GO TO 290 |
|
|
217 C BLANK LINE |
|
|
218 STR = DEXP(-ZR)*CSSR(KFLAG) |
|
|
219 STI = -STR*DSIN(ZI) |
|
|
220 STR = STR*DCOS(ZI) |
|
|
221 CALL ZMLT(COEFR, COEFI, STR, STI, COEFR, COEFI) |
|
|
222 120 CONTINUE |
|
|
223 IF (DABS(DNU).EQ.0.5D0) GO TO 300 |
|
|
224 C----------------------------------------------------------------------- |
|
|
225 C MILLER ALGORITHM FOR CABS(Z).GT.R1 |
|
|
226 C----------------------------------------------------------------------- |
|
|
227 AK = DCOS(DPI*DNU) |
|
|
228 AK = DABS(AK) |
|
|
229 IF (AK.EQ.CZEROR) GO TO 300 |
|
|
230 FHS = DABS(0.25D0-DNU2) |
|
|
231 IF (FHS.EQ.CZEROR) GO TO 300 |
|
|
232 C----------------------------------------------------------------------- |
|
|
233 C COMPUTE R2=F(E). IF CABS(Z).GE.R2, USE FORWARD RECURRENCE TO |
|
|
234 C DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON |
|
|
235 C 12.LE.E.LE.60. E IS COMPUTED FROM 2**(-E)=B**(1-I1MACH(14))= |
|
|
236 C TOL WHERE B IS THE BASE OF THE ARITHMETIC. |
|
|
237 C----------------------------------------------------------------------- |
|
|
238 T1 = DBLE(FLOAT(I1MACH(14)-1)) |
|
|
239 T1 = T1*D1MACH(5)*3.321928094D0 |
|
|
240 T1 = DMAX1(T1,12.0D0) |
|
|
241 T1 = DMIN1(T1,60.0D0) |
|
|
242 T2 = TTH*T1 - 6.0D0 |
|
|
243 IF (ZR.NE.0.0D0) GO TO 130 |
|
|
244 T1 = HPI |
|
|
245 GO TO 140 |
|
|
246 130 CONTINUE |
|
|
247 T1 = DATAN(ZI/ZR) |
|
|
248 T1 = DABS(T1) |
|
|
249 140 CONTINUE |
|
|
250 IF (T2.GT.CAZ) GO TO 170 |
|
|
251 C----------------------------------------------------------------------- |
|
|
252 C FORWARD RECURRENCE LOOP WHEN CABS(Z).GE.R2 |
|
|
253 C----------------------------------------------------------------------- |
|
|
254 ETEST = AK/(DPI*CAZ*TOL) |
|
|
255 FK = CONER |
|
|
256 IF (ETEST.LT.CONER) GO TO 180 |
|
|
257 FKS = CTWOR |
|
|
258 CKR = CAZ + CAZ + CTWOR |
|
|
259 P1R = CZEROR |
|
|
260 P2R = CONER |
|
|
261 DO 150 I=1,KMAX |
|
|
262 AK = FHS/FKS |
|
|
263 CBR = CKR/(FK+CONER) |
|
|
264 PTR = P2R |
|
|
265 P2R = CBR*P2R - P1R*AK |
|
|
266 P1R = PTR |
|
|
267 CKR = CKR + CTWOR |
|
|
268 FKS = FKS + FK + FK + CTWOR |
|
|
269 FHS = FHS + FK + FK |
|
|
270 FK = FK + CONER |
|
|
271 STR = DABS(P2R)*FK |
|
|
272 IF (ETEST.LT.STR) GO TO 160 |
|
|
273 150 CONTINUE |
|
|
274 GO TO 310 |
|
|
275 160 CONTINUE |
|
|
276 FK = FK + SPI*T1*DSQRT(T2/CAZ) |
|
|
277 FHS = DABS(0.25D0-DNU2) |
|
|
278 GO TO 180 |
|
|
279 170 CONTINUE |
|
|
280 C----------------------------------------------------------------------- |
|
|
281 C COMPUTE BACKWARD INDEX K FOR CABS(Z).LT.R2 |
|
|
282 C----------------------------------------------------------------------- |
|
|
283 A2 = DSQRT(CAZ) |
|
|
284 AK = FPI*AK/(TOL*DSQRT(A2)) |
|
|
285 AA = 3.0D0*T1/(1.0D0+CAZ) |
|
|
286 BB = 14.7D0*T1/(28.0D0+CAZ) |
|
|
287 AK = (DLOG(AK)+CAZ*DCOS(AA)/(1.0D0+0.008D0*CAZ))/DCOS(BB) |
|
|
288 FK = 0.12125D0*AK*AK/CAZ + 1.5D0 |
|
|
289 180 CONTINUE |
|
|
290 C----------------------------------------------------------------------- |
|
|
291 C BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM |
|
|
292 C----------------------------------------------------------------------- |
|
|
293 K = INT(SNGL(FK)) |
|
|
294 FK = DBLE(FLOAT(K)) |
|
|
295 FKS = FK*FK |
|
|
296 P1R = CZEROR |
|
|
297 P1I = CZEROI |
|
|
298 P2R = TOL |
|
|
299 P2I = CZEROI |
|
|
300 CSR = P2R |
|
|
301 CSI = P2I |
|
|
302 DO 190 I=1,K |
|
|
303 A1 = FKS - FK |
|
|
304 AK = (FKS+FK)/(A1+FHS) |
|
|
305 RAK = 2.0D0/(FK+CONER) |
|
|
306 CBR = (FK+ZR)*RAK |
|
|
307 CBI = ZI*RAK |
|
|
308 PTR = P2R |
|
|
309 PTI = P2I |
|
|
310 P2R = (PTR*CBR-PTI*CBI-P1R)*AK |
|
|
311 P2I = (PTI*CBR+PTR*CBI-P1I)*AK |
|
|
312 P1R = PTR |
|
|
313 P1I = PTI |
|
|
314 CSR = CSR + P2R |
|
|
315 CSI = CSI + P2I |
|
|
316 FKS = A1 - FK + CONER |
|
|
317 FK = FK - CONER |
|
|
318 190 CONTINUE |
|
|
319 C----------------------------------------------------------------------- |
|
|
320 C COMPUTE (P2/CS)=(P2/CABS(CS))*(CONJG(CS)/CABS(CS)) FOR BETTER |
|
|
321 C SCALING |
|
|
322 C----------------------------------------------------------------------- |
|
|
323 TM = XZABS(CSR,CSI) |
|
|
324 PTR = 1.0D0/TM |
|
|
325 S1R = P2R*PTR |
|
|
326 S1I = P2I*PTR |
|
|
327 CSR = CSR*PTR |
|
|
328 CSI = -CSI*PTR |
|
|
329 CALL ZMLT(COEFR, COEFI, S1R, S1I, STR, STI) |
|
|
330 CALL ZMLT(STR, STI, CSR, CSI, S1R, S1I) |
|
|
331 IF (INU.GT.0 .OR. N.GT.1) GO TO 200 |
|
|
332 ZDR = ZR |
|
|
333 ZDI = ZI |
|
|
334 IF(IFLAG.EQ.1) GO TO 270 |
|
|
335 GO TO 240 |
|
|
336 200 CONTINUE |
|
|
337 C----------------------------------------------------------------------- |
|
|
338 C COMPUTE P1/P2=(P1/CABS(P2)*CONJG(P2)/CABS(P2) FOR SCALING |
|
|
339 C----------------------------------------------------------------------- |
|
|
340 TM = XZABS(P2R,P2I) |
|
|
341 PTR = 1.0D0/TM |
|
|
342 P1R = P1R*PTR |
|
|
343 P1I = P1I*PTR |
|
|
344 P2R = P2R*PTR |
|
|
345 P2I = -P2I*PTR |
|
|
346 CALL ZMLT(P1R, P1I, P2R, P2I, PTR, PTI) |
|
|
347 STR = DNU + 0.5D0 - PTR |
|
|
348 STI = -PTI |
|
|
349 CALL ZDIV(STR, STI, ZR, ZI, STR, STI) |
|
|
350 STR = STR + 1.0D0 |
|
|
351 CALL ZMLT(STR, STI, S1R, S1I, S2R, S2I) |
|
|
352 C----------------------------------------------------------------------- |
|
|
353 C FORWARD RECURSION ON THE THREE TERM RECURSION WITH RELATION WITH |
|
|
354 C SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3 |
|
|
355 C----------------------------------------------------------------------- |
|
|
356 210 CONTINUE |
|
|
357 STR = DNU + 1.0D0 |
|
|
358 CKR = STR*RZR |
|
|
359 CKI = STR*RZI |
|
|
360 IF (N.EQ.1) INU = INU - 1 |
|
|
361 IF (INU.GT.0) GO TO 220 |
|
|
362 IF (N.GT.1) GO TO 215 |
|
|
363 S1R = S2R |
|
|
364 S1I = S2I |
|
|
365 215 CONTINUE |
|
|
366 ZDR = ZR |
|
|
367 ZDI = ZI |
|
|
368 IF(IFLAG.EQ.1) GO TO 270 |
|
|
369 GO TO 240 |
|
|
370 220 CONTINUE |
|
|
371 INUB = 1 |
|
|
372 IF(IFLAG.EQ.1) GO TO 261 |
|
|
373 225 CONTINUE |
|
|
374 P1R = CSRR(KFLAG) |
|
|
375 ASCLE = BRY(KFLAG) |
|
|
376 DO 230 I=INUB,INU |
|
|
377 STR = S2R |
|
|
378 STI = S2I |
|
|
379 S2R = CKR*STR - CKI*STI + S1R |
|
|
380 S2I = CKR*STI + CKI*STR + S1I |
|
|
381 S1R = STR |
|
|
382 S1I = STI |
|
|
383 CKR = CKR + RZR |
|
|
384 CKI = CKI + RZI |
|
|
385 IF (KFLAG.GE.3) GO TO 230 |
|
|
386 P2R = S2R*P1R |
|
|
387 P2I = S2I*P1R |
|
|
388 STR = DABS(P2R) |
|
|
389 STI = DABS(P2I) |
|
|
390 P2M = DMAX1(STR,STI) |
|
|
391 IF (P2M.LE.ASCLE) GO TO 230 |
|
|
392 KFLAG = KFLAG + 1 |
|
|
393 ASCLE = BRY(KFLAG) |
|
|
394 S1R = S1R*P1R |
|
|
395 S1I = S1I*P1R |
|
|
396 S2R = P2R |
|
|
397 S2I = P2I |
|
|
398 STR = CSSR(KFLAG) |
|
|
399 S1R = S1R*STR |
|
|
400 S1I = S1I*STR |
|
|
401 S2R = S2R*STR |
|
|
402 S2I = S2I*STR |
|
|
403 P1R = CSRR(KFLAG) |
|
|
404 230 CONTINUE |
|
|
405 IF (N.NE.1) GO TO 240 |
|
|
406 S1R = S2R |
|
|
407 S1I = S2I |
|
|
408 240 CONTINUE |
|
|
409 STR = CSRR(KFLAG) |
|
|
410 YR(1) = S1R*STR |
|
|
411 YI(1) = S1I*STR |
|
|
412 IF (N.EQ.1) RETURN |
|
|
413 YR(2) = S2R*STR |
|
|
414 YI(2) = S2I*STR |
|
|
415 IF (N.EQ.2) RETURN |
|
|
416 KK = 2 |
|
|
417 250 CONTINUE |
|
|
418 KK = KK + 1 |
|
|
419 IF (KK.GT.N) RETURN |
|
|
420 P1R = CSRR(KFLAG) |
|
|
421 ASCLE = BRY(KFLAG) |
|
|
422 DO 260 I=KK,N |
|
|
423 P2R = S2R |
|
|
424 P2I = S2I |
|
|
425 S2R = CKR*P2R - CKI*P2I + S1R |
|
|
426 S2I = CKI*P2R + CKR*P2I + S1I |
|
|
427 S1R = P2R |
|
|
428 S1I = P2I |
|
|
429 CKR = CKR + RZR |
|
|
430 CKI = CKI + RZI |
|
|
431 P2R = S2R*P1R |
|
|
432 P2I = S2I*P1R |
|
|
433 YR(I) = P2R |
|
|
434 YI(I) = P2I |
|
|
435 IF (KFLAG.GE.3) GO TO 260 |
|
|
436 STR = DABS(P2R) |
|
|
437 STI = DABS(P2I) |
|
|
438 P2M = DMAX1(STR,STI) |
|
|
439 IF (P2M.LE.ASCLE) GO TO 260 |
|
|
440 KFLAG = KFLAG + 1 |
|
|
441 ASCLE = BRY(KFLAG) |
|
|
442 S1R = S1R*P1R |
|
|
443 S1I = S1I*P1R |
|
|
444 S2R = P2R |
|
|
445 S2I = P2I |
|
|
446 STR = CSSR(KFLAG) |
|
|
447 S1R = S1R*STR |
|
|
448 S1I = S1I*STR |
|
|
449 S2R = S2R*STR |
|
|
450 S2I = S2I*STR |
|
|
451 P1R = CSRR(KFLAG) |
|
|
452 260 CONTINUE |
|
|
453 RETURN |
|
|
454 C----------------------------------------------------------------------- |
|
|
455 C IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW |
|
|
456 C----------------------------------------------------------------------- |
|
|
457 261 CONTINUE |
|
|
458 HELIM = 0.5D0*ELIM |
|
|
459 ELM = DEXP(-ELIM) |
|
|
460 CELMR = ELM |
|
|
461 ASCLE = BRY(1) |
|
|
462 ZDR = ZR |
|
|
463 ZDI = ZI |
|
|
464 IC = -1 |
|
|
465 J = 2 |
|
|
466 DO 262 I=1,INU |
|
|
467 STR = S2R |
|
|
468 STI = S2I |
|
|
469 S2R = STR*CKR-STI*CKI+S1R |
|
|
470 S2I = STI*CKR+STR*CKI+S1I |
|
|
471 S1R = STR |
|
|
472 S1I = STI |
|
|
473 CKR = CKR+RZR |
|
|
474 CKI = CKI+RZI |
|
|
475 AS = XZABS(S2R,S2I) |
|
|
476 ALAS = DLOG(AS) |
|
|
477 P2R = -ZDR+ALAS |
|
|
478 IF(P2R.LT.(-ELIM)) GO TO 263 |
|
|
479 CALL XZLOG(S2R,S2I,STR,STI,IDUM) |
|
|
480 P2R = -ZDR+STR |
|
|
481 P2I = -ZDI+STI |
|
|
482 P2M = DEXP(P2R)/TOL |
|
|
483 P1R = P2M*DCOS(P2I) |
|
|
484 P1I = P2M*DSIN(P2I) |
|
|
485 CALL ZUCHK(P1R,P1I,NW,ASCLE,TOL) |
|
|
486 IF(NW.NE.0) GO TO 263 |
|
|
487 J = 3 - J |
|
|
488 CYR(J) = P1R |
|
|
489 CYI(J) = P1I |
|
|
490 IF(IC.EQ.(I-1)) GO TO 264 |
|
|
491 IC = I |
|
|
492 GO TO 262 |
|
|
493 263 CONTINUE |
|
|
494 IF(ALAS.LT.HELIM) GO TO 262 |
|
|
495 ZDR = ZDR-ELIM |
|
|
496 S1R = S1R*CELMR |
|
|
497 S1I = S1I*CELMR |
|
|
498 S2R = S2R*CELMR |
|
|
499 S2I = S2I*CELMR |
|
|
500 262 CONTINUE |
|
|
501 IF(N.NE.1) GO TO 270 |
|
|
502 S1R = S2R |
|
|
503 S1I = S2I |
|
|
504 GO TO 270 |
|
|
505 264 CONTINUE |
|
|
506 KFLAG = 1 |
|
|
507 INUB = I+1 |
|
|
508 S2R = CYR(J) |
|
|
509 S2I = CYI(J) |
|
|
510 J = 3 - J |
|
|
511 S1R = CYR(J) |
|
|
512 S1I = CYI(J) |
|
|
513 IF(INUB.LE.INU) GO TO 225 |
|
|
514 IF(N.NE.1) GO TO 240 |
|
|
515 S1R = S2R |
|
|
516 S1I = S2I |
|
|
517 GO TO 240 |
|
|
518 270 CONTINUE |
|
|
519 YR(1) = S1R |
|
|
520 YI(1) = S1I |
|
|
521 IF(N.EQ.1) GO TO 280 |
|
|
522 YR(2) = S2R |
|
|
523 YI(2) = S2I |
|
|
524 280 CONTINUE |
|
|
525 ASCLE = BRY(1) |
|
|
526 CALL ZKSCL(ZDR,ZDI,FNU,N,YR,YI,NZ,RZR,RZI,ASCLE,TOL,ELIM) |
|
|
527 INU = N - NZ |
|
|
528 IF (INU.LE.0) RETURN |
|
|
529 KK = NZ + 1 |
|
|
530 S1R = YR(KK) |
|
|
531 S1I = YI(KK) |
|
|
532 YR(KK) = S1R*CSRR(1) |
|
|
533 YI(KK) = S1I*CSRR(1) |
|
|
534 IF (INU.EQ.1) RETURN |
|
|
535 KK = NZ + 2 |
|
|
536 S2R = YR(KK) |
|
|
537 S2I = YI(KK) |
|
|
538 YR(KK) = S2R*CSRR(1) |
|
|
539 YI(KK) = S2I*CSRR(1) |
|
|
540 IF (INU.EQ.2) RETURN |
|
|
541 T2 = FNU + DBLE(FLOAT(KK-1)) |
|
|
542 CKR = T2*RZR |
|
|
543 CKI = T2*RZI |
|
|
544 KFLAG = 1 |
|
|
545 GO TO 250 |
|
|
546 290 CONTINUE |
|
|
547 C----------------------------------------------------------------------- |
|
|
548 C SCALE BY DEXP(Z), IFLAG = 1 CASES |
|
|
549 C----------------------------------------------------------------------- |
|
|
550 KODED = 2 |
|
|
551 IFLAG = 1 |
|
|
552 KFLAG = 2 |
|
|
553 GO TO 120 |
|
|
554 C----------------------------------------------------------------------- |
|
|
555 C FNU=HALF ODD INTEGER CASE, DNU=-0.5 |
|
|
556 C----------------------------------------------------------------------- |
|
|
557 300 CONTINUE |
|
|
558 S1R = COEFR |
|
|
559 S1I = COEFI |
|
|
560 S2R = COEFR |
|
|
561 S2I = COEFI |
|
|
562 GO TO 210 |
|
|
563 C |
|
|
564 C |
|
|
565 310 CONTINUE |
|
|
566 NZ=-2 |
|
|
567 RETURN |
|
|
568 END |
