github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/mathext/internal/amos/amoslib/zbknu.f (about) 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,ZABS,ZDIV, 9 C ZEXP,ZLOG,ZMLT,ZSQRT 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, ZABS, 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 = ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0))) 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 ZLOG(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 ZEXP(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 ZEXP(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 ZEXP(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 ZSQRT(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 = ZABS(CMPLX(CSR,CSI,kind=KIND(1.0D0))) 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 = ZABS(CMPLX(P2R,P2I,kind=KIND(1.0D0))) 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 = ZABS(CMPLX(S2R,S2I,kind=KIND(1.0D0))) 476 ALAS = DLOG(AS) 477 P2R = -ZDR+ALAS 478 IF(P2R.LT.(-ELIM)) GO TO 263 479 CALL ZLOG(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