cmlri.f

Go to the documentation of this file.

      SUBROUTINE cmlri(Z, FNU, KODE, N, Y, NZ, TOL)
C***BEGIN PROLOGUE  CMLRI
C***REFER TO  CBESI,CBESK
C
C     CMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE
C     MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES.
C
C***ROUTINES CALLED  GAMLN,R1MACH
C***END PROLOGUE  CMLRI
      COMPLEX CK, CNORM, CONE, CTWO, CZERO, PT, P1, P2, RZ, SUM, Y, Z
      REAL ACK, AK, AP, AT, AZ, BK, FKAP, FKK, FLAM, FNF, FNU, RHO,
     * RHO2, SCLE, TFNF, TOL, TST, X, GAMLN, R1MACH
      INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N
      dimension y(n)
      DATA czero,cone,ctwo /(0.0e0,0.0e0),(1.0e0,0.0e0),(2.0e0,0.0e0)/
      scle = 1.0e+3*r1mach(1)/tol
      nz=0
      az = cabs(z)
      x = real(z)
      iaz = int(az)
      ifnu = int(fnu)
      inu = ifnu + n - 1
      at = float(iaz) + 1.0e0
      ck = cmplx(at,0.0e0)/z
      rz = ctwo/z
      p1 = czero
      p2 = cone
      ack = (at+1.0e0)/az
      rho = ack + sqrt(ack*ack-1.0e0)
      rho2 = rho*rho
      tst = (rho2+rho2)/((rho2-1.0e0)*(rho-1.0e0))
      tst = tst/tol
C-----------------------------------------------------------------------
C     COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
C-----------------------------------------------------------------------
      ak = at
      DO 10 i=1,80
        pt = p2
        p2 = p1 - ck*p2
        p1 = pt
        ck = ck + rz
        ap = cabs(p2)
        IF (ap.GT.tst*ak*ak) GO TO 20
        ak = ak + 1.0e0
   10 CONTINUE
      GO TO 110
   20 CONTINUE
      i = i + 1
      k = 0
      IF (inu.LT.iaz) GO TO 40
C-----------------------------------------------------------------------
C     COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
C-----------------------------------------------------------------------
      p1 = czero
      p2 = cone
      at = float(inu) + 1.0e0
      ck = cmplx(at,0.0e0)/z
      ack = at/az
      tst = sqrt(ack/tol)
      itime = 1
      DO 30 k=1,80
        pt = p2
        p2 = p1 - ck*p2
        p1 = pt
        ck = ck + rz
        ap = cabs(p2)
        IF (ap.LT.tst) GO TO 30
        IF (itime.EQ.2) GO TO 40
        ack = cabs(ck)
        flam = ack + sqrt(ack*ack-1.0e0)
        fkap = ap/cabs(p1)
        rho = amin1(flam,fkap)
        tst = tst*sqrt(rho/(rho*rho-1.0e0))
        itime = 2
   30 CONTINUE
      GO TO 110
   40 CONTINUE
C-----------------------------------------------------------------------
C     BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
C-----------------------------------------------------------------------
      k = k + 1
      kk = max0(i+iaz,k+inu)
      fkk = float(kk)
      p1 = czero
C-----------------------------------------------------------------------
C     SCALE P2 AND SUM BY SCLE
C-----------------------------------------------------------------------
      p2 = cmplx(scle,0.0e0)
      fnf = fnu - float(ifnu)
      tfnf = fnf + fnf
      bk = gamln(fkk+tfnf+1.0e0,idum) - gamln(fkk+1.0e0,idum)
     *     -gamln(tfnf+1.0e0,idum)
      bk = exp(bk)
      sum = czero
      km = kk - inu
      DO 50 i=1,km
        pt = p2
        p2 = p1 + cmplx(fkk+fnf,0.0e0)*rz*p2
        p1 = pt
        ak = 1.0e0 - tfnf/(fkk+tfnf)
        ack = bk*ak
        sum = sum + cmplx(ack+bk,0.0e0)*p1
        bk = ack
        fkk = fkk - 1.0e0
   50 CONTINUE
      y(n) = p2
      IF (n.EQ.1) GO TO 70
      DO 60 i=2,n
        pt = p2
        p2 = p1 + cmplx(fkk+fnf,0.0e0)*rz*p2
        p1 = pt
        ak = 1.0e0 - tfnf/(fkk+tfnf)
        ack = bk*ak
        sum = sum + cmplx(ack+bk,0.0e0)*p1
        bk = ack
        fkk = fkk - 1.0e0
        m = n - i + 1
        y(m) = p2
   60 CONTINUE
   70 CONTINUE
      IF (ifnu.LE.0) GO TO 90
      DO 80 i=1,ifnu
        pt = p2
        p2 = p1 + cmplx(fkk+fnf,0.0e0)*rz*p2
        p1 = pt
        ak = 1.0e0 - tfnf/(fkk+tfnf)
        ack = bk*ak
        sum = sum + cmplx(ack+bk,0.0e0)*p1
        bk = ack
        fkk = fkk - 1.0e0
   80 CONTINUE
   90 CONTINUE
      pt = z
      IF (kode.EQ.2) pt = pt - cmplx(x,0.0e0)
      p1 = -cmplx(fnf,0.0e0)*clog(rz) + pt
      ap = gamln(1.0e0+fnf,idum)
      pt = p1 - cmplx(ap,0.0e0)
C-----------------------------------------------------------------------
C     THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
C     IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
C-----------------------------------------------------------------------
      p2 = p2 + sum
      ap = cabs(p2)
      p1 = cmplx(1.0e0/ap,0.0e0)
      ck = cexp(pt)*p1
      pt = conjg(p2)*p1
      cnorm = ck*pt
      DO 100 i=1,n
        y(i) = y(i)*cnorm
  100 CONTINUE
      RETURN
  110 CONTINUE
      nz=-2
      RETURN
      END