amos/zunk1.f

   1       SUBROUTINE ZUNK1(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM,
   2      * ALIM)
   3 C***BEGIN PROLOGUE  ZUNK1
   4 C***REFER TO  ZBESK
   5 C
   6 C     ZUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
   7 C     RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
   8 C     UNIFORM ASYMPTOTIC EXPANSION.
   9 C     MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
  10 C     NZ=-1 MEANS AN OVERFLOW WILL OCCUR
  11 C
  12 C***ROUTINES CALLED  ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,AZABS
  13 C***END PROLOGUE  ZUNK1
  14 C     COMPLEX CFN,CK,CONE,CRSC,CS,CSCL,CSGN,CSPN,CSR,CSS,CWRK,CY,CZERO,
  15 C    *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, AZABS
  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 /
  30 C
  31       KDFLG = 1
  32       NZ = 0
  33 C-----------------------------------------------------------------------
  34 C     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
  35 C     THE UNDERFLOW LIMIT
  36 C-----------------------------------------------------------------------
  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
  56 C-----------------------------------------------------------------------
  57 C     J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
  58 C-----------------------------------------------------------------------
  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/AZABS(STR,STI)
  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
  79 C-----------------------------------------------------------------------
  80 C     TEST FOR UNDERFLOW AND OVERFLOW
  81 C-----------------------------------------------------------------------
  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
  85 C-----------------------------------------------------------------------
  86 C     REFINE  TEST AND SCALE
  87 C-----------------------------------------------------------------------
  88         APHI = AZABS(PHIR(J),PHII(J))
  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
  95 C-----------------------------------------------------------------------
  96 C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
  97 C     EXPONENT EXTREMES
  98 C-----------------------------------------------------------------------
  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
 120 C-----------------------------------------------------------------------
 121 C     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
 122 C-----------------------------------------------------------------------
 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/AZABS(ZRR,ZRI)
 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
 145 C-----------------------------------------------------------------------
 146 C     TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
 147 C     ON UNDERFLOW.
 148 C-----------------------------------------------------------------------
 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/AZABS(STR,STI)
 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
 172 C----------------------------------------------------------------------------
 173 C     REFINE ESTIMATE AND TEST
 174 C-------------------------------------------------------------------------
 175       APHI = AZABS(PHIDR,PHIDI)
 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
 180 C-----------------------------------------------------------------------
 181 C     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
 182 C-----------------------------------------------------------------------
 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
 190 C---------------------------------------------------------------------------
 191 C     FORWARD RECUR FOR REMAINDER OF THE SEQUENCE
 192 C----------------------------------------------------------------------------
 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
 232 C-----------------------------------------------------------------------
 233 C     ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
 234 C-----------------------------------------------------------------------
 235       NZ = 0
 236       FMR = DBLE(FLOAT(MR))
 237       SGN = -DSIGN(PI,FMR)
 238 C-----------------------------------------------------------------------
 239 C     CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
 240 C-----------------------------------------------------------------------
 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))
 260 C-----------------------------------------------------------------------
 261 C     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
 262 C     FUNCTION ABOVE
 263 C-----------------------------------------------------------------------
 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/AZABS(STR,STI)
 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
 300 C-----------------------------------------------------------------------
 301 C     TEST FOR UNDERFLOW AND OVERFLOW
 302 C-----------------------------------------------------------------------
 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
 307 C-----------------------------------------------------------------------
 308 C     REFINE  TEST AND SCALE
 309 C-----------------------------------------------------------------------
 310         APHI = AZABS(PHIDR,PHIDI)
 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)
 339 C-----------------------------------------------------------------------
 340 C     ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
 341 C-----------------------------------------------------------------------
 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
 370 C-----------------------------------------------------------------------
 371 C     RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
 372 C     K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
 373 C     INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
 374 C-----------------------------------------------------------------------
 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