crati.f

Go to the documentation of this file.

      SUBROUTINE crati(Z, FNU, N, CY, TOL)
C***BEGIN PROLOGUE  CRATI
C***REFER TO  CBESI,CBESK,CBESH
C
C     CRATI 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  (NONE)
C***END PROLOGUE  CRATI
      COMPLEX CDFNU, CONE, CY, CZERO, PT, P1, P2, RZ, T1, Z
      REAL AK, AMAGZ, AP1, AP2, ARG, AZ, DFNU, FDNU, FLAM, FNU, FNUP,
     * RAP1, RHO, TEST, TEST1, TOL
      INTEGER I, ID, IDNU, INU, ITIME, K, KK, MAGZ, N
      dimension cy(n)
      DATA czero, cone / (0.0e0,0.0e0), (1.0e0,0.0e0) /
      az = cabs(z)
      inu = int(fnu)
      idnu = inu + n - 1
      fdnu = float(idnu)
      magz = int(az)
      amagz = float(magz+1)
      fnup = amax1(amagz,fdnu)
      id = idnu - magz - 1
      itime = 1
      k = 1
      rz = (cone+cone)/z
      t1 = cmplx(fnup,0.0e0)*rz
      p2 = -t1
      p1 = cone
      t1 = t1 + rz
      IF (id.GT.0) id = 0
      ap2 = cabs(p2)
      ap1 = cabs(p1)
C-----------------------------------------------------------------------
C     THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNX
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 = sqrt(arg)
      test = test1
      rap1 = 1.0e0/ap1
      p1 = p1*cmplx(rap1,0.0e0)
      p2 = p2*cmplx(rap1,0.0e0)
      ap2 = ap2*rap1
   10 CONTINUE
      k = k + 1
      ap1 = ap2
      pt = p2
      p2 = p1 - t1*p2
      p1 = pt
      t1 = t1 + rz
      ap2 = cabs(p2)
      IF (ap1.LE.test) GO TO 10
      IF (itime.EQ.2) GO TO 20
      ak = cabs(t1)*0.5e0
      flam = ak + sqrt(ak*ak-1.0e0)
      rho = amin1(ap2/ap1,flam)
      test = test1*sqrt(rho/(rho*rho-1.0e0))
      itime = 2
      GO TO 10
   20 CONTINUE
      kk = k + 1 - id
      ak = float(kk)
      dfnu = fnu + float(n-1)
      cdfnu = cmplx(dfnu,0.0e0)
      t1 = cmplx(ak,0.0e0)
      p1 = cmplx(1.0e0/ap2,0.0e0)
      p2 = czero
      DO 30 i=1,kk
        pt = p1
        p1 = rz*(cdfnu+t1)*p1 + p2
        p2 = pt
        t1 = t1 - cone
   30 CONTINUE
      IF (real(p1).NE.0.0e0 .OR. aimag(p1).NE.0.0e0) GO TO 40
      p1 = cmplx(tol,tol)
   40 CONTINUE
      cy(n) = p2/p1
      IF (n.EQ.1) RETURN
      k = n - 1
      ak = float(k)
      t1 = cmplx(ak,0.0e0)
      cdfnu = cmplx(fnu,0.0e0)*rz
      DO 60 i=2,n
        pt = cdfnu + t1*rz + cy(k+1)
        IF (real(pt).NE.0.0e0 .OR. aimag(pt).NE.0.0e0) GO TO 50
        pt = cmplx(tol,tol)
   50   CONTINUE
        cy(k) = cone/pt
        t1 = t1 - cone
        k = k - 1
   60 CONTINUE
      RETURN
      END