GNU Octave: libcruft/amos/zunk2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZUNK2(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM,
00002      * ALIM)
00003 C***BEGIN PROLOGUE  ZUNK2
00004 C***REFER TO  ZBESK
00005 C
00006 C     ZUNK2 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
00007 C     RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
00008 C     UNIFORM ASYMPTOTIC EXPANSIONS FOR H(KIND,FNU,ZN) AND J(FNU,ZN)
00009 C     WHERE ZN IS IN THE RIGHT HALF PLANE, KIND=(3-MR)/2, MR=+1 OR
00010 C     -1. HERE ZN=ZR*I OR -ZR*I WHERE ZR=Z IF Z IS IN THE RIGHT
00011 C     HALF PLANE OR ZR=-Z IF Z IS IN THE LEFT HALF PLANE. MR INDIC-
00012 C     ATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
00013 C     NZ=-1 MEANS AN OVERFLOW WILL OCCUR
00014 C
00015 C***ROUTINES CALLED  ZAIRY,ZKSCL,ZS1S2,ZUCHK,ZUNHJ,D1MACH,XZABS
00016 C***END PROLOGUE  ZUNK2
00017 C     COMPLEX AI,ARG,ARGD,ASUM,ASUMD,BSUM,BSUMD,CFN,CI,CIP,CK,CONE,CRSC,
00018 C    *CR1,CR2,CS,CSCL,CSGN,CSPN,CSR,CSS,CY,CZERO,C1,C2,DAI,PHI,PHID,RZ,
00019 C    *S1,S2,Y,Z,ZB,ZETA1,ZETA1D,ZETA2,ZETA2D,ZN,ZR
00020       DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGDI,
00021      * ARGDR, ARGI, ARGR, ASC, ASCLE, ASUMDI, ASUMDR, ASUMI, ASUMR,
00022      * BRY, BSUMDI, BSUMDR, BSUMI, BSUMR, CAR, CIPI, CIPR, CKI, CKR,
00023      * CONER, CRSC, CR1I, CR1R, CR2I, CR2R, CSCL, CSGNI, CSI,
00024      * CSPNI, CSPNR, CSR, CSRR, CSSR, CYI, CYR, C1I, C1R, C2I, C2M,
00025      * C2R, DAII, DAIR, ELIM, FMR, FN, FNF, FNU, HPI, PHIDI, PHIDR,
00026      * PHII, PHIR, PI, PTI, PTR, RAST, RAZR, RS1, RZI, RZR, SAR, SGN,
00027      * STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, YY, ZBI, ZBR, ZEROI,
00028      * ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZET1DI, ZET1DR, ZET2DI,
00029      * ZET2DR, ZI, ZNI, ZNR, ZR, ZRI, ZRR, D1MACH, XZABS
00030       INTEGER I, IB, IFLAG, IFN, IL, IN, INU, IUF, K, KDFLG, KFLAG, KK,
00031      * KODE, MR, N, NAI, NDAI, NW, NZ, IDUM, J, IPARD, IC
00032       DIMENSION BRY(3), YR(N), YI(N), ASUMR(2), ASUMI(2), BSUMR(2),
00033      * BSUMI(2), PHIR(2), PHII(2), ARGR(2), ARGI(2), ZETA1R(2),
00034      * ZETA1I(2), ZETA2R(2), ZETA2I(2), CYR(2), CYI(2), CIPR(4),
00035      * CIPI(4), CSSR(3), CSRR(3)
00036       DATA ZEROR,ZEROI,CONER,CR1R,CR1I,CR2R,CR2I /
00037      1         0.0D0, 0.0D0, 1.0D0,
00038      1 1.0D0,1.73205080756887729D0 , -0.5D0,-8.66025403784438647D-01 /
00039       DATA HPI, PI, AIC /
00040      1     1.57079632679489662D+00,     3.14159265358979324D+00,
00041      1     1.26551212348464539D+00/
00042       DATA CIPR(1),CIPI(1),CIPR(2),CIPI(2),CIPR(3),CIPI(3),CIPR(4),
00043      * CIPI(4) /
00044      1  1.0D0,0.0D0 ,  0.0D0,-1.0D0 ,  -1.0D0,0.0D0 ,  0.0D0,1.0D0 /
00045 C
00046       KDFLG = 1
00047       NZ = 0
00048 C-----------------------------------------------------------------------
00049 C     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
00050 C     THE UNDERFLOW LIMIT
00051 C-----------------------------------------------------------------------
00052       CSCL = 1.0D0/TOL
00053       CRSC = TOL
00054       CSSR(1) = CSCL
00055       CSSR(2) = CONER
00056       CSSR(3) = CRSC
00057       CSRR(1) = CRSC
00058       CSRR(2) = CONER
00059       CSRR(3) = CSCL
00060       BRY(1) = 1.0D+3*D1MACH(1)/TOL
00061       BRY(2) = 1.0D0/BRY(1)
00062       BRY(3) = D1MACH(2)
00063       ZRR = ZR
00064       ZRI = ZI
00065       IF (ZR.GE.0.0D0) GO TO 10
00066       ZRR = -ZR
00067       ZRI = -ZI
00068    10 CONTINUE
00069       YY = ZRI
00070       ZNR = ZRI
00071       ZNI = -ZRR
00072       ZBR = ZRR
00073       ZBI = ZRI
00074       INU = INT(SNGL(FNU))
00075       FNF = FNU - DBLE(FLOAT(INU))
00076       ANG = -HPI*FNF
00077       CAR = DCOS(ANG)
00078       SAR = DSIN(ANG)
00079       C2R = HPI*SAR
00080       C2I = -HPI*CAR
00081       KK = MOD(INU,4) + 1
00082       STR = C2R*CIPR(KK) - C2I*CIPI(KK)
00083       STI = C2R*CIPI(KK) + C2I*CIPR(KK)
00084       CSR = CR1R*STR - CR1I*STI
00085       CSI = CR1R*STI + CR1I*STR
00086       IF (YY.GT.0.0D0) GO TO 20
00087       ZNR = -ZNR
00088       ZBI = -ZBI
00089    20 CONTINUE
00090 C-----------------------------------------------------------------------
00091 C     K(FNU,Z) IS COMPUTED FROM H(2,FNU,-I*Z) WHERE Z IS IN THE FIRST
00092 C     QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
00093 C     CONJUGATION SINCE THE K FUNCTION IS REAL ON THE POSITIVE REAL AXIS
00094 C-----------------------------------------------------------------------
00095       J = 2
00096       DO 80 I=1,N
00097 C-----------------------------------------------------------------------
00098 C     J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
00099 C-----------------------------------------------------------------------
00100         J = 3 - J
00101         FN = FNU + DBLE(FLOAT(I-1))
00102         CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIR(J), PHII(J), ARGR(J),
00103      *   ARGI(J), ZETA1R(J), ZETA1I(J), ZETA2R(J), ZETA2I(J), ASUMR(J),
00104      *   ASUMI(J), BSUMR(J), BSUMI(J))
00105         IF (KODE.EQ.1) GO TO 30
00106         STR = ZBR + ZETA2R(J)
00107         STI = ZBI + ZETA2I(J)
00108         RAST = FN/XZABS(STR,STI)
00109         STR = STR*RAST*RAST
00110         STI = -STI*RAST*RAST
00111         S1R = ZETA1R(J) - STR
00112         S1I = ZETA1I(J) - STI
00113         GO TO 40
00114    30   CONTINUE
00115         S1R = ZETA1R(J) - ZETA2R(J)
00116         S1I = ZETA1I(J) - ZETA2I(J)
00117    40   CONTINUE
00118 C-----------------------------------------------------------------------
00119 C     TEST FOR UNDERFLOW AND OVERFLOW
00120 C-----------------------------------------------------------------------
00121         RS1 = S1R
00122         IF (DABS(RS1).GT.ELIM) GO TO 70
00123         IF (KDFLG.EQ.1) KFLAG = 2
00124         IF (DABS(RS1).LT.ALIM) GO TO 50
00125 C-----------------------------------------------------------------------
00126 C     REFINE  TEST AND SCALE
00127 C-----------------------------------------------------------------------
00128         APHI = XZABS(PHIR(J),PHII(J))
00129         AARG = XZABS(ARGR(J),ARGI(J))
00130         RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC
00131         IF (DABS(RS1).GT.ELIM) GO TO 70
00132         IF (KDFLG.EQ.1) KFLAG = 1
00133         IF (RS1.LT.0.0D0) GO TO 50
00134         IF (KDFLG.EQ.1) KFLAG = 3
00135    50   CONTINUE
00136 C-----------------------------------------------------------------------
00137 C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
00138 C     EXPONENT EXTREMES
00139 C-----------------------------------------------------------------------
00140         C2R = ARGR(J)*CR2R - ARGI(J)*CR2I
00141         C2I = ARGR(J)*CR2I + ARGI(J)*CR2R
00142         CALL ZAIRY(C2R, C2I, 0, 2, AIR, AII, NAI, IDUM)
00143         CALL ZAIRY(C2R, C2I, 1, 2, DAIR, DAII, NDAI, IDUM)
00144         STR = DAIR*BSUMR(J) - DAII*BSUMI(J)
00145         STI = DAIR*BSUMI(J) + DAII*BSUMR(J)
00146         PTR = STR*CR2R - STI*CR2I
00147         PTI = STR*CR2I + STI*CR2R
00148         STR = PTR + (AIR*ASUMR(J)-AII*ASUMI(J))
00149         STI = PTI + (AIR*ASUMI(J)+AII*ASUMR(J))
00150         PTR = STR*PHIR(J) - STI*PHII(J)
00151         PTI = STR*PHII(J) + STI*PHIR(J)
00152         S2R = PTR*CSR - PTI*CSI
00153         S2I = PTR*CSI + PTI*CSR
00154         STR = DEXP(S1R)*CSSR(KFLAG)
00155         S1R = STR*DCOS(S1I)
00156         S1I = STR*DSIN(S1I)
00157         STR = S2R*S1R - S2I*S1I
00158         S2I = S1R*S2I + S2R*S1I
00159         S2R = STR
00160         IF (KFLAG.NE.1) GO TO 60
00161         CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
00162         IF (NW.NE.0) GO TO 70
00163    60   CONTINUE
00164         IF (YY.LE.0.0D0) S2I = -S2I
00165         CYR(KDFLG) = S2R
00166         CYI(KDFLG) = S2I
00167         YR(I) = S2R*CSRR(KFLAG)
00168         YI(I) = S2I*CSRR(KFLAG)
00169         STR = CSI
00170         CSI = -CSR
00171         CSR = STR
00172         IF (KDFLG.EQ.2) GO TO 85
00173         KDFLG = 2
00174         GO TO 80
00175    70   CONTINUE
00176         IF (RS1.GT.0.0D0) GO TO 320
00177 C-----------------------------------------------------------------------
00178 C     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
00179 C-----------------------------------------------------------------------
00180         IF (ZR.LT.0.0D0) GO TO 320
00181         KDFLG = 1
00182         YR(I)=ZEROR
00183         YI(I)=ZEROI
00184         NZ=NZ+1
00185         STR = CSI
00186         CSI =-CSR
00187         CSR = STR
00188         IF (I.EQ.1) GO TO 80
00189         IF ((YR(I-1).EQ.ZEROR).AND.(YI(I-1).EQ.ZEROI)) GO TO 80
00190         YR(I-1)=ZEROR
00191         YI(I-1)=ZEROI
00192         NZ=NZ+1
00193    80 CONTINUE
00194       I = N
00195    85 CONTINUE
00196       RAZR = 1.0D0/XZABS(ZRR,ZRI)
00197       STR = ZRR*RAZR
00198       STI = -ZRI*RAZR
00199       RZR = (STR+STR)*RAZR
00200       RZI = (STI+STI)*RAZR
00201       CKR = FN*RZR
00202       CKI = FN*RZI
00203       IB = I + 1
00204       IF (N.LT.IB) GO TO 180
00205 C-----------------------------------------------------------------------
00206 C     TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
00207 C     ON UNDERFLOW.
00208 C-----------------------------------------------------------------------
00209       FN = FNU + DBLE(FLOAT(N-1))
00210       IPARD = 1
00211       IF (MR.NE.0) IPARD = 0
00212       CALL ZUNHJ(ZNR, ZNI, FN, IPARD, TOL, PHIDR, PHIDI, ARGDR, ARGDI,
00213      * ZET1DR, ZET1DI, ZET2DR, ZET2DI, ASUMDR, ASUMDI, BSUMDR, BSUMDI)
00214       IF (KODE.EQ.1) GO TO 90
00215       STR = ZBR + ZET2DR
00216       STI = ZBI + ZET2DI
00217       RAST = FN/XZABS(STR,STI)
00218       STR = STR*RAST*RAST
00219       STI = -STI*RAST*RAST
00220       S1R = ZET1DR - STR
00221       S1I = ZET1DI - STI
00222       GO TO 100
00223    90 CONTINUE
00224       S1R = ZET1DR - ZET2DR
00225       S1I = ZET1DI - ZET2DI
00226   100 CONTINUE
00227       RS1 = S1R
00228       IF (DABS(RS1).GT.ELIM) GO TO 105
00229       IF (DABS(RS1).LT.ALIM) GO TO 120
00230 C----------------------------------------------------------------------------
00231 C     REFINE ESTIMATE AND TEST
00232 C-------------------------------------------------------------------------
00233       APHI = XZABS(PHIDR,PHIDI)
00234       RS1 = RS1+DLOG(APHI)
00235       IF (DABS(RS1).LT.ELIM) GO TO 120
00236   105 CONTINUE
00237       IF (RS1.GT.0.0D0) GO TO 320
00238 C-----------------------------------------------------------------------
00239 C     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
00240 C-----------------------------------------------------------------------
00241       IF (ZR.LT.0.0D0) GO TO 320
00242       NZ = N
00243       DO 106 I=1,N
00244         YR(I) = ZEROR
00245         YI(I) = ZEROI
00246   106 CONTINUE
00247       RETURN
00248   120 CONTINUE
00249       S1R = CYR(1)
00250       S1I = CYI(1)
00251       S2R = CYR(2)
00252       S2I = CYI(2)
00253       C1R = CSRR(KFLAG)
00254       ASCLE = BRY(KFLAG)
00255       DO 130 I=IB,N
00256         C2R = S2R
00257         C2I = S2I
00258         S2R = CKR*C2R - CKI*C2I + S1R
00259         S2I = CKR*C2I + CKI*C2R + S1I
00260         S1R = C2R
00261         S1I = C2I
00262         CKR = CKR + RZR
00263         CKI = CKI + RZI
00264         C2R = S2R*C1R
00265         C2I = S2I*C1R
00266         YR(I) = C2R
00267         YI(I) = C2I
00268         IF (KFLAG.GE.3) GO TO 130
00269         STR = DABS(C2R)
00270         STI = DABS(C2I)
00271         C2M = DMAX1(STR,STI)
00272         IF (C2M.LE.ASCLE) GO TO 130
00273         KFLAG = KFLAG + 1
00274         ASCLE = BRY(KFLAG)
00275         S1R = S1R*C1R
00276         S1I = S1I*C1R
00277         S2R = C2R
00278         S2I = C2I
00279         S1R = S1R*CSSR(KFLAG)
00280         S1I = S1I*CSSR(KFLAG)
00281         S2R = S2R*CSSR(KFLAG)
00282         S2I = S2I*CSSR(KFLAG)
00283         C1R = CSRR(KFLAG)
00284   130 CONTINUE
00285   180 CONTINUE
00286       IF (MR.EQ.0) RETURN
00287 C-----------------------------------------------------------------------
00288 C     ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
00289 C-----------------------------------------------------------------------
00290       NZ = 0
00291       FMR = DBLE(FLOAT(MR))
00292       SGN = -DSIGN(PI,FMR)
00293 C-----------------------------------------------------------------------
00294 C     CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP.
00295 C-----------------------------------------------------------------------
00296       CSGNI = SGN
00297       IF (YY.LE.0.0D0) CSGNI = -CSGNI
00298       IFN = INU + N - 1
00299       ANG = FNF*SGN
00300       CSPNR = DCOS(ANG)
00301       CSPNI = DSIN(ANG)
00302       IF (MOD(IFN,2).EQ.0) GO TO 190
00303       CSPNR = -CSPNR
00304       CSPNI = -CSPNI
00305   190 CONTINUE
00306 C-----------------------------------------------------------------------
00307 C     CS=COEFF OF THE J FUNCTION TO GET THE I FUNCTION. I(FNU,Z) IS
00308 C     COMPUTED FROM EXP(I*FNU*HPI)*J(FNU,-I*Z) WHERE Z IS IN THE FIRST
00309 C     QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
00310 C     CONJUGATION SINCE THE I FUNCTION IS REAL ON THE POSITIVE REAL AXIS
00311 C-----------------------------------------------------------------------
00312       CSR = SAR*CSGNI
00313       CSI = CAR*CSGNI
00314       IN = MOD(IFN,4) + 1
00315       C2R = CIPR(IN)
00316       C2I = CIPI(IN)
00317       STR = CSR*C2R + CSI*C2I
00318       CSI = -CSR*C2I + CSI*C2R
00319       CSR = STR
00320       ASC = BRY(1)
00321       IUF = 0
00322       KK = N
00323       KDFLG = 1
00324       IB = IB - 1
00325       IC = IB - 1
00326       DO 290 K=1,N
00327         FN = FNU + DBLE(FLOAT(KK-1))
00328 C-----------------------------------------------------------------------
00329 C     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
00330 C     FUNCTION ABOVE
00331 C-----------------------------------------------------------------------
00332         IF (N.GT.2) GO TO 175
00333   172   CONTINUE
00334         PHIDR = PHIR(J)
00335         PHIDI = PHII(J)
00336         ARGDR = ARGR(J)
00337         ARGDI = ARGI(J)
00338         ZET1DR = ZETA1R(J)
00339         ZET1DI = ZETA1I(J)
00340         ZET2DR = ZETA2R(J)
00341         ZET2DI = ZETA2I(J)
00342         ASUMDR = ASUMR(J)
00343         ASUMDI = ASUMI(J)
00344         BSUMDR = BSUMR(J)
00345         BSUMDI = BSUMI(J)
00346         J = 3 - J
00347         GO TO 210
00348   175   CONTINUE
00349         IF ((KK.EQ.N).AND.(IB.LT.N)) GO TO 210
00350         IF ((KK.EQ.IB).OR.(KK.EQ.IC)) GO TO 172
00351         CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIDR, PHIDI, ARGDR,
00352      *   ARGDI, ZET1DR, ZET1DI, ZET2DR, ZET2DI, ASUMDR,
00353      *   ASUMDI, BSUMDR, BSUMDI)
00354   210   CONTINUE
00355         IF (KODE.EQ.1) GO TO 220
00356         STR = ZBR + ZET2DR
00357         STI = ZBI + ZET2DI
00358         RAST = FN/XZABS(STR,STI)
00359         STR = STR*RAST*RAST
00360         STI = -STI*RAST*RAST
00361         S1R = -ZET1DR + STR
00362         S1I = -ZET1DI + STI
00363         GO TO 230
00364   220   CONTINUE
00365         S1R = -ZET1DR + ZET2DR
00366         S1I = -ZET1DI + ZET2DI
00367   230   CONTINUE
00368 C-----------------------------------------------------------------------
00369 C     TEST FOR UNDERFLOW AND OVERFLOW
00370 C-----------------------------------------------------------------------
00371         RS1 = S1R
00372         IF (DABS(RS1).GT.ELIM) GO TO 280
00373         IF (KDFLG.EQ.1) IFLAG = 2
00374         IF (DABS(RS1).LT.ALIM) GO TO 240
00375 C-----------------------------------------------------------------------
00376 C     REFINE  TEST AND SCALE
00377 C-----------------------------------------------------------------------
00378         APHI = XZABS(PHIDR,PHIDI)
00379         AARG = XZABS(ARGDR,ARGDI)
00380         RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC
00381         IF (DABS(RS1).GT.ELIM) GO TO 280
00382         IF (KDFLG.EQ.1) IFLAG = 1
00383         IF (RS1.LT.0.0D0) GO TO 240
00384         IF (KDFLG.EQ.1) IFLAG = 3
00385   240   CONTINUE
00386         CALL ZAIRY(ARGDR, ARGDI, 0, 2, AIR, AII, NAI, IDUM)
00387         CALL ZAIRY(ARGDR, ARGDI, 1, 2, DAIR, DAII, NDAI, IDUM)
00388         STR = DAIR*BSUMDR - DAII*BSUMDI
00389         STI = DAIR*BSUMDI + DAII*BSUMDR
00390         STR = STR + (AIR*ASUMDR-AII*ASUMDI)
00391         STI = STI + (AIR*ASUMDI+AII*ASUMDR)
00392         PTR = STR*PHIDR - STI*PHIDI
00393         PTI = STR*PHIDI + STI*PHIDR
00394         S2R = PTR*CSR - PTI*CSI
00395         S2I = PTR*CSI + PTI*CSR
00396         STR = DEXP(S1R)*CSSR(IFLAG)
00397         S1R = STR*DCOS(S1I)
00398         S1I = STR*DSIN(S1I)
00399         STR = S2R*S1R - S2I*S1I
00400         S2I = S2R*S1I + S2I*S1R
00401         S2R = STR
00402         IF (IFLAG.NE.1) GO TO 250
00403         CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
00404         IF (NW.EQ.0) GO TO 250
00405         S2R = ZEROR
00406         S2I = ZEROI
00407   250   CONTINUE
00408         IF (YY.LE.0.0D0) S2I = -S2I
00409         CYR(KDFLG) = S2R
00410         CYI(KDFLG) = S2I
00411         C2R = S2R
00412         C2I = S2I
00413         S2R = S2R*CSRR(IFLAG)
00414         S2I = S2I*CSRR(IFLAG)
00415 C-----------------------------------------------------------------------
00416 C     ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
00417 C-----------------------------------------------------------------------
00418         S1R = YR(KK)
00419         S1I = YI(KK)
00420         IF (KODE.EQ.1) GO TO 270
00421         CALL ZS1S2(ZRR, ZRI, S1R, S1I, S2R, S2I, NW, ASC, ALIM, IUF)
00422         NZ = NZ + NW
00423   270   CONTINUE
00424         YR(KK) = S1R*CSPNR - S1I*CSPNI + S2R
00425         YI(KK) = S1R*CSPNI + S1I*CSPNR + S2I
00426         KK = KK - 1
00427         CSPNR = -CSPNR
00428         CSPNI = -CSPNI
00429         STR = CSI
00430         CSI = -CSR
00431         CSR = STR
00432         IF (C2R.NE.0.0D0 .OR. C2I.NE.0.0D0) GO TO 255
00433         KDFLG = 1
00434         GO TO 290
00435   255   CONTINUE
00436         IF (KDFLG.EQ.2) GO TO 295
00437         KDFLG = 2
00438         GO TO 290
00439   280   CONTINUE
00440         IF (RS1.GT.0.0D0) GO TO 320
00441         S2R = ZEROR
00442         S2I = ZEROI
00443         GO TO 250
00444   290 CONTINUE
00445       K = N
00446   295 CONTINUE
00447       IL = N - K
00448       IF (IL.EQ.0) RETURN
00449 C-----------------------------------------------------------------------
00450 C     RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
00451 C     K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
00452 C     INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
00453 C-----------------------------------------------------------------------
00454       S1R = CYR(1)
00455       S1I = CYI(1)
00456       S2R = CYR(2)
00457       S2I = CYI(2)
00458       CSR = CSRR(IFLAG)
00459       ASCLE = BRY(IFLAG)
00460       FN = DBLE(FLOAT(INU+IL))
00461       DO 310 I=1,IL
00462         C2R = S2R
00463         C2I = S2I
00464         S2R = S1R + (FN+FNF)*(RZR*C2R-RZI*C2I)
00465         S2I = S1I + (FN+FNF)*(RZR*C2I+RZI*C2R)
00466         S1R = C2R
00467         S1I = C2I
00468         FN = FN - 1.0D0
00469         C2R = S2R*CSR
00470         C2I = S2I*CSR
00471         CKR = C2R
00472         CKI = C2I
00473         C1R = YR(KK)
00474         C1I = YI(KK)
00475         IF (KODE.EQ.1) GO TO 300
00476         CALL ZS1S2(ZRR, ZRI, C1R, C1I, C2R, C2I, NW, ASC, ALIM, IUF)
00477         NZ = NZ + NW
00478   300   CONTINUE
00479         YR(KK) = C1R*CSPNR - C1I*CSPNI + C2R
00480         YI(KK) = C1R*CSPNI + C1I*CSPNR + C2I
00481         KK = KK - 1
00482         CSPNR = -CSPNR
00483         CSPNI = -CSPNI
00484         IF (IFLAG.GE.3) GO TO 310
00485         C2R = DABS(CKR)
00486         C2I = DABS(CKI)
00487         C2M = DMAX1(C2R,C2I)
00488         IF (C2M.LE.ASCLE) GO TO 310
00489         IFLAG = IFLAG + 1
00490         ASCLE = BRY(IFLAG)
00491         S1R = S1R*CSR
00492         S1I = S1I*CSR
00493         S2R = CKR
00494         S2I = CKI
00495         S1R = S1R*CSSR(IFLAG)
00496         S1I = S1I*CSSR(IFLAG)
00497         S2R = S2R*CSSR(IFLAG)
00498         S2I = S2I*CSSR(IFLAG)
00499         CSR = CSRR(IFLAG)
00500   310 CONTINUE
00501       RETURN
00502   320 CONTINUE
00503       NZ = -1
00504       RETURN
00505       END