amos/zunk2.f

   1       SUBROUTINE ZUNK2(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM,
   2      * ALIM)
   3 C***BEGIN PROLOGUE  ZUNK2
   4 C***REFER TO  ZBESK
   5 C
   6 C     ZUNK2 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 EXPANSIONS FOR H(KIND,FNU,ZN) AND J(FNU,ZN)
   9 C     WHERE ZN IS IN THE RIGHT HALF PLANE, KIND=(3-MR)/2, MR=+1 OR
  10 C     -1. HERE ZN=ZR*I OR -ZR*I WHERE ZR=Z IF Z IS IN THE RIGHT
  11 C     HALF PLANE OR ZR=-Z IF Z IS IN THE LEFT HALF PLANE. MR INDIC-
  12 C     ATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
  13 C     NZ=-1 MEANS AN OVERFLOW WILL OCCUR
  14 C
  15 C***ROUTINES CALLED  ZAIRY,ZKSCL,ZS1S2,ZUCHK,ZUNHJ,D1MACH,AZABS
  16 C***END PROLOGUE  ZUNK2
  17 C     COMPLEX AI,ARG,ARGD,ASUM,ASUMD,BSUM,BSUMD,CFN,CI,CIP,CK,CONE,CRSC,
  18 C    *CR1,CR2,CS,CSCL,CSGN,CSPN,CSR,CSS,CY,CZERO,C1,C2,DAI,PHI,PHID,RZ,
  19 C    *S1,S2,Y,Z,ZB,ZETA1,ZETA1D,ZETA2,ZETA2D,ZN,ZR
  20       DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGDI,
  21      * ARGDR, ARGI, ARGR, ASC, ASCLE, ASUMDI, ASUMDR, ASUMI, ASUMR,
  22      * BRY, BSUMDI, BSUMDR, BSUMI, BSUMR, CAR, CIPI, CIPR, CKI, CKR,
  23      * CONER, CRSC, CR1I, CR1R, CR2I, CR2R, CSCL, CSGNI, CSI,
  24      * CSPNI, CSPNR, CSR, CSRR, CSSR, CYI, CYR, C1I, C1R, C2I, C2M,
  25      * C2R, DAII, DAIR, ELIM, FMR, FN, FNF, FNU, HPI, PHIDI, PHIDR,
  26      * PHII, PHIR, PI, PTI, PTR, RAST, RAZR, RS1, RZI, RZR, SAR, SGN,
  27      * STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, YY, ZBI, ZBR, ZEROI,
  28      * ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZET1DI, ZET1DR, ZET2DI,
  29      * ZET2DR, ZI, ZNI, ZNR, ZR, ZRI, ZRR, D1MACH, AZABS
  30       INTEGER I, IB, IFLAG, IFN, IL, IN, INU, IUF, K, KDFLG, KFLAG, KK,
  31      * KODE, MR, N, NAI, NDAI, NW, NZ, IDUM, J, IPARD, IC
  32       DIMENSION BRY(3), YR(N), YI(N), ASUMR(2), ASUMI(2), BSUMR(2),
  33      * BSUMI(2), PHIR(2), PHII(2), ARGR(2), ARGI(2), ZETA1R(2),
  34      * ZETA1I(2), ZETA2R(2), ZETA2I(2), CYR(2), CYI(2), CIPR(4),
  35      * CIPI(4), CSSR(3), CSRR(3)
  36       DATA ZEROR,ZEROI,CONER,CR1R,CR1I,CR2R,CR2I /
  37      1         0.0D0, 0.0D0, 1.0D0,
  38      1 1.0D0,1.73205080756887729D0 , -0.5D0,-8.66025403784438647D-01 /
  39       DATA HPI, PI, AIC /
  40      1     1.57079632679489662D+00,     3.14159265358979324D+00,
  41      1     1.26551212348464539D+00/
  42       DATA CIPR(1),CIPI(1),CIPR(2),CIPI(2),CIPR(3),CIPI(3),CIPR(4),
  43      * CIPI(4) /
  44      1  1.0D0,0.0D0 ,  0.0D0,-1.0D0 ,  -1.0D0,0.0D0 ,  0.0D0,1.0D0 /
  45 C
  46       KDFLG = 1
  47       NZ = 0
  48 C-----------------------------------------------------------------------
  49 C     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
  50 C     THE UNDERFLOW LIMIT
  51 C-----------------------------------------------------------------------
  52       CSCL = 1.0D0/TOL
  53       CRSC = TOL
  54       CSSR(1) = CSCL
  55       CSSR(2) = CONER
  56       CSSR(3) = CRSC
  57       CSRR(1) = CRSC
  58       CSRR(2) = CONER
  59       CSRR(3) = CSCL
  60       BRY(1) = 1.0D+3*D1MACH(1)/TOL
  61       BRY(2) = 1.0D0/BRY(1)
  62       BRY(3) = D1MACH(2)
  63       ZRR = ZR
  64       ZRI = ZI
  65       IF (ZR.GE.0.0D0) GO TO 10
  66       ZRR = -ZR
  67       ZRI = -ZI
  68    10 CONTINUE
  69       YY = ZRI
  70       ZNR = ZRI
  71       ZNI = -ZRR
  72       ZBR = ZRR
  73       ZBI = ZRI
  74       INU = INT(SNGL(FNU))
  75       FNF = FNU - DBLE(FLOAT(INU))
  76       ANG = -HPI*FNF
  77       CAR = DCOS(ANG)
  78       SAR = DSIN(ANG)
  79       C2R = HPI*SAR
  80       C2I = -HPI*CAR
  81       KK = MOD(INU,4) + 1
  82       STR = C2R*CIPR(KK) - C2I*CIPI(KK)
  83       STI = C2R*CIPI(KK) + C2I*CIPR(KK)
  84       CSR = CR1R*STR - CR1I*STI
  85       CSI = CR1R*STI + CR1I*STR
  86       IF (YY.GT.0.0D0) GO TO 20
  87       ZNR = -ZNR
  88       ZBI = -ZBI
  89    20 CONTINUE
  90 C-----------------------------------------------------------------------
  91 C     K(FNU,Z) IS COMPUTED FROM H(2,FNU,-I*Z) WHERE Z IS IN THE FIRST
  92 C     QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
  93 C     CONJUGATION SINCE THE K FUNCTION IS REAL ON THE POSITIVE REAL AXIS
  94 C-----------------------------------------------------------------------
  95       J = 2
  96       DO 80 I=1,N
  97 C-----------------------------------------------------------------------
  98 C     J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
  99 C-----------------------------------------------------------------------
 100         J = 3 - J
 101         FN = FNU + DBLE(FLOAT(I-1))
 102         CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIR(J), PHII(J), ARGR(J),
 103      *   ARGI(J), ZETA1R(J), ZETA1I(J), ZETA2R(J), ZETA2I(J), ASUMR(J),
 104      *   ASUMI(J), BSUMR(J), BSUMI(J))
 105         IF (KODE.EQ.1) GO TO 30
 106         STR = ZBR + ZETA2R(J)
 107         STI = ZBI + ZETA2I(J)
 108         RAST = FN/AZABS(STR,STI)
 109         STR = STR*RAST*RAST
 110         STI = -STI*RAST*RAST
 111         S1R = ZETA1R(J) - STR
 112         S1I = ZETA1I(J) - STI
 113         GO TO 40
 114    30   CONTINUE
 115         S1R = ZETA1R(J) - ZETA2R(J)
 116         S1I = ZETA1I(J) - ZETA2I(J)
 117    40   CONTINUE
 118 C-----------------------------------------------------------------------
 119 C     TEST FOR UNDERFLOW AND OVERFLOW
 120 C-----------------------------------------------------------------------
 121         RS1 = S1R
 122         IF (DABS(RS1).GT.ELIM) GO TO 70
 123         IF (KDFLG.EQ.1) KFLAG = 2
 124         IF (DABS(RS1).LT.ALIM) GO TO 50
 125 C-----------------------------------------------------------------------
 126 C     REFINE  TEST AND SCALE
 127 C-----------------------------------------------------------------------
 128         APHI = AZABS(PHIR(J),PHII(J))
 129         AARG = AZABS(ARGR(J),ARGI(J))
 130         RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC
 131         IF (DABS(RS1).GT.ELIM) GO TO 70
 132         IF (KDFLG.EQ.1) KFLAG = 1
 133         IF (RS1.LT.0.0D0) GO TO 50
 134         IF (KDFLG.EQ.1) KFLAG = 3
 135    50   CONTINUE
 136 C-----------------------------------------------------------------------
 137 C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
 138 C     EXPONENT EXTREMES
 139 C-----------------------------------------------------------------------
 140         C2R = ARGR(J)*CR2R - ARGI(J)*CR2I
 141         C2I = ARGR(J)*CR2I + ARGI(J)*CR2R
 142         CALL ZAIRY(C2R, C2I, 0, 2, AIR, AII, NAI, IDUM)
 143         CALL ZAIRY(C2R, C2I, 1, 2, DAIR, DAII, NDAI, IDUM)
 144         STR = DAIR*BSUMR(J) - DAII*BSUMI(J)
 145         STI = DAIR*BSUMI(J) + DAII*BSUMR(J)
 146         PTR = STR*CR2R - STI*CR2I
 147         PTI = STR*CR2I + STI*CR2R
 148         STR = PTR + (AIR*ASUMR(J)-AII*ASUMI(J))
 149         STI = PTI + (AIR*ASUMI(J)+AII*ASUMR(J))
 150         PTR = STR*PHIR(J) - STI*PHII(J)
 151         PTI = STR*PHII(J) + STI*PHIR(J)
 152         S2R = PTR*CSR - PTI*CSI
 153         S2I = PTR*CSI + PTI*CSR
 154         STR = DEXP(S1R)*CSSR(KFLAG)
 155         S1R = STR*DCOS(S1I)
 156         S1I = STR*DSIN(S1I)
 157         STR = S2R*S1R - S2I*S1I
 158         S2I = S1R*S2I + S2R*S1I
 159         S2R = STR
 160         IF (KFLAG.NE.1) GO TO 60
 161         CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
 162         IF (NW.NE.0) GO TO 70
 163    60   CONTINUE
 164         IF (YY.LE.0.0D0) S2I = -S2I
 165         CYR(KDFLG) = S2R
 166         CYI(KDFLG) = S2I
 167         YR(I) = S2R*CSRR(KFLAG)
 168         YI(I) = S2I*CSRR(KFLAG)
 169         STR = CSI
 170         CSI = -CSR
 171         CSR = STR
 172         IF (KDFLG.EQ.2) GO TO 85
 173         KDFLG = 2
 174         GO TO 80
 175    70   CONTINUE
 176         IF (RS1.GT.0.0D0) GO TO 320
 177 C-----------------------------------------------------------------------
 178 C     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
 179 C-----------------------------------------------------------------------
 180         IF (ZR.LT.0.0D0) GO TO 320
 181         KDFLG = 1
 182         YR(I)=ZEROR
 183         YI(I)=ZEROI
 184         NZ=NZ+1
 185         STR = CSI
 186         CSI =-CSR
 187         CSR = STR
 188         IF (I.EQ.1) GO TO 80
 189         IF ((YR(I-1).EQ.ZEROR).AND.(YI(I-1).EQ.ZEROI)) GO TO 80
 190         YR(I-1)=ZEROR
 191         YI(I-1)=ZEROI
 192         NZ=NZ+1
 193    80 CONTINUE
 194       I = N
 195    85 CONTINUE
 196       RAZR = 1.0D0/AZABS(ZRR,ZRI)
 197       STR = ZRR*RAZR
 198       STI = -ZRI*RAZR
 199       RZR = (STR+STR)*RAZR
 200       RZI = (STI+STI)*RAZR
 201       CKR = FN*RZR
 202       CKI = FN*RZI
 203       IB = I + 1
 204       IF (N.LT.IB) GO TO 180
 205 C-----------------------------------------------------------------------
 206 C     TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
 207 C     ON UNDERFLOW.
 208 C-----------------------------------------------------------------------
 209       FN = FNU + DBLE(FLOAT(N-1))
 210       IPARD = 1
 211       IF (MR.NE.0) IPARD = 0
 212       CALL ZUNHJ(ZNR, ZNI, FN, IPARD, TOL, PHIDR, PHIDI, ARGDR, ARGDI,
 213      * ZET1DR, ZET1DI, ZET2DR, ZET2DI, ASUMDR, ASUMDI, BSUMDR, BSUMDI)
 214       IF (KODE.EQ.1) GO TO 90
 215       STR = ZBR + ZET2DR
 216       STI = ZBI + ZET2DI
 217       RAST = FN/AZABS(STR,STI)
 218       STR = STR*RAST*RAST
 219       STI = -STI*RAST*RAST
 220       S1R = ZET1DR - STR
 221       S1I = ZET1DI - STI
 222       GO TO 100
 223    90 CONTINUE
 224       S1R = ZET1DR - ZET2DR
 225       S1I = ZET1DI - ZET2DI
 226   100 CONTINUE
 227       RS1 = S1R
 228       IF (DABS(RS1).GT.ELIM) GO TO 105
 229       IF (DABS(RS1).LT.ALIM) GO TO 120
 230 C----------------------------------------------------------------------------
 231 C     REFINE ESTIMATE AND TEST
 232 C-------------------------------------------------------------------------
 233       APHI = AZABS(PHIDR,PHIDI)
 234       RS1 = RS1+DLOG(APHI)
 235       IF (DABS(RS1).LT.ELIM) GO TO 120
 236   105 CONTINUE
 237       IF (RS1.GT.0.0D0) GO TO 320
 238 C-----------------------------------------------------------------------
 239 C     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
 240 C-----------------------------------------------------------------------
 241       IF (ZR.LT.0.0D0) GO TO 320
 242       NZ = N
 243       DO 106 I=1,N
 244         YR(I) = ZEROR
 245         YI(I) = ZEROI
 246   106 CONTINUE
 247       RETURN
 248   120 CONTINUE
 249       S1R = CYR(1)
 250       S1I = CYI(1)
 251       S2R = CYR(2)
 252       S2I = CYI(2)
 253       C1R = CSRR(KFLAG)
 254       ASCLE = BRY(KFLAG)
 255       DO 130 I=IB,N
 256         C2R = S2R
 257         C2I = S2I
 258         S2R = CKR*C2R - CKI*C2I + S1R
 259         S2I = CKR*C2I + CKI*C2R + S1I
 260         S1R = C2R
 261         S1I = C2I
 262         CKR = CKR + RZR
 263         CKI = CKI + RZI
 264         C2R = S2R*C1R
 265         C2I = S2I*C1R
 266         YR(I) = C2R
 267         YI(I) = C2I
 268         IF (KFLAG.GE.3) GO TO 130
 269         STR = DABS(C2R)
 270         STI = DABS(C2I)
 271         C2M = DMAX1(STR,STI)
 272         IF (C2M.LE.ASCLE) GO TO 130
 273         KFLAG = KFLAG + 1
 274         ASCLE = BRY(KFLAG)
 275         S1R = S1R*C1R
 276         S1I = S1I*C1R
 277         S2R = C2R
 278         S2I = C2I
 279         S1R = S1R*CSSR(KFLAG)
 280         S1I = S1I*CSSR(KFLAG)
 281         S2R = S2R*CSSR(KFLAG)
 282         S2I = S2I*CSSR(KFLAG)
 283         C1R = CSRR(KFLAG)
 284   130 CONTINUE
 285   180 CONTINUE
 286       IF (MR.EQ.0) RETURN
 287 C-----------------------------------------------------------------------
 288 C     ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
 289 C-----------------------------------------------------------------------
 290       NZ = 0
 291       FMR = DBLE(FLOAT(MR))
 292       SGN = -DSIGN(PI,FMR)
 293 C-----------------------------------------------------------------------
 294 C     CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
 295 C-----------------------------------------------------------------------
 296       CSGNI = SGN
 297       IF (YY.LE.0.0D0) CSGNI = -CSGNI
 298       IFN = INU + N - 1
 299       ANG = FNF*SGN
 300       CSPNR = DCOS(ANG)
 301       CSPNI = DSIN(ANG)
 302       IF (MOD(IFN,2).EQ.0) GO TO 190
 303       CSPNR = -CSPNR
 304       CSPNI = -CSPNI
 305   190 CONTINUE
 306 C-----------------------------------------------------------------------
 307 C     CS=COEFF OF THE J FUNCTION TO GET THE I FUNCTION. I(FNU,Z) IS
 308 C     COMPUTED FROM EXP(I*FNU*HPI)*J(FNU,-I*Z) WHERE Z IS IN THE FIRST
 309 C     QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
 310 C     CONJUGATION SINCE THE I FUNCTION IS REAL ON THE POSITIVE REAL AXIS
 311 C-----------------------------------------------------------------------
 312       CSR = SAR*CSGNI
 313       CSI = CAR*CSGNI
 314       IN = MOD(IFN,4) + 1
 315       C2R = CIPR(IN)
 316       C2I = CIPI(IN)
 317       STR = CSR*C2R + CSI*C2I
 318       CSI = -CSR*C2I + CSI*C2R
 319       CSR = STR
 320       ASC = BRY(1)
 321       IUF = 0
 322       KK = N
 323       KDFLG = 1
 324       IB = IB - 1
 325       IC = IB - 1
 326       DO 290 K=1,N
 327         FN = FNU + DBLE(FLOAT(KK-1))
 328 C-----------------------------------------------------------------------
 329 C     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
 330 C     FUNCTION ABOVE
 331 C-----------------------------------------------------------------------
 332         IF (N.GT.2) GO TO 175
 333   172   CONTINUE
 334         PHIDR = PHIR(J)
 335         PHIDI = PHII(J)
 336         ARGDR = ARGR(J)
 337         ARGDI = ARGI(J)
 338         ZET1DR = ZETA1R(J)
 339         ZET1DI = ZETA1I(J)
 340         ZET2DR = ZETA2R(J)
 341         ZET2DI = ZETA2I(J)
 342         ASUMDR = ASUMR(J)
 343         ASUMDI = ASUMI(J)
 344         BSUMDR = BSUMR(J)
 345         BSUMDI = BSUMI(J)
 346         J = 3 - J
 347         GO TO 210
 348   175   CONTINUE
 349         IF ((KK.EQ.N).AND.(IB.LT.N)) GO TO 210
 350         IF ((KK.EQ.IB).OR.(KK.EQ.IC)) GO TO 172
 351         CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIDR, PHIDI, ARGDR,
 352      *   ARGDI, ZET1DR, ZET1DI, ZET2DR, ZET2DI, ASUMDR,
 353      *   ASUMDI, BSUMDR, BSUMDI)
 354   210   CONTINUE
 355         IF (KODE.EQ.1) GO TO 220
 356         STR = ZBR + ZET2DR
 357         STI = ZBI + ZET2DI
 358         RAST = FN/AZABS(STR,STI)
 359         STR = STR*RAST*RAST
 360         STI = -STI*RAST*RAST
 361         S1R = -ZET1DR + STR
 362         S1I = -ZET1DI + STI
 363         GO TO 230
 364   220   CONTINUE
 365         S1R = -ZET1DR + ZET2DR
 366         S1I = -ZET1DI + ZET2DI
 367   230   CONTINUE
 368 C-----------------------------------------------------------------------
 369 C     TEST FOR UNDERFLOW AND OVERFLOW
 370 C-----------------------------------------------------------------------
 371         RS1 = S1R
 372         IF (DABS(RS1).GT.ELIM) GO TO 280
 373         IF (KDFLG.EQ.1) IFLAG = 2
 374         IF (DABS(RS1).LT.ALIM) GO TO 240
 375 C-----------------------------------------------------------------------
 376 C     REFINE  TEST AND SCALE
 377 C-----------------------------------------------------------------------
 378         APHI = AZABS(PHIDR,PHIDI)
 379         AARG = AZABS(ARGDR,ARGDI)
 380         RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC
 381         IF (DABS(RS1).GT.ELIM) GO TO 280
 382         IF (KDFLG.EQ.1) IFLAG = 1
 383         IF (RS1.LT.0.0D0) GO TO 240
 384         IF (KDFLG.EQ.1) IFLAG = 3
 385   240   CONTINUE
 386         CALL ZAIRY(ARGDR, ARGDI, 0, 2, AIR, AII, NAI, IDUM)
 387         CALL ZAIRY(ARGDR, ARGDI, 1, 2, DAIR, DAII, NDAI, IDUM)
 388         STR = DAIR*BSUMDR - DAII*BSUMDI
 389         STI = DAIR*BSUMDI + DAII*BSUMDR
 390         STR = STR + (AIR*ASUMDR-AII*ASUMDI)
 391         STI = STI + (AIR*ASUMDI+AII*ASUMDR)
 392         PTR = STR*PHIDR - STI*PHIDI
 393         PTI = STR*PHIDI + STI*PHIDR
 394         S2R = PTR*CSR - PTI*CSI
 395         S2I = PTR*CSI + PTI*CSR
 396         STR = DEXP(S1R)*CSSR(IFLAG)
 397         S1R = STR*DCOS(S1I)
 398         S1I = STR*DSIN(S1I)
 399         STR = S2R*S1R - S2I*S1I
 400         S2I = S2R*S1I + S2I*S1R
 401         S2R = STR
 402         IF (IFLAG.NE.1) GO TO 250
 403         CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
 404         IF (NW.EQ.0) GO TO 250
 405         S2R = ZEROR
 406         S2I = ZEROI
 407   250   CONTINUE
 408         IF (YY.LE.0.0D0) S2I = -S2I
 409         CYR(KDFLG) = S2R
 410         CYI(KDFLG) = S2I
 411         C2R = S2R
 412         C2I = S2I
 413         S2R = S2R*CSRR(IFLAG)
 414         S2I = S2I*CSRR(IFLAG)
 415 C-----------------------------------------------------------------------
 416 C     ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
 417 C-----------------------------------------------------------------------
 418         S1R = YR(KK)
 419         S1I = YI(KK)
 420         IF (KODE.EQ.1) GO TO 270
 421         CALL ZS1S2(ZRR, ZRI, S1R, S1I, S2R, S2I, NW, ASC, ALIM, IUF)
 422         NZ = NZ + NW
 423   270   CONTINUE
 424         YR(KK) = S1R*CSPNR - S1I*CSPNI + S2R
 425         YI(KK) = S1R*CSPNI + S1I*CSPNR + S2I
 426         KK = KK - 1
 427         CSPNR = -CSPNR
 428         CSPNI = -CSPNI
 429         STR = CSI
 430         CSI = -CSR
 431         CSR = STR
 432         IF (C2R.NE.0.0D0 .OR. C2I.NE.0.0D0) GO TO 255
 433         KDFLG = 1
 434         GO TO 290
 435   255   CONTINUE
 436         IF (KDFLG.EQ.2) GO TO 295
 437         KDFLG = 2
 438         GO TO 290
 439   280   CONTINUE
 440         IF (RS1.GT.0.0D0) GO TO 320
 441         S2R = ZEROR
 442         S2I = ZEROI
 443         GO TO 250
 444   290 CONTINUE
 445       K = N
 446   295 CONTINUE
 447       IL = N - K
 448       IF (IL.EQ.0) RETURN
 449 C-----------------------------------------------------------------------
 450 C     RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
 451 C     K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
 452 C     INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
 453 C-----------------------------------------------------------------------
 454       S1R = CYR(1)
 455       S1I = CYI(1)
 456       S2R = CYR(2)
 457       S2I = CYI(2)
 458       CSR = CSRR(IFLAG)
 459       ASCLE = BRY(IFLAG)
 460       FN = DBLE(FLOAT(INU+IL))
 461       DO 310 I=1,IL
 462         C2R = S2R
 463         C2I = S2I
 464         S2R = S1R + (FN+FNF)*(RZR*C2R-RZI*C2I)
 465         S2I = S1I + (FN+FNF)*(RZR*C2I+RZI*C2R)
 466         S1R = C2R
 467         S1I = C2I
 468         FN = FN - 1.0D0
 469         C2R = S2R*CSR
 470         C2I = S2I*CSR
 471         CKR = C2R
 472         CKI = C2I
 473         C1R = YR(KK)
 474         C1I = YI(KK)
 475         IF (KODE.EQ.1) GO TO 300
 476         CALL ZS1S2(ZRR, ZRI, C1R, C1I, C2R, C2I, NW, ASC, ALIM, IUF)
 477         NZ = NZ + NW
 478   300   CONTINUE
 479         YR(KK) = C1R*CSPNR - C1I*CSPNI + C2R
 480         YI(KK) = C1R*CSPNI + C1I*CSPNR + C2I
 481         KK = KK - 1
 482         CSPNR = -CSPNR
 483         CSPNI = -CSPNI
 484         IF (IFLAG.GE.3) GO TO 310
 485         C2R = DABS(CKR)
 486         C2I = DABS(CKI)
 487         C2M = DMAX1(C2R,C2I)
 488         IF (C2M.LE.ASCLE) GO TO 310
 489         IFLAG = IFLAG + 1
 490         ASCLE = BRY(IFLAG)
 491         S1R = S1R*CSR
 492         S1I = S1I*CSR
 493         S2R = CKR
 494         S2I = CKI
 495         S1R = S1R*CSSR(IFLAG)
 496         S1I = S1I*CSSR(IFLAG)
 497         S2R = S2R*CSSR(IFLAG)
 498         S2I = S2I*CSSR(IFLAG)
 499         CSR = CSRR(IFLAG)
 500   310 CONTINUE
 501       RETURN
 502   320 CONTINUE
 503       NZ = -1
 504       RETURN
 505       END