cseri.f

Go to the documentation of this file.

      SUBROUTINE cseri(Z, FNU, KODE, N, Y, NZ, TOL, ELIM, ALIM)
C***BEGIN PROLOGUE  CSERI
C***REFER TO  CBESI,CBESK
C
C     CSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
C     MEANS OF THE POWER SERIES FOR LARGE CABS(Z) IN THE
C     REGION CABS(Z).LE.2*SQRT(FNU+1). NZ=0 IS A NORMAL RETURN.
C     NZ.GT.0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO
C     DUE TO UNDERFLOW. NZ.LT.0 MEANS UNDERFLOW OCCURRED, BUT THE
C     CONDITION CABS(Z).LE.2*SQRT(FNU+1) WAS VIOLATED AND THE
C     COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ).
C
C***ROUTINES CALLED  CUCHK,GAMLN,R1MACH
C***END PROLOGUE  CSERI
      COMPLEX AK1, CK, COEF, CONE, CRSC, CZ, CZERO, HZ, RZ, S1, S2, W,
     * y, z
      REAL AA, ACZ, AK, ALIM, ARM, ASCLE, ATOL, AZ, DFNU, ELIM, FNU,
     * fnup, rak1, rs, rtr1, s, ss, tol, x, gamln, r1mach
      INTEGER I, IB, IDUM, IFLAG, IL, K, KODE, L, M, N, NN, NW, NZ
      dimension y(n), w(2)
      DATA czero, cone / (0.0e0,0.0e0), (1.0e0,0.0e0) /
C
      nz = 0
      az = cabs(z)
      IF (az.EQ.0.0e0) go to 150
      x = REAL(z)
      arm = 1.0e+3*r1mach(1)
      rtr1 = sqrt(arm)
      crsc = cmplx(1.0e0,0.0e0)
      iflag = 0
      IF (az.LT.arm) go to 140
      hz = z*cmplx(0.5e0,0.0e0)
      cz = czero
      IF (az.GT.rtr1) cz = hz*hz
      acz = cabs(cz)
      nn = n
      ck = clog(hz)
   10 CONTINUE
      dfnu = fnu + float(nn-1)
      fnup = dfnu + 1.0e0
C-----------------------------------------------------------------------
C     UNDERFLOW TEST
C-----------------------------------------------------------------------
      ak1 = ck*cmplx(dfnu,0.0e0)
      ak = gamln(fnup,idum)
      ak1 = ak1 - cmplx(ak,0.0e0)
      IF (kode.EQ.2) ak1 = ak1 - cmplx(x,0.0e0)
      rak1 = REAL(ak1)
      IF (rak1.GT.(-elim)) go to 30
   20 CONTINUE
      nz = nz + 1
      y(nn) = czero
      IF (acz.GT.dfnu) go to 170
      nn = nn - 1
      IF (nn.EQ.0) RETURN
      go to 10
   30 CONTINUE
      IF (rak1.GT.(-alim)) go to 40
      iflag = 1
      ss = 1.0e0/tol
      crsc = cmplx(tol,0.0e0)
      ascle = arm*ss
   40 CONTINUE
      ak = aimag(ak1)
      aa = exp(rak1)
      IF (iflag.EQ.1) aa = aa*ss
      coef = cmplx(aa,0.0e0)*cmplx(cos(ak),sin(ak))
      atol = tol*acz/fnup
      il = min0(2,nn)
      DO 80 i=1,il
        dfnu = fnu + float(nn-i)
        fnup = dfnu + 1.0e0
        s1 = cone
        IF (acz.LT.tol*fnup) go to 60
        ak1 = cone
        ak = fnup + 2.0e0
        s = fnup
        aa = 2.0e0
   50   CONTINUE
        rs = 1.0e0/s
        ak1 = ak1*cz*cmplx(rs,0.0e0)
        s1 = s1 + ak1
        s = s + ak
        ak = ak + 2.0e0
        aa = aa*acz*rs
        IF (aa.GT.atol) go to 50
   60   CONTINUE
        m = nn - i + 1
        s2 = s1*coef
        w(i) = s2
        IF (iflag.EQ.0) go to 70
        CALL cuchk(s2, nw, ascle, tol)
        IF (nw.NE.0) go to 20
   70   CONTINUE
        y(m) = s2*crsc
        IF (i.NE.il) coef = coef*cmplx(dfnu,0.0e0)/hz
   80 CONTINUE
      IF (nn.LE.2) RETURN
      k = nn - 2
      ak = float(k)
      rz = (cone+cone)/z
      IF (iflag.EQ.1) go to 110
      ib = 3
   90 CONTINUE
      DO 100 i=ib,nn
        y(k) = cmplx(ak+fnu,0.0e0)*rz*y(k+1) + y(k+2)
        ak = ak - 1.0e0
        k = k - 1
  100 CONTINUE
      RETURN
C-----------------------------------------------------------------------
C     RECUR BACKWARD WITH SCALED VALUES
C-----------------------------------------------------------------------
  110 CONTINUE
C-----------------------------------------------------------------------
C     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE
C     UNDERFLOW LIMIT = ASCLE = R1MACH(1)*CSCL*1.0E+3
C-----------------------------------------------------------------------
      s1 = w(1)
      s2 = w(2)
      DO 120 l=3,nn
        ck = s2
        s2 = s1 + cmplx(ak+fnu,0.0e0)*rz*s2
        s1 = ck
        ck = s2*crsc
        y(k) = ck
        ak = ak - 1.0e0
        k = k - 1
        IF (cabs(ck).GT.ascle) go to 130
  120 CONTINUE
      RETURN
  130 CONTINUE
      ib = l + 1
      IF (ib.GT.nn) RETURN
      go to 90
  140 CONTINUE
      nz = n
      IF (fnu.EQ.0.0e0) nz = nz - 1
  150 CONTINUE
      y(1) = czero
      IF (fnu.EQ.0.0e0) y(1) = cone
      IF (n.EQ.1) RETURN
      DO 160 i=2,n
        y(i) = czero
  160 CONTINUE
      RETURN
C-----------------------------------------------------------------------
C     RETURN WITH NZ.LT.0 IF CABS(Z*Z/4).GT.FNU+N-NZ-1 COMPLETE
C     THE CALCULATION IN CBINU WITH N=N-IABS(NZ)
C-----------------------------------------------------------------------
  170 CONTINUE
      nz = -nz
      RETURN
      END