zasyi.f

Go to the documentation of this file.

      SUBROUTINE zasyi(ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM,
     * alim)
C***BEGIN PROLOGUE  ZASYI
C***REFER TO  ZBESI,ZBESK
C
C     ZASYI 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  D1MACH,XZABS,ZDIV,XZEXP,ZMLT,XZSQRT
C***END PROLOGUE  ZASYI
C     COMPLEX AK1,CK,CONE,CS1,CS2,CZ,CZERO,DK,EZ,P1,RZ,S2,Y,Z
      DOUBLE PRECISION aa, aez, ak, ak1i, ak1r, alim, arg, arm, atol,
     * az, bb, bk, cki, ckr, conei, coner, cs1i, cs1r, cs2i, cs2r, czi,
     * czr, dfnu, dki, dkr, dnu2, elim, ezi, ezr, fdn, fnu, pi, p1i,
     * p1r, raz, rl, rtpi, rtr1, rzi, rzr, s, sgn, sqk, sti, str, s2i,
     * s2r, tol, tzi, tzr, yi, yr, zeroi, zeror, zi, zr, d1mach, xzabs
      INTEGER i, ib, il, inu, j, jl, k, kode, koded, m, n, nn, nz
      dimension yr(n), yi(n)
      DATA pi, rtpi  /3.14159265358979324d0 , 0.159154943091895336d0 /
      DATA zeror,zeroi,coner,conei / 0.0d0, 0.0d0, 1.0d0, 0.0d0 /
C
      nz = 0
      az = xzabs(zr,zi)
      arm = 1.0d+3*d1mach(1)
      rtr1 = dsqrt(arm)
      il = min0(2,n)
      dfnu = fnu + dble(float(n-il))
C-----------------------------------------------------------------------
C     OVERFLOW TEST
C-----------------------------------------------------------------------
      raz = 1.0d0/az
      str = zr*raz
      sti = -zi*raz
      ak1r = rtpi*str*raz
      ak1i = rtpi*sti*raz
      CALL xzsqrt(ak1r, ak1i, ak1r, ak1i)
      czr = zr
      czi = zi
      IF (kode.NE.2) go to 10
      czr = zeror
      czi = zi
   10 CONTINUE
      IF (dabs(czr).GT.elim) go to 100
      dnu2 = dfnu + dfnu
      koded = 1
      IF ((dabs(czr).GT.alim) .AND. (n.GT.2)) go to 20
      koded = 0
      CALL xzexp(czr, czi, str, sti)
      CALL zmlt(ak1r, ak1i, str, sti, ak1r, ak1i)
   20 CONTINUE
      fdn = 0.0d0
      IF (dnu2.GT.rtr1) fdn = dnu2*dnu2
      ezr = zr*8.0d0
      ezi = zi*8.0d0
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.0d0*az
      s = tol/aez
      jl = int(sngl(rl+rl)) + 2
      p1r = zeror
      p1i = zeroi
      IF (zi.EQ.0.0d0) go to 30
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(sngl(fnu))
      arg = (fnu-dble(float(inu)))*pi
      inu = inu + n - il
      ak = -dsin(arg)
      bk = dcos(arg)
      IF (zi.LT.0.0d0) bk = -bk
      p1r = ak
      p1i = bk
      IF (mod(inu,2).EQ.0) go to 30
      p1r = -p1r
      p1i = -p1i
   30 CONTINUE
      DO 70 k=1,il
        sqk = fdn - 1.0d0
        atol = s*dabs(sqk)
        sgn = 1.0d0
        cs1r = coner
        cs1i = conei
        cs2r = coner
        cs2i = conei
        ckr = coner
        cki = conei
        ak = 0.0d0
        aa = 1.0d0
        bb = aez
        dkr = ezr
        dki = ezi
        DO 40 j=1,jl
          CALL zdiv(ckr, cki, dkr, dki, str, sti)
          ckr = str*sqk
          cki = sti*sqk
          cs2r = cs2r + ckr
          cs2i = cs2i + cki
          sgn = -sgn
          cs1r = cs1r + ckr*sgn
          cs1i = cs1i + cki*sgn
          dkr = dkr + ezr
          dki = dki + ezi
          aa = aa*dabs(sqk)/bb
          bb = bb + aez
          ak = ak + 8.0d0
          sqk = sqk - ak
          IF (aa.LE.atol) go to 50
   40   CONTINUE
        go to 110
   50   CONTINUE
        s2r = cs1r
        s2i = cs1i
        IF (zr+zr.GE.elim) go to 60
        tzr = zr + zr
        tzi = zi + zi
        CALL xzexp(-tzr, -tzi, str, sti)
        CALL zmlt(str, sti, p1r, p1i, str, sti)
        CALL zmlt(str, sti, cs2r, cs2i, str, sti)
        s2r = s2r + str
        s2i = s2i + sti
   60   CONTINUE
        fdn = fdn + 8.0d0*dfnu + 4.0d0
        p1r = -p1r
        p1i = -p1i
        m = n - il + k
        yr(m) = s2r*ak1r - s2i*ak1i
        yi(m) = s2r*ak1i + s2i*ak1r
   70 CONTINUE
      IF (n.LE.2) RETURN
      nn = n
      k = nn - 2
      ak = dble(float(k))
      str = zr*raz
      sti = -zi*raz
      rzr = (str+str)*raz
      rzi = (sti+sti)*raz
      ib = 3
      DO 80 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
   80 CONTINUE
      IF (koded.EQ.0) RETURN
      CALL xzexp(czr, czi, ckr, cki)
      DO 90 i=1,nn
        str = yr(i)*ckr - yi(i)*cki
        yi(i) = yr(i)*cki + yi(i)*ckr
        yr(i) = str
   90 CONTINUE
      RETURN
  100 CONTINUE
      nz = -1
      RETURN
  110 CONTINUE
      nz=-2
      RETURN
      END