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

Go to the documentation of this file.
00001       SUBROUTINE CUNK1(Z, FNU, KODE, MR, N, Y, NZ, TOL, ELIM, ALIM)
00002 C***BEGIN PROLOGUE  CUNK1
00003 C***REFER TO  CBESK
00004 C
00005 C     CUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
00006 C     RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
00007 C     UNIFORM ASYMPTOTIC EXPANSION.
00008 C     MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
00009 C     NZ=-1 MEANS AN OVERFLOW WILL OCCUR
00010 C
00011 C***ROUTINES CALLED  CS1S2,CUCHK,CUNIK,R1MACH
00012 C***END PROLOGUE  CUNK1
00013       COMPLEX CFN, CK, CONE, CRSC, CS, CSCL, CSGN, CSPN, CSR, CSS,
00014      * CWRK, CY, CZERO, C1, C2, PHI,  RZ, SUM,  S1, S2, Y, Z,
00015      * ZETA1,  ZETA2,  ZR, PHID, ZETA1D, ZETA2D, SUMD
00016       REAL ALIM, ANG, APHI, ASC, ASCLE, BRY, CPN, C2I, C2M, C2R, ELIM,
00017      * FMR, FN, FNF, FNU, PI, RS1, SGN, SPN, TOL, X, R1MACH
00018       INTEGER I, IB, IFLAG, IFN, IL, INIT, INU, IUF, K, KDFLG, KFLAG,
00019      * KK, KODE, MR, N, NW, NZ, J, IPARD, INITD, IC
00020       DIMENSION BRY(3), INIT(2), Y(N), SUM(2), PHI(2), ZETA1(2),
00021      * ZETA2(2), CY(2), CWRK(16,3), CSS(3), CSR(3)
00022       DATA CZERO, CONE / (0.0E0,0.0E0) , (1.0E0,0.0E0) /
00023       DATA PI / 3.14159265358979324E0 /
00024 C
00025       KDFLG = 1
00026       NZ = 0
00027 C-----------------------------------------------------------------------
00028 C     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
00029 C     THE UNDERFLOW LIMIT
00030 C-----------------------------------------------------------------------
00031       CSCL = CMPLX(1.0E0/TOL,0.0E0)
00032       CRSC = CMPLX(TOL,0.0E0)
00033       CSS(1) = CSCL
00034       CSS(2) = CONE
00035       CSS(3) = CRSC
00036       CSR(1) = CRSC
00037       CSR(2) = CONE
00038       CSR(3) = CSCL
00039       BRY(1) = 1.0E+3*R1MACH(1)/TOL
00040       BRY(2) = 1.0E0/BRY(1)
00041       BRY(3) = R1MACH(2)
00042       X = REAL(Z)
00043       ZR = Z
00044       IF (X.LT.0.0E0) ZR = -Z
00045       J=2
00046       DO 70 I=1,N
00047 C-----------------------------------------------------------------------
00048 C     J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
00049 C-----------------------------------------------------------------------
00050         J = 3 - J
00051         FN = FNU + FLOAT(I-1)
00052         INIT(J) = 0
00053         CALL CUNIK(ZR, FN, 2, 0, TOL, INIT(J), PHI(J), ZETA1(J),
00054      *   ZETA2(J), SUM(J), CWRK(1,J))
00055         IF (KODE.EQ.1) GO TO 20
00056         CFN = CMPLX(FN,0.0E0)
00057         S1 = ZETA1(J) - CFN*(CFN/(ZR+ZETA2(J)))
00058         GO TO 30
00059    20   CONTINUE
00060         S1 = ZETA1(J) - ZETA2(J)
00061    30   CONTINUE
00062 C-----------------------------------------------------------------------
00063 C     TEST FOR UNDERFLOW AND OVERFLOW
00064 C-----------------------------------------------------------------------
00065         RS1 = REAL(S1)
00066         IF (ABS(RS1).GT.ELIM) GO TO 60
00067         IF (KDFLG.EQ.1) KFLAG = 2
00068         IF (ABS(RS1).LT.ALIM) GO TO 40
00069 C-----------------------------------------------------------------------
00070 C     REFINE  TEST AND SCALE
00071 C-----------------------------------------------------------------------
00072         APHI = CABS(PHI(J))
00073         RS1 = RS1 + ALOG(APHI)
00074         IF (ABS(RS1).GT.ELIM) GO TO 60
00075         IF (KDFLG.EQ.1) KFLAG = 1
00076         IF (RS1.LT.0.0E0) GO TO 40
00077         IF (KDFLG.EQ.1) KFLAG = 3
00078    40   CONTINUE
00079 C-----------------------------------------------------------------------
00080 C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
00081 C     EXPONENT EXTREMES
00082 C-----------------------------------------------------------------------
00083         S2 = PHI(J)*SUM(J)
00084         C2R = REAL(S1)
00085         C2I = AIMAG(S1)
00086         C2M = EXP(C2R)*REAL(CSS(KFLAG))
00087         S1 = CMPLX(C2M,0.0E0)*CMPLX(COS(C2I),SIN(C2I))
00088         S2 = S2*S1
00089         IF (KFLAG.NE.1) GO TO 50
00090         CALL CUCHK(S2, NW, BRY(1), TOL)
00091         IF (NW.NE.0) GO TO 60
00092    50   CONTINUE
00093         CY(KDFLG) = S2
00094         Y(I) = S2*CSR(KFLAG)
00095         IF (KDFLG.EQ.2) GO TO 75
00096         KDFLG = 2
00097         GO TO 70
00098    60   CONTINUE
00099         IF (RS1.GT.0.0E0) GO TO 290
00100 C-----------------------------------------------------------------------
00101 C     FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
00102 C-----------------------------------------------------------------------
00103         IF (X.LT.0.0E0) GO TO 290
00104         KDFLG = 1
00105         Y(I) = CZERO
00106         NZ=NZ+1
00107         IF (I.EQ.1) GO TO 70
00108         IF (Y(I-1).EQ.CZERO) GO TO 70
00109         Y(I-1) = CZERO
00110         NZ=NZ+1
00111    70 CONTINUE
00112       I=N
00113    75 CONTINUE
00114       RZ = CMPLX(2.0E0,0.0E0)/ZR
00115       CK = CMPLX(FN,0.0E0)*RZ
00116       IB = I+1
00117       IF (N.LT.IB) GO TO 160
00118 C-----------------------------------------------------------------------
00119 C     TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW, SET SEQUENCE TO ZERO
00120 C     ON UNDERFLOW
00121 C-----------------------------------------------------------------------
00122       FN = FNU+FLOAT(N-1)
00123       IPARD = 1
00124       IF (MR.NE.0) IPARD = 0
00125       INITD = 0
00126       CALL CUNIK(ZR,FN,2,IPARD,TOL,INITD,PHID,ZETA1D,ZETA2D,SUMD,
00127      *CWRK(1,3))
00128       IF (KODE.EQ.1) GO TO 80
00129       CFN=CMPLX(FN,0.0E0)
00130       S1=ZETA1D-CFN*(CFN/(ZR+ZETA2D))
00131       GO TO 90
00132    80 CONTINUE
00133       S1=ZETA1D-ZETA2D
00134    90 CONTINUE
00135       RS1=REAL(S1)
00136       IF (ABS(RS1).GT.ELIM) GO TO 95
00137       IF (ABS(RS1).LT.ALIM) GO TO 100
00138 C-----------------------------------------------------------------------
00139 C     REFINE ESTIMATE AND TEST
00140 C-----------------------------------------------------------------------
00141       APHI=CABS(PHID)
00142       RS1=RS1+ALOG(APHI)
00143       IF (ABS(RS1).LT.ELIM) GO TO 100
00144    95 CONTINUE
00145       IF (RS1.GT.0.0E0) GO TO 290
00146 C-----------------------------------------------------------------------
00147 C     FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
00148 C-----------------------------------------------------------------------
00149       IF (X.LT.0.0E0) GO TO 290
00150       NZ=N
00151       DO 96 I=1,N
00152         Y(I) = CZERO
00153    96 CONTINUE
00154       RETURN
00155   100 CONTINUE
00156 C-----------------------------------------------------------------------
00157 C     RECUR FORWARD FOR REMAINDER OF THE SEQUENCE
00158 C-----------------------------------------------------------------------
00159       S1 = CY(1)
00160       S2 = CY(2)
00161       C1 = CSR(KFLAG)
00162       ASCLE = BRY(KFLAG)
00163       DO 120 I=IB,N
00164         C2 = S2
00165         S2 = CK*S2 + S1
00166         S1 = C2
00167         CK = CK + RZ
00168         C2 = S2*C1
00169         Y(I) = C2
00170         IF (KFLAG.GE.3) GO TO 120
00171         C2R = REAL(C2)
00172         C2I = AIMAG(C2)
00173         C2R = ABS(C2R)
00174         C2I = ABS(C2I)
00175         C2M = AMAX1(C2R,C2I)
00176         IF (C2M.LE.ASCLE) GO TO 120
00177         KFLAG = KFLAG + 1
00178         ASCLE = BRY(KFLAG)
00179         S1 = S1*C1
00180         S2 = C2
00181         S1 = S1*CSS(KFLAG)
00182         S2 = S2*CSS(KFLAG)
00183         C1 = CSR(KFLAG)
00184   120 CONTINUE
00185   160 CONTINUE
00186       IF (MR.EQ.0) RETURN
00187 C-----------------------------------------------------------------------
00188 C     ANALYTIC CONTINUATION FOR RE(Z).LT.0.0E0
00189 C-----------------------------------------------------------------------
00190       NZ = 0
00191       FMR = FLOAT(MR)
00192       SGN = -SIGN(PI,FMR)
00193 C-----------------------------------------------------------------------
00194 C     CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP.
00195 C-----------------------------------------------------------------------
00196       CSGN = CMPLX(0.0E0,SGN)
00197       INU = INT(FNU)
00198       FNF = FNU - FLOAT(INU)
00199       IFN = INU + N - 1
00200       ANG = FNF*SGN
00201       CPN = COS(ANG)
00202       SPN = SIN(ANG)
00203       CSPN = CMPLX(CPN,SPN)
00204       IF (MOD(IFN,2).EQ.1) CSPN = -CSPN
00205       ASC = BRY(1)
00206       KK = N
00207       IUF = 0
00208       KDFLG = 1
00209       IB = IB-1
00210       IC = IB-1
00211       DO 260 K=1,N
00212         FN = FNU + FLOAT(KK-1)
00213 C-----------------------------------------------------------------------
00214 C     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
00215 C     FUNCTION ABOVE
00216 C-----------------------------------------------------------------------
00217         M=3
00218         IF (N.GT.2) GO TO 175
00219   170   CONTINUE
00220         INITD = INIT(J)
00221         PHID = PHI(J)
00222         ZETA1D = ZETA1(J)
00223         ZETA2D = ZETA2(J)
00224         SUMD = SUM(J)
00225         M = J
00226         J = 3 - J
00227         GO TO 180
00228   175   CONTINUE
00229         IF ((KK.EQ.N).AND.(IB.LT.N)) GO TO 180
00230         IF ((KK.EQ.IB).OR.(KK.EQ.IC)) GO TO 170
00231         INITD = 0
00232   180   CONTINUE
00233         CALL CUNIK(ZR, FN, 1, 0, TOL, INITD, PHID, ZETA1D,
00234      *   ZETA2D, SUMD, CWRK(1,M))
00235         IF (KODE.EQ.1) GO TO 190
00236         CFN = CMPLX(FN,0.0E0)
00237         S1 = -ZETA1D + CFN*(CFN/(ZR+ZETA2D))
00238         GO TO 200
00239   190   CONTINUE
00240         S1 = -ZETA1D + ZETA2D
00241   200   CONTINUE
00242 C-----------------------------------------------------------------------
00243 C     TEST FOR UNDERFLOW AND OVERFLOW
00244 C-----------------------------------------------------------------------
00245         RS1 = REAL(S1)
00246         IF (ABS(RS1).GT.ELIM) GO TO 250
00247         IF (KDFLG.EQ.1) IFLAG = 2
00248         IF (ABS(RS1).LT.ALIM) GO TO 210
00249 C-----------------------------------------------------------------------
00250 C     REFINE  TEST AND SCALE
00251 C-----------------------------------------------------------------------
00252         APHI = CABS(PHID)
00253         RS1 = RS1 + ALOG(APHI)
00254         IF (ABS(RS1).GT.ELIM) GO TO 250
00255         IF (KDFLG.EQ.1) IFLAG = 1
00256         IF (RS1.LT.0.0E0) GO TO 210
00257         IF (KDFLG.EQ.1) IFLAG = 3
00258   210   CONTINUE
00259         S2 = CSGN*PHID*SUMD
00260         C2R = REAL(S1)
00261         C2I = AIMAG(S1)
00262         C2M = EXP(C2R)*REAL(CSS(IFLAG))
00263         S1 = CMPLX(C2M,0.0E0)*CMPLX(COS(C2I),SIN(C2I))
00264         S2 = S2*S1
00265         IF (IFLAG.NE.1) GO TO 220
00266         CALL CUCHK(S2, NW, BRY(1), TOL)
00267         IF (NW.NE.0) S2 = CMPLX(0.0E0,0.0E0)
00268   220   CONTINUE
00269         CY(KDFLG) = S2
00270         C2 = S2
00271         S2 = S2*CSR(IFLAG)
00272 C-----------------------------------------------------------------------
00273 C     ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
00274 C-----------------------------------------------------------------------
00275         S1 = Y(KK)
00276         IF (KODE.EQ.1) GO TO 240
00277         CALL CS1S2(ZR, S1, S2, NW, ASC, ALIM, IUF)
00278         NZ = NZ + NW
00279   240   CONTINUE
00280         Y(KK) = S1*CSPN + S2
00281         KK = KK - 1
00282         CSPN = -CSPN
00283         IF (C2.NE.CZERO) GO TO 245
00284         KDFLG = 1
00285         GO TO 260
00286   245   CONTINUE
00287         IF (KDFLG.EQ.2) GO TO 265
00288         KDFLG = 2
00289         GO TO 260
00290   250   CONTINUE
00291         IF (RS1.GT.0.0E0) GO TO 290
00292         S2 = CZERO
00293         GO TO 220
00294   260 CONTINUE
00295       K = N
00296   265 CONTINUE
00297       IL = N - K
00298       IF (IL.EQ.0) RETURN
00299 C-----------------------------------------------------------------------
00300 C     RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
00301 C     K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
00302 C     INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
00303 C-----------------------------------------------------------------------
00304       S1 = CY(1)
00305       S2 = CY(2)
00306       CS = CSR(IFLAG)
00307       ASCLE = BRY(IFLAG)
00308       FN = FLOAT(INU+IL)
00309       DO 280 I=1,IL
00310         C2 = S2
00311         S2 = S1 + CMPLX(FN+FNF,0.0E0)*RZ*S2
00312         S1 = C2
00313         FN = FN - 1.0E0
00314         C2 = S2*CS
00315         CK = C2
00316         C1 = Y(KK)
00317         IF (KODE.EQ.1) GO TO 270
00318         CALL CS1S2(ZR, C1, C2, NW, ASC, ALIM, IUF)
00319         NZ = NZ + NW
00320   270   CONTINUE
00321         Y(KK) = C1*CSPN + C2
00322         KK = KK - 1
00323         CSPN = -CSPN
00324         IF (IFLAG.GE.3) GO TO 280
00325         C2R = REAL(CK)
00326         C2I = AIMAG(CK)
00327         C2R = ABS(C2R)
00328         C2I = ABS(C2I)
00329         C2M = AMAX1(C2R,C2I)
00330         IF (C2M.LE.ASCLE) GO TO 280
00331         IFLAG = IFLAG + 1
00332         ASCLE = BRY(IFLAG)
00333         S1 = S1*CS
00334         S2 = CK
00335         S1 = S1*CSS(IFLAG)
00336         S2 = S2*CSS(IFLAG)
00337         CS = CSR(IFLAG)
00338   280 CONTINUE
00339       RETURN
00340   290 CONTINUE
00341       NZ = -1
00342       RETURN
00343       END