zbuni.f

Go to the documentation of this file.

      SUBROUTINE zbuni(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NUI, NLAST,
     * FNUL, TOL, ELIM, ALIM)
C***BEGIN PROLOGUE  ZBUNI
C***REFER TO  ZBESI,ZBESK
C
C     ZBUNI 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  ZUNI1,ZUNI2,XZABS,D1MACH
C***END PROLOGUE  ZBUNI
C     COMPLEX CSCL,CSCR,CY,RZ,ST,S1,S2,Y,Z
      DOUBLE PRECISION ALIM, AX, AY, CSCLR, CSCRR, CYI, CYR, DFNU,
     * elim, fnu, fnui, fnul, gnu, raz, rzi, rzr, sti, str, s1i, s1r,
     * s2i, s2r, tol, yi, yr, zi, zr, xzabs, ascle, bry, c1r, c1i, c1m,
     * d1mach
      INTEGER I, IFLAG, IFORM, K, KODE, N, NL, NLAST, NUI, NW, NZ
      dimension yr(n), yi(n), cyr(2), cyi(2), bry(3)
      nz = 0
      ax = dabs(zr)*1.7321d0
      ay = dabs(zi)
      iform = 1
      IF (ay.GT.ax) iform = 2
      IF (nui.EQ.0) GO TO 60
      fnui = dble(float(nui))
      dfnu = fnu + dble(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 zuni1(zr, zi, gnu, kode, 2, cyr, cyi, 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 zuni2(zr, zi, gnu, kode, 2, cyr, cyi, nw, nlast, fnul, tol,
     * elim, alim)
   20 CONTINUE
      IF (nw.LT.0) GO TO 50
      IF (nw.NE.0) GO TO 90
      str = xzabs(cyr(1),cyi(1))
C----------------------------------------------------------------------
C     SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED
C----------------------------------------------------------------------
      bry(1)=1.0d+3*d1mach(1)/tol
      bry(2) = 1.0d0/bry(1)
      bry(3) = bry(2)
      iflag = 2
      ascle = bry(2)
      csclr = 1.0d0
      IF (str.GT.bry(1)) GO TO 21
      iflag = 1
      ascle = bry(1)
      csclr = 1.0d0/tol
      GO TO 25
   21 CONTINUE
      IF (str.LT.bry(2)) GO TO 25
      iflag = 3
      ascle=bry(3)
      csclr = tol
   25 CONTINUE
      cscrr = 1.0d0/csclr
      s1r = cyr(2)*csclr
      s1i = cyi(2)*csclr
      s2r = cyr(1)*csclr
      s2i = cyi(1)*csclr
      raz = 1.0d0/xzabs(zr,zi)
      str = zr*raz
      sti = -zi*raz
      rzr = (str+str)*raz
      rzi = (sti+sti)*raz
      DO 30 i=1,nui
        str = s2r
        sti = s2i
        s2r = (dfnu+fnui)*(rzr*str-rzi*sti) + s1r
        s2i = (dfnu+fnui)*(rzr*sti+rzi*str) + s1i
        s1r = str
        s1i = sti
        fnui = fnui - 1.0d0
        IF (iflag.GE.3) GO TO 30
        str = s2r*cscrr
        sti = s2i*cscrr
        c1r = dabs(str)
        c1i = dabs(sti)
        c1m = dmax1(c1r,c1i)
        IF (c1m.LE.ascle) GO TO 30
        iflag = iflag+1
        ascle = bry(iflag)
        s1r = s1r*cscrr
        s1i = s1i*cscrr
        s2r = str
        s2i = sti
        csclr = csclr*tol
        cscrr = 1.0d0/csclr
        s1r = s1r*csclr
        s1i = s1i*csclr
        s2r = s2r*csclr
        s2i = s2i*csclr
   30 CONTINUE
      yr(n) = s2r*cscrr
      yi(n) = s2i*cscrr
      IF (n.EQ.1) RETURN
      nl = n - 1
      fnui = dble(float(nl))
      k = nl
      DO 40 i=1,nl
        str = s2r
        sti = s2i
        s2r = (fnu+fnui)*(rzr*str-rzi*sti) + s1r
        s2i = (fnu+fnui)*(rzr*sti+rzi*str) + s1i
        s1r = str
        s1i = sti
        str = s2r*cscrr
        sti = s2i*cscrr
        yr(k) = str
        yi(k) = sti
        fnui = fnui - 1.0d0
        k = k - 1
        IF (iflag.GE.3) GO TO 40
        c1r = dabs(str)
        c1i = dabs(sti)
        c1m = dmax1(c1r,c1i)
        IF (c1m.LE.ascle) GO TO 40
        iflag = iflag+1
        ascle = bry(iflag)
        s1r = s1r*cscrr
        s1i = s1i*cscrr
        s2r = str
        s2i = sti
        csclr = csclr*tol
        cscrr = 1.0d0/csclr
        s1r = s1r*csclr
        s1i = s1i*csclr
        s2r = s2r*csclr
        s2i = s2i*csclr
   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 zuni1(zr, zi, fnu, kode, n, yr, yi, 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 zuni2(zr, zi, fnu, kode, n, yr, yi, nw, nlast, fnul, tol,
     * elim, alim)
   80 CONTINUE
      IF (nw.LT.0) GO TO 50
      nz = nw
      RETURN
   90 CONTINUE
      nlast = n
      RETURN
      END