zseri.f

Go to the documentation of this file.

      SUBROUTINE zseri(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM,
     * alim)
C***BEGIN PROLOGUE  ZSERI
C***REFER TO  ZBESI,ZBESK
C
C     ZSERI 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  DGAMLN,D1MACH,ZUCHK,XZABS,ZDIV,XZLOG,ZMLT
C***END PROLOGUE  ZSERI
C     COMPLEX AK1,CK,COEF,CONE,CRSC,CSCL,CZ,CZERO,HZ,RZ,S1,S2,Y,Z
      DOUBLE PRECISION aa, acz, ak, ak1i, ak1r, alim, arm, ascle, atol,
     * az, cki, ckr, coefi, coefr, conei, coner, crscr, czi, czr, dfnu,
     * elim, fnu, fnup, hzi, hzr, raz, rs, rtr1, rzi, rzr, s, ss, sti,
     * str, s1i, s1r, s2i, s2r, tol, yi, yr, wi, wr, zeroi, zeror, zi,
     * zr, dgamln, d1mach, xzabs
      INTEGER i, ib, idum, iflag, il, k, kode, l, m, n, nn, nz, nw
      dimension yr(n), yi(n), wr(2), wi(2)
      DATA zeror,zeroi,coner,conei / 0.0d0, 0.0d0, 1.0d0, 0.0d0 /
C
      nz = 0
      az = xzabs(zr,zi)
      IF (az.EQ.0.0d0) go to 160
      arm = 1.0d+3*d1mach(1)
      rtr1 = dsqrt(arm)
      crscr = 1.0d0
      iflag = 0
      IF (az.LT.arm) go to 150
      hzr = 0.5d0*zr
      hzi = 0.5d0*zi
      czr = zeror
      czi = zeroi
      IF (az.LE.rtr1) go to 10
      CALL zmlt(hzr, hzi, hzr, hzi, czr, czi)
   10 CONTINUE
      acz = xzabs(czr,czi)
      nn = n
      CALL xzlog(hzr, hzi, ckr, cki, idum)
   20 CONTINUE
      dfnu = fnu + dble(float(nn-1))
      fnup = dfnu + 1.0d0
C-----------------------------------------------------------------------
C     UNDERFLOW TEST
C-----------------------------------------------------------------------
      ak1r = ckr*dfnu
      ak1i = cki*dfnu
      ak = dgamln(fnup,idum)
      ak1r = ak1r - ak
      IF (kode.EQ.2) ak1r = ak1r - zr
      IF (ak1r.GT.(-elim)) go to 40
   30 CONTINUE
      nz = nz + 1
      yr(nn) = zeror
      yi(nn) = zeroi
      IF (acz.GT.dfnu) go to 190
      nn = nn - 1
      IF (nn.EQ.0) RETURN
      go to 20
   40 CONTINUE
      IF (ak1r.GT.(-alim)) go to 50
      iflag = 1
      ss = 1.0d0/tol
      crscr = tol
      ascle = arm*ss
   50 CONTINUE
      aa = dexp(ak1r)
      IF (iflag.EQ.1) aa = aa*ss
      coefr = aa*dcos(ak1i)
      coefi = aa*dsin(ak1i)
      atol = tol*acz/fnup
      il = min0(2,nn)
      DO 90 i=1,il
        dfnu = fnu + dble(float(nn-i))
        fnup = dfnu + 1.0d0
        s1r = coner
        s1i = conei
        IF (acz.LT.tol*fnup) go to 70
        ak1r = coner
        ak1i = conei
        ak = fnup + 2.0d0
        s = fnup
        aa = 2.0d0
   60   CONTINUE
        rs = 1.0d0/s
        str = ak1r*czr - ak1i*czi
        sti = ak1r*czi + ak1i*czr
        ak1r = str*rs
        ak1i = sti*rs
        s1r = s1r + ak1r
        s1i = s1i + ak1i
        s = s + ak
        ak = ak + 2.0d0
        aa = aa*acz*rs
        IF (aa.GT.atol) go to 60
   70   CONTINUE
        s2r = s1r*coefr - s1i*coefi
        s2i = s1r*coefi + s1i*coefr
        wr(i) = s2r
        wi(i) = s2i
        IF (iflag.EQ.0) go to 80
        CALL zuchk(s2r, s2i, nw, ascle, tol)
        IF (nw.NE.0) go to 30
   80   CONTINUE
        m = nn - i + 1
        yr(m) = s2r*crscr
        yi(m) = s2i*crscr
        IF (i.EQ.il) go to 90
        CALL zdiv(coefr, coefi, hzr, hzi, str, sti)
        coefr = str*dfnu
        coefi = sti*dfnu
   90 CONTINUE
      IF (nn.LE.2) RETURN
      k = nn - 2
      ak = dble(float(k))
      raz = 1.0d0/az
      str = zr*raz
      sti = -zi*raz
      rzr = (str+str)*raz
      rzi = (sti+sti)*raz
      IF (iflag.EQ.1) go to 120
      ib = 3
  100 CONTINUE
      DO 110 i=ib,nn
        yr(k) = (ak+fnu)*(rzr*yr(k+1)-rzi*yi(k+1)) + yr(k+2)
        yi(k) = (ak+fnu)*(rzr*yi(k+1)+rzi*yr(k+1)) + yi(k+2)
        ak = ak - 1.0d0
        k = k - 1
  110 CONTINUE
      RETURN
C-----------------------------------------------------------------------
C     RECUR BACKWARD WITH SCALED VALUES
C-----------------------------------------------------------------------
  120 CONTINUE
C-----------------------------------------------------------------------
C     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE
C     UNDERFLOW LIMIT = ASCLE = D1MACH(1)*SS*1.0D+3
C-----------------------------------------------------------------------
      s1r = wr(1)
      s1i = wi(1)
      s2r = wr(2)
      s2i = wi(2)
      DO 130 l=3,nn
        ckr = s2r
        cki = s2i
        s2r = s1r + (ak+fnu)*(rzr*ckr-rzi*cki)
        s2i = s1i + (ak+fnu)*(rzr*cki+rzi*ckr)
        s1r = ckr
        s1i = cki
        ckr = s2r*crscr
        cki = s2i*crscr
        yr(k) = ckr
        yi(k) = cki
        ak = ak - 1.0d0
        k = k - 1
        IF (xzabs(ckr,cki).GT.ascle) go to 140
  130 CONTINUE
      RETURN
  140 CONTINUE
      ib = l + 1
      IF (ib.GT.nn) RETURN
      go to 100
  150 CONTINUE
      nz = n
      IF (fnu.EQ.0.0d0) nz = nz - 1
  160 CONTINUE
      yr(1) = zeror
      yi(1) = zeroi
      IF (fnu.NE.0.0d0) go to 170
      yr(1) = coner
      yi(1) = conei
  170 CONTINUE
      IF (n.EQ.1) RETURN
      DO 180 i=2,n
        yr(i) = zeror
        yi(i) = zeroi
  180 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-----------------------------------------------------------------------
  190 CONTINUE
      nz = -nz
      RETURN
      END