cuoik.f

Go to the documentation of this file.

      SUBROUTINE cuoik(Z, FNU, KODE, IKFLG, N, Y, NUF, TOL, ELIM, ALIM)
C***BEGIN PROLOGUE  CUOIK
C***REFER TO  CBESI,CBESK,CBESH
C
C     CUOIK 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  CUCHK,CUNHJ,CUNIK,R1MACH
C***END PROLOGUE  CUOIK
      COMPLEX arg, asum, bsum, cwrk, cz, czero, phi, sum, y, z, zb,
     * zeta1, zeta2, zn, zr
      REAL aarg, aic, alim, aphi, ascle, ax, ay, elim, fnn, fnu, gnn,
     * gnu, rcz, tol, x, yy
      INTEGER i, iform, ikflg, init, kode, n, nn, nuf, nw
      dimension y(n), cwrk(16)
      DATA czero / (0.0e0,0.0e0) /
      DATA aic / 1.265512123484645396e+00 /
      nuf = 0
      nn = n
      x = REAL(z)
      zr = z
      IF (x.LT.0.0e0) zr = -z
      zb = zr
      yy = aimag(zr)
      ax = abs(x)*1.7321e0
      ay = abs(yy)
      iform = 1
      IF (ay.GT.ax) iform = 2
      gnu = amax1(fnu,1.0e0)
      IF (ikflg.EQ.1) go to 10
      fnn = float(nn)
      gnn = fnu + fnn - 1.0e0
      gnu = amax1(gnn,fnn)
   10 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 20
      init = 0
      CALL cunik(zr, gnu, ikflg, 1, tol, init, phi, zeta1, zeta2, sum,
     * cwrk)
      cz = -zeta1 + zeta2
      go to 40
   20 CONTINUE
      zn = -zr*cmplx(0.0e0,1.0e0)
      IF (yy.GT.0.0e0) go to 30
      zn = conjg(-zn)
   30 CONTINUE
      CALL cunhj(zn, gnu, 1, tol, phi, arg, zeta1, zeta2, asum, bsum)
      cz = -zeta1 + zeta2
      aarg = cabs(arg)
   40 CONTINUE
      IF (kode.EQ.2) cz = cz - zb
      IF (ikflg.EQ.2) cz = -cz
      aphi = cabs(phi)
      rcz = REAL(cz)
C-----------------------------------------------------------------------
C     OVERFLOW TEST
C-----------------------------------------------------------------------
      IF (rcz.GT.elim) go to 170
      IF (rcz.LT.alim) go to 50
      rcz = rcz + alog(aphi)
      IF (iform.EQ.2) rcz = rcz - 0.25e0*alog(aarg) - aic
      IF (rcz.GT.elim) go to 170
      go to 100
   50 CONTINUE
C-----------------------------------------------------------------------
C     UNDERFLOW TEST
C-----------------------------------------------------------------------
      IF (rcz.LT.(-elim)) go to 60
      IF (rcz.GT.(-alim)) go to 100
      rcz = rcz + alog(aphi)
      IF (iform.EQ.2) rcz = rcz - 0.25e0*alog(aarg) - aic
      IF (rcz.GT.(-elim)) go to 80
   60 CONTINUE
      DO 70 i=1,nn
        y(i) = czero
   70 CONTINUE
      nuf = nn
      RETURN
   80 CONTINUE
      ascle = 1.0e+3*r1mach(1)/tol
      cz = cz + clog(phi)
      IF (iform.EQ.1) go to 90
      cz = cz - cmplx(0.25e0,0.0e0)*clog(arg) - cmplx(aic,0.0e0)
   90 CONTINUE
      ax = exp(rcz)/tol
      ay = aimag(cz)
      cz = cmplx(ax,0.0e0)*cmplx(cos(ay),sin(ay))
      CALL cuchk(cz, nw, ascle, tol)
      IF (nw.EQ.1) go to 60
  100 CONTINUE
      IF (ikflg.EQ.2) RETURN
      IF (n.EQ.1) RETURN
C-----------------------------------------------------------------------
C     SET UNDERFLOWS ON I SEQUENCE
C-----------------------------------------------------------------------
  110 CONTINUE
      gnu = fnu + float(nn-1)
      IF (iform.EQ.2) go to 120
      init = 0
      CALL cunik(zr, gnu, ikflg, 1, tol, init, phi, zeta1, zeta2, sum,
     * cwrk)
      cz = -zeta1 + zeta2
      go to 130
  120 CONTINUE
      CALL cunhj(zn, gnu, 1, tol, phi, arg, zeta1, zeta2, asum, bsum)
      cz = -zeta1 + zeta2
      aarg = cabs(arg)
  130 CONTINUE
      IF (kode.EQ.2) cz = cz - zb
      aphi = cabs(phi)
      rcz = REAL(cz)
      IF (rcz.LT.(-elim)) go to 140
      IF (rcz.GT.(-alim)) RETURN
      rcz = rcz + alog(aphi)
      IF (iform.EQ.2) rcz = rcz - 0.25e0*alog(aarg) - aic
      IF (rcz.GT.(-elim)) go to 150
  140 CONTINUE
      y(nn) = czero
      nn = nn - 1
      nuf = nuf + 1
      IF (nn.EQ.0) RETURN
      go to 110
  150 CONTINUE
      ascle = 1.0e+3*r1mach(1)/tol
      cz = cz + clog(phi)
      IF (iform.EQ.1) go to 160
      cz = cz - cmplx(0.25e0,0.0e0)*clog(arg) - cmplx(aic,0.0e0)
  160 CONTINUE
      ax = exp(rcz)/tol
      ay = aimag(cz)
      cz = cmplx(ax,0.0e0)*cmplx(cos(ay),sin(ay))
      CALL cuchk(cz, nw, ascle, tol)
      IF (nw.EQ.1) go to 140
      RETURN
  170 CONTINUE
      nuf = -1
      RETURN
      END