1 SUBROUTINE ZUNK1(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM, 2 * ALIM) 3C***BEGIN PROLOGUE ZUNK1 4C***REFER TO ZBESK 5C 6C ZUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE 7C RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE 8C UNIFORM ASYMPTOTIC EXPANSION. 9C MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION. 10C NZ=-1 MEANS AN OVERFLOW WILL OCCUR 11C 12C***ROUTINES CALLED ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,ZABS 13C***END PROLOGUE ZUNK1 14C COMPLEX CFN,CK,CONE,CRSC,CS,CSCL,CSGN,CSPN,CSR,CSS,CWRK,CY,CZERO, 15C *C1,C2,PHI,PHID,RZ,SUM,SUMD,S1,S2,Y,Z,ZETA1,ZETA1D,ZETA2,ZETA2D,ZR 16 DOUBLE PRECISION ALIM, ANG, APHI, ASC, ASCLE, BRY, CKI, CKR, 17 * CONER, CRSC, CSCL, CSGNI, CSPNI, CSPNR, CSR, CSRR, CSSR, 18 * CWRKI, CWRKR, CYI, CYR, C1I, C1R, C2I, C2M, C2R, ELIM, FMR, FN, 19 * FNF, FNU, PHIDI, PHIDR, PHII, PHIR, PI, RAST, RAZR, RS1, RZI, 20 * RZR, SGN, STI, STR, SUMDI, SUMDR, SUMI, SUMR, S1I, S1R, S2I, 21 * S2R, TOL, YI, YR, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, 22 * ZET1DI, ZET1DR, ZET2DI, ZET2DR, ZI, ZR, ZRI, ZRR, D1MACH, ZABS 23 INTEGER I, IB, IFLAG, IFN, IL, INIT, INU, IUF, K, KDFLG, KFLAG, 24 * KK, KODE, MR, N, NW, NZ, INITD, IC, IPARD, J 25 DIMENSION BRY(3), INIT(2), YR(N), YI(N), SUMR(2), SUMI(2), 26 * ZETA1R(2), ZETA1I(2), ZETA2R(2), ZETA2I(2), CYR(2), CYI(2), 27 * CWRKR(16,3), CWRKI(16,3), CSSR(3), CSRR(3), PHIR(2), PHII(2) 28 DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 / 29 DATA PI / 3.14159265358979324D0 / 30C 31 KDFLG = 1 32 NZ = 0 33C----------------------------------------------------------------------- 34C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN 35C THE UNDERFLOW LIMIT 36C----------------------------------------------------------------------- 37 CSCL = 1.0D0/TOL 38 CRSC = TOL 39 CSSR(1) = CSCL 40 CSSR(2) = CONER 41 CSSR(3) = CRSC 42 CSRR(1) = CRSC 43 CSRR(2) = CONER 44 CSRR(3) = CSCL 45 BRY(1) = 1.0D+3*D1MACH(1)/TOL 46 BRY(2) = 1.0D0/BRY(1) 47 BRY(3) = D1MACH(2) 48 ZRR = ZR 49 ZRI = ZI 50 IF (ZR.GE.0.0D0) GO TO 10 51 ZRR = -ZR 52 ZRI = -ZI 53 10 CONTINUE 54 J = 2 55 DO 70 I=1,N 56C----------------------------------------------------------------------- 57C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J 58C----------------------------------------------------------------------- 59 J = 3 - J 60 FN = FNU + DBLE(FLOAT(I-1)) 61 INIT(J) = 0 62 CALL ZUNIK(ZRR, ZRI, FN, 2, 0, TOL, INIT(J), PHIR(J), PHII(J), 63 * ZETA1R(J), ZETA1I(J), ZETA2R(J), ZETA2I(J), SUMR(J), SUMI(J), 64 * CWRKR(1,J), CWRKI(1,J)) 65 IF (KODE.EQ.1) GO TO 20 66 STR = ZRR + ZETA2R(J) 67 STI = ZRI + ZETA2I(J) 68 RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0))) 69 STR = STR*RAST*RAST 70 STI = -STI*RAST*RAST 71 S1R = ZETA1R(J) - STR 72 S1I = ZETA1I(J) - STI 73 GO TO 30 74 20 CONTINUE 75 S1R = ZETA1R(J) - ZETA2R(J) 76 S1I = ZETA1I(J) - ZETA2I(J) 77 30 CONTINUE 78 RS1 = S1R 79C----------------------------------------------------------------------- 80C TEST FOR UNDERFLOW AND OVERFLOW 81C----------------------------------------------------------------------- 82 IF (DABS(RS1).GT.ELIM) GO TO 60 83 IF (KDFLG.EQ.1) KFLAG = 2 84 IF (DABS(RS1).LT.ALIM) GO TO 40 85C----------------------------------------------------------------------- 86C REFINE TEST AND SCALE 87C----------------------------------------------------------------------- 88 APHI = ZABS(CMPLX(PHIR(J),PHII(J),kind=KIND(1.0D0))) 89 RS1 = RS1 + DLOG(APHI) 90 IF (DABS(RS1).GT.ELIM) GO TO 60 91 IF (KDFLG.EQ.1) KFLAG = 1 92 IF (RS1.LT.0.0D0) GO TO 40 93 IF (KDFLG.EQ.1) KFLAG = 3 94 40 CONTINUE 95C----------------------------------------------------------------------- 96C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR 97C EXPONENT EXTREMES 98C----------------------------------------------------------------------- 99 S2R = PHIR(J)*SUMR(J) - PHII(J)*SUMI(J) 100 S2I = PHIR(J)*SUMI(J) + PHII(J)*SUMR(J) 101 STR = DEXP(S1R)*CSSR(KFLAG) 102 S1R = STR*DCOS(S1I) 103 S1I = STR*DSIN(S1I) 104 STR = S2R*S1R - S2I*S1I 105 S2I = S1R*S2I + S2R*S1I 106 S2R = STR 107 IF (KFLAG.NE.1) GO TO 50 108 CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL) 109 IF (NW.NE.0) GO TO 60 110 50 CONTINUE 111 CYR(KDFLG) = S2R 112 CYI(KDFLG) = S2I 113 YR(I) = S2R*CSRR(KFLAG) 114 YI(I) = S2I*CSRR(KFLAG) 115 IF (KDFLG.EQ.2) GO TO 75 116 KDFLG = 2 117 GO TO 70 118 60 CONTINUE 119 IF (RS1.GT.0.0D0) GO TO 300 120C----------------------------------------------------------------------- 121C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW 122C----------------------------------------------------------------------- 123 IF (ZR.LT.0.0D0) GO TO 300 124 KDFLG = 1 125 YR(I)=ZEROR 126 YI(I)=ZEROI 127 NZ=NZ+1 128 IF (I.EQ.1) GO TO 70 129 IF ((YR(I-1).EQ.ZEROR).AND.(YI(I-1).EQ.ZEROI)) GO TO 70 130 YR(I-1)=ZEROR 131 YI(I-1)=ZEROI 132 NZ=NZ+1 133 70 CONTINUE 134 I = N 135 75 CONTINUE 136 RAZR = 1.0D0/ZABS(CMPLX(ZRR,ZRI,kind=KIND(1.0D0))) 137 STR = ZRR*RAZR 138 STI = -ZRI*RAZR 139 RZR = (STR+STR)*RAZR 140 RZI = (STI+STI)*RAZR 141 CKR = FN*RZR 142 CKI = FN*RZI 143 IB = I + 1 144 IF (N.LT.IB) GO TO 160 145C----------------------------------------------------------------------- 146C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO 147C ON UNDERFLOW. 148C----------------------------------------------------------------------- 149 FN = FNU + DBLE(FLOAT(N-1)) 150 IPARD = 1 151 IF (MR.NE.0) IPARD = 0 152 INITD = 0 153 CALL ZUNIK(ZRR, ZRI, FN, 2, IPARD, TOL, INITD, PHIDR, PHIDI, 154 * ZET1DR, ZET1DI, ZET2DR, ZET2DI, SUMDR, SUMDI, CWRKR(1,3), 155 * CWRKI(1,3)) 156 IF (KODE.EQ.1) GO TO 80 157 STR = ZRR + ZET2DR 158 STI = ZRI + ZET2DI 159 RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0))) 160 STR = STR*RAST*RAST 161 STI = -STI*RAST*RAST 162 S1R = ZET1DR - STR 163 S1I = ZET1DI - STI 164 GO TO 90 165 80 CONTINUE 166 S1R = ZET1DR - ZET2DR 167 S1I = ZET1DI - ZET2DI 168 90 CONTINUE 169 RS1 = S1R 170 IF (DABS(RS1).GT.ELIM) GO TO 95 171 IF (DABS(RS1).LT.ALIM) GO TO 100 172C---------------------------------------------------------------------------- 173C REFINE ESTIMATE AND TEST 174C------------------------------------------------------------------------- 175 APHI = ZABS(CMPLX(PHIDR,PHIDI,kind=KIND(1.0D0))) 176 RS1 = RS1+DLOG(APHI) 177 IF (DABS(RS1).LT.ELIM) GO TO 100 178 95 CONTINUE 179 IF (DABS(RS1).GT.0.0D0) GO TO 300 180C----------------------------------------------------------------------- 181C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW 182C----------------------------------------------------------------------- 183 IF (ZR.LT.0.0D0) GO TO 300 184 NZ = N 185 DO 96 I=1,N 186 YR(I) = ZEROR 187 YI(I) = ZEROI 188 96 CONTINUE 189 RETURN 190C--------------------------------------------------------------------------- 191C FORWARD RECUR FOR REMAINDER OF THE SEQUENCE 192C---------------------------------------------------------------------------- 193 100 CONTINUE 194 S1R = CYR(1) 195 S1I = CYI(1) 196 S2R = CYR(2) 197 S2I = CYI(2) 198 C1R = CSRR(KFLAG) 199 ASCLE = BRY(KFLAG) 200 DO 120 I=IB,N 201 C2R = S2R 202 C2I = S2I 203 S2R = CKR*C2R - CKI*C2I + S1R 204 S2I = CKR*C2I + CKI*C2R + S1I 205 S1R = C2R 206 S1I = C2I 207 CKR = CKR + RZR 208 CKI = CKI + RZI 209 C2R = S2R*C1R 210 C2I = S2I*C1R 211 YR(I) = C2R 212 YI(I) = C2I 213 IF (KFLAG.GE.3) GO TO 120 214 STR = DABS(C2R) 215 STI = DABS(C2I) 216 C2M = DMAX1(STR,STI) 217 IF (C2M.LE.ASCLE) GO TO 120 218 KFLAG = KFLAG + 1 219 ASCLE = BRY(KFLAG) 220 S1R = S1R*C1R 221 S1I = S1I*C1R 222 S2R = C2R 223 S2I = C2I 224 S1R = S1R*CSSR(KFLAG) 225 S1I = S1I*CSSR(KFLAG) 226 S2R = S2R*CSSR(KFLAG) 227 S2I = S2I*CSSR(KFLAG) 228 C1R = CSRR(KFLAG) 229 120 CONTINUE 230 160 CONTINUE 231 IF (MR.EQ.0) RETURN 232C----------------------------------------------------------------------- 233C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0 234C----------------------------------------------------------------------- 235 NZ = 0 236 FMR = DBLE(FLOAT(MR)) 237 SGN = -DSIGN(PI,FMR) 238C----------------------------------------------------------------------- 239C CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP. 240C----------------------------------------------------------------------- 241 CSGNI = SGN 242 INU = INT(SNGL(FNU)) 243 FNF = FNU - DBLE(FLOAT(INU)) 244 IFN = INU + N - 1 245 ANG = FNF*SGN 246 CSPNR = DCOS(ANG) 247 CSPNI = DSIN(ANG) 248 IF (MOD(IFN,2).EQ.0) GO TO 170 249 CSPNR = -CSPNR 250 CSPNI = -CSPNI 251 170 CONTINUE 252 ASC = BRY(1) 253 IUF = 0 254 KK = N 255 KDFLG = 1 256 IB = IB - 1 257 IC = IB - 1 258 DO 270 K=1,N 259 FN = FNU + DBLE(FLOAT(KK-1)) 260C----------------------------------------------------------------------- 261C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K 262C FUNCTION ABOVE 263C----------------------------------------------------------------------- 264 M=3 265 IF (N.GT.2) GO TO 175 266 172 CONTINUE 267 INITD = INIT(J) 268 PHIDR = PHIR(J) 269 PHIDI = PHII(J) 270 ZET1DR = ZETA1R(J) 271 ZET1DI = ZETA1I(J) 272 ZET2DR = ZETA2R(J) 273 ZET2DI = ZETA2I(J) 274 SUMDR = SUMR(J) 275 SUMDI = SUMI(J) 276 M = J 277 J = 3 - J 278 GO TO 180 279 175 CONTINUE 280 IF ((KK.EQ.N).AND.(IB.LT.N)) GO TO 180 281 IF ((KK.EQ.IB).OR.(KK.EQ.IC)) GO TO 172 282 INITD = 0 283 180 CONTINUE 284 CALL ZUNIK(ZRR, ZRI, FN, 1, 0, TOL, INITD, PHIDR, PHIDI, 285 * ZET1DR, ZET1DI, ZET2DR, ZET2DI, SUMDR, SUMDI, 286 * CWRKR(1,M), CWRKI(1,M)) 287 IF (KODE.EQ.1) GO TO 200 288 STR = ZRR + ZET2DR 289 STI = ZRI + ZET2DI 290 RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0))) 291 STR = STR*RAST*RAST 292 STI = -STI*RAST*RAST 293 S1R = -ZET1DR + STR 294 S1I = -ZET1DI + STI 295 GO TO 210 296 200 CONTINUE 297 S1R = -ZET1DR + ZET2DR 298 S1I = -ZET1DI + ZET2DI 299 210 CONTINUE 300C----------------------------------------------------------------------- 301C TEST FOR UNDERFLOW AND OVERFLOW 302C----------------------------------------------------------------------- 303 RS1 = S1R 304 IF (DABS(RS1).GT.ELIM) GO TO 260 305 IF (KDFLG.EQ.1) IFLAG = 2 306 IF (DABS(RS1).LT.ALIM) GO TO 220 307C----------------------------------------------------------------------- 308C REFINE TEST AND SCALE 309C----------------------------------------------------------------------- 310 APHI = ZABS(CMPLX(PHIDR,PHIDI,kind=KIND(1.0D0))) 311 RS1 = RS1 + DLOG(APHI) 312 IF (DABS(RS1).GT.ELIM) GO TO 260 313 IF (KDFLG.EQ.1) IFLAG = 1 314 IF (RS1.LT.0.0D0) GO TO 220 315 IF (KDFLG.EQ.1) IFLAG = 3 316 220 CONTINUE 317 STR = PHIDR*SUMDR - PHIDI*SUMDI 318 STI = PHIDR*SUMDI + PHIDI*SUMDR 319 S2R = -CSGNI*STI 320 S2I = CSGNI*STR 321 STR = DEXP(S1R)*CSSR(IFLAG) 322 S1R = STR*DCOS(S1I) 323 S1I = STR*DSIN(S1I) 324 STR = S2R*S1R - S2I*S1I 325 S2I = S2R*S1I + S2I*S1R 326 S2R = STR 327 IF (IFLAG.NE.1) GO TO 230 328 CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL) 329 IF (NW.EQ.0) GO TO 230 330 S2R = ZEROR 331 S2I = ZEROI 332 230 CONTINUE 333 CYR(KDFLG) = S2R 334 CYI(KDFLG) = S2I 335 C2R = S2R 336 C2I = S2I 337 S2R = S2R*CSRR(IFLAG) 338 S2I = S2I*CSRR(IFLAG) 339C----------------------------------------------------------------------- 340C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N 341C----------------------------------------------------------------------- 342 S1R = YR(KK) 343 S1I = YI(KK) 344 IF (KODE.EQ.1) GO TO 250 345 CALL ZS1S2(ZRR, ZRI, S1R, S1I, S2R, S2I, NW, ASC, ALIM, IUF) 346 NZ = NZ + NW 347 250 CONTINUE 348 YR(KK) = S1R*CSPNR - S1I*CSPNI + S2R 349 YI(KK) = CSPNR*S1I + CSPNI*S1R + S2I 350 KK = KK - 1 351 CSPNR = -CSPNR 352 CSPNI = -CSPNI 353 IF (C2R.NE.0.0D0 .OR. C2I.NE.0.0D0) GO TO 255 354 KDFLG = 1 355 GO TO 270 356 255 CONTINUE 357 IF (KDFLG.EQ.2) GO TO 275 358 KDFLG = 2 359 GO TO 270 360 260 CONTINUE 361 IF (RS1.GT.0.0D0) GO TO 300 362 S2R = ZEROR 363 S2I = ZEROI 364 GO TO 230 365 270 CONTINUE 366 K = N 367 275 CONTINUE 368 IL = N - K 369 IF (IL.EQ.0) RETURN 370C----------------------------------------------------------------------- 371C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE 372C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP 373C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES. 374C----------------------------------------------------------------------- 375 S1R = CYR(1) 376 S1I = CYI(1) 377 S2R = CYR(2) 378 S2I = CYI(2) 379 CSR = CSRR(IFLAG) 380 ASCLE = BRY(IFLAG) 381 FN = DBLE(FLOAT(INU+IL)) 382 DO 290 I=1,IL 383 C2R = S2R 384 C2I = S2I 385 S2R = S1R + (FN+FNF)*(RZR*C2R-RZI*C2I) 386 S2I = S1I + (FN+FNF)*(RZR*C2I+RZI*C2R) 387 S1R = C2R 388 S1I = C2I 389 FN = FN - 1.0D0 390 C2R = S2R*CSR 391 C2I = S2I*CSR 392 CKR = C2R 393 CKI = C2I 394 C1R = YR(KK) 395 C1I = YI(KK) 396 IF (KODE.EQ.1) GO TO 280 397 CALL ZS1S2(ZRR, ZRI, C1R, C1I, C2R, C2I, NW, ASC, ALIM, IUF) 398 NZ = NZ + NW 399 280 CONTINUE 400 YR(KK) = C1R*CSPNR - C1I*CSPNI + C2R 401 YI(KK) = C1R*CSPNI + C1I*CSPNR + C2I 402 KK = KK - 1 403 CSPNR = -CSPNR 404 CSPNI = -CSPNI 405 IF (IFLAG.GE.3) GO TO 290 406 C2R = DABS(CKR) 407 C2I = DABS(CKI) 408 C2M = DMAX1(C2R,C2I) 409 IF (C2M.LE.ASCLE) GO TO 290 410 IFLAG = IFLAG + 1 411 ASCLE = BRY(IFLAG) 412 S1R = S1R*CSR 413 S1I = S1I*CSR 414 S2R = CKR 415 S2I = CKI 416 S1R = S1R*CSSR(IFLAG) 417 S1I = S1I*CSSR(IFLAG) 418 S2R = S2R*CSSR(IFLAG) 419 S2I = S2I*CSSR(IFLAG) 420 CSR = CSRR(IFLAG) 421 290 CONTINUE 422 RETURN 423 300 CONTINUE 424 NZ = -1 425 RETURN 426 END 427