00001 SUBROUTINE ZBKNU(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM,
00002 * ALIM)
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 DOUBLE PRECISION AA, AK, ALIM, ASCLE, A1, A2, BB, BK, BRY, CAZ,
00013 * CBI, CBR, CC, CCHI, CCHR, CKI, CKR, COEFI, COEFR, CONEI, CONER,
00014 * CRSCR, CSCLR, CSHI, CSHR, CSI, CSR, CSRR, CSSR, CTWOR,
00015 * CZEROI, CZEROR, CZI, CZR, DNU, DNU2, DPI, ELIM, ETEST, FC, FHS,
00016 * FI, FK, FKS, FMUI, FMUR, FNU, FPI, FR, G1, G2, HPI, PI, PR, PTI,
00017 * PTR, P1I, P1R, P2I, P2M, P2R, QI, QR, RAK, RCAZ, RTHPI, RZI,
00018 * RZR, R1, S, SMUI, SMUR, SPI, STI, STR, S1I, S1R, S2I, S2R, TM,
00019 * TOL, TTH, T1, T2, YI, YR, ZI, ZR, DGAMLN, D1MACH, XZABS, ELM,
00020 * CELMR, ZDR, ZDI, AS, ALAS, HELIM, CYR, CYI
00021 INTEGER I, IFLAG, INU, K, KFLAG, KK, KMAX, KODE, KODED, N, NZ,
00022 * IDUM, I1MACH, J, IC, INUB, NW
00023 DIMENSION YR(N), YI(N), CC(8), CSSR(3), CSRR(3), BRY(3), CYR(2),
00024 * CYI(2)
00025
00026
00027
00028 DATA KMAX / 30 /
00029 DATA CZEROR,CZEROI,CONER,CONEI,CTWOR,R1/
00030 1 0.0D0 , 0.0D0 , 1.0D0 , 0.0D0 , 2.0D0 , 2.0D0 /
00031 DATA DPI, RTHPI, SPI ,HPI, FPI, TTH /
00032 1 3.14159265358979324D0, 1.25331413731550025D0,
00033 2 1.90985931710274403D0, 1.57079632679489662D0,
00034 3 1.89769999331517738D0, 6.66666666666666666D-01/
00035 DATA CC(1), CC(2), CC(3), CC(4), CC(5), CC(6), CC(7), CC(8)/
00036 1 5.77215664901532861D-01, -4.20026350340952355D-02,
00037 2 -4.21977345555443367D-02, 7.21894324666309954D-03,
00038 3 -2.15241674114950973D-04, -2.01348547807882387D-05,
00039 4 1.13302723198169588D-06, 6.11609510448141582D-09/
00040
00041 CAZ = XZABS(ZR,ZI)
00042 CSCLR = 1.0D0/TOL
00043 CRSCR = TOL
00044 CSSR(1) = CSCLR
00045 CSSR(2) = 1.0D0
00046 CSSR(3) = CRSCR
00047 CSRR(1) = CRSCR
00048 CSRR(2) = 1.0D0
00049 CSRR(3) = CSCLR
00050 BRY(1) = 1.0D+3*D1MACH(1)/TOL
00051 BRY(2) = 1.0D0/BRY(1)
00052 BRY(3) = D1MACH(2)
00053 NZ = 0
00054 IFLAG = 0
00055 KODED = KODE
00056 RCAZ = 1.0D0/CAZ
00057 STR = ZR*RCAZ
00058 STI = -ZI*RCAZ
00059 RZR = (STR+STR)*RCAZ
00060 RZI = (STI+STI)*RCAZ
00061 INU = INT(SNGL(FNU+0.5D0))
00062 DNU = FNU - DBLE(FLOAT(INU))
00063 IF (DABS(DNU).EQ.0.5D0) GO TO 110
00064 DNU2 = 0.0D0
00065 IF (DABS(DNU).GT.TOL) DNU2 = DNU*DNU
00066 IF (CAZ.GT.R1) GO TO 110
00067
00068
00069
00070 FC = 1.0D0
00071 CALL XZLOG(RZR, RZI, SMUR, SMUI, IDUM)
00072 FMUR = SMUR*DNU
00073 FMUI = SMUI*DNU
00074 CALL ZSHCH(FMUR, FMUI, CSHR, CSHI, CCHR, CCHI)
00075 IF (DNU.EQ.0.0D0) GO TO 10
00076 FC = DNU*DPI
00077 FC = FC/DSIN(FC)
00078 SMUR = CSHR/DNU
00079 SMUI = CSHI/DNU
00080 10 CONTINUE
00081 A2 = 1.0D0 + DNU
00082
00083
00084
00085 T2 = DEXP(-DGAMLN(A2,IDUM))
00086 T1 = 1.0D0/(T2*FC)
00087 IF (DABS(DNU).GT.0.1D0) GO TO 40
00088
00089
00090
00091 AK = 1.0D0
00092 S = CC(1)
00093 DO 20 K=2,8
00094 AK = AK*DNU2
00095 TM = CC(K)*AK
00096 S = S + TM
00097 IF (DABS(TM).LT.TOL) GO TO 30
00098 20 CONTINUE
00099 30 G1 = -S
00100 GO TO 50
00101 40 CONTINUE
00102 G1 = (T1-T2)/(DNU+DNU)
00103 50 CONTINUE
00104 G2 = (T1+T2)*0.5D0
00105 FR = FC*(CCHR*G1+SMUR*G2)
00106 FI = FC*(CCHI*G1+SMUI*G2)
00107 CALL XZEXP(FMUR, FMUI, STR, STI)
00108 PR = 0.5D0*STR/T2
00109 PI = 0.5D0*STI/T2
00110 CALL ZDIV(0.5D0, 0.0D0, STR, STI, PTR, PTI)
00111 QR = PTR/T1
00112 QI = PTI/T1
00113 S1R = FR
00114 S1I = FI
00115 S2R = PR
00116 S2I = PI
00117 AK = 1.0D0
00118 A1 = 1.0D0
00119 CKR = CONER
00120 CKI = CONEI
00121 BK = 1.0D0 - DNU2
00122 IF (INU.GT.0 .OR. N.GT.1) GO TO 80
00123
00124
00125
00126 IF (CAZ.LT.TOL) GO TO 70
00127 CALL ZMLT(ZR, ZI, ZR, ZI, CZR, CZI)
00128 CZR = 0.25D0*CZR
00129 CZI = 0.25D0*CZI
00130 T1 = 0.25D0*CAZ*CAZ
00131 60 CONTINUE
00132 FR = (FR*AK+PR+QR)/BK
00133 FI = (FI*AK+PI+QI)/BK
00134 STR = 1.0D0/(AK-DNU)
00135 PR = PR*STR
00136 PI = PI*STR
00137 STR = 1.0D0/(AK+DNU)
00138 QR = QR*STR
00139 QI = QI*STR
00140 STR = CKR*CZR - CKI*CZI
00141 RAK = 1.0D0/AK
00142 CKI = (CKR*CZI+CKI*CZR)*RAK
00143 CKR = STR*RAK
00144 S1R = CKR*FR - CKI*FI + S1R
00145 S1I = CKR*FI + CKI*FR + S1I
00146 A1 = A1*T1*RAK
00147 BK = BK + AK + AK + 1.0D0
00148 AK = AK + 1.0D0
00149 IF (A1.GT.TOL) GO TO 60
00150 70 CONTINUE
00151 YR(1) = S1R
00152 YI(1) = S1I
00153 IF (KODED.EQ.1) RETURN
00154 CALL XZEXP(ZR, ZI, STR, STI)
00155 CALL ZMLT(S1R, S1I, STR, STI, YR(1), YI(1))
00156 RETURN
00157
00158
00159
00160 80 CONTINUE
00161 IF (CAZ.LT.TOL) GO TO 100
00162 CALL ZMLT(ZR, ZI, ZR, ZI, CZR, CZI)
00163 CZR = 0.25D0*CZR
00164 CZI = 0.25D0*CZI
00165 T1 = 0.25D0*CAZ*CAZ
00166 90 CONTINUE
00167 FR = (FR*AK+PR+QR)/BK
00168 FI = (FI*AK+PI+QI)/BK
00169 STR = 1.0D0/(AK-DNU)
00170 PR = PR*STR
00171 PI = PI*STR
00172 STR = 1.0D0/(AK+DNU)
00173 QR = QR*STR
00174 QI = QI*STR
00175 STR = CKR*CZR - CKI*CZI
00176 RAK = 1.0D0/AK
00177 CKI = (CKR*CZI+CKI*CZR)*RAK
00178 CKR = STR*RAK
00179 S1R = CKR*FR - CKI*FI + S1R
00180 S1I = CKR*FI + CKI*FR + S1I
00181 STR = PR - FR*AK
00182 STI = PI - FI*AK
00183 S2R = CKR*STR - CKI*STI + S2R
00184 S2I = CKR*STI + CKI*STR + S2I
00185 A1 = A1*T1*RAK
00186 BK = BK + AK + AK + 1.0D0
00187 AK = AK + 1.0D0
00188 IF (A1.GT.TOL) GO TO 90
00189 100 CONTINUE
00190 KFLAG = 2
00191 A1 = FNU + 1.0D0
00192 AK = A1*DABS(SMUR)
00193 IF (AK.GT.ALIM) KFLAG = 3
00194 STR = CSSR(KFLAG)
00195 P2R = S2R*STR
00196 P2I = S2I*STR
00197 CALL ZMLT(P2R, P2I, RZR, RZI, S2R, S2I)
00198 S1R = S1R*STR
00199 S1I = S1I*STR
00200 IF (KODED.EQ.1) GO TO 210
00201 CALL XZEXP(ZR, ZI, FR, FI)
00202 CALL ZMLT(S1R, S1I, FR, FI, S1R, S1I)
00203 CALL ZMLT(S2R, S2I, FR, FI, S2R, S2I)
00204 GO TO 210
00205
00206
00207
00208
00209
00210
00211 110 CONTINUE
00212 CALL XZSQRT(ZR, ZI, STR, STI)
00213 CALL ZDIV(RTHPI, CZEROI, STR, STI, COEFR, COEFI)
00214 KFLAG = 2
00215 IF (KODED.EQ.2) GO TO 120
00216 IF (ZR.GT.ALIM) GO TO 290
00217
00218 STR = DEXP(-ZR)*CSSR(KFLAG)
00219 STI = -STR*DSIN(ZI)
00220 STR = STR*DCOS(ZI)
00221 CALL ZMLT(COEFR, COEFI, STR, STI, COEFR, COEFI)
00222 120 CONTINUE
00223 IF (DABS(DNU).EQ.0.5D0) GO TO 300
00224
00225
00226
00227 AK = DCOS(DPI*DNU)
00228 AK = DABS(AK)
00229 IF (AK.EQ.CZEROR) GO TO 300
00230 FHS = DABS(0.25D0-DNU2)
00231 IF (FHS.EQ.CZEROR) GO TO 300
00232
00233
00234
00235
00236
00237
00238 T1 = DBLE(FLOAT(I1MACH(14)-1))
00239 T1 = T1*D1MACH(5)*3.321928094D0
00240 T1 = DMAX1(T1,12.0D0)
00241 T1 = DMIN1(T1,60.0D0)
00242 T2 = TTH*T1 - 6.0D0
00243 IF (ZR.NE.0.0D0) GO TO 130
00244 T1 = HPI
00245 GO TO 140
00246 130 CONTINUE
00247 T1 = DATAN(ZI/ZR)
00248 T1 = DABS(T1)
00249 140 CONTINUE
00250 IF (T2.GT.CAZ) GO TO 170
00251
00252
00253
00254 ETEST = AK/(DPI*CAZ*TOL)
00255 FK = CONER
00256 IF (ETEST.LT.CONER) GO TO 180
00257 FKS = CTWOR
00258 CKR = CAZ + CAZ + CTWOR
00259 P1R = CZEROR
00260 P2R = CONER
00261 DO 150 I=1,KMAX
00262 AK = FHS/FKS
00263 CBR = CKR/(FK+CONER)
00264 PTR = P2R
00265 P2R = CBR*P2R - P1R*AK
00266 P1R = PTR
00267 CKR = CKR + CTWOR
00268 FKS = FKS + FK + FK + CTWOR
00269 FHS = FHS + FK + FK
00270 FK = FK + CONER
00271 STR = DABS(P2R)*FK
00272 IF (ETEST.LT.STR) GO TO 160
00273 150 CONTINUE
00274 GO TO 310
00275 160 CONTINUE
00276 FK = FK + SPI*T1*DSQRT(T2/CAZ)
00277 FHS = DABS(0.25D0-DNU2)
00278 GO TO 180
00279 170 CONTINUE
00280
00281
00282
00283 A2 = DSQRT(CAZ)
00284 AK = FPI*AK/(TOL*DSQRT(A2))
00285 AA = 3.0D0*T1/(1.0D0+CAZ)
00286 BB = 14.7D0*T1/(28.0D0+CAZ)
00287 AK = (DLOG(AK)+CAZ*DCOS(AA)/(1.0D0+0.008D0*CAZ))/DCOS(BB)
00288 FK = 0.12125D0*AK*AK/CAZ + 1.5D0
00289 180 CONTINUE
00290
00291
00292
00293 K = INT(SNGL(FK))
00294 FK = DBLE(FLOAT(K))
00295 FKS = FK*FK
00296 P1R = CZEROR
00297 P1I = CZEROI
00298 P2R = TOL
00299 P2I = CZEROI
00300 CSR = P2R
00301 CSI = P2I
00302 DO 190 I=1,K
00303 A1 = FKS - FK
00304 AK = (FKS+FK)/(A1+FHS)
00305 RAK = 2.0D0/(FK+CONER)
00306 CBR = (FK+ZR)*RAK
00307 CBI = ZI*RAK
00308 PTR = P2R
00309 PTI = P2I
00310 P2R = (PTR*CBR-PTI*CBI-P1R)*AK
00311 P2I = (PTI*CBR+PTR*CBI-P1I)*AK
00312 P1R = PTR
00313 P1I = PTI
00314 CSR = CSR + P2R
00315 CSI = CSI + P2I
00316 FKS = A1 - FK + CONER
00317 FK = FK - CONER
00318 190 CONTINUE
00319
00320
00321
00322
00323 TM = XZABS(CSR,CSI)
00324 PTR = 1.0D0/TM
00325 S1R = P2R*PTR
00326 S1I = P2I*PTR
00327 CSR = CSR*PTR
00328 CSI = -CSI*PTR
00329 CALL ZMLT(COEFR, COEFI, S1R, S1I, STR, STI)
00330 CALL ZMLT(STR, STI, CSR, CSI, S1R, S1I)
00331 IF (INU.GT.0 .OR. N.GT.1) GO TO 200
00332 ZDR = ZR
00333 ZDI = ZI
00334 IF(IFLAG.EQ.1) GO TO 270
00335 GO TO 240
00336 200 CONTINUE
00337
00338
00339
00340 TM = XZABS(P2R,P2I)
00341 PTR = 1.0D0/TM
00342 P1R = P1R*PTR
00343 P1I = P1I*PTR
00344 P2R = P2R*PTR
00345 P2I = -P2I*PTR
00346 CALL ZMLT(P1R, P1I, P2R, P2I, PTR, PTI)
00347 STR = DNU + 0.5D0 - PTR
00348 STI = -PTI
00349 CALL ZDIV(STR, STI, ZR, ZI, STR, STI)
00350 STR = STR + 1.0D0
00351 CALL ZMLT(STR, STI, S1R, S1I, S2R, S2I)
00352
00353
00354
00355
00356 210 CONTINUE
00357 STR = DNU + 1.0D0
00358 CKR = STR*RZR
00359 CKI = STR*RZI
00360 IF (N.EQ.1) INU = INU - 1
00361 IF (INU.GT.0) GO TO 220
00362 IF (N.GT.1) GO TO 215
00363 S1R = S2R
00364 S1I = S2I
00365 215 CONTINUE
00366 ZDR = ZR
00367 ZDI = ZI
00368 IF(IFLAG.EQ.1) GO TO 270
00369 GO TO 240
00370 220 CONTINUE
00371 INUB = 1
00372 IF(IFLAG.EQ.1) GO TO 261
00373 225 CONTINUE
00374 P1R = CSRR(KFLAG)
00375 ASCLE = BRY(KFLAG)
00376 DO 230 I=INUB,INU
00377 STR = S2R
00378 STI = S2I
00379 S2R = CKR*STR - CKI*STI + S1R
00380 S2I = CKR*STI + CKI*STR + S1I
00381 S1R = STR
00382 S1I = STI
00383 CKR = CKR + RZR
00384 CKI = CKI + RZI
00385 IF (KFLAG.GE.3) GO TO 230
00386 P2R = S2R*P1R
00387 P2I = S2I*P1R
00388 STR = DABS(P2R)
00389 STI = DABS(P2I)
00390 P2M = DMAX1(STR,STI)
00391 IF (P2M.LE.ASCLE) GO TO 230
00392 KFLAG = KFLAG + 1
00393 ASCLE = BRY(KFLAG)
00394 S1R = S1R*P1R
00395 S1I = S1I*P1R
00396 S2R = P2R
00397 S2I = P2I
00398 STR = CSSR(KFLAG)
00399 S1R = S1R*STR
00400 S1I = S1I*STR
00401 S2R = S2R*STR
00402 S2I = S2I*STR
00403 P1R = CSRR(KFLAG)
00404 230 CONTINUE
00405 IF (N.NE.1) GO TO 240
00406 S1R = S2R
00407 S1I = S2I
00408 240 CONTINUE
00409 STR = CSRR(KFLAG)
00410 YR(1) = S1R*STR
00411 YI(1) = S1I*STR
00412 IF (N.EQ.1) RETURN
00413 YR(2) = S2R*STR
00414 YI(2) = S2I*STR
00415 IF (N.EQ.2) RETURN
00416 KK = 2
00417 250 CONTINUE
00418 KK = KK + 1
00419 IF (KK.GT.N) RETURN
00420 P1R = CSRR(KFLAG)
00421 ASCLE = BRY(KFLAG)
00422 DO 260 I=KK,N
00423 P2R = S2R
00424 P2I = S2I
00425 S2R = CKR*P2R - CKI*P2I + S1R
00426 S2I = CKI*P2R + CKR*P2I + S1I
00427 S1R = P2R
00428 S1I = P2I
00429 CKR = CKR + RZR
00430 CKI = CKI + RZI
00431 P2R = S2R*P1R
00432 P2I = S2I*P1R
00433 YR(I) = P2R
00434 YI(I) = P2I
00435 IF (KFLAG.GE.3) GO TO 260
00436 STR = DABS(P2R)
00437 STI = DABS(P2I)
00438 P2M = DMAX1(STR,STI)
00439 IF (P2M.LE.ASCLE) GO TO 260
00440 KFLAG = KFLAG + 1
00441 ASCLE = BRY(KFLAG)
00442 S1R = S1R*P1R
00443 S1I = S1I*P1R
00444 S2R = P2R
00445 S2I = P2I
00446 STR = CSSR(KFLAG)
00447 S1R = S1R*STR
00448 S1I = S1I*STR
00449 S2R = S2R*STR
00450 S2I = S2I*STR
00451 P1R = CSRR(KFLAG)
00452 260 CONTINUE
00453 RETURN
00454
00455
00456
00457 261 CONTINUE
00458 HELIM = 0.5D0*ELIM
00459 ELM = DEXP(-ELIM)
00460 CELMR = ELM
00461 ASCLE = BRY(1)
00462 ZDR = ZR
00463 ZDI = ZI
00464 IC = -1
00465 J = 2
00466 DO 262 I=1,INU
00467 STR = S2R
00468 STI = S2I
00469 S2R = STR*CKR-STI*CKI+S1R
00470 S2I = STI*CKR+STR*CKI+S1I
00471 S1R = STR
00472 S1I = STI
00473 CKR = CKR+RZR
00474 CKI = CKI+RZI
00475 AS = XZABS(S2R,S2I)
00476 ALAS = DLOG(AS)
00477 P2R = -ZDR+ALAS
00478 IF(P2R.LT.(-ELIM)) GO TO 263
00479 CALL XZLOG(S2R,S2I,STR,STI,IDUM)
00480 P2R = -ZDR+STR
00481 P2I = -ZDI+STI
00482 P2M = DEXP(P2R)/TOL
00483 P1R = P2M*DCOS(P2I)
00484 P1I = P2M*DSIN(P2I)
00485 CALL ZUCHK(P1R,P1I,NW,ASCLE,TOL)
00486 IF(NW.NE.0) GO TO 263
00487 J = 3 - J
00488 CYR(J) = P1R
00489 CYI(J) = P1I
00490 IF(IC.EQ.(I-1)) GO TO 264
00491 IC = I
00492 GO TO 262
00493 263 CONTINUE
00494 IF(ALAS.LT.HELIM) GO TO 262
00495 ZDR = ZDR-ELIM
00496 S1R = S1R*CELMR
00497 S1I = S1I*CELMR
00498 S2R = S2R*CELMR
00499 S2I = S2I*CELMR
00500 262 CONTINUE
00501 IF(N.NE.1) GO TO 270
00502 S1R = S2R
00503 S1I = S2I
00504 GO TO 270
00505 264 CONTINUE
00506 KFLAG = 1
00507 INUB = I+1
00508 S2R = CYR(J)
00509 S2I = CYI(J)
00510 J = 3 - J
00511 S1R = CYR(J)
00512 S1I = CYI(J)
00513 IF(INUB.LE.INU) GO TO 225
00514 IF(N.NE.1) GO TO 240
00515 S1R = S2R
00516 S1I = S2I
00517 GO TO 240
00518 270 CONTINUE
00519 YR(1) = S1R
00520 YI(1) = S1I
00521 IF(N.EQ.1) GO TO 280
00522 YR(2) = S2R
00523 YI(2) = S2I
00524 280 CONTINUE
00525 ASCLE = BRY(1)
00526 CALL ZKSCL(ZDR,ZDI,FNU,N,YR,YI,NZ,RZR,RZI,ASCLE,TOL,ELIM)
00527 INU = N - NZ
00528 IF (INU.LE.0) RETURN
00529 KK = NZ + 1
00530 S1R = YR(KK)
00531 S1I = YI(KK)
00532 YR(KK) = S1R*CSRR(1)
00533 YI(KK) = S1I*CSRR(1)
00534 IF (INU.EQ.1) RETURN
00535 KK = NZ + 2
00536 S2R = YR(KK)
00537 S2I = YI(KK)
00538 YR(KK) = S2R*CSRR(1)
00539 YI(KK) = S2I*CSRR(1)
00540 IF (INU.EQ.2) RETURN
00541 T2 = FNU + DBLE(FLOAT(KK-1))
00542 CKR = T2*RZR
00543 CKI = T2*RZI
00544 KFLAG = 1
00545 GO TO 250
00546 290 CONTINUE
00547
00548
00549
00550 KODED = 2
00551 IFLAG = 1
00552 KFLAG = 2
00553 GO TO 120
00554
00555
00556
00557 300 CONTINUE
00558 S1R = COEFR
00559 S1I = COEFI
00560 S2R = COEFR
00561 S2I = COEFI
00562 GO TO 210
00563
00564
00565 310 CONTINUE
00566 NZ=-2
00567 RETURN
00568 END