cbuni.f

Go to the documentation of this file.

      SUBROUTINE cbuni(Z, FNU, KODE, N, Y, NZ, NUI, NLAST, FNUL, TOL,
     * elim, alim)
C***BEGIN PROLOGUE  CBUNI
C***REFER TO  CBESI,CBESK
C
C     CBUNI COMPUTES THE I BESSEL FUNCTION FOR LARGE CABS(Z).GT.
C     FNUL AND FNU+N-1.LT.FNUL. THE ORDER IS INCREASED FROM
C     FNU+N-1 GREATER THAN FNUL BY ADDING NUI AND COMPUTING
C     ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR I(FNU,Z)
C     ON IFORM=1 AND THE EXPANSION FOR J(FNU,Z) ON IFORM=2
C
C***ROUTINES CALLED  CUNI1,CUNI2,R1MACH
C***END PROLOGUE  CBUNI
      COMPLEX cscl, cscr, cy, rz, st, s1, s2, y, z
      REAL alim, ax, ay, dfnu, elim, fnu, fnui, fnul, gnu, tol, xx, yy,
     * ascle, bry, str, sti, stm, r1mach
      INTEGER i, iflag, iform, k, kode, n, nl, nlast, nui, nw, nz
      dimension y(n), cy(2), bry(3)
      nz = 0
      xx = REAL(z)
      yy = aimag(z)
      ax = abs(xx)*1.7321e0
      ay = abs(yy)
      iform = 1
      IF (ay.GT.ax) iform = 2
      IF (nui.EQ.0) go to 60
      fnui = float(nui)
      dfnu = fnu + float(n-1)
      gnu = dfnu + fnui
      IF (iform.EQ.2) go to 10
C-----------------------------------------------------------------------
C     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
C     -PI/3.LE.ARG(Z).LE.PI/3
C-----------------------------------------------------------------------
      CALL cuni1(z, gnu, kode, 2, cy, nw, nlast, fnul, tol, elim, alim)
      go to 20
   10 CONTINUE
C-----------------------------------------------------------------------
C     ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
C     APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
C     AND HPI=PI/2
C-----------------------------------------------------------------------
      CALL cuni2(z, gnu, kode, 2, cy, nw, nlast, fnul, tol, elim, alim)
   20 CONTINUE
      IF (nw.LT.0) go to 50
      IF (nw.NE.0) go to 90
      ay = cabs(cy(1))
C----------------------------------------------------------------------
C     SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED
C----------------------------------------------------------------------
      bry(1) = 1.0e+3*r1mach(1)/tol
      bry(2) = 1.0e0/bry(1)
      bry(3) = bry(2)
      iflag = 2
      ascle = bry(2)
      ax = 1.0e0
      cscl = cmplx(ax,0.0e0)
      IF (ay.GT.bry(1)) go to 21
      iflag = 1
      ascle = bry(1)
      ax = 1.0e0/tol
      cscl = cmplx(ax,0.0e0)
      go to 25
   21 CONTINUE
      IF (ay.LT.bry(2)) go to 25
      iflag = 3
      ascle = bry(3)
      ax = tol
      cscl = cmplx(ax,0.0e0)
   25 CONTINUE
      ay = 1.0e0/ax
      cscr = cmplx(ay,0.0e0)
      s1 = cy(2)*cscl
      s2 = cy(1)*cscl
      rz = cmplx(2.0e0,0.0e0)/z
      DO 30 i=1,nui
        st = s2
        s2 = cmplx(dfnu+fnui,0.0e0)*rz*s2 + s1
        s1 = st
        fnui = fnui - 1.0e0
        IF (iflag.GE.3) go to 30
        st = s2*cscr
        str = REAL(st)
        sti = aimag(st)
        str = abs(str)
        sti = abs(sti)
        stm = amax1(str,sti)
        IF (stm.LE.ascle) go to 30
        iflag = iflag+1
        ascle = bry(iflag)
        s1 = s1*cscr
        s2 = st
        ax = ax*tol
        ay = 1.0e0/ax
        cscl = cmplx(ax,0.0e0)
        cscr = cmplx(ay,0.0e0)
        s1 = s1*cscl
        s2 = s2*cscl
   30 CONTINUE
      y(n) = s2*cscr
      IF (n.EQ.1) RETURN
      nl = n - 1
      fnui = float(nl)
      k = nl
      DO 40 i=1,nl
        st = s2
        s2 = cmplx(fnu+fnui,0.0e0)*rz*s2 + s1
        s1 = st
        st = s2*cscr
        y(k) = st
        fnui = fnui - 1.0e0
        k = k - 1
        IF (iflag.GE.3) go to 40
        str = REAL(st)
        sti = aimag(st)
        str = abs(str)
        sti = abs(sti)
        stm = amax1(str,sti)
        IF (stm.LE.ascle) go to 40
        iflag = iflag+1
        ascle = bry(iflag)
        s1 = s1*cscr
        s2 = st
        ax = ax*tol
        ay = 1.0e0/ax
        cscl = cmplx(ax,0.0e0)
        cscr = cmplx(ay,0.0e0)
        s1 = s1*cscl
        s2 = s2*cscl
   40 CONTINUE
      RETURN
   50 CONTINUE
      nz = -1
      IF(nw.EQ.(-2)) nz=-2
      RETURN
   60 CONTINUE
      IF (iform.EQ.2) go to 70
C-----------------------------------------------------------------------
C     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
C     -PI/3.LE.ARG(Z).LE.PI/3
C-----------------------------------------------------------------------
      CALL cuni1(z, fnu, kode, n, y, nw, nlast, fnul, tol, elim, alim)
      go to 80
   70 CONTINUE
C-----------------------------------------------------------------------
C     ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
C     APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
C     AND HPI=PI/2
C-----------------------------------------------------------------------
      CALL cuni2(z, fnu, kode, n, y, nw, nlast, fnul, tol, elim, alim)
   80 CONTINUE
      IF (nw.LT.0) go to 50
      nz = nw
      RETURN
   90 CONTINUE
      nlast = n
      RETURN
      END