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

Go to the documentation of this file.
00001       SUBROUTINE ZUNI2(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL,
00002      * TOL, ELIM, ALIM)
00003 C***BEGIN PROLOGUE  ZUNI2
00004 C***REFER TO  ZBESI,ZBESK
00005 C
00006 C     ZUNI2 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  ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,XZABS
00017 C***END PROLOGUE  ZUNI2
00018 C     COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS,
00019 C    *CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN
00020       DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGI,
00021      * ARGR, ASCLE, ASUMI, ASUMR, BRY, BSUMI, BSUMR, CIDI, CIPI, CIPR,
00022      * CONER, CRSC, CSCL, CSRR, CSSR, C1R, C2I, C2M, C2R, DAII,
00023      * DAIR, ELIM, FN, FNU, FNUL, HPI, PHII, PHIR, RAST, RAZ, RS1, RZI,
00024      * RZR, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZBI, ZBR, ZEROI,
00025      * ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI, ZNI, ZNR, ZR, CYR,
00026      * CYI, D1MACH, XZABS, CAR, SAR
00027       INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST,
00028      * NN, NUF, NW, NZ, IDUM
00029       DIMENSION BRY(3), YR(N), YI(N), CIPR(4), CIPI(4), CSSR(3),
00030      * CSRR(3), CYR(2), CYI(2)
00031       DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 /
00032       DATA CIPR(1),CIPI(1),CIPR(2),CIPI(2),CIPR(3),CIPI(3),CIPR(4),
00033      * CIPI(4)/ 1.0D0,0.0D0, 0.0D0,1.0D0, -1.0D0,0.0D0, 0.0D0,-1.0D0/
00034       DATA HPI, AIC  /
00035      1      1.57079632679489662D+00,     1.265512123484645396D+00/
00036 C
00037       NZ = 0
00038       ND = N
00039       NLAST = 0
00040 C-----------------------------------------------------------------------
00041 C     COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
00042 C     NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
00043 C     EXP(ALIM)=EXP(ELIM)*TOL
00044 C-----------------------------------------------------------------------
00045       CSCL = 1.0D0/TOL
00046       CRSC = TOL
00047       CSSR(1) = CSCL
00048       CSSR(2) = CONER
00049       CSSR(3) = CRSC
00050       CSRR(1) = CRSC
00051       CSRR(2) = CONER
00052       CSRR(3) = CSCL
00053       BRY(1) = 1.0D+3*D1MACH(1)/TOL
00054 C-----------------------------------------------------------------------
00055 C     ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI
00056 C-----------------------------------------------------------------------
00057       ZNR = ZI
00058       ZNI = -ZR
00059       ZBR = ZR
00060       ZBI = ZI
00061       CIDI = -CONER
00062       INU = INT(SNGL(FNU))
00063       ANG = HPI*(FNU-DBLE(FLOAT(INU)))
00064       C2R = DCOS(ANG)
00065       C2I = DSIN(ANG)
00066       CAR = C2R
00067       SAR = C2I
00068       IN = INU + N - 1
00069       IN = MOD(IN,4) + 1
00070       STR = C2R*CIPR(IN) - C2I*CIPI(IN)
00071       C2I = C2R*CIPI(IN) + C2I*CIPR(IN)
00072       C2R = STR
00073       IF (ZI.GT.0.0D0) GO TO 10
00074       ZNR = -ZNR
00075       ZBI = -ZBI
00076       CIDI = -CIDI
00077       C2I = -C2I
00078    10 CONTINUE
00079 C-----------------------------------------------------------------------
00080 C     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
00081 C-----------------------------------------------------------------------
00082       FN = DMAX1(FNU,1.0D0)
00083       CALL ZUNHJ(ZNR, ZNI, FN, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R,
00084      * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
00085       IF (KODE.EQ.1) GO TO 20
00086       STR = ZBR + ZETA2R
00087       STI = ZBI + ZETA2I
00088       RAST = FN/XZABS(STR,STI)
00089       STR = STR*RAST*RAST
00090       STI = -STI*RAST*RAST
00091       S1R = -ZETA1R + STR
00092       S1I = -ZETA1I + STI
00093       GO TO 30
00094    20 CONTINUE
00095       S1R = -ZETA1R + ZETA2R
00096       S1I = -ZETA1I + ZETA2I
00097    30 CONTINUE
00098       RS1 = S1R
00099       IF (DABS(RS1).GT.ELIM) GO TO 150
00100    40 CONTINUE
00101       NN = MIN0(2,ND)
00102       DO 90 I=1,NN
00103         FN = FNU + DBLE(FLOAT(ND-I))
00104         CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIR, PHII, ARGR, ARGI,
00105      *   ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
00106         IF (KODE.EQ.1) GO TO 50
00107         STR = ZBR + ZETA2R
00108         STI = ZBI + ZETA2I
00109         RAST = FN/XZABS(STR,STI)
00110         STR = STR*RAST*RAST
00111         STI = -STI*RAST*RAST
00112         S1R = -ZETA1R + STR
00113         S1I = -ZETA1I + STI + DABS(ZI)
00114         GO TO 60
00115    50   CONTINUE
00116         S1R = -ZETA1R + ZETA2R
00117         S1I = -ZETA1I + ZETA2I
00118    60   CONTINUE
00119 C-----------------------------------------------------------------------
00120 C     TEST FOR UNDERFLOW AND OVERFLOW
00121 C-----------------------------------------------------------------------
00122         RS1 = S1R
00123         IF (DABS(RS1).GT.ELIM) GO TO 120
00124         IF (I.EQ.1) IFLAG = 2
00125         IF (DABS(RS1).LT.ALIM) GO TO 70
00126 C-----------------------------------------------------------------------
00127 C     REFINE  TEST AND SCALE
00128 C-----------------------------------------------------------------------
00129 C-----------------------------------------------------------------------
00130         APHI = XZABS(PHIR,PHII)
00131         AARG = XZABS(ARGR,ARGI)
00132         RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC
00133         IF (DABS(RS1).GT.ELIM) GO TO 120
00134         IF (I.EQ.1) IFLAG = 1
00135         IF (RS1.LT.0.0D0) GO TO 70
00136         IF (I.EQ.1) IFLAG = 3
00137    70   CONTINUE
00138 C-----------------------------------------------------------------------
00139 C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
00140 C     EXPONENT EXTREMES
00141 C-----------------------------------------------------------------------
00142         CALL ZAIRY(ARGR, ARGI, 0, 2, AIR, AII, NAI, IDUM)
00143         CALL ZAIRY(ARGR, ARGI, 1, 2, DAIR, DAII, NDAI, IDUM)
00144         STR = DAIR*BSUMR - DAII*BSUMI
00145         STI = DAIR*BSUMI + DAII*BSUMR
00146         STR = STR + (AIR*ASUMR-AII*ASUMI)
00147         STI = STI + (AIR*ASUMI+AII*ASUMR)
00148         S2R = PHIR*STR - PHII*STI
00149         S2I = PHIR*STI + PHII*STR
00150         STR = DEXP(S1R)*CSSR(IFLAG)
00151         S1R = STR*DCOS(S1I)
00152         S1I = STR*DSIN(S1I)
00153         STR = S2R*S1R - S2I*S1I
00154         S2I = S2R*S1I + S2I*S1R
00155         S2R = STR
00156         IF (IFLAG.NE.1) GO TO 80
00157         CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
00158         IF (NW.NE.0) GO TO 120
00159    80   CONTINUE
00160         IF (ZI.LE.0.0D0) S2I = -S2I
00161         STR = S2R*C2R - S2I*C2I
00162         S2I = S2R*C2I + S2I*C2R
00163         S2R = STR
00164         CYR(I) = S2R
00165         CYI(I) = S2I
00166         J = ND - I + 1
00167         YR(J) = S2R*CSRR(IFLAG)
00168         YI(J) = S2I*CSRR(IFLAG)
00169         STR = -C2I*CIDI
00170         C2I = C2R*CIDI
00171         C2R = STR
00172    90 CONTINUE
00173       IF (ND.LE.2) GO TO 110
00174       RAZ = 1.0D0/XZABS(ZR,ZI)
00175       STR = ZR*RAZ
00176       STI = -ZI*RAZ
00177       RZR = (STR+STR)*RAZ
00178       RZI = (STI+STI)*RAZ
00179       BRY(2) = 1.0D0/BRY(1)
00180       BRY(3) = D1MACH(2)
00181       S1R = CYR(1)
00182       S1I = CYI(1)
00183       S2R = CYR(2)
00184       S2I = CYI(2)
00185       C1R = CSRR(IFLAG)
00186       ASCLE = BRY(IFLAG)
00187       K = ND - 2
00188       FN = DBLE(FLOAT(K))
00189       DO 100 I=3,ND
00190         C2R = S2R
00191         C2I = S2I
00192         S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I)
00193         S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R)
00194         S1R = C2R
00195         S1I = C2I
00196         C2R = S2R*C1R
00197         C2I = S2I*C1R
00198         YR(K) = C2R
00199         YI(K) = C2I
00200         K = K - 1
00201         FN = FN - 1.0D0
00202         IF (IFLAG.GE.3) GO TO 100
00203         STR = DABS(C2R)
00204         STI = DABS(C2I)
00205         C2M = DMAX1(STR,STI)
00206         IF (C2M.LE.ASCLE) GO TO 100
00207         IFLAG = IFLAG + 1
00208         ASCLE = BRY(IFLAG)
00209         S1R = S1R*C1R
00210         S1I = S1I*C1R
00211         S2R = C2R
00212         S2I = C2I
00213         S1R = S1R*CSSR(IFLAG)
00214         S1I = S1I*CSSR(IFLAG)
00215         S2R = S2R*CSSR(IFLAG)
00216         S2I = S2I*CSSR(IFLAG)
00217         C1R = CSRR(IFLAG)
00218   100 CONTINUE
00219   110 CONTINUE
00220       RETURN
00221   120 CONTINUE
00222       IF (RS1.GT.0.0D0) GO TO 140
00223 C-----------------------------------------------------------------------
00224 C     SET UNDERFLOW AND UPDATE PARAMETERS
00225 C-----------------------------------------------------------------------
00226       YR(ND) = ZEROR
00227       YI(ND) = ZEROI
00228       NZ = NZ + 1
00229       ND = ND - 1
00230       IF (ND.EQ.0) GO TO 110
00231       CALL ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, NUF, TOL, ELIM, ALIM)
00232       IF (NUF.LT.0) GO TO 140
00233       ND = ND - NUF
00234       NZ = NZ + NUF
00235       IF (ND.EQ.0) GO TO 110
00236       FN = FNU + DBLE(FLOAT(ND-1))
00237       IF (FN.LT.FNUL) GO TO 130
00238 C      FN = CIDI
00239 C      J = NUF + 1
00240 C      K = MOD(J,4) + 1
00241 C      S1R = CIPR(K)
00242 C      S1I = CIPI(K)
00243 C      IF (FN.LT.0.0D0) S1I = -S1I
00244 C      STR = C2R*S1R - C2I*S1I
00245 C      C2I = C2R*S1I + C2I*S1R
00246 C      C2R = STR
00247       IN = INU + ND - 1
00248       IN = MOD(IN,4) + 1
00249       C2R = CAR*CIPR(IN) - SAR*CIPI(IN)
00250       C2I = CAR*CIPI(IN) + SAR*CIPR(IN)
00251       IF (ZI.LE.0.0D0) C2I = -C2I
00252       GO TO 40
00253   130 CONTINUE
00254       NLAST = ND
00255       RETURN
00256   140 CONTINUE
00257       NZ = -1
00258       RETURN
00259   150 CONTINUE
00260       IF (RS1.GT.0.0D0) GO TO 140
00261       NZ = N
00262       DO 160 I=1,N
00263         YR(I) = ZEROR
00264         YI(I) = ZEROI
00265   160 CONTINUE
00266       RETURN
00267       END