zacon.f

Go to the documentation of this file.

      SUBROUTINE zacon(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, FNUL,
     * tol, elim, alim)
C***BEGIN PROLOGUE  ZACON
C***REFER TO  ZBESK,ZBESH
C
C     ZACON 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  ZBINU,ZBKNU,ZS1S2,D1MACH,XZABS,ZMLT
C***END PROLOGUE  ZACON
C     COMPLEX CK,CONE,CSCL,CSCR,CSGN,CSPN,CY,CZERO,C1,C2,RZ,SC1,SC2,ST,
C    *S1,S2,Y,Z,ZN
      DOUBLE PRECISION ALIM, ARG, ASCLE, AS2, AZN, BRY, BSCLE, CKI,
     * ckr, coner, cpn, cscl, cscr, csgni, csgnr, cspni, cspnr,
     * csr, csrr, cssr, cyi, cyr, c1i, c1m, c1r, c2i, c2r, elim, fmr,
     * fn, fnu, fnul, pi, pti, ptr, razn, rl, rzi, rzr, sc1i, sc1r,
     * sc2i, sc2r, sgn, spn, sti, str, s1i, s1r, s2i, s2r, tol, yi, yr,
     * yy, zeror, zi, zni, znr, zr, d1mach, xzabs
      INTEGER I, INU, IUF, KFLAG, KODE, MR, N, NN, NW, NZ
      dimension yr(n), yi(n), cyr(2), cyi(2), cssr(3), csrr(3), bry(3)
      DATA pi / 3.14159265358979324d0 /
      DATA zeror,coner / 0.0d0,1.0d0 /
      nz = 0
      znr = -zr
      zni = -zi
      nn = n
      CALL zbinu(znr, zni, fnu, kode, nn, yr, yi, nw, rl, fnul, tol,
     * elim, alim)
      IF (nw.LT.0) go to 90
C-----------------------------------------------------------------------
C     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
C-----------------------------------------------------------------------
      nn = min0(2,n)
      CALL zbknu(znr, zni, fnu, kode, nn, cyr, cyi, nw, tol, elim, alim)
      IF (nw.NE.0) go to 90
      s1r = cyr(1)
      s1i = cyi(1)
      fmr = dble(float(mr))
      sgn = -dsign(pi,fmr)
      csgnr = zeror
      csgni = sgn
      IF (kode.EQ.1) go to 10
      yy = -zni
      cpn = dcos(yy)
      spn = dsin(yy)
      CALL zmlt(csgnr, csgni, cpn, spn, csgnr, csgni)
   10 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
      cpn = dcos(arg)
      spn = dsin(arg)
      cspnr = cpn
      cspni = spn
      IF (mod(inu,2).EQ.0) go to 20
      cspnr = -cspnr
      cspni = -cspni
   20 CONTINUE
      iuf = 0
      c1r = s1r
      c1i = s1i
      c2r = yr(1)
      c2i = yi(1)
      ascle = 1.0d+3*d1mach(1)/tol
      IF (kode.EQ.1) go to 30
      CALL zs1s2(znr, zni, c1r, c1i, c2r, c2i, nw, ascle, alim, iuf)
      nz = nz + nw
      sc1r = c1r
      sc1i = c1i
   30 CONTINUE
      CALL zmlt(cspnr, cspni, c1r, c1i, str, sti)
      CALL zmlt(csgnr, csgni, c2r, c2i, ptr, pti)
      yr(1) = str + ptr
      yi(1) = sti + pti
      IF (n.EQ.1) RETURN
      cspnr = -cspnr
      cspni = -cspni
      s2r = cyr(2)
      s2i = cyi(2)
      c1r = s2r
      c1i = s2i
      c2r = yr(2)
      c2i = yi(2)
      IF (kode.EQ.1) go to 40
      CALL zs1s2(znr, zni, c1r, c1i, c2r, c2i, nw, ascle, alim, iuf)
      nz = nz + nw
      sc2r = c1r
      sc2i = c1i
   40 CONTINUE
      CALL zmlt(cspnr, cspni, c1r, c1i, str, sti)
      CALL zmlt(csgnr, csgni, c2r, c2i, ptr, pti)
      yr(2) = str + ptr
      yi(2) = sti + pti
      IF (n.EQ.2) RETURN
      cspnr = -cspnr
      cspni = -cspni
      azn = xzabs(znr,zni)
      razn = 1.0d0/azn
      str = znr*razn
      sti = -zni*razn
      rzr = (str+str)*razn
      rzi = (sti+sti)*razn
      fn = fnu + 1.0d0
      ckr = fn*rzr
      cki = fn*rzi
C-----------------------------------------------------------------------
C     SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS
C-----------------------------------------------------------------------
      cscl = 1.0d0/tol
      cscr = tol
      cssr(1) = cscl
      cssr(2) = coner
      cssr(3) = cscr
      csrr(1) = cscr
      csrr(2) = coner
      csrr(3) = cscl
      bry(1) = ascle
      bry(2) = 1.0d0/ascle
      bry(3) = d1mach(2)
      as2 = xzabs(s2r,s2i)
      kflag = 2
      IF (as2.GT.bry(1)) go to 50
      kflag = 1
      go to 60
   50 CONTINUE
      IF (as2.LT.bry(2)) go to 60
      kflag = 3
   60 CONTINUE
      bscle = bry(kflag)
      s1r = s1r*cssr(kflag)
      s1i = s1i*cssr(kflag)
      s2r = s2r*cssr(kflag)
      s2i = s2i*cssr(kflag)
      csr = csrr(kflag)
      DO 80 i=3,n
        str = s2r
        sti = s2i
        s2r = ckr*str - cki*sti + s1r
        s2i = ckr*sti + cki*str + s1i
        s1r = str
        s1i = sti
        c1r = s2r*csr
        c1i = s2i*csr
        str = c1r
        sti = c1i
        c2r = yr(i)
        c2i = yi(i)
        IF (kode.EQ.1) go to 70
        IF (iuf.LT.0) go to 70
        CALL zs1s2(znr, zni, c1r, c1i, c2r, c2i, nw, ascle, alim, iuf)
        nz = nz + nw
        sc1r = sc2r
        sc1i = sc2i
        sc2r = c1r
        sc2i = c1i
        IF (iuf.NE.3) go to 70
        iuf = -4
        s1r = sc1r*cssr(kflag)
        s1i = sc1i*cssr(kflag)
        s2r = sc2r*cssr(kflag)
        s2i = sc2i*cssr(kflag)
        str = sc2r
        sti = sc2i
   70   CONTINUE
        ptr = cspnr*c1r - cspni*c1i
        pti = cspnr*c1i + cspni*c1r
        yr(i) = ptr + csgnr*c2r - csgni*c2i
        yi(i) = pti + csgnr*c2i + csgni*c2r
        ckr = ckr + rzr
        cki = cki + rzi
        cspnr = -cspnr
        cspni = -cspni
        IF (kflag.GE.3) go to 80
        ptr = dabs(c1r)
        pti = dabs(c1i)
        c1m = dmax1(ptr,pti)
        IF (c1m.LE.bscle) go to 80
        kflag = kflag + 1
        bscle = bry(kflag)
        s1r = s1r*csr
        s1i = s1i*csr
        s2r = str
        s2i = sti
        s1r = s1r*cssr(kflag)
        s1i = s1i*cssr(kflag)
        s2r = s2r*cssr(kflag)
        s2i = s2i*cssr(kflag)
        csr = csrr(kflag)
   80 CONTINUE
      RETURN
   90 CONTINUE
      nz = -1
      IF(nw.EQ.(-2)) nz=-2
      RETURN
      END