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