zmlri.f

Go to the documentation of this file.

      SUBROUTINE zmlri(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL)
C***BEGIN PROLOGUE  ZMLRI
C***REFER TO  ZBESI,ZBESK
C
C     ZMLRI 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  DGAMLN,D1MACH,XZABS,XZEXP,XZLOG,ZMLT
C***END PROLOGUE  ZMLRI
C     COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z
      DOUBLE PRECISION ack, ak, ap, at, az, bk, cki, ckr, cnormi,
     * cnormr, conei, coner, fkap, fkk, flam, fnf, fnu, pti, ptr, p1i,
     * p1r, p2i, p2r, raz, rho, rho2, rzi, rzr, scle, sti, str, sumi,
     * sumr, tfnf, tol, tst, yi, yr, zeroi, zeror, zi, zr, dgamln,
     * d1mach, xzabs
      INTEGER i, iaz, idum, ifnu, inu, itime, k, kk, km, kode, m, n, nz
      dimension yr(n), yi(n)
      DATA zeror,zeroi,coner,conei / 0.0d0, 0.0d0, 1.0d0, 0.0d0 /
      scle = d1mach(1)/tol
      nz=0
      az = xzabs(zr,zi)
      iaz = int(sngl(az))
      ifnu = int(sngl(fnu))
      inu = ifnu + n - 1
      at = dble(float(iaz)) + 1.0d0
      raz = 1.0d0/az
      str = zr*raz
      sti = -zi*raz
      ckr = str*at*raz
      cki = sti*at*raz
      rzr = (str+str)*raz
      rzi = (sti+sti)*raz
      p1r = zeror
      p1i = zeroi
      p2r = coner
      p2i = conei
      ack = (at+1.0d0)*raz
      rho = ack + dsqrt(ack*ack-1.0d0)
      rho2 = rho*rho
      tst = (rho2+rho2)/((rho2-1.0d0)*(rho-1.0d0))
      tst = tst/tol
C-----------------------------------------------------------------------
C     COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
C-----------------------------------------------------------------------
      ak = at
      DO 10 i=1,80
        ptr = p2r
        pti = p2i
        p2r = p1r - (ckr*ptr-cki*pti)
        p2i = p1i - (cki*ptr+ckr*pti)
        p1r = ptr
        p1i = pti
        ckr = ckr + rzr
        cki = cki + rzi
        ap = xzabs(p2r,p2i)
        IF (ap.GT.tst*ak*ak) go to 20
        ak = ak + 1.0d0
   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-----------------------------------------------------------------------
      p1r = zeror
      p1i = zeroi
      p2r = coner
      p2i = conei
      at = dble(float(inu)) + 1.0d0
      str = zr*raz
      sti = -zi*raz
      ckr = str*at*raz
      cki = sti*at*raz
      ack = at*raz
      tst = dsqrt(ack/tol)
      itime = 1
      DO 30 k=1,80
        ptr = p2r
        pti = p2i
        p2r = p1r - (ckr*ptr-cki*pti)
        p2i = p1i - (ckr*pti+cki*ptr)
        p1r = ptr
        p1i = pti
        ckr = ckr + rzr
        cki = cki + rzi
        ap = xzabs(p2r,p2i)
        IF (ap.LT.tst) go to 30
        IF (itime.EQ.2) go to 40
        ack = xzabs(ckr,cki)
        flam = ack + dsqrt(ack*ack-1.0d0)
        fkap = ap/xzabs(p1r,p1i)
        rho = dmin1(flam,fkap)
        tst = tst*dsqrt(rho/(rho*rho-1.0d0))
        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 = dble(float(kk))
      p1r = zeror
      p1i = zeroi
C-----------------------------------------------------------------------
C     SCALE P2 AND SUM BY SCLE
C-----------------------------------------------------------------------
      p2r = scle
      p2i = zeroi
      fnf = fnu - dble(float(ifnu))
      tfnf = fnf + fnf
      bk = dgamln(fkk+tfnf+1.0d0,idum) - dgamln(fkk+1.0d0,idum) -
     * dgamln(tfnf+1.0d0,idum)
      bk = dexp(bk)
      sumr = zeror
      sumi = zeroi
      km = kk - inu
      DO 50 i=1,km
        ptr = p2r
        pti = p2i
        p2r = p1r + (fkk+fnf)*(rzr*ptr-rzi*pti)
        p2i = p1i + (fkk+fnf)*(rzi*ptr+rzr*pti)
        p1r = ptr
        p1i = pti
        ak = 1.0d0 - tfnf/(fkk+tfnf)
        ack = bk*ak
        sumr = sumr + (ack+bk)*p1r
        sumi = sumi + (ack+bk)*p1i
        bk = ack
        fkk = fkk - 1.0d0
   50 CONTINUE
      yr(n) = p2r
      yi(n) = p2i
      IF (n.EQ.1) go to 70
      DO 60 i=2,n
        ptr = p2r
        pti = p2i
        p2r = p1r + (fkk+fnf)*(rzr*ptr-rzi*pti)
        p2i = p1i + (fkk+fnf)*(rzi*ptr+rzr*pti)
        p1r = ptr
        p1i = pti
        ak = 1.0d0 - tfnf/(fkk+tfnf)
        ack = bk*ak
        sumr = sumr + (ack+bk)*p1r
        sumi = sumi + (ack+bk)*p1i
        bk = ack
        fkk = fkk - 1.0d0
        m = n - i + 1
        yr(m) = p2r
        yi(m) = p2i
   60 CONTINUE
   70 CONTINUE
      IF (ifnu.LE.0) go to 90
      DO 80 i=1,ifnu
        ptr = p2r
        pti = p2i
        p2r = p1r + (fkk+fnf)*(rzr*ptr-rzi*pti)
        p2i = p1i + (fkk+fnf)*(rzr*pti+rzi*ptr)
        p1r = ptr
        p1i = pti
        ak = 1.0d0 - tfnf/(fkk+tfnf)
        ack = bk*ak
        sumr = sumr + (ack+bk)*p1r
        sumi = sumi + (ack+bk)*p1i
        bk = ack
        fkk = fkk - 1.0d0
   80 CONTINUE
   90 CONTINUE
      ptr = zr
      pti = zi
      IF (kode.EQ.2) ptr = zeror
      CALL xzlog(rzr, rzi, str, sti, idum)
      p1r = -fnf*str + ptr
      p1i = -fnf*sti + pti
      ap = dgamln(1.0d0+fnf,idum)
      ptr = p1r - ap
      pti = p1i
C-----------------------------------------------------------------------
C     THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
C     IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
C-----------------------------------------------------------------------
      p2r = p2r + sumr
      p2i = p2i + sumi
      ap = xzabs(p2r,p2i)
      p1r = 1.0d0/ap
      CALL xzexp(ptr, pti, str, sti)
      ckr = str*p1r
      cki = sti*p1r
      ptr = p2r*p1r
      pti = -p2i*p1r
      CALL zmlt(ckr, cki, ptr, pti, cnormr, cnormi)
      DO 100 i=1,n
        str = yr(i)*cnormr - yi(i)*cnormi
        yi(i) = yr(i)*cnormi + yi(i)*cnormr
        yr(i) = str
  100 CONTINUE
      RETURN
  110 CONTINUE
      nz=-2
      RETURN
      END