zrati.f

Go to the documentation of this file.

      SUBROUTINE zrati(ZR, ZI, FNU, N, CYR, CYI, TOL)
C***BEGIN PROLOGUE  ZRATI
C***REFER TO  ZBESI,ZBESK,ZBESH
C
C     ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD
C     RECURRENCE.  THE STARTING INDEX IS DETERMINED BY FORWARD
C     RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B,
C     MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973,
C     BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER,
C     BY D. J. SOOKNE.
C
C***ROUTINES CALLED  XZABS,ZDIV
C***END PROLOGUE  ZRATI
C     COMPLEX Z,CY(1),CONE,CZERO,P1,P2,T1,RZ,PT,CDFNU
      DOUBLE PRECISION AK, AMAGZ, AP1, AP2, ARG, AZ, CDFNUI, CDFNUR,
     * CONEI, CONER, CYI, CYR, CZEROI, CZEROR, DFNU, FDNU, FLAM, FNU,
     * FNUP, PTI, PTR, P1I, P1R, P2I, P2R, RAK, RAP1, RHO, RT2, RZI,
     * RZR, TEST, TEST1, TOL, TTI, TTR, T1I, T1R, ZI, ZR, XZABS
      INTEGER I, ID, IDNU, INU, ITIME, K, KK, MAGZ, N
      dimension cyr(n), cyi(n)
      DATA czeror,czeroi,coner,conei,rt2/
     1 0.0d0, 0.0d0, 1.0d0, 0.0d0, 1.41421356237309505d0 /
      az = xzabs(zr,zi)
      inu = int(sngl(fnu))
      idnu = inu + n - 1
      magz = int(sngl(az))
      amagz = dble(float(magz+1))
      fdnu = dble(float(idnu))
      fnup = dmax1(amagz,fdnu)
      id = idnu - magz - 1
      itime = 1
      k = 1
      ptr = 1.0d0/az
      rzr = ptr*(zr+zr)*ptr
      rzi = -ptr*(zi+zi)*ptr
      t1r = rzr*fnup
      t1i = rzi*fnup
      p2r = -t1r
      p2i = -t1i
      p1r = coner
      p1i = conei
      t1r = t1r + rzr
      t1i = t1i + rzi
      IF (id.GT.0) id = 0
      ap2 = xzabs(p2r,p2i)
      ap1 = xzabs(p1r,p1i)
C-----------------------------------------------------------------------
C     THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU
C     GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT
C     P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR
C     PREMATURELY.
C-----------------------------------------------------------------------
      arg = (ap2+ap2)/(ap1*tol)
      test1 = dsqrt(arg)
      test = test1
      rap1 = 1.0d0/ap1
      p1r = p1r*rap1
      p1i = p1i*rap1
      p2r = p2r*rap1
      p2i = p2i*rap1
      ap2 = ap2*rap1
   10 CONTINUE
      k = k + 1
      ap1 = ap2
      ptr = p2r
      pti = p2i
      p2r = p1r - (t1r*ptr-t1i*pti)
      p2i = p1i - (t1r*pti+t1i*ptr)
      p1r = ptr
      p1i = pti
      t1r = t1r + rzr
      t1i = t1i + rzi
      ap2 = xzabs(p2r,p2i)
      IF (ap1.LE.test) GO TO 10
      IF (itime.EQ.2) GO TO 20
      ak = xzabs(t1r,t1i)*0.5d0
      flam = ak + dsqrt(ak*ak-1.0d0)
      rho = dmin1(ap2/ap1,flam)
      test = test1*dsqrt(rho/(rho*rho-1.0d0))
      itime = 2
      GO TO 10
   20 CONTINUE
      kk = k + 1 - id
      ak = dble(float(kk))
      t1r = ak
      t1i = czeroi
      dfnu = fnu + dble(float(n-1))
      p1r = 1.0d0/ap2
      p1i = czeroi
      p2r = czeror
      p2i = czeroi
      DO 30 i=1,kk
        ptr = p1r
        pti = p1i
        rap1 = dfnu + t1r
        ttr = rzr*rap1
        tti = rzi*rap1
        p1r = (ptr*ttr-pti*tti) + p2r
        p1i = (ptr*tti+pti*ttr) + p2i
        p2r = ptr
        p2i = pti
        t1r = t1r - coner
   30 CONTINUE
      IF (p1r.NE.czeror .OR. p1i.NE.czeroi) GO TO 40
      p1r = tol
      p1i = tol
   40 CONTINUE
      CALL zdiv(p2r, p2i, p1r, p1i, cyr(n), cyi(n))
      IF (n.EQ.1) RETURN
      k = n - 1
      ak = dble(float(k))
      t1r = ak
      t1i = czeroi
      cdfnur = fnu*rzr
      cdfnui = fnu*rzi
      DO 60 i=2,n
        ptr = cdfnur + (t1r*rzr-t1i*rzi) + cyr(k+1)
        pti = cdfnui + (t1r*rzi+t1i*rzr) + cyi(k+1)
        ak = xzabs(ptr,pti)
        IF (ak.NE.czeror) GO TO 50
        ptr = tol
        pti = tol
        ak = tol*rt2
   50   CONTINUE
        rak = coner/ak
        cyr(k) = rak*ptr*rak
        cyi(k) = -rak*pti*rak
        t1r = t1r - coner
        k = k - 1
   60 CONTINUE
      RETURN
      END