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

Go to the documentation of this file.
00001       SUBROUTINE ZUOIK(ZR, ZI, FNU, KODE, IKFLG, N, YR, YI, NUF, TOL,
00002      * ELIM, ALIM)
00003 C***BEGIN PROLOGUE  ZUOIK
00004 C***REFER TO  ZBESI,ZBESK,ZBESH
00005 C
00006 C     ZUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC
00007 C     EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM
00008 C     (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW
00009 C     WHERE ALIM.LT.ELIM. IF THE MAGNITUDE, BASED ON THE LEADING
00010 C     EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN
00011 C     THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER
00012 C     MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE
00013 C     EXP(-ELIM)=SMALLEST MACHINE NUMBER*1.0E+3 AND EXP(-ALIM)=
00014 C     EXP(-ELIM)/TOL
00015 C
00016 C     IKFLG=1 MEANS THE I SEQUENCE IS TESTED
00017 C          =2 MEANS THE K SEQUENCE IS TESTED
00018 C     NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE
00019 C         =-1 MEANS AN OVERFLOW WOULD OCCUR
00020 C     IKFLG=1 AND NUF.GT.0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO
00021 C             THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE
00022 C     IKFLG=2 AND NUF.EQ.N MEANS ALL Y VALUES WERE SET TO ZERO
00023 C     IKFLG=2 AND 0.LT.NUF.LT.N NOT CONSIDERED. Y MUST BE SET BY
00024 C             ANOTHER ROUTINE
00025 C
00026 C***ROUTINES CALLED  ZUCHK,ZUNHJ,ZUNIK,D1MACH,XZABS,XZLOG
00027 C***END PROLOGUE  ZUOIK
00028 C     COMPLEX ARG,ASUM,BSUM,CWRK,CZ,CZERO,PHI,SUM,Y,Z,ZB,ZETA1,ZETA2,ZN,
00029 C    *ZR
00030       DOUBLE PRECISION AARG, AIC, ALIM, APHI, ARGI, ARGR, ASUMI, ASUMR,
00031      * ASCLE, AX, AY, BSUMI, BSUMR, CWRKI, CWRKR, CZI, CZR, ELIM, FNN,
00032      * FNU, GNN, GNU, PHII, PHIR, RCZ, STR, STI, SUMI, SUMR, TOL, YI,
00033      * YR, ZBI, ZBR, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI,
00034      * ZNI, ZNR, ZR, ZRI, ZRR, D1MACH, XZABS
00035       INTEGER I, IDUM, IFORM, IKFLG, INIT, KODE, N, NN, NUF, NW
00036       DIMENSION YR(N), YI(N), CWRKR(16), CWRKI(16)
00037       DATA ZEROR,ZEROI / 0.0D0, 0.0D0 /
00038       DATA AIC / 1.265512123484645396D+00 /
00039       NUF = 0
00040       NN = N
00041       ZRR = ZR
00042       ZRI = ZI
00043       IF (ZR.GE.0.0D0) GO TO 10
00044       ZRR = -ZR
00045       ZRI = -ZI
00046    10 CONTINUE
00047       ZBR = ZRR
00048       ZBI = ZRI
00049       AX = DABS(ZR)*1.7321D0
00050       AY = DABS(ZI)
00051       IFORM = 1
00052       IF (AY.GT.AX) IFORM = 2
00053       GNU = DMAX1(FNU,1.0D0)
00054       IF (IKFLG.EQ.1) GO TO 20
00055       FNN = DBLE(FLOAT(NN))
00056       GNN = FNU + FNN - 1.0D0
00057       GNU = DMAX1(GNN,FNN)
00058    20 CONTINUE
00059 C-----------------------------------------------------------------------
00060 C     ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE
00061 C     REAL PARTS OF ZETA1, ZETA2 AND ZB. NO ATTEMPT IS MADE TO GET
00062 C     THE SIGN OF THE IMAGINARY PART CORRECT.
00063 C-----------------------------------------------------------------------
00064       IF (IFORM.EQ.2) GO TO 30
00065       INIT = 0
00066       CALL ZUNIK(ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, PHIR, PHII,
00067      * ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
00068       CZR = -ZETA1R + ZETA2R
00069       CZI = -ZETA1I + ZETA2I
00070       GO TO 50
00071    30 CONTINUE
00072       ZNR = ZRI
00073       ZNI = -ZRR
00074       IF (ZI.GT.0.0D0) GO TO 40
00075       ZNR = -ZNR
00076    40 CONTINUE
00077       CALL ZUNHJ(ZNR, ZNI, GNU, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R,
00078      * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
00079       CZR = -ZETA1R + ZETA2R
00080       CZI = -ZETA1I + ZETA2I
00081       AARG = XZABS(ARGR,ARGI)
00082    50 CONTINUE
00083       IF (KODE.EQ.1) GO TO 60
00084       CZR = CZR - ZBR
00085       CZI = CZI - ZBI
00086    60 CONTINUE
00087       IF (IKFLG.EQ.1) GO TO 70
00088       CZR = -CZR
00089       CZI = -CZI
00090    70 CONTINUE
00091       APHI = XZABS(PHIR,PHII)
00092       RCZ = CZR
00093 C-----------------------------------------------------------------------
00094 C     OVERFLOW TEST
00095 C-----------------------------------------------------------------------
00096       IF (RCZ.GT.ELIM) GO TO 210
00097       IF (RCZ.LT.ALIM) GO TO 80
00098       RCZ = RCZ + DLOG(APHI)
00099       IF (IFORM.EQ.2) RCZ = RCZ - 0.25D0*DLOG(AARG) - AIC
00100       IF (RCZ.GT.ELIM) GO TO 210
00101       GO TO 130
00102    80 CONTINUE
00103 C-----------------------------------------------------------------------
00104 C     UNDERFLOW TEST
00105 C-----------------------------------------------------------------------
00106       IF (RCZ.LT.(-ELIM)) GO TO 90
00107       IF (RCZ.GT.(-ALIM)) GO TO 130
00108       RCZ = RCZ + DLOG(APHI)
00109       IF (IFORM.EQ.2) RCZ = RCZ - 0.25D0*DLOG(AARG) - AIC
00110       IF (RCZ.GT.(-ELIM)) GO TO 110
00111    90 CONTINUE
00112       DO 100 I=1,NN
00113         YR(I) = ZEROR
00114         YI(I) = ZEROI
00115   100 CONTINUE
00116       NUF = NN
00117       RETURN
00118   110 CONTINUE
00119       ASCLE = 1.0D+3*D1MACH(1)/TOL
00120       CALL XZLOG(PHIR, PHII, STR, STI, IDUM)
00121       CZR = CZR + STR
00122       CZI = CZI + STI
00123       IF (IFORM.EQ.1) GO TO 120
00124       CALL XZLOG(ARGR, ARGI, STR, STI, IDUM)
00125       CZR = CZR - 0.25D0*STR - AIC
00126       CZI = CZI - 0.25D0*STI
00127   120 CONTINUE
00128       AX = DEXP(RCZ)/TOL
00129       AY = CZI
00130       CZR = AX*DCOS(AY)
00131       CZI = AX*DSIN(AY)
00132       CALL ZUCHK(CZR, CZI, NW, ASCLE, TOL)
00133       IF (NW.NE.0) GO TO 90
00134   130 CONTINUE
00135       IF (IKFLG.EQ.2) RETURN
00136       IF (N.EQ.1) RETURN
00137 C-----------------------------------------------------------------------
00138 C     SET UNDERFLOWS ON I SEQUENCE
00139 C-----------------------------------------------------------------------
00140   140 CONTINUE
00141       GNU = FNU + DBLE(FLOAT(NN-1))
00142       IF (IFORM.EQ.2) GO TO 150
00143       INIT = 0
00144       CALL ZUNIK(ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, PHIR, PHII,
00145      * ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
00146       CZR = -ZETA1R + ZETA2R
00147       CZI = -ZETA1I + ZETA2I
00148       GO TO 160
00149   150 CONTINUE
00150       CALL ZUNHJ(ZNR, ZNI, GNU, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R,
00151      * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
00152       CZR = -ZETA1R + ZETA2R
00153       CZI = -ZETA1I + ZETA2I
00154       AARG = XZABS(ARGR,ARGI)
00155   160 CONTINUE
00156       IF (KODE.EQ.1) GO TO 170
00157       CZR = CZR - ZBR
00158       CZI = CZI - ZBI
00159   170 CONTINUE
00160       APHI = XZABS(PHIR,PHII)
00161       RCZ = CZR
00162       IF (RCZ.LT.(-ELIM)) GO TO 180
00163       IF (RCZ.GT.(-ALIM)) RETURN
00164       RCZ = RCZ + DLOG(APHI)
00165       IF (IFORM.EQ.2) RCZ = RCZ - 0.25D0*DLOG(AARG) - AIC
00166       IF (RCZ.GT.(-ELIM)) GO TO 190
00167   180 CONTINUE
00168       YR(NN) = ZEROR
00169       YI(NN) = ZEROI
00170       NN = NN - 1
00171       NUF = NUF + 1
00172       IF (NN.EQ.0) RETURN
00173       GO TO 140
00174   190 CONTINUE
00175       ASCLE = 1.0D+3*D1MACH(1)/TOL
00176       CALL XZLOG(PHIR, PHII, STR, STI, IDUM)
00177       CZR = CZR + STR
00178       CZI = CZI + STI
00179       IF (IFORM.EQ.1) GO TO 200
00180       CALL XZLOG(ARGR, ARGI, STR, STI, IDUM)
00181       CZR = CZR - 0.25D0*STR - AIC
00182       CZI = CZI - 0.25D0*STI
00183   200 CONTINUE
00184       AX = DEXP(RCZ)/TOL
00185       AY = CZI
00186       CZR = AX*DCOS(AY)
00187       CZI = AX*DSIN(AY)
00188       CALL ZUCHK(CZR, CZI, NW, ASCLE, TOL)
00189       IF (NW.NE.0) GO TO 180
00190       RETURN
00191   210 CONTINUE
00192       NUF = -1
00193       RETURN
00194       END