cacon.f

Go to the documentation of this file.

      SUBROUTINE cacon(Z, FNU, KODE, MR, N, Y, NZ, RL, FNUL, TOL, ELIM,
     * ALIM)
C***BEGIN PROLOGUE  CACON
C***REFER TO  CBESK,CBESH
C
C     CACON 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
C
C***ROUTINES CALLED  CBINU,CBKNU,CS1S2,R1MACH
C***END PROLOGUE  CACON
      COMPLEX CK, CONE, CS, CSCL, CSCR, CSGN, CSPN, CSS, CSR, C1, C2,
     * rz, sc1, sc2, st, s1, s2, y, z, zn, cy
      REAL ALIM, ARG, ASCLE, AS2, BSCLE, BRY, CPN, C1I, C1M, C1R, ELIM,
     * fmr, fnu, fnul, pi, rl, sgn, spn, tol, yy, r1mach
      INTEGER I, INU, IUF, KFLAG, KODE, MR, N, NN, NW, NZ
      dimension y(n), cy(2), css(3), csr(3), bry(3)
      DATA pi / 3.14159265358979324e0 /
      DATA cone / (1.0e0,0.0e0) /
      nz = 0
      zn = -z
      nn = n
      CALL cbinu(zn, fnu, kode, nn, y, nw, rl, fnul, tol, elim, alim)
      IF (nw.LT.0) GO TO 80
C-----------------------------------------------------------------------
C     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
C-----------------------------------------------------------------------
      nn = min0(2,n)
      CALL cbknu(zn, fnu, kode, nn, cy, nw, tol, elim, alim)
      IF (nw.NE.0) GO TO 80
      s1 = cy(1)
      fmr = float(mr)
      sgn = -sign(pi,fmr)
      csgn = cmplx(0.0e0,sgn)
      IF (kode.EQ.1) GO TO 10
      yy = -aimag(zn)
      cpn = cos(yy)
      spn = sin(yy)
      csgn = csgn*cmplx(cpn,spn)
   10 CONTINUE
C-----------------------------------------------------------------------
C     CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
C     WHEN FNU IS LARGE
C-----------------------------------------------------------------------
      inu = int(fnu)
      arg = (fnu-float(inu))*sgn
      cpn = cos(arg)
      spn = sin(arg)
      cspn = cmplx(cpn,spn)
      IF (mod(inu,2).EQ.1) cspn = -cspn
      iuf = 0
      c1 = s1
      c2 = y(1)
      ascle = 1.0e+3*r1mach(1)/tol
      IF (kode.EQ.1) GO TO 20
      CALL cs1s2(zn, c1, c2, nw, ascle, alim, iuf)
      nz = nz + nw
      sc1 = c1
   20 CONTINUE
      y(1) = cspn*c1 + csgn*c2
      IF (n.EQ.1) RETURN
      cspn = -cspn
      s2 = cy(2)
      c1 = s2
      c2 = y(2)
      IF (kode.EQ.1) GO TO 30
      CALL cs1s2(zn, c1, c2, nw, ascle, alim, iuf)
      nz = nz + nw
      sc2 = c1
   30 CONTINUE
      y(2) = cspn*c1 + csgn*c2
      IF (n.EQ.2) RETURN
      cspn = -cspn
      rz = cmplx(2.0e0,0.0e0)/zn
      ck = cmplx(fnu+1.0e0,0.0e0)*rz
C-----------------------------------------------------------------------
C     SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS
C-----------------------------------------------------------------------
      cscl = cmplx(1.0e0/tol,0.0e0)
      cscr = cmplx(tol,0.0e0)
      css(1) = cscl
      css(2) = cone
      css(3) = cscr
      csr(1) = cscr
      csr(2) = cone
      csr(3) = cscl
      bry(1) = ascle
      bry(2) = 1.0e0/ascle
      bry(3) = r1mach(2)
      as2 = cabs(s2)
      kflag = 2
      IF (as2.GT.bry(1)) GO TO 40
      kflag = 1
      GO TO 50
   40 CONTINUE
      IF (as2.LT.bry(2)) GO TO 50
      kflag = 3
   50 CONTINUE
      bscle = bry(kflag)
      s1 = s1*css(kflag)
      s2 = s2*css(kflag)
      cs = csr(kflag)
      DO 70 i=3,n
        st = s2
        s2 = ck*s2 + s1
        s1 = st
        c1 = s2*cs
        st = c1
        c2 = y(i)
        IF (kode.EQ.1) GO TO 60
        IF (iuf.LT.0) GO TO 60
        CALL cs1s2(zn, c1, c2, nw, ascle, alim, iuf)
        nz = nz + nw
        sc1 = sc2
        sc2 = c1
        IF (iuf.NE.3) GO TO 60
        iuf = -4
        s1 = sc1*css(kflag)
        s2 = sc2*css(kflag)
        st = sc2
   60   CONTINUE
        y(i) = cspn*c1 + csgn*c2
        ck = ck + rz
        cspn = -cspn
        IF (kflag.GE.3) GO TO 70
        c1r = REAL(C1)
        c1i = aimag(c1)
        c1r = abs(c1r)
        c1i = abs(c1i)
        c1m = amax1(c1r,c1i)
        IF (c1m.LE.bscle) GO TO 70
        kflag = kflag + 1
        bscle = bry(kflag)
        s1 = s1*cs
        s2 = st
        s1 = s1*css(kflag)
        s2 = s2*css(kflag)
        cs = csr(kflag)
   70 CONTINUE
      RETURN
   80 CONTINUE
      nz = -1
      IF(nw.EQ.(-2)) nz=-2
      RETURN
      END