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

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