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

Go to the documentation of this file.
00001       SUBROUTINE CUNI2(Z, FNU, KODE, N, Y, NZ, NLAST, FNUL, TOL, ELIM,
00002      * ALIM)
00003 C***BEGIN PROLOGUE  CUNI2
00004 C***REFER TO  CBESI,CBESK
00005 C
00006 C     CUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF
00007 C     UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I
00008 C     OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO.
00009 C
00010 C     FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
00011 C     EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
00012 C     NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
00013 C     FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
00014 C     Y(I)=CZERO FOR I=NLAST+1,N
00015 C
00016 C***ROUTINES CALLED  CAIRY,CUCHK,CUNHJ,CUOIK,R1MACH
00017 C***END PROLOGUE  CUNI2
00018       COMPLEX AI, ARG, ASUM, BSUM, CFN, CI, CID, CIP, CONE, CRSC, CSCL,
00019      * CSR, CSS, CY, CZERO, C1, C2, DAI, PHI, RZ, S1, S2, Y, Z, ZB,
00020      * ZETA1, ZETA2, ZN, ZAR
00021       REAL AARG, AIC, ALIM, ANG, APHI, ASCLE, AY, BRY, CAR, C2I, C2M,
00022      * C2R, ELIM, FN, FNU, FNUL, HPI, RS1, SAR, TOL, YY, R1MACH
00023       INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST,
00024      * NN, NUF, NW, NZ, IDUM
00025       DIMENSION BRY(3), Y(N), CIP(4), CSS(3), CSR(3), CY(2)
00026       DATA CZERO,CONE,CI/(0.0E0,0.0E0),(1.0E0,0.0E0),(0.0E0,1.0E0)/
00027       DATA CIP(1),CIP(2),CIP(3),CIP(4)/
00028      1 (1.0E0,0.0E0), (0.0E0,1.0E0), (-1.0E0,0.0E0), (0.0E0,-1.0E0)/
00029       DATA HPI, AIC  /
00030      1      1.57079632679489662E+00,     1.265512123484645396E+00/
00031 C
00032       NZ = 0
00033       ND = N
00034       NLAST = 0
00035 C-----------------------------------------------------------------------
00036 C     COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
00037 C     NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
00038 C     EXP(ALIM)=EXP(ELIM)*TOL
00039 C-----------------------------------------------------------------------
00040       CSCL = CMPLX(1.0E0/TOL,0.0E0)
00041       CRSC = CMPLX(TOL,0.0E0)
00042       CSS(1) = CSCL
00043       CSS(2) = CONE
00044       CSS(3) = CRSC
00045       CSR(1) = CRSC
00046       CSR(2) = CONE
00047       CSR(3) = CSCL
00048       BRY(1) = 1.0E+3*R1MACH(1)/TOL
00049       YY = AIMAG(Z)
00050 C-----------------------------------------------------------------------
00051 C     ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI
00052 C-----------------------------------------------------------------------
00053       ZN = -Z*CI
00054       ZB = Z
00055       CID = -CI
00056       INU = INT(FNU)
00057       ANG = HPI*(FNU-FLOAT(INU))
00058       CAR = COS(ANG)
00059       SAR = SIN(ANG)
00060       C2 = CMPLX(CAR,SAR)
00061       ZAR = C2
00062       IN = INU + N - 1
00063       IN = MOD(IN,4)
00064       C2 = C2*CIP(IN+1)
00065       IF (YY.GT.0.0E0) GO TO 10
00066       ZN = CONJG(-ZN)
00067       ZB = CONJG(ZB)
00068       CID = -CID
00069       C2 = CONJG(C2)
00070    10 CONTINUE
00071 C-----------------------------------------------------------------------
00072 C     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
00073 C-----------------------------------------------------------------------
00074       FN = AMAX1(FNU,1.0E0)
00075       CALL CUNHJ(ZN, FN, 1, TOL, PHI, ARG, ZETA1, ZETA2, ASUM, BSUM)
00076       IF (KODE.EQ.1) GO TO 20
00077       CFN = CMPLX(FNU,0.0E0)
00078       S1 = -ZETA1 + CFN*(CFN/(ZB+ZETA2))
00079       GO TO 30
00080    20 CONTINUE
00081       S1 = -ZETA1 + ZETA2
00082    30 CONTINUE
00083       RS1 = REAL(S1)
00084       IF (ABS(RS1).GT.ELIM) GO TO 150
00085    40 CONTINUE
00086       NN = MIN0(2,ND)
00087       DO 90 I=1,NN
00088         FN = FNU + FLOAT(ND-I)
00089         CALL CUNHJ(ZN, FN, 0, TOL, PHI, ARG, ZETA1, ZETA2, ASUM, BSUM)
00090         IF (KODE.EQ.1) GO TO 50
00091         CFN = CMPLX(FN,0.0E0)
00092         AY = ABS(YY)
00093         S1 = -ZETA1 + CFN*(CFN/(ZB+ZETA2)) + CMPLX(0.0E0,AY)
00094         GO TO 60
00095    50   CONTINUE
00096         S1 = -ZETA1 + ZETA2
00097    60   CONTINUE
00098 C-----------------------------------------------------------------------
00099 C     TEST FOR UNDERFLOW AND OVERFLOW
00100 C-----------------------------------------------------------------------
00101         RS1 = REAL(S1)
00102         IF (ABS(RS1).GT.ELIM) GO TO 120
00103         IF (I.EQ.1) IFLAG = 2
00104         IF (ABS(RS1).LT.ALIM) GO TO 70
00105 C-----------------------------------------------------------------------
00106 C     REFINE  TEST AND SCALE
00107 C-----------------------------------------------------------------------
00108 C-----------------------------------------------------------------------
00109         APHI = CABS(PHI)
00110         AARG = CABS(ARG)
00111         RS1 = RS1 + ALOG(APHI) - 0.25E0*ALOG(AARG) - AIC
00112         IF (ABS(RS1).GT.ELIM) GO TO 120
00113         IF (I.EQ.1) IFLAG = 1
00114         IF (RS1.LT.0.0E0) GO TO 70
00115         IF (I.EQ.1) IFLAG = 3
00116    70   CONTINUE
00117 C-----------------------------------------------------------------------
00118 C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
00119 C     EXPONENT EXTREMES
00120 C-----------------------------------------------------------------------
00121         CALL CAIRY(ARG, 0, 2, AI, NAI, IDUM)
00122         CALL CAIRY(ARG, 1, 2, DAI, NDAI, IDUM)
00123         S2 = PHI*(AI*ASUM+DAI*BSUM)
00124         C2R = REAL(S1)
00125         C2I = AIMAG(S1)
00126         C2M = EXP(C2R)*REAL(CSS(IFLAG))
00127         S1 = CMPLX(C2M,0.0E0)*CMPLX(COS(C2I),SIN(C2I))
00128         S2 = S2*S1
00129         IF (IFLAG.NE.1) GO TO 80
00130         CALL CUCHK(S2, NW, BRY(1), TOL)
00131         IF (NW.NE.0) GO TO 120
00132    80   CONTINUE
00133         IF (YY.LE.0.0E0) S2 = CONJG(S2)
00134         J = ND - I + 1
00135         S2 = S2*C2
00136         CY(I) = S2
00137         Y(J) = S2*CSR(IFLAG)
00138         C2 = C2*CID
00139    90 CONTINUE
00140       IF (ND.LE.2) GO TO 110
00141       RZ = CMPLX(2.0E0,0.0E0)/Z
00142       BRY(2) = 1.0E0/BRY(1)
00143       BRY(3) = R1MACH(2)
00144       S1 = CY(1)
00145       S2 = CY(2)
00146       C1 = CSR(IFLAG)
00147       ASCLE = BRY(IFLAG)
00148       K = ND - 2
00149       FN = FLOAT(K)
00150       DO 100 I=3,ND
00151         C2 = S2
00152         S2 = S1 + CMPLX(FNU+FN,0.0E0)*RZ*S2
00153         S1 = C2
00154         C2 = S2*C1
00155         Y(K) = C2
00156         K = K - 1
00157         FN = FN - 1.0E0
00158         IF (IFLAG.GE.3) GO TO 100
00159         C2R = REAL(C2)
00160         C2I = AIMAG(C2)
00161         C2R = ABS(C2R)
00162         C2I = ABS(C2I)
00163         C2M = AMAX1(C2R,C2I)
00164         IF (C2M.LE.ASCLE) GO TO 100
00165         IFLAG = IFLAG + 1
00166         ASCLE = BRY(IFLAG)
00167         S1 = S1*C1
00168         S2 = C2
00169         S1 = S1*CSS(IFLAG)
00170         S2 = S2*CSS(IFLAG)
00171         C1 = CSR(IFLAG)
00172   100 CONTINUE
00173   110 CONTINUE
00174       RETURN
00175   120 CONTINUE
00176       IF (RS1.GT.0.0E0) GO TO 140
00177 C-----------------------------------------------------------------------
00178 C     SET UNDERFLOW AND UPDATE PARAMETERS
00179 C-----------------------------------------------------------------------
00180       Y(ND) = CZERO
00181       NZ = NZ + 1
00182       ND = ND - 1
00183       IF (ND.EQ.0) GO TO 110
00184       CALL CUOIK(Z, FNU, KODE, 1, ND, Y, NUF, TOL, ELIM, ALIM)
00185       IF (NUF.LT.0) GO TO 140
00186       ND = ND - NUF
00187       NZ = NZ + NUF
00188       IF (ND.EQ.0) GO TO 110
00189       FN = FNU + FLOAT(ND-1)
00190       IF (FN.LT.FNUL) GO TO 130
00191 C      FN = AIMAG(CID)
00192 C      J = NUF + 1
00193 C      K = MOD(J,4) + 1
00194 C      S1 = CIP(K)
00195 C      IF (FN.LT.0.0E0) S1 = CONJG(S1)
00196 C      C2 = C2*S1
00197       IN = INU + ND - 1
00198       IN = MOD(IN,4) + 1
00199       C2 = ZAR*CIP(IN)
00200       IF (YY.LE.0.0E0)C2=CONJG(C2)
00201       GO TO 40
00202   130 CONTINUE
00203       NLAST = ND
00204       RETURN
00205   140 CONTINUE
00206       NZ = -1
00207       RETURN
00208   150 CONTINUE
00209       IF (RS1.GT.0.0E0) GO TO 140
00210       NZ = N
00211       DO 160 I=1,N
00212         Y(I) = CZERO
00213   160 CONTINUE
00214       RETURN
00215       END