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

Go to the documentation of this file.
00001       SUBROUTINE ZUNI1(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL,
00002      * TOL, ELIM, ALIM)
00003 C***BEGIN PROLOGUE  ZUNI1
00004 C***REFER TO  ZBESI,ZBESK
00005 C
00006 C     ZUNI1 COMPUTES I(FNU,Z)  BY MEANS OF THE UNIFORM ASYMPTOTIC
00007 C     EXPANSION FOR I(FNU,Z) IN -PI/3.LE.ARG Z.LE.PI/3.
00008 C
00009 C     FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
00010 C     EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
00011 C     NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
00012 C     FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
00013 C     Y(I)=CZERO FOR I=NLAST+1,N
00014 C
00015 C***ROUTINES CALLED  ZUCHK,ZUNIK,ZUOIK,D1MACH,XZABS
00016 C***END PROLOGUE  ZUNI1
00017 C     COMPLEX CFN,CONE,CRSC,CSCL,CSR,CSS,CWRK,CZERO,C1,C2,PHI,RZ,SUM,S1,
00018 C    *S2,Y,Z,ZETA1,ZETA2
00019       DOUBLE PRECISION ALIM, APHI, ASCLE, BRY, CONER, CRSC,
00020      * CSCL, CSRR, CSSR, CWRKI, CWRKR, C1R, C2I, C2M, C2R, ELIM, FN,
00021      * FNU, FNUL, PHII, PHIR, RAST, RS1, RZI, RZR, STI, STR, SUMI,
00022      * SUMR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZEROI, ZEROR, ZETA1I,
00023      * ZETA1R, ZETA2I, ZETA2R, ZI, ZR, CYR, CYI, D1MACH, XZABS
00024       INTEGER I, IFLAG, INIT, K, KODE, M, N, ND, NLAST, NN, NUF, NW, NZ
00025       DIMENSION BRY(3), YR(N), YI(N), CWRKR(16), CWRKI(16), CSSR(3),
00026      * CSRR(3), CYR(2), CYI(2)
00027       DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 /
00028 C
00029       NZ = 0
00030       ND = N
00031       NLAST = 0
00032 C-----------------------------------------------------------------------
00033 C     COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
00034 C     NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
00035 C     EXP(ALIM)=EXP(ELIM)*TOL
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 C-----------------------------------------------------------------------
00047 C     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
00048 C-----------------------------------------------------------------------
00049       FN = DMAX1(FNU,1.0D0)
00050       INIT = 0
00051       CALL ZUNIK(ZR, ZI, FN, 1, 1, TOL, INIT, PHIR, PHII, ZETA1R,
00052      * ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
00053       IF (KODE.EQ.1) GO TO 10
00054       STR = ZR + ZETA2R
00055       STI = ZI + ZETA2I
00056       RAST = FN/XZABS(STR,STI)
00057       STR = STR*RAST*RAST
00058       STI = -STI*RAST*RAST
00059       S1R = -ZETA1R + STR
00060       S1I = -ZETA1I + STI
00061       GO TO 20
00062    10 CONTINUE
00063       S1R = -ZETA1R + ZETA2R
00064       S1I = -ZETA1I + ZETA2I
00065    20 CONTINUE
00066       RS1 = S1R
00067       IF (DABS(RS1).GT.ELIM) GO TO 130
00068    30 CONTINUE
00069       NN = MIN0(2,ND)
00070       DO 80 I=1,NN
00071         FN = FNU + DBLE(FLOAT(ND-I))
00072         INIT = 0
00073         CALL ZUNIK(ZR, ZI, FN, 1, 0, TOL, INIT, PHIR, PHII, ZETA1R,
00074      *   ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
00075         IF (KODE.EQ.1) GO TO 40
00076         STR = ZR + ZETA2R
00077         STI = ZI + ZETA2I
00078         RAST = FN/XZABS(STR,STI)
00079         STR = STR*RAST*RAST
00080         STI = -STI*RAST*RAST
00081         S1R = -ZETA1R + STR
00082         S1I = -ZETA1I + STI + ZI
00083         GO TO 50
00084    40   CONTINUE
00085         S1R = -ZETA1R + ZETA2R
00086         S1I = -ZETA1I + ZETA2I
00087    50   CONTINUE
00088 C-----------------------------------------------------------------------
00089 C     TEST FOR UNDERFLOW AND OVERFLOW
00090 C-----------------------------------------------------------------------
00091         RS1 = S1R
00092         IF (DABS(RS1).GT.ELIM) GO TO 110
00093         IF (I.EQ.1) IFLAG = 2
00094         IF (DABS(RS1).LT.ALIM) GO TO 60
00095 C-----------------------------------------------------------------------
00096 C     REFINE  TEST AND SCALE
00097 C-----------------------------------------------------------------------
00098         APHI = XZABS(PHIR,PHII)
00099         RS1 = RS1 + DLOG(APHI)
00100         IF (DABS(RS1).GT.ELIM) GO TO 110
00101         IF (I.EQ.1) IFLAG = 1
00102         IF (RS1.LT.0.0D0) GO TO 60
00103         IF (I.EQ.1) IFLAG = 3
00104    60   CONTINUE
00105 C-----------------------------------------------------------------------
00106 C     SCALE S1 IF CABS(S1).LT.ASCLE
00107 C-----------------------------------------------------------------------
00108         S2R = PHIR*SUMR - PHII*SUMI
00109         S2I = PHIR*SUMI + PHII*SUMR
00110         STR = DEXP(S1R)*CSSR(IFLAG)
00111         S1R = STR*DCOS(S1I)
00112         S1I = STR*DSIN(S1I)
00113         STR = S2R*S1R - S2I*S1I
00114         S2I = S2R*S1I + S2I*S1R
00115         S2R = STR
00116         IF (IFLAG.NE.1) GO TO 70
00117         CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
00118         IF (NW.NE.0) GO TO 110
00119    70   CONTINUE
00120         CYR(I) = S2R
00121         CYI(I) = S2I
00122         M = ND - I + 1
00123         YR(M) = S2R*CSRR(IFLAG)
00124         YI(M) = S2I*CSRR(IFLAG)
00125    80 CONTINUE
00126       IF (ND.LE.2) GO TO 100
00127       RAST = 1.0D0/XZABS(ZR,ZI)
00128       STR = ZR*RAST
00129       STI = -ZI*RAST
00130       RZR = (STR+STR)*RAST
00131       RZI = (STI+STI)*RAST
00132       BRY(2) = 1.0D0/BRY(1)
00133       BRY(3) = D1MACH(2)
00134       S1R = CYR(1)
00135       S1I = CYI(1)
00136       S2R = CYR(2)
00137       S2I = CYI(2)
00138       C1R = CSRR(IFLAG)
00139       ASCLE = BRY(IFLAG)
00140       K = ND - 2
00141       FN = DBLE(FLOAT(K))
00142       DO 90 I=3,ND
00143         C2R = S2R
00144         C2I = S2I
00145         S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I)
00146         S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R)
00147         S1R = C2R
00148         S1I = C2I
00149         C2R = S2R*C1R
00150         C2I = S2I*C1R
00151         YR(K) = C2R
00152         YI(K) = C2I
00153         K = K - 1
00154         FN = FN - 1.0D0
00155         IF (IFLAG.GE.3) GO TO 90
00156         STR = DABS(C2R)
00157         STI = DABS(C2I)
00158         C2M = DMAX1(STR,STI)
00159         IF (C2M.LE.ASCLE) GO TO 90
00160         IFLAG = IFLAG + 1
00161         ASCLE = BRY(IFLAG)
00162         S1R = S1R*C1R
00163         S1I = S1I*C1R
00164         S2R = C2R
00165         S2I = C2I
00166         S1R = S1R*CSSR(IFLAG)
00167         S1I = S1I*CSSR(IFLAG)
00168         S2R = S2R*CSSR(IFLAG)
00169         S2I = S2I*CSSR(IFLAG)
00170         C1R = CSRR(IFLAG)
00171    90 CONTINUE
00172   100 CONTINUE
00173       RETURN
00174 C-----------------------------------------------------------------------
00175 C     SET UNDERFLOW AND UPDATE PARAMETERS
00176 C-----------------------------------------------------------------------
00177   110 CONTINUE
00178       IF (RS1.GT.0.0D0) GO TO 120
00179       YR(ND) = ZEROR
00180       YI(ND) = ZEROI
00181       NZ = NZ + 1
00182       ND = ND - 1
00183       IF (ND.EQ.0) GO TO 100
00184       CALL ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, NUF, TOL, ELIM, ALIM)
00185       IF (NUF.LT.0) GO TO 120
00186       ND = ND - NUF
00187       NZ = NZ + NUF
00188       IF (ND.EQ.0) GO TO 100
00189       FN = FNU + DBLE(FLOAT(ND-1))
00190       IF (FN.GE.FNUL) GO TO 30
00191       NLAST = ND
00192       RETURN
00193   120 CONTINUE
00194       NZ = -1
00195       RETURN
00196   130 CONTINUE
00197       IF (RS1.GT.0.0D0) GO TO 120
00198       NZ = N
00199       DO 140 I=1,N
00200         YR(I) = ZEROR
00201         YI(I) = ZEROI
00202   140 CONTINUE
00203       RETURN
00204       END