zuoik.f

Go to the documentation of this file.

      SUBROUTINE zuoik(ZR, ZI, FNU, KODE, IKFLG, N, YR, YI, NUF, TOL,
     * elim, alim)
C***BEGIN PROLOGUE  ZUOIK
C***REFER TO  ZBESI,ZBESK,ZBESH
C
C     ZUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC
C     EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM
C     (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW
C     WHERE ALIM.LT.ELIM. IF THE MAGNITUDE, BASED ON THE LEADING
C     EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN
C     THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER
C     MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE
C     EXP(-ELIM)=SMALLEST MACHINE NUMBER*1.0E+3 AND EXP(-ALIM)=
C     EXP(-ELIM)/TOL
C
C     IKFLG=1 MEANS THE I SEQUENCE IS TESTED
C          =2 MEANS THE K SEQUENCE IS TESTED
C     NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE
C         =-1 MEANS AN OVERFLOW WOULD OCCUR
C     IKFLG=1 AND NUF.GT.0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO
C             THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE
C     IKFLG=2 AND NUF.EQ.N MEANS ALL Y VALUES WERE SET TO ZERO
C     IKFLG=2 AND 0.LT.NUF.LT.N NOT CONSIDERED. Y MUST BE SET BY
C             ANOTHER ROUTINE
C
C***ROUTINES CALLED  ZUCHK,ZUNHJ,ZUNIK,D1MACH,XZABS,XZLOG
C***END PROLOGUE  ZUOIK
C     COMPLEX ARG,ASUM,BSUM,CWRK,CZ,CZERO,PHI,SUM,Y,Z,ZB,ZETA1,ZETA2,ZN,
C    *ZR
      DOUBLE PRECISION aarg, aic, alim, aphi, argi, argr, asumi, asumr,
     * ascle, ax, ay, bsumi, bsumr, cwrki, cwrkr, czi, czr, elim, fnn,
     * fnu, gnn, gnu, phii, phir, rcz, str, sti, sumi, sumr, tol, yi,
     * yr, zbi, zbr, zeroi, zeror, zeta1i, zeta1r, zeta2i, zeta2r, zi,
     * zni, znr, zr, zri, zrr, d1mach, xzabs
      INTEGER i, idum, iform, ikflg, init, kode, n, nn, nuf, nw
      dimension yr(n), yi(n), cwrkr(16), cwrki(16)
      DATA zeror,zeroi / 0.0d0, 0.0d0 /
      DATA aic / 1.265512123484645396d+00 /
      nuf = 0
      nn = n
      zrr = zr
      zri = zi
      IF (zr.GE.0.0d0) go to 10
      zrr = -zr
      zri = -zi
   10 CONTINUE
      zbr = zrr
      zbi = zri
      ax = dabs(zr)*1.7321d0
      ay = dabs(zi)
      iform = 1
      IF (ay.GT.ax) iform = 2
      gnu = dmax1(fnu,1.0d0)
      IF (ikflg.EQ.1) go to 20
      fnn = dble(float(nn))
      gnn = fnu + fnn - 1.0d0
      gnu = dmax1(gnn,fnn)
   20 CONTINUE
C-----------------------------------------------------------------------
C     ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE
C     REAL PARTS OF ZETA1, ZETA2 AND ZB. NO ATTEMPT IS MADE TO GET
C     THE SIGN OF THE IMAGINARY PART CORRECT.
C-----------------------------------------------------------------------
      IF (iform.EQ.2) go to 30
      init = 0
      CALL zunik(zrr, zri, gnu, ikflg, 1, tol, init, phir, phii,
     * zeta1r, zeta1i, zeta2r, zeta2i, sumr, sumi, cwrkr, cwrki)
      czr = -zeta1r + zeta2r
      czi = -zeta1i + zeta2i
      go to 50
   30 CONTINUE
      znr = zri
      zni = -zrr
      IF (zi.GT.0.0d0) go to 40
      znr = -znr
   40 CONTINUE
      CALL zunhj(znr, zni, gnu, 1, tol, phir, phii, argr, argi, zeta1r,
     * zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
      czr = -zeta1r + zeta2r
      czi = -zeta1i + zeta2i
      aarg = xzabs(argr,argi)
   50 CONTINUE
      IF (kode.EQ.1) go to 60
      czr = czr - zbr
      czi = czi - zbi
   60 CONTINUE
      IF (ikflg.EQ.1) go to 70
      czr = -czr
      czi = -czi
   70 CONTINUE
      aphi = xzabs(phir,phii)
      rcz = czr
C-----------------------------------------------------------------------
C     OVERFLOW TEST
C-----------------------------------------------------------------------
      IF (rcz.GT.elim) go to 210
      IF (rcz.LT.alim) go to 80
      rcz = rcz + dlog(aphi)
      IF (iform.EQ.2) rcz = rcz - 0.25d0*dlog(aarg) - aic
      IF (rcz.GT.elim) go to 210
      go to 130
   80 CONTINUE
C-----------------------------------------------------------------------
C     UNDERFLOW TEST
C-----------------------------------------------------------------------
      IF (rcz.LT.(-elim)) go to 90
      IF (rcz.GT.(-alim)) go to 130
      rcz = rcz + dlog(aphi)
      IF (iform.EQ.2) rcz = rcz - 0.25d0*dlog(aarg) - aic
      IF (rcz.GT.(-elim)) go to 110
   90 CONTINUE
      DO 100 i=1,nn
        yr(i) = zeror
        yi(i) = zeroi
  100 CONTINUE
      nuf = nn
      RETURN
  110 CONTINUE
      ascle = 1.0d+3*d1mach(1)/tol
      CALL xzlog(phir, phii, str, sti, idum)
      czr = czr + str
      czi = czi + sti
      IF (iform.EQ.1) go to 120
      CALL xzlog(argr, argi, str, sti, idum)
      czr = czr - 0.25d0*str - aic
      czi = czi - 0.25d0*sti
  120 CONTINUE
      ax = dexp(rcz)/tol
      ay = czi
      czr = ax*dcos(ay)
      czi = ax*dsin(ay)
      CALL zuchk(czr, czi, nw, ascle, tol)
      IF (nw.NE.0) go to 90
  130 CONTINUE
      IF (ikflg.EQ.2) RETURN
      IF (n.EQ.1) RETURN
C-----------------------------------------------------------------------
C     SET UNDERFLOWS ON I SEQUENCE
C-----------------------------------------------------------------------
  140 CONTINUE
      gnu = fnu + dble(float(nn-1))
      IF (iform.EQ.2) go to 150
      init = 0
      CALL zunik(zrr, zri, gnu, ikflg, 1, tol, init, phir, phii,
     * zeta1r, zeta1i, zeta2r, zeta2i, sumr, sumi, cwrkr, cwrki)
      czr = -zeta1r + zeta2r
      czi = -zeta1i + zeta2i
      go to 160
  150 CONTINUE
      CALL zunhj(znr, zni, gnu, 1, tol, phir, phii, argr, argi, zeta1r,
     * zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
      czr = -zeta1r + zeta2r
      czi = -zeta1i + zeta2i
      aarg = xzabs(argr,argi)
  160 CONTINUE
      IF (kode.EQ.1) go to 170
      czr = czr - zbr
      czi = czi - zbi
  170 CONTINUE
      aphi = xzabs(phir,phii)
      rcz = czr
      IF (rcz.LT.(-elim)) go to 180
      IF (rcz.GT.(-alim)) RETURN
      rcz = rcz + dlog(aphi)
      IF (iform.EQ.2) rcz = rcz - 0.25d0*dlog(aarg) - aic
      IF (rcz.GT.(-elim)) go to 190
  180 CONTINUE
      yr(nn) = zeror
      yi(nn) = zeroi
      nn = nn - 1
      nuf = nuf + 1
      IF (nn.EQ.0) RETURN
      go to 140
  190 CONTINUE
      ascle = 1.0d+3*d1mach(1)/tol
      CALL xzlog(phir, phii, str, sti, idum)
      czr = czr + str
      czi = czi + sti
      IF (iform.EQ.1) go to 200
      CALL xzlog(argr, argi, str, sti, idum)
      czr = czr - 0.25d0*str - aic
      czi = czi - 0.25d0*sti
  200 CONTINUE
      ax = dexp(rcz)/tol
      ay = czi
      czr = ax*dcos(ay)
      czi = ax*dsin(ay)
      CALL zuchk(czr, czi, nw, ascle, tol)
      IF (nw.NE.0) go to 180
      RETURN
  210 CONTINUE
      nuf = -1
      RETURN
      END