|
3217
|
1 SUBROUTINE ZUNI2(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL, |
|
|
2 * TOL, ELIM, ALIM) |
|
|
3 C***BEGIN PROLOGUE ZUNI2 |
|
|
4 C***REFER TO ZBESI,ZBESK |
|
|
5 C |
|
|
6 C ZUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF |
|
|
7 C UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I |
|
|
8 C OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO. |
|
|
9 C |
|
|
10 C FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC |
|
|
11 C EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET. |
|
|
12 C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER |
|
|
13 C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL. |
|
|
14 C Y(I)=CZERO FOR I=NLAST+1,N |
|
|
15 C |
|
|
16 C***ROUTINES CALLED ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,XZABS |
|
|
17 C***END PROLOGUE ZUNI2 |
|
|
18 C COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS, |
|
|
19 C *CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN |
|
|
20 DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGI, |
|
|
21 * ARGR, ASCLE, ASUMI, ASUMR, BRY, BSUMI, BSUMR, CIDI, CIPI, CIPR, |
|
|
22 * CONER, CRSC, CSCL, CSRR, CSSR, C1R, C2I, C2M, C2R, DAII, |
|
|
23 * DAIR, ELIM, FN, FNU, FNUL, HPI, PHII, PHIR, RAST, RAZ, RS1, RZI, |
|
|
24 * RZR, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZBI, ZBR, ZEROI, |
|
|
25 * ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI, ZNI, ZNR, ZR, CYR, |
|
|
26 * CYI, D1MACH, XZABS, CAR, SAR |
|
|
27 INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST, |
|
|
28 * NN, NUF, NW, NZ, IDUM |
|
|
29 DIMENSION BRY(3), YR(N), YI(N), CIPR(4), CIPI(4), CSSR(3), |
|
|
30 * CSRR(3), CYR(2), CYI(2) |
|
|
31 DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 / |
|
|
32 DATA CIPR(1),CIPI(1),CIPR(2),CIPI(2),CIPR(3),CIPI(3),CIPR(4), |
|
|
33 * CIPI(4)/ 1.0D0,0.0D0, 0.0D0,1.0D0, -1.0D0,0.0D0, 0.0D0,-1.0D0/ |
|
|
34 DATA HPI, AIC / |
|
|
35 1 1.57079632679489662D+00, 1.265512123484645396D+00/ |
|
|
36 C |
|
|
37 NZ = 0 |
|
|
38 ND = N |
|
|
39 NLAST = 0 |
|
|
40 C----------------------------------------------------------------------- |
|
|
41 C COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG- |
|
|
42 C NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE, |
|
|
43 C EXP(ALIM)=EXP(ELIM)*TOL |
|
|
44 C----------------------------------------------------------------------- |
|
|
45 CSCL = 1.0D0/TOL |
|
|
46 CRSC = TOL |
|
|
47 CSSR(1) = CSCL |
|
|
48 CSSR(2) = CONER |
|
|
49 CSSR(3) = CRSC |
|
|
50 CSRR(1) = CRSC |
|
|
51 CSRR(2) = CONER |
|
|
52 CSRR(3) = CSCL |
|
|
53 BRY(1) = 1.0D+3*D1MACH(1)/TOL |
|
|
54 C----------------------------------------------------------------------- |
|
|
55 C ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI |
|
|
56 C----------------------------------------------------------------------- |
|
|
57 ZNR = ZI |
|
|
58 ZNI = -ZR |
|
|
59 ZBR = ZR |
|
|
60 ZBI = ZI |
|
|
61 CIDI = -CONER |
|
|
62 INU = INT(SNGL(FNU)) |
|
|
63 ANG = HPI*(FNU-DBLE(FLOAT(INU))) |
|
|
64 C2R = DCOS(ANG) |
|
|
65 C2I = DSIN(ANG) |
|
|
66 CAR = C2R |
|
|
67 SAR = C2I |
|
|
68 IN = INU + N - 1 |
|
|
69 IN = MOD(IN,4) + 1 |
|
|
70 STR = C2R*CIPR(IN) - C2I*CIPI(IN) |
|
|
71 C2I = C2R*CIPI(IN) + C2I*CIPR(IN) |
|
|
72 C2R = STR |
|
|
73 IF (ZI.GT.0.0D0) GO TO 10 |
|
|
74 ZNR = -ZNR |
|
|
75 ZBI = -ZBI |
|
|
76 CIDI = -CIDI |
|
|
77 C2I = -C2I |
|
|
78 10 CONTINUE |
|
|
79 C----------------------------------------------------------------------- |
|
|
80 C CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER |
|
|
81 C----------------------------------------------------------------------- |
|
|
82 FN = DMAX1(FNU,1.0D0) |
|
|
83 CALL ZUNHJ(ZNR, ZNI, FN, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, |
|
|
84 * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) |
|
|
85 IF (KODE.EQ.1) GO TO 20 |
|
|
86 STR = ZBR + ZETA2R |
|
|
87 STI = ZBI + ZETA2I |
|
|
88 RAST = FN/XZABS(STR,STI) |
|
|
89 STR = STR*RAST*RAST |
|
|
90 STI = -STI*RAST*RAST |
|
|
91 S1R = -ZETA1R + STR |
|
|
92 S1I = -ZETA1I + STI |
|
|
93 GO TO 30 |
|
|
94 20 CONTINUE |
|
|
95 S1R = -ZETA1R + ZETA2R |
|
|
96 S1I = -ZETA1I + ZETA2I |
|
|
97 30 CONTINUE |
|
|
98 RS1 = S1R |
|
|
99 IF (DABS(RS1).GT.ELIM) GO TO 150 |
|
|
100 40 CONTINUE |
|
|
101 NN = MIN0(2,ND) |
|
|
102 DO 90 I=1,NN |
|
|
103 FN = FNU + DBLE(FLOAT(ND-I)) |
|
|
104 CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIR, PHII, ARGR, ARGI, |
|
|
105 * ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) |
|
|
106 IF (KODE.EQ.1) GO TO 50 |
|
|
107 STR = ZBR + ZETA2R |
|
|
108 STI = ZBI + ZETA2I |
|
|
109 RAST = FN/XZABS(STR,STI) |
|
|
110 STR = STR*RAST*RAST |
|
|
111 STI = -STI*RAST*RAST |
|
|
112 S1R = -ZETA1R + STR |
|
|
113 S1I = -ZETA1I + STI + DABS(ZI) |
|
|
114 GO TO 60 |
|
|
115 50 CONTINUE |
|
|
116 S1R = -ZETA1R + ZETA2R |
|
|
117 S1I = -ZETA1I + ZETA2I |
|
|
118 60 CONTINUE |
|
|
119 C----------------------------------------------------------------------- |
|
|
120 C TEST FOR UNDERFLOW AND OVERFLOW |
|
|
121 C----------------------------------------------------------------------- |
|
|
122 RS1 = S1R |
|
|
123 IF (DABS(RS1).GT.ELIM) GO TO 120 |
|
|
124 IF (I.EQ.1) IFLAG = 2 |
|
|
125 IF (DABS(RS1).LT.ALIM) GO TO 70 |
|
|
126 C----------------------------------------------------------------------- |
|
|
127 C REFINE TEST AND SCALE |
|
|
128 C----------------------------------------------------------------------- |
|
|
129 C----------------------------------------------------------------------- |
|
|
130 APHI = XZABS(PHIR,PHII) |
|
|
131 AARG = XZABS(ARGR,ARGI) |
|
|
132 RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC |
|
|
133 IF (DABS(RS1).GT.ELIM) GO TO 120 |
|
|
134 IF (I.EQ.1) IFLAG = 1 |
|
|
135 IF (RS1.LT.0.0D0) GO TO 70 |
|
|
136 IF (I.EQ.1) IFLAG = 3 |
|
|
137 70 CONTINUE |
|
|
138 C----------------------------------------------------------------------- |
|
|
139 C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR |
|
|
140 C EXPONENT EXTREMES |
|
|
141 C----------------------------------------------------------------------- |
|
|
142 CALL ZAIRY(ARGR, ARGI, 0, 2, AIR, AII, NAI, IDUM) |
|
|
143 CALL ZAIRY(ARGR, ARGI, 1, 2, DAIR, DAII, NDAI, IDUM) |
|
|
144 STR = DAIR*BSUMR - DAII*BSUMI |
|
|
145 STI = DAIR*BSUMI + DAII*BSUMR |
|
|
146 STR = STR + (AIR*ASUMR-AII*ASUMI) |
|
|
147 STI = STI + (AIR*ASUMI+AII*ASUMR) |
|
|
148 S2R = PHIR*STR - PHII*STI |
|
|
149 S2I = PHIR*STI + PHII*STR |
|
|
150 STR = DEXP(S1R)*CSSR(IFLAG) |
|
|
151 S1R = STR*DCOS(S1I) |
|
|
152 S1I = STR*DSIN(S1I) |
|
|
153 STR = S2R*S1R - S2I*S1I |
|
|
154 S2I = S2R*S1I + S2I*S1R |
|
|
155 S2R = STR |
|
|
156 IF (IFLAG.NE.1) GO TO 80 |
|
|
157 CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL) |
|
|
158 IF (NW.NE.0) GO TO 120 |
|
|
159 80 CONTINUE |
|
|
160 IF (ZI.LE.0.0D0) S2I = -S2I |
|
|
161 STR = S2R*C2R - S2I*C2I |
|
|
162 S2I = S2R*C2I + S2I*C2R |
|
|
163 S2R = STR |
|
|
164 CYR(I) = S2R |
|
|
165 CYI(I) = S2I |
|
|
166 J = ND - I + 1 |
|
|
167 YR(J) = S2R*CSRR(IFLAG) |
|
|
168 YI(J) = S2I*CSRR(IFLAG) |
|
|
169 STR = -C2I*CIDI |
|
|
170 C2I = C2R*CIDI |
|
|
171 C2R = STR |
|
|
172 90 CONTINUE |
|
|
173 IF (ND.LE.2) GO TO 110 |
|
|
174 RAZ = 1.0D0/XZABS(ZR,ZI) |
|
|
175 STR = ZR*RAZ |
|
|
176 STI = -ZI*RAZ |
|
|
177 RZR = (STR+STR)*RAZ |
|
|
178 RZI = (STI+STI)*RAZ |
|
|
179 BRY(2) = 1.0D0/BRY(1) |
|
|
180 BRY(3) = D1MACH(2) |
|
|
181 S1R = CYR(1) |
|
|
182 S1I = CYI(1) |
|
|
183 S2R = CYR(2) |
|
|
184 S2I = CYI(2) |
|
|
185 C1R = CSRR(IFLAG) |
|
|
186 ASCLE = BRY(IFLAG) |
|
|
187 K = ND - 2 |
|
|
188 FN = DBLE(FLOAT(K)) |
|
|
189 DO 100 I=3,ND |
|
|
190 C2R = S2R |
|
|
191 C2I = S2I |
|
|
192 S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I) |
|
|
193 S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R) |
|
|
194 S1R = C2R |
|
|
195 S1I = C2I |
|
|
196 C2R = S2R*C1R |
|
|
197 C2I = S2I*C1R |
|
|
198 YR(K) = C2R |
|
|
199 YI(K) = C2I |
|
|
200 K = K - 1 |
|
|
201 FN = FN - 1.0D0 |
|
|
202 IF (IFLAG.GE.3) GO TO 100 |
|
|
203 STR = DABS(C2R) |
|
|
204 STI = DABS(C2I) |
|
|
205 C2M = DMAX1(STR,STI) |
|
|
206 IF (C2M.LE.ASCLE) GO TO 100 |
|
|
207 IFLAG = IFLAG + 1 |
|
|
208 ASCLE = BRY(IFLAG) |
|
|
209 S1R = S1R*C1R |
|
|
210 S1I = S1I*C1R |
|
|
211 S2R = C2R |
|
|
212 S2I = C2I |
|
|
213 S1R = S1R*CSSR(IFLAG) |
|
|
214 S1I = S1I*CSSR(IFLAG) |
|
|
215 S2R = S2R*CSSR(IFLAG) |
|
|
216 S2I = S2I*CSSR(IFLAG) |
|
|
217 C1R = CSRR(IFLAG) |
|
|
218 100 CONTINUE |
|
|
219 110 CONTINUE |
|
|
220 RETURN |
|
|
221 120 CONTINUE |
|
|
222 IF (RS1.GT.0.0D0) GO TO 140 |
|
|
223 C----------------------------------------------------------------------- |
|
|
224 C SET UNDERFLOW AND UPDATE PARAMETERS |
|
|
225 C----------------------------------------------------------------------- |
|
|
226 YR(ND) = ZEROR |
|
|
227 YI(ND) = ZEROI |
|
|
228 NZ = NZ + 1 |
|
|
229 ND = ND - 1 |
|
|
230 IF (ND.EQ.0) GO TO 110 |
|
|
231 CALL ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, NUF, TOL, ELIM, ALIM) |
|
|
232 IF (NUF.LT.0) GO TO 140 |
|
|
233 ND = ND - NUF |
|
|
234 NZ = NZ + NUF |
|
|
235 IF (ND.EQ.0) GO TO 110 |
|
|
236 FN = FNU + DBLE(FLOAT(ND-1)) |
|
|
237 IF (FN.LT.FNUL) GO TO 130 |
|
|
238 C FN = CIDI |
|
|
239 C J = NUF + 1 |
|
|
240 C K = MOD(J,4) + 1 |
|
|
241 C S1R = CIPR(K) |
|
|
242 C S1I = CIPI(K) |
|
|
243 C IF (FN.LT.0.0D0) S1I = -S1I |
|
|
244 C STR = C2R*S1R - C2I*S1I |
|
|
245 C C2I = C2R*S1I + C2I*S1R |
|
|
246 C C2R = STR |
|
|
247 IN = INU + ND - 1 |
|
|
248 IN = MOD(IN,4) + 1 |
|
|
249 C2R = CAR*CIPR(IN) - SAR*CIPI(IN) |
|
|
250 C2I = CAR*CIPI(IN) + SAR*CIPR(IN) |
|
|
251 IF (ZI.LE.0.0D0) C2I = -C2I |
|
|
252 GO TO 40 |
|
|
253 130 CONTINUE |
|
|
254 NLAST = ND |
|
|
255 RETURN |
|
|
256 140 CONTINUE |
|
|
257 NZ = -1 |
|
|
258 RETURN |
|
|
259 150 CONTINUE |
|
|
260 IF (RS1.GT.0.0D0) GO TO 140 |
|
|
261 NZ = N |
|
|
262 DO 160 I=1,N |
|
|
263 YR(I) = ZEROR |
|
|
264 YI(I) = ZEROI |
|
|
265 160 CONTINUE |
|
|
266 RETURN |
|
|
267 END |
