zacai.f

Go to the documentation of this file.

      SUBROUTINE zacai(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, TOL,
     * elim, alim)
C***BEGIN PROLOGUE  ZACAI
C***REFER TO  ZAIRY
C
C     ZACAI APPLIES THE ANALYTIC CONTINUATION FORMULA
C
C         K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
C                 MP=PI*MR*CMPLX(0.0,1.0)
C
C     TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
C     HALF Z PLANE FOR USE WITH ZAIRY WHERE FNU=1/3 OR 2/3 AND N=1.
C     ZACAI IS THE SAME AS ZACON WITH THE PARTS FOR LARGER ORDERS AND
C     RECURRENCE REMOVED. A RECURSIVE CALL TO ZACON CAN RESULT IF ZACON
C     IS CALLED FROM ZAIRY.
C
C***ROUTINES CALLED  ZASYI,ZBKNU,ZMLRI,ZSERI,ZS1S2,D1MACH,XZABS
C***END PROLOGUE  ZACAI
C     COMPLEX CSGN,CSPN,C1,C2,Y,Z,ZN,CY
      DOUBLE PRECISION ALIM, ARG, ASCLE, AZ, CSGNR, CSGNI, CSPNR,
     * cspni, c1r, c1i, c2r, c2i, cyr, cyi, dfnu, elim, fmr, fnu, pi,
     * rl, sgn, tol, yy, yr, yi, zr, zi, znr, zni, d1mach, xzabs
      INTEGER INU, IUF, KODE, MR, N, NN, NW, NZ
      dimension yr(n), yi(n), cyr(2), cyi(2)
      DATA pi / 3.14159265358979324d0 /
      nz = 0
      znr = -zr
      zni = -zi
      az = xzabs(zr,zi)
      nn = n
      dfnu = fnu + dble(float(n-1))
      IF (az.LE.2.0d0) go to 10
      IF (az*az*0.25d0.GT.dfnu+1.0d0) go to 20
   10 CONTINUE
C-----------------------------------------------------------------------
C     POWER SERIES FOR THE I FUNCTION
C-----------------------------------------------------------------------
      CALL zseri(znr, zni, fnu, kode, nn, yr, yi, nw, tol, elim, alim)
      go to 40
   20 CONTINUE
      IF (az.LT.rl) go to 30
C-----------------------------------------------------------------------
C     ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION
C-----------------------------------------------------------------------
      CALL zasyi(znr, zni, fnu, kode, nn, yr, yi, nw, rl, tol, elim,
     * alim)
      IF (nw.LT.0) go to 80
      go to 40
   30 CONTINUE
C-----------------------------------------------------------------------
C     MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION
C-----------------------------------------------------------------------
      CALL zmlri(znr, zni, fnu, kode, nn, yr, yi, nw, tol)
      IF(nw.LT.0) go to 80
   40 CONTINUE
C-----------------------------------------------------------------------
C     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
C-----------------------------------------------------------------------
      CALL zbknu(znr, zni, fnu, kode, 1, cyr, cyi, nw, tol, elim, alim)
      IF (nw.NE.0) go to 80
      fmr = dble(float(mr))
      sgn = -dsign(pi,fmr)
      csgnr = 0.0d0
      csgni = sgn
      IF (kode.EQ.1) go to 50
      yy = -zni
      csgnr = -csgni*dsin(yy)
      csgni = csgni*dcos(yy)
   50 CONTINUE
C-----------------------------------------------------------------------
C     CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
C     WHEN FNU IS LARGE
C-----------------------------------------------------------------------
      inu = int(sngl(fnu))
      arg = (fnu-dble(float(inu)))*sgn
      cspnr = dcos(arg)
      cspni = dsin(arg)
      IF (mod(inu,2).EQ.0) go to 60
      cspnr = -cspnr
      cspni = -cspni
   60 CONTINUE
      c1r = cyr(1)
      c1i = cyi(1)
      c2r = yr(1)
      c2i = yi(1)
      IF (kode.EQ.1) go to 70
      iuf = 0
      ascle = 1.0d+3*d1mach(1)/tol
      CALL zs1s2(znr, zni, c1r, c1i, c2r, c2i, nw, ascle, alim, iuf)
      nz = nz + nw
   70 CONTINUE
      yr(1) = cspnr*c1r - cspni*c1i + csgnr*c2r - csgni*c2i
      yi(1) = cspnr*c1i + cspni*c1r + csgnr*c2i + csgni*c2r
      RETURN
   80 CONTINUE
      nz = -1
      IF(nw.EQ.(-2)) nz=-2
      RETURN
      END