casyi.f

Go to the documentation of this file.

      SUBROUTINE casyi(Z, FNU, KODE, N, Y, NZ, RL, TOL, ELIM, ALIM)
C***BEGIN PROLOGUE  CASYI
C***REFER TO  CBESI,CBESK
C
C     CASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
C     MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE
C     REGION CABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN.
C     NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1.
C
C***ROUTINES CALLED  R1MACH
C***END PROLOGUE  CASYI
      COMPLEX AK1, CK, CONE, CS1, CS2, CZ, CZERO, DK, EZ, P1, RZ, S2,
     * Y, Z
      REAL AA, ACZ, AEZ, AK, ALIM, ARG, ARM, ATOL, AZ, BB, BK, DFNU,
     * DNU2, ELIM, FDN, FNU, PI, RL, RTPI, RTR1, S, SGN, SQK, TOL, X,
     * YY, R1MACH
      INTEGER I, IB, IL, INU, J, JL, K, KODE, KODED, M, N, NN, NZ
      dimension y(n)
      DATA pi, rtpi  /3.14159265358979324e0 , 0.159154943091895336e0 /
      DATA czero, cone / (0.0e0,0.0e0), (1.0e0,0.0e0) /
C
      nz = 0
      az = cabs(z)
      x = REAL(Z)
      arm = 1.0e+3*r1mach(1)
      rtr1 = sqrt(arm)
      il = min0(2,n)
      dfnu = fnu + float(n-il)
C-----------------------------------------------------------------------
C     OVERFLOW TEST
C-----------------------------------------------------------------------
      ak1 = cmplx(rtpi,0.0e0)/z
      ak1 = csqrt(ak1)
      cz = z
      IF (kode.EQ.2) cz = z - cmplx(x,0.0e0)
      acz = REAL(CZ)
      IF (abs(acz).GT.elim) GO TO 80
      dnu2 = dfnu + dfnu
      koded = 1
      IF ((abs(acz).GT.alim) .AND. (n.GT.2)) GO TO 10
      koded = 0
      ak1 = ak1*cexp(cz)
   10 CONTINUE
      fdn = 0.0e0
      IF (dnu2.GT.rtr1) fdn = dnu2*dnu2
      ez = z*cmplx(8.0e0,0.0e0)
C-----------------------------------------------------------------------
C     WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE
C     FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE
C     EXPANSION FOR THE IMAGINARY PART.
C-----------------------------------------------------------------------
      aez = 8.0e0*az
      s = tol/aez
      jl = int(rl+rl) + 2
      yy = aimag(z)
      p1 = czero
      IF (yy.EQ.0.0e0) GO TO 20
C-----------------------------------------------------------------------
C     CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF
C     SIGNIFICANCE WHEN FNU OR N IS LARGE
C-----------------------------------------------------------------------
      inu = int(fnu)
      arg = (fnu-float(inu))*pi
      inu = inu + n - il
      ak = -sin(arg)
      bk = cos(arg)
      IF (yy.LT.0.0e0) bk = -bk
      p1 = cmplx(ak,bk)
      IF (mod(inu,2).EQ.1) p1 = -p1
   20 CONTINUE
      DO 50 k=1,il
        sqk = fdn - 1.0e0
        atol = s*abs(sqk)
        sgn = 1.0e0
        cs1 = cone
        cs2 = cone
        ck = cone
        ak = 0.0e0
        aa = 1.0e0
        bb = aez
        dk = ez
        DO 30 j=1,jl
          ck = ck*cmplx(sqk,0.0e0)/dk
          cs2 = cs2 + ck
          sgn = -sgn
          cs1 = cs1 + ck*cmplx(sgn,0.0e0)
          dk = dk + ez
          aa = aa*abs(sqk)/bb
          bb = bb + aez
          ak = ak + 8.0e0
          sqk = sqk - ak
          IF (aa.LE.atol) GO TO 40
   30   CONTINUE
        GO TO 90
   40   CONTINUE
        s2 = cs1
        IF (x+x.LT.elim) s2 = s2 + p1*cs2*cexp(-z-z)
        fdn = fdn + 8.0e0*dfnu + 4.0e0
        p1 = -p1
        m = n - il + k
        y(m) = s2*ak1
   50 CONTINUE
      IF (n.LE.2) RETURN
      nn = n
      k = nn - 2
      ak = float(k)
      rz = (cone+cone)/z
      ib = 3
      DO 60 i=ib,nn
        y(k) = cmplx(ak+fnu,0.0e0)*rz*y(k+1) + y(k+2)
        ak = ak - 1.0e0
        k = k - 1
   60 CONTINUE
      IF (koded.EQ.0) RETURN
      ck = cexp(cz)
      DO 70 i=1,nn
        y(i) = y(i)*ck
   70 CONTINUE
      RETURN
   80 CONTINUE
      nz = -1
      RETURN
   90 CONTINUE
      nz=-2
      RETURN
      END