zunk2.f

Go to the documentation of this file.

      SUBROUTINE zunk2(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM,
     * ALIM)
C***BEGIN PROLOGUE  ZUNK2
C***REFER TO  ZBESK
C
C     ZUNK2 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
C     RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
C     UNIFORM ASYMPTOTIC EXPANSIONS FOR H(KIND,FNU,ZN) AND J(FNU,ZN)
C     WHERE ZN IS IN THE RIGHT HALF PLANE, KIND=(3-MR)/2, MR=+1 OR
C     -1. HERE ZN=ZR*I OR -ZR*I WHERE ZR=Z IF Z IS IN THE RIGHT
C     HALF PLANE OR ZR=-Z IF Z IS IN THE LEFT HALF PLANE. MR INDIC-
C     ATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
C     NZ=-1 MEANS AN OVERFLOW WILL OCCUR
C
C***ROUTINES CALLED  ZAIRY,ZKSCL,ZS1S2,ZUCHK,ZUNHJ,D1MACH,XZABS
C***END PROLOGUE  ZUNK2
C     COMPLEX AI,ARG,ARGD,ASUM,ASUMD,BSUM,BSUMD,CFN,CI,CIP,CK,CONE,CRSC,
C    *CR1,CR2,CS,CSCL,CSGN,CSPN,CSR,CSS,CY,CZERO,C1,C2,DAI,PHI,PHID,RZ,
C    *S1,S2,Y,Z,ZB,ZETA1,ZETA1D,ZETA2,ZETA2D,ZN,ZR
      DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGDI,
     * argdr, argi, argr, asc, ascle, asumdi, asumdr, asumi, asumr,
     * bry, bsumdi, bsumdr, bsumi, bsumr, car, cipi, cipr, cki, ckr,
     * coner, crsc, cr1i, cr1r, cr2i, cr2r, cscl, csgni, csi,
     * cspni, cspnr, csr, csrr, cssr, cyi, cyr, c1i, c1r, c2i, c2m,
     * c2r, daii, dair, elim, fmr, fn, fnf, fnu, hpi, phidi, phidr,
     * phii, phir, pi, pti, ptr, rast, razr, rs1, rzi, rzr, sar, sgn,
     * sti, str, s1i, s1r, s2i, s2r, tol, yi, yr, yy, zbi, zbr, zeroi,
     * zeror, zeta1i, zeta1r, zeta2i, zeta2r, zet1di, zet1dr, zet2di,
     * zet2dr, zi, zni, znr, zr, zri, zrr, d1mach, xzabs
      INTEGER I, IB, IFLAG, IFN, IL, IN, INU, IUF, K, KDFLG, KFLAG, KK,
     * kode, mr, n, nai, ndai, nw, nz, idum, j, ipard, ic
      dimension bry(3), yr(n), yi(n), asumr(2), asumi(2), bsumr(2),
     * bsumi(2), phir(2), phii(2), argr(2), argi(2), zeta1r(2),
     * zeta1i(2), zeta2r(2), zeta2i(2), cyr(2), cyi(2), cipr(4),
     * cipi(4), cssr(3), csrr(3)
      DATA zeror,zeroi,coner,cr1r,cr1i,cr2r,cr2i /
     1         0.0d0, 0.0d0, 1.0d0,
     1 1.0d0,1.73205080756887729d0 , -0.5d0,-8.66025403784438647d-01 /
      DATA hpi, pi, aic /
     1     1.57079632679489662d+00,     3.14159265358979324d+00,
     1     1.26551212348464539d+00/
      DATA cipr(1),cipi(1),cipr(2),cipi(2),cipr(3),cipi(3),cipr(4),
     * cipi(4) /
     1  1.0d0,0.0d0 ,  0.0d0,-1.0d0 ,  -1.0d0,0.0d0 ,  0.0d0,1.0d0 /
C
      kdflg = 1
      nz = 0
C-----------------------------------------------------------------------
C     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
C     THE UNDERFLOW LIMIT
C-----------------------------------------------------------------------
      cscl = 1.0d0/tol
      crsc = tol
      cssr(1) = cscl
      cssr(2) = coner
      cssr(3) = crsc
      csrr(1) = crsc
      csrr(2) = coner
      csrr(3) = cscl
      bry(1) = 1.0d+3*d1mach(1)/tol
      bry(2) = 1.0d0/bry(1)
      bry(3) = d1mach(2)
      zrr = zr
      zri = zi
      IF (zr.GE.0.0d0) GO TO 10
      zrr = -zr
      zri = -zi
   10 CONTINUE
      yy = zri
      znr = zri
      zni = -zrr
      zbr = zrr
      zbi = zri
      inu = int(sngl(fnu))
      fnf = fnu - dble(float(inu))
      ang = -hpi*fnf
      car = dcos(ang)
      sar = dsin(ang)
      c2r = hpi*sar
      c2i = -hpi*car
      kk = mod(inu,4) + 1
      str = c2r*cipr(kk) - c2i*cipi(kk)
      sti = c2r*cipi(kk) + c2i*cipr(kk)
      csr = cr1r*str - cr1i*sti
      csi = cr1r*sti + cr1i*str
      IF (yy.GT.0.0d0) GO TO 20
      znr = -znr
      zbi = -zbi
   20 CONTINUE
C-----------------------------------------------------------------------
C     K(FNU,Z) IS COMPUTED FROM H(2,FNU,-I*Z) WHERE Z IS IN THE FIRST
C     QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
C     CONJUGATION SINCE THE K FUNCTION IS REAL ON THE POSITIVE REAL AXIS
C-----------------------------------------------------------------------
      j = 2
      DO 80 i=1,n
C-----------------------------------------------------------------------
C     J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
C-----------------------------------------------------------------------
        j = 3 - j
        fn = fnu + dble(float(i-1))
        CALL zunhj(znr, zni, fn, 0, tol, phir(j), phii(j), argr(j),
     *   argi(j), zeta1r(j), zeta1i(j), zeta2r(j), zeta2i(j), asumr(j),
     *   asumi(j), bsumr(j), bsumi(j))
        IF (kode.EQ.1) GO TO 30
        str = zbr + zeta2r(j)
        sti = zbi + zeta2i(j)
        rast = fn/xzabs(str,sti)
        str = str*rast*rast
        sti = -sti*rast*rast
        s1r = zeta1r(j) - str
        s1i = zeta1i(j) - sti
        GO TO 40
   30   CONTINUE
        s1r = zeta1r(j) - zeta2r(j)
        s1i = zeta1i(j) - zeta2i(j)
   40   CONTINUE
C-----------------------------------------------------------------------
C     TEST FOR UNDERFLOW AND OVERFLOW
C-----------------------------------------------------------------------
        rs1 = s1r
        IF (dabs(rs1).GT.elim) GO TO 70
        IF (kdflg.EQ.1) kflag = 2
        IF (dabs(rs1).LT.alim) GO TO 50
C-----------------------------------------------------------------------
C     REFINE  TEST AND SCALE
C-----------------------------------------------------------------------
        aphi = xzabs(phir(j),phii(j))
        aarg = xzabs(argr(j),argi(j))
        rs1 = rs1 + dlog(aphi) - 0.25d0*dlog(aarg) - aic
        IF (dabs(rs1).GT.elim) GO TO 70
        IF (kdflg.EQ.1) kflag = 1
        IF (rs1.LT.0.0d0) GO TO 50
        IF (kdflg.EQ.1) kflag = 3
   50   CONTINUE
C-----------------------------------------------------------------------
C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
C     EXPONENT EXTREMES
C-----------------------------------------------------------------------
        c2r = argr(j)*cr2r - argi(j)*cr2i
        c2i = argr(j)*cr2i + argi(j)*cr2r
        CALL zairy(c2r, c2i, 0, 2, air, aii, nai, idum)
        CALL zairy(c2r, c2i, 1, 2, dair, daii, ndai, idum)
        str = dair*bsumr(j) - daii*bsumi(j)
        sti = dair*bsumi(j) + daii*bsumr(j)
        ptr = str*cr2r - sti*cr2i
        pti = str*cr2i + sti*cr2r
        str = ptr + (air*asumr(j)-aii*asumi(j))
        sti = pti + (air*asumi(j)+aii*asumr(j))
        ptr = str*phir(j) - sti*phii(j)
        pti = str*phii(j) + sti*phir(j)
        s2r = ptr*csr - pti*csi
        s2i = ptr*csi + pti*csr
        str = dexp(s1r)*cssr(kflag)
        s1r = str*dcos(s1i)
        s1i = str*dsin(s1i)
        str = s2r*s1r - s2i*s1i
        s2i = s1r*s2i + s2r*s1i
        s2r = str
        IF (kflag.NE.1) GO TO 60
        CALL zuchk(s2r, s2i, nw, bry(1), tol)
        IF (nw.NE.0) GO TO 70
   60   CONTINUE
        IF (yy.LE.0.0d0) s2i = -s2i
        cyr(kdflg) = s2r
        cyi(kdflg) = s2i
        yr(i) = s2r*csrr(kflag)
        yi(i) = s2i*csrr(kflag)
        str = csi
        csi = -csr
        csr = str
        IF (kdflg.EQ.2) GO TO 85
        kdflg = 2
        GO TO 80
   70   CONTINUE
        IF (rs1.GT.0.0d0) GO TO 320
C-----------------------------------------------------------------------
C     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
C-----------------------------------------------------------------------
        IF (zr.LT.0.0d0) GO TO 320
        kdflg = 1
        yr(i)=zeror
        yi(i)=zeroi
        nz=nz+1
        str = csi
        csi =-csr
        csr = str
        IF (i.EQ.1) GO TO 80
        IF ((yr(i-1).EQ.zeror).AND.(yi(i-1).EQ.zeroi)) GO TO 80
        yr(i-1)=zeror
        yi(i-1)=zeroi
        nz=nz+1
   80 CONTINUE
      i = n
   85 CONTINUE
      razr = 1.0d0/xzabs(zrr,zri)
      str = zrr*razr
      sti = -zri*razr
      rzr = (str+str)*razr
      rzi = (sti+sti)*razr
      ckr = fn*rzr
      cki = fn*rzi
      ib = i + 1
      IF (n.LT.ib) GO TO 180
C-----------------------------------------------------------------------
C     TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
C     ON UNDERFLOW.
C-----------------------------------------------------------------------
      fn = fnu + dble(float(n-1))
      ipard = 1
      IF (mr.NE.0) ipard = 0
      CALL zunhj(znr, zni, fn, ipard, tol, phidr, phidi, argdr, argdi,
     * zet1dr, zet1di, zet2dr, zet2di, asumdr, asumdi, bsumdr, bsumdi)
      IF (kode.EQ.1) GO TO 90
      str = zbr + zet2dr
      sti = zbi + zet2di
      rast = fn/xzabs(str,sti)
      str = str*rast*rast
      sti = -sti*rast*rast
      s1r = zet1dr - str
      s1i = zet1di - sti
      GO TO 100
   90 CONTINUE
      s1r = zet1dr - zet2dr
      s1i = zet1di - zet2di
  100 CONTINUE
      rs1 = s1r
      IF (dabs(rs1).GT.elim) GO TO 105
      IF (dabs(rs1).LT.alim) GO TO 120
C----------------------------------------------------------------------------
C     REFINE ESTIMATE AND TEST
C-------------------------------------------------------------------------
      aphi = xzabs(phidr,phidi)
      rs1 = rs1+dlog(aphi)
      IF (dabs(rs1).LT.elim) GO TO 120
  105 CONTINUE
      IF (rs1.GT.0.0d0) GO TO 320
C-----------------------------------------------------------------------
C     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
C-----------------------------------------------------------------------
      IF (zr.LT.0.0d0) GO TO 320
      nz = n
      DO 106 i=1,n
        yr(i) = zeror
        yi(i) = zeroi
  106 CONTINUE
      RETURN
  120 CONTINUE
      s1r = cyr(1)
      s1i = cyi(1)
      s2r = cyr(2)
      s2i = cyi(2)
      c1r = csrr(kflag)
      ascle = bry(kflag)
      DO 130 i=ib,n
        c2r = s2r
        c2i = s2i
        s2r = ckr*c2r - cki*c2i + s1r
        s2i = ckr*c2i + cki*c2r + s1i
        s1r = c2r
        s1i = c2i
        ckr = ckr + rzr
        cki = cki + rzi
        c2r = s2r*c1r
        c2i = s2i*c1r
        yr(i) = c2r
        yi(i) = c2i
        IF (kflag.GE.3) GO TO 130
        str = dabs(c2r)
        sti = dabs(c2i)
        c2m = dmax1(str,sti)
        IF (c2m.LE.ascle) GO TO 130
        kflag = kflag + 1
        ascle = bry(kflag)
        s1r = s1r*c1r
        s1i = s1i*c1r
        s2r = c2r
        s2i = c2i
        s1r = s1r*cssr(kflag)
        s1i = s1i*cssr(kflag)
        s2r = s2r*cssr(kflag)
        s2i = s2i*cssr(kflag)
        c1r = csrr(kflag)
  130 CONTINUE
  180 CONTINUE
      IF (mr.EQ.0) RETURN
C-----------------------------------------------------------------------
C     ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
C-----------------------------------------------------------------------
      nz = 0
      fmr = dble(float(mr))
      sgn = -dsign(pi,fmr)
C-----------------------------------------------------------------------
C     CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP.
C-----------------------------------------------------------------------
      csgni = sgn
      IF (yy.LE.0.0d0) csgni = -csgni
      ifn = inu + n - 1
      ang = fnf*sgn
      cspnr = dcos(ang)
      cspni = dsin(ang)
      IF (mod(ifn,2).EQ.0) GO TO 190
      cspnr = -cspnr
      cspni = -cspni
  190 CONTINUE
C-----------------------------------------------------------------------
C     CS=COEFF OF THE J FUNCTION TO GET THE I FUNCTION. I(FNU,Z) IS
C     COMPUTED FROM EXP(I*FNU*HPI)*J(FNU,-I*Z) WHERE Z IS IN THE FIRST
C     QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
C     CONJUGATION SINCE THE I FUNCTION IS REAL ON THE POSITIVE REAL AXIS
C-----------------------------------------------------------------------
      csr = sar*csgni
      csi = car*csgni
      in = mod(ifn,4) + 1
      c2r = cipr(in)
      c2i = cipi(in)
      str = csr*c2r + csi*c2i
      csi = -csr*c2i + csi*c2r
      csr = str
      asc = bry(1)
      iuf = 0
      kk = n
      kdflg = 1
      ib = ib - 1
      ic = ib - 1
      DO 290 k=1,n
        fn = fnu + dble(float(kk-1))
C-----------------------------------------------------------------------
C     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
C     FUNCTION ABOVE
C-----------------------------------------------------------------------
        IF (n.GT.2) GO TO 175
  172   CONTINUE
        phidr = phir(j)
        phidi = phii(j)
        argdr = argr(j)
        argdi = argi(j)
        zet1dr = zeta1r(j)
        zet1di = zeta1i(j)
        zet2dr = zeta2r(j)
        zet2di = zeta2i(j)
        asumdr = asumr(j)
        asumdi = asumi(j)
        bsumdr = bsumr(j)
        bsumdi = bsumi(j)
        j = 3 - j
        GO TO 210
  175   CONTINUE
        IF ((kk.EQ.n).AND.(ib.LT.n)) GO TO 210
        IF ((kk.EQ.ib).OR.(kk.EQ.ic)) GO TO 172
        CALL zunhj(znr, zni, fn, 0, tol, phidr, phidi, argdr,
     *   argdi, zet1dr, zet1di, zet2dr, zet2di, asumdr,
     *   asumdi, bsumdr, bsumdi)
  210   CONTINUE
        IF (kode.EQ.1) GO TO 220
        str = zbr + zet2dr
        sti = zbi + zet2di
        rast = fn/xzabs(str,sti)
        str = str*rast*rast
        sti = -sti*rast*rast
        s1r = -zet1dr + str
        s1i = -zet1di + sti
        GO TO 230
  220   CONTINUE
        s1r = -zet1dr + zet2dr
        s1i = -zet1di + zet2di
  230   CONTINUE
C-----------------------------------------------------------------------
C     TEST FOR UNDERFLOW AND OVERFLOW
C-----------------------------------------------------------------------
        rs1 = s1r
        IF (dabs(rs1).GT.elim) GO TO 280
        IF (kdflg.EQ.1) iflag = 2
        IF (dabs(rs1).LT.alim) GO TO 240
C-----------------------------------------------------------------------
C     REFINE  TEST AND SCALE
C-----------------------------------------------------------------------
        aphi = xzabs(phidr,phidi)
        aarg = xzabs(argdr,argdi)
        rs1 = rs1 + dlog(aphi) - 0.25d0*dlog(aarg) - aic
        IF (dabs(rs1).GT.elim) GO TO 280
        IF (kdflg.EQ.1) iflag = 1
        IF (rs1.LT.0.0d0) GO TO 240
        IF (kdflg.EQ.1) iflag = 3
  240   CONTINUE
        CALL zairy(argdr, argdi, 0, 2, air, aii, nai, idum)
        CALL zairy(argdr, argdi, 1, 2, dair, daii, ndai, idum)
        str = dair*bsumdr - daii*bsumdi
        sti = dair*bsumdi + daii*bsumdr
        str = str + (air*asumdr-aii*asumdi)
        sti = sti + (air*asumdi+aii*asumdr)
        ptr = str*phidr - sti*phidi
        pti = str*phidi + sti*phidr
        s2r = ptr*csr - pti*csi
        s2i = ptr*csi + pti*csr
        str = dexp(s1r)*cssr(iflag)
        s1r = str*dcos(s1i)
        s1i = str*dsin(s1i)
        str = s2r*s1r - s2i*s1i
        s2i = s2r*s1i + s2i*s1r
        s2r = str
        IF (iflag.NE.1) GO TO 250
        CALL zuchk(s2r, s2i, nw, bry(1), tol)
        IF (nw.EQ.0) GO TO 250
        s2r = zeror
        s2i = zeroi
  250   CONTINUE
        IF (yy.LE.0.0d0) s2i = -s2i
        cyr(kdflg) = s2r
        cyi(kdflg) = s2i
        c2r = s2r
        c2i = s2i
        s2r = s2r*csrr(iflag)
        s2i = s2i*csrr(iflag)
C-----------------------------------------------------------------------
C     ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
C-----------------------------------------------------------------------
        s1r = yr(kk)
        s1i = yi(kk)
        IF (kode.EQ.1) GO TO 270
        CALL zs1s2(zrr, zri, s1r, s1i, s2r, s2i, nw, asc, alim, iuf)
        nz = nz + nw
  270   CONTINUE
        yr(kk) = s1r*cspnr - s1i*cspni + s2r
        yi(kk) = s1r*cspni + s1i*cspnr + s2i
        kk = kk - 1
        cspnr = -cspnr
        cspni = -cspni
        str = csi
        csi = -csr
        csr = str
        IF (c2r.NE.0.0d0 .OR. c2i.NE.0.0d0) GO TO 255
        kdflg = 1
        GO TO 290
  255   CONTINUE
        IF (kdflg.EQ.2) GO TO 295
        kdflg = 2
        GO TO 290
  280   CONTINUE
        IF (rs1.GT.0.0d0) GO TO 320
        s2r = zeror
        s2i = zeroi
        GO TO 250
  290 CONTINUE
      k = n
  295 CONTINUE
      il = n - k
      IF (il.EQ.0) RETURN
C-----------------------------------------------------------------------
C     RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
C     K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
C     INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
C-----------------------------------------------------------------------
      s1r = cyr(1)
      s1i = cyi(1)
      s2r = cyr(2)
      s2i = cyi(2)
      csr = csrr(iflag)
      ascle = bry(iflag)
      fn = dble(float(inu+il))
      DO 310 i=1,il
        c2r = s2r
        c2i = s2i
        s2r = s1r + (fn+fnf)*(rzr*c2r-rzi*c2i)
        s2i = s1i + (fn+fnf)*(rzr*c2i+rzi*c2r)
        s1r = c2r
        s1i = c2i
        fn = fn - 1.0d0
        c2r = s2r*csr
        c2i = s2i*csr
        ckr = c2r
        cki = c2i
        c1r = yr(kk)
        c1i = yi(kk)
        IF (kode.EQ.1) GO TO 300
        CALL zs1s2(zrr, zri, c1r, c1i, c2r, c2i, nw, asc, alim, iuf)
        nz = nz + nw
  300   CONTINUE
        yr(kk) = c1r*cspnr - c1i*cspni + c2r
        yi(kk) = c1r*cspni + c1i*cspnr + c2i
        kk = kk - 1
        cspnr = -cspnr
        cspni = -cspni
        IF (iflag.GE.3) GO TO 310
        c2r = dabs(ckr)
        c2i = dabs(cki)
        c2m = dmax1(c2r,c2i)
        IF (c2m.LE.ascle) GO TO 310
        iflag = iflag + 1
        ascle = bry(iflag)
        s1r = s1r*csr
        s1i = s1i*csr
        s2r = ckr
        s2i = cki
        s1r = s1r*cssr(iflag)
        s1i = s1i*cssr(iflag)
        s2r = s2r*cssr(iflag)
        s2i = s2i*cssr(iflag)
        csr = csrr(iflag)
  310 CONTINUE
      RETURN
  320 CONTINUE
      nz = -1
      RETURN
      END