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

Go to the documentation of this file.
00001       SUBROUTINE ZUNK1(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM,
00002      * ALIM)
00003 C***BEGIN PROLOGUE  ZUNK1
00004 C***REFER TO  ZBESK
00005 C
00006 C     ZUNK1 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 EXPANSION.
00009 C     MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
00010 C     NZ=-1 MEANS AN OVERFLOW WILL OCCUR
00011 C
00012 C***ROUTINES CALLED  ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,XZABS
00013 C***END PROLOGUE  ZUNK1
00014 C     COMPLEX CFN,CK,CONE,CRSC,CS,CSCL,CSGN,CSPN,CSR,CSS,CWRK,CY,CZERO,
00015 C    *C1,C2,PHI,PHID,RZ,SUM,SUMD,S1,S2,Y,Z,ZETA1,ZETA1D,ZETA2,ZETA2D,ZR
00016       DOUBLE PRECISION ALIM, ANG, APHI, ASC, ASCLE, BRY, CKI, CKR,
00017      * CONER, CRSC, CSCL, CSGNI, CSPNI, CSPNR, CSR, CSRR, CSSR,
00018      * CWRKI, CWRKR, CYI, CYR, C1I, C1R, C2I, C2M, C2R, ELIM, FMR, FN,
00019      * FNF, FNU, PHIDI, PHIDR, PHII, PHIR, PI, RAST, RAZR, RS1, RZI,
00020      * RZR, SGN, STI, STR, SUMDI, SUMDR, SUMI, SUMR, S1I, S1R, S2I,
00021      * S2R, TOL, YI, YR, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R,
00022      * ZET1DI, ZET1DR, ZET2DI, ZET2DR, ZI, ZR, ZRI, ZRR, D1MACH, XZABS
00023       INTEGER I, IB, IFLAG, IFN, IL, INIT, INU, IUF, K, KDFLG, KFLAG,
00024      * KK, KODE, MR, N, NW, NZ, INITD, IC, IPARD, J
00025       DIMENSION BRY(3), INIT(2), YR(N), YI(N), SUMR(2), SUMI(2),
00026      * ZETA1R(2), ZETA1I(2), ZETA2R(2), ZETA2I(2), CYR(2), CYI(2),
00027      * CWRKR(16,3), CWRKI(16,3), CSSR(3), CSRR(3), PHIR(2), PHII(2)
00028       DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 /
00029       DATA PI / 3.14159265358979324D0 /
00030 C
00031       KDFLG = 1
00032       NZ = 0
00033 C-----------------------------------------------------------------------
00034 C     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
00035 C     THE UNDERFLOW LIMIT
00036 C-----------------------------------------------------------------------
00037       CSCL = 1.0D0/TOL
00038       CRSC = TOL
00039       CSSR(1) = CSCL
00040       CSSR(2) = CONER
00041       CSSR(3) = CRSC
00042       CSRR(1) = CRSC
00043       CSRR(2) = CONER
00044       CSRR(3) = CSCL
00045       BRY(1) = 1.0D+3*D1MACH(1)/TOL
00046       BRY(2) = 1.0D0/BRY(1)
00047       BRY(3) = D1MACH(2)
00048       ZRR = ZR
00049       ZRI = ZI
00050       IF (ZR.GE.0.0D0) GO TO 10
00051       ZRR = -ZR
00052       ZRI = -ZI
00053    10 CONTINUE
00054       J = 2
00055       DO 70 I=1,N
00056 C-----------------------------------------------------------------------
00057 C     J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
00058 C-----------------------------------------------------------------------
00059         J = 3 - J
00060         FN = FNU + DBLE(FLOAT(I-1))
00061         INIT(J) = 0
00062         CALL ZUNIK(ZRR, ZRI, FN, 2, 0, TOL, INIT(J), PHIR(J), PHII(J),
00063      *   ZETA1R(J), ZETA1I(J), ZETA2R(J), ZETA2I(J), SUMR(J), SUMI(J),
00064      *   CWRKR(1,J), CWRKI(1,J))
00065         IF (KODE.EQ.1) GO TO 20
00066         STR = ZRR + ZETA2R(J)
00067         STI = ZRI + ZETA2I(J)
00068         RAST = FN/XZABS(STR,STI)
00069         STR = STR*RAST*RAST
00070         STI = -STI*RAST*RAST
00071         S1R = ZETA1R(J) - STR
00072         S1I = ZETA1I(J) - STI
00073         GO TO 30
00074    20   CONTINUE
00075         S1R = ZETA1R(J) - ZETA2R(J)
00076         S1I = ZETA1I(J) - ZETA2I(J)
00077    30   CONTINUE
00078         RS1 = S1R
00079 C-----------------------------------------------------------------------
00080 C     TEST FOR UNDERFLOW AND OVERFLOW
00081 C-----------------------------------------------------------------------
00082         IF (DABS(RS1).GT.ELIM) GO TO 60
00083         IF (KDFLG.EQ.1) KFLAG = 2
00084         IF (DABS(RS1).LT.ALIM) GO TO 40
00085 C-----------------------------------------------------------------------
00086 C     REFINE  TEST AND SCALE
00087 C-----------------------------------------------------------------------
00088         APHI = XZABS(PHIR(J),PHII(J))
00089         RS1 = RS1 + DLOG(APHI)
00090         IF (DABS(RS1).GT.ELIM) GO TO 60
00091         IF (KDFLG.EQ.1) KFLAG = 1
00092         IF (RS1.LT.0.0D0) GO TO 40
00093         IF (KDFLG.EQ.1) KFLAG = 3
00094    40   CONTINUE
00095 C-----------------------------------------------------------------------
00096 C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
00097 C     EXPONENT EXTREMES
00098 C-----------------------------------------------------------------------
00099         S2R = PHIR(J)*SUMR(J) - PHII(J)*SUMI(J)
00100         S2I = PHIR(J)*SUMI(J) + PHII(J)*SUMR(J)
00101         STR = DEXP(S1R)*CSSR(KFLAG)
00102         S1R = STR*DCOS(S1I)
00103         S1I = STR*DSIN(S1I)
00104         STR = S2R*S1R - S2I*S1I
00105         S2I = S1R*S2I + S2R*S1I
00106         S2R = STR
00107         IF (KFLAG.NE.1) GO TO 50
00108         CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
00109         IF (NW.NE.0) GO TO 60
00110    50   CONTINUE
00111         CYR(KDFLG) = S2R
00112         CYI(KDFLG) = S2I
00113         YR(I) = S2R*CSRR(KFLAG)
00114         YI(I) = S2I*CSRR(KFLAG)
00115         IF (KDFLG.EQ.2) GO TO 75
00116         KDFLG = 2
00117         GO TO 70
00118    60   CONTINUE
00119         IF (RS1.GT.0.0D0) GO TO 300
00120 C-----------------------------------------------------------------------
00121 C     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
00122 C-----------------------------------------------------------------------
00123         IF (ZR.LT.0.0D0) GO TO 300
00124         KDFLG = 1
00125         YR(I)=ZEROR
00126         YI(I)=ZEROI
00127         NZ=NZ+1
00128         IF (I.EQ.1) GO TO 70
00129         IF ((YR(I-1).EQ.ZEROR).AND.(YI(I-1).EQ.ZEROI)) GO TO 70
00130         YR(I-1)=ZEROR
00131         YI(I-1)=ZEROI
00132         NZ=NZ+1
00133    70 CONTINUE
00134       I = N
00135    75 CONTINUE
00136       RAZR = 1.0D0/XZABS(ZRR,ZRI)
00137       STR = ZRR*RAZR
00138       STI = -ZRI*RAZR
00139       RZR = (STR+STR)*RAZR
00140       RZI = (STI+STI)*RAZR
00141       CKR = FN*RZR
00142       CKI = FN*RZI
00143       IB = I + 1
00144       IF (N.LT.IB) GO TO 160
00145 C-----------------------------------------------------------------------
00146 C     TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
00147 C     ON UNDERFLOW.
00148 C-----------------------------------------------------------------------
00149       FN = FNU + DBLE(FLOAT(N-1))
00150       IPARD = 1
00151       IF (MR.NE.0) IPARD = 0
00152       INITD = 0
00153       CALL ZUNIK(ZRR, ZRI, FN, 2, IPARD, TOL, INITD, PHIDR, PHIDI,
00154      * ZET1DR, ZET1DI, ZET2DR, ZET2DI, SUMDR, SUMDI, CWRKR(1,3),
00155      * CWRKI(1,3))
00156       IF (KODE.EQ.1) GO TO 80
00157       STR = ZRR + ZET2DR
00158       STI = ZRI + ZET2DI
00159       RAST = FN/XZABS(STR,STI)
00160       STR = STR*RAST*RAST
00161       STI = -STI*RAST*RAST
00162       S1R = ZET1DR - STR
00163       S1I = ZET1DI - STI
00164       GO TO 90
00165    80 CONTINUE
00166       S1R = ZET1DR - ZET2DR
00167       S1I = ZET1DI - ZET2DI
00168    90 CONTINUE
00169       RS1 = S1R
00170       IF (DABS(RS1).GT.ELIM) GO TO 95
00171       IF (DABS(RS1).LT.ALIM) GO TO 100
00172 C----------------------------------------------------------------------------
00173 C     REFINE ESTIMATE AND TEST
00174 C-------------------------------------------------------------------------
00175       APHI = XZABS(PHIDR,PHIDI)
00176       RS1 = RS1+DLOG(APHI)
00177       IF (DABS(RS1).LT.ELIM) GO TO 100
00178    95 CONTINUE
00179       IF (DABS(RS1).GT.0.0D0) GO TO 300
00180 C-----------------------------------------------------------------------
00181 C     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
00182 C-----------------------------------------------------------------------
00183       IF (ZR.LT.0.0D0) GO TO 300
00184       NZ = N
00185       DO 96 I=1,N
00186         YR(I) = ZEROR
00187         YI(I) = ZEROI
00188    96 CONTINUE
00189       RETURN
00190 C---------------------------------------------------------------------------
00191 C     FORWARD RECUR FOR REMAINDER OF THE SEQUENCE
00192 C----------------------------------------------------------------------------
00193   100 CONTINUE
00194       S1R = CYR(1)
00195       S1I = CYI(1)
00196       S2R = CYR(2)
00197       S2I = CYI(2)
00198       C1R = CSRR(KFLAG)
00199       ASCLE = BRY(KFLAG)
00200       DO 120 I=IB,N
00201         C2R = S2R
00202         C2I = S2I
00203         S2R = CKR*C2R - CKI*C2I + S1R
00204         S2I = CKR*C2I + CKI*C2R + S1I
00205         S1R = C2R
00206         S1I = C2I
00207         CKR = CKR + RZR
00208         CKI = CKI + RZI
00209         C2R = S2R*C1R
00210         C2I = S2I*C1R
00211         YR(I) = C2R
00212         YI(I) = C2I
00213         IF (KFLAG.GE.3) GO TO 120
00214         STR = DABS(C2R)
00215         STI = DABS(C2I)
00216         C2M = DMAX1(STR,STI)
00217         IF (C2M.LE.ASCLE) GO TO 120
00218         KFLAG = KFLAG + 1
00219         ASCLE = BRY(KFLAG)
00220         S1R = S1R*C1R
00221         S1I = S1I*C1R
00222         S2R = C2R
00223         S2I = C2I
00224         S1R = S1R*CSSR(KFLAG)
00225         S1I = S1I*CSSR(KFLAG)
00226         S2R = S2R*CSSR(KFLAG)
00227         S2I = S2I*CSSR(KFLAG)
00228         C1R = CSRR(KFLAG)
00229   120 CONTINUE
00230   160 CONTINUE
00231       IF (MR.EQ.0) RETURN
00232 C-----------------------------------------------------------------------
00233 C     ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
00234 C-----------------------------------------------------------------------
00235       NZ = 0
00236       FMR = DBLE(FLOAT(MR))
00237       SGN = -DSIGN(PI,FMR)
00238 C-----------------------------------------------------------------------
00239 C     CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
00240 C-----------------------------------------------------------------------
00241       CSGNI = SGN
00242       INU = INT(SNGL(FNU))
00243       FNF = FNU - DBLE(FLOAT(INU))
00244       IFN = INU + N - 1
00245       ANG = FNF*SGN
00246       CSPNR = DCOS(ANG)
00247       CSPNI = DSIN(ANG)
00248       IF (MOD(IFN,2).EQ.0) GO TO 170
00249       CSPNR = -CSPNR
00250       CSPNI = -CSPNI
00251   170 CONTINUE
00252       ASC = BRY(1)
00253       IUF = 0
00254       KK = N
00255       KDFLG = 1
00256       IB = IB - 1
00257       IC = IB - 1
00258       DO 270 K=1,N
00259         FN = FNU + DBLE(FLOAT(KK-1))
00260 C-----------------------------------------------------------------------
00261 C     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
00262 C     FUNCTION ABOVE
00263 C-----------------------------------------------------------------------
00264         M=3
00265         IF (N.GT.2) GO TO 175
00266   172   CONTINUE
00267         INITD = INIT(J)
00268         PHIDR = PHIR(J)
00269         PHIDI = PHII(J)
00270         ZET1DR = ZETA1R(J)
00271         ZET1DI = ZETA1I(J)
00272         ZET2DR = ZETA2R(J)
00273         ZET2DI = ZETA2I(J)
00274         SUMDR = SUMR(J)
00275         SUMDI = SUMI(J)
00276         M = J
00277         J = 3 - J
00278         GO TO 180
00279   175   CONTINUE
00280         IF ((KK.EQ.N).AND.(IB.LT.N)) GO TO 180
00281         IF ((KK.EQ.IB).OR.(KK.EQ.IC)) GO TO 172
00282         INITD = 0
00283   180   CONTINUE
00284         CALL ZUNIK(ZRR, ZRI, FN, 1, 0, TOL, INITD, PHIDR, PHIDI,
00285      *   ZET1DR, ZET1DI, ZET2DR, ZET2DI, SUMDR, SUMDI,
00286      *   CWRKR(1,M), CWRKI(1,M))
00287         IF (KODE.EQ.1) GO TO 200
00288         STR = ZRR + ZET2DR
00289         STI = ZRI + ZET2DI
00290         RAST = FN/XZABS(STR,STI)
00291         STR = STR*RAST*RAST
00292         STI = -STI*RAST*RAST
00293         S1R = -ZET1DR + STR
00294         S1I = -ZET1DI + STI
00295         GO TO 210
00296   200   CONTINUE
00297         S1R = -ZET1DR + ZET2DR
00298         S1I = -ZET1DI + ZET2DI
00299   210   CONTINUE
00300 C-----------------------------------------------------------------------
00301 C     TEST FOR UNDERFLOW AND OVERFLOW
00302 C-----------------------------------------------------------------------
00303         RS1 = S1R
00304         IF (DABS(RS1).GT.ELIM) GO TO 260
00305         IF (KDFLG.EQ.1) IFLAG = 2
00306         IF (DABS(RS1).LT.ALIM) GO TO 220
00307 C-----------------------------------------------------------------------
00308 C     REFINE  TEST AND SCALE
00309 C-----------------------------------------------------------------------
00310         APHI = XZABS(PHIDR,PHIDI)
00311         RS1 = RS1 + DLOG(APHI)
00312         IF (DABS(RS1).GT.ELIM) GO TO 260
00313         IF (KDFLG.EQ.1) IFLAG = 1
00314         IF (RS1.LT.0.0D0) GO TO 220
00315         IF (KDFLG.EQ.1) IFLAG = 3
00316   220   CONTINUE
00317         STR = PHIDR*SUMDR - PHIDI*SUMDI
00318         STI = PHIDR*SUMDI + PHIDI*SUMDR
00319         S2R = -CSGNI*STI
00320         S2I = CSGNI*STR
00321         STR = DEXP(S1R)*CSSR(IFLAG)
00322         S1R = STR*DCOS(S1I)
00323         S1I = STR*DSIN(S1I)
00324         STR = S2R*S1R - S2I*S1I
00325         S2I = S2R*S1I + S2I*S1R
00326         S2R = STR
00327         IF (IFLAG.NE.1) GO TO 230
00328         CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
00329         IF (NW.EQ.0) GO TO 230
00330         S2R = ZEROR
00331         S2I = ZEROI
00332   230   CONTINUE
00333         CYR(KDFLG) = S2R
00334         CYI(KDFLG) = S2I
00335         C2R = S2R
00336         C2I = S2I
00337         S2R = S2R*CSRR(IFLAG)
00338         S2I = S2I*CSRR(IFLAG)
00339 C-----------------------------------------------------------------------
00340 C     ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
00341 C-----------------------------------------------------------------------
00342         S1R = YR(KK)
00343         S1I = YI(KK)
00344         IF (KODE.EQ.1) GO TO 250
00345         CALL ZS1S2(ZRR, ZRI, S1R, S1I, S2R, S2I, NW, ASC, ALIM, IUF)
00346         NZ = NZ + NW
00347   250   CONTINUE
00348         YR(KK) = S1R*CSPNR - S1I*CSPNI + S2R
00349         YI(KK) = CSPNR*S1I + CSPNI*S1R + S2I
00350         KK = KK - 1
00351         CSPNR = -CSPNR
00352         CSPNI = -CSPNI
00353         IF (C2R.NE.0.0D0 .OR. C2I.NE.0.0D0) GO TO 255
00354         KDFLG = 1
00355         GO TO 270
00356   255   CONTINUE
00357         IF (KDFLG.EQ.2) GO TO 275
00358         KDFLG = 2
00359         GO TO 270
00360   260   CONTINUE
00361         IF (RS1.GT.0.0D0) GO TO 300
00362         S2R = ZEROR
00363         S2I = ZEROI
00364         GO TO 230
00365   270 CONTINUE
00366       K = N
00367   275 CONTINUE
00368       IL = N - K
00369       IF (IL.EQ.0) RETURN
00370 C-----------------------------------------------------------------------
00371 C     RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
00372 C     K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
00373 C     INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
00374 C-----------------------------------------------------------------------
00375       S1R = CYR(1)
00376       S1I = CYI(1)
00377       S2R = CYR(2)
00378       S2I = CYI(2)
00379       CSR = CSRR(IFLAG)
00380       ASCLE = BRY(IFLAG)
00381       FN = DBLE(FLOAT(INU+IL))
00382       DO 290 I=1,IL
00383         C2R = S2R
00384         C2I = S2I
00385         S2R = S1R + (FN+FNF)*(RZR*C2R-RZI*C2I)
00386         S2I = S1I + (FN+FNF)*(RZR*C2I+RZI*C2R)
00387         S1R = C2R
00388         S1I = C2I
00389         FN = FN - 1.0D0
00390         C2R = S2R*CSR
00391         C2I = S2I*CSR
00392         CKR = C2R
00393         CKI = C2I
00394         C1R = YR(KK)
00395         C1I = YI(KK)
00396         IF (KODE.EQ.1) GO TO 280
00397         CALL ZS1S2(ZRR, ZRI, C1R, C1I, C2R, C2I, NW, ASC, ALIM, IUF)
00398         NZ = NZ + NW
00399   280   CONTINUE
00400         YR(KK) = C1R*CSPNR - C1I*CSPNI + C2R
00401         YI(KK) = C1R*CSPNI + C1I*CSPNR + C2I
00402         KK = KK - 1
00403         CSPNR = -CSPNR
00404         CSPNI = -CSPNI
00405         IF (IFLAG.GE.3) GO TO 290
00406         C2R = DABS(CKR)
00407         C2I = DABS(CKI)
00408         C2M = DMAX1(C2R,C2I)
00409         IF (C2M.LE.ASCLE) GO TO 290
00410         IFLAG = IFLAG + 1
00411         ASCLE = BRY(IFLAG)
00412         S1R = S1R*CSR
00413         S1I = S1I*CSR
00414         S2R = CKR
00415         S2I = CKI
00416         S1R = S1R*CSSR(IFLAG)
00417         S1I = S1I*CSSR(IFLAG)
00418         S2R = S2R*CSSR(IFLAG)
00419         S2I = S2I*CSSR(IFLAG)
00420         CSR = CSRR(IFLAG)
00421   290 CONTINUE
00422       RETURN
00423   300 CONTINUE
00424       NZ = -1
00425       RETURN
00426       END