src/numerical/slatec/fortran/zunk1.f

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