zuni2.f

Go to the documentation of this file.

      SUBROUTINE zuni2(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL,
     * TOL, ELIM, ALIM)
C***BEGIN PROLOGUE  ZUNI2
C***REFER TO  ZBESI,ZBESK
C
C     ZUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF
C     UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I
C     OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO.
C
C     FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
C     EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
C     NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
C     FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
C     Y(I)=CZERO FOR I=NLAST+1,N
C
C***ROUTINES CALLED  ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,XZABS
C***END PROLOGUE  ZUNI2
C     COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS,
C    *CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN
      DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGI,
     * argr, ascle, asumi, asumr, bry, bsumi, bsumr, cidi, cipi, cipr,
     * coner, crsc, cscl, csrr, cssr, c1r, c2i, c2m, c2r, daii,
     * dair, elim, fn, fnu, fnul, hpi, phii, phir, rast, raz, rs1, rzi,
     * rzr, sti, str, s1i, s1r, s2i, s2r, tol, yi, yr, zbi, zbr, zeroi,
     * zeror, zeta1i, zeta1r, zeta2i, zeta2r, zi, zni, znr, zr, cyr,
     * cyi, d1mach, xzabs, car, sar
      INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST,
     * nn, nuf, nw, nz, idum
      dimension bry(3), yr(n), yi(n), cipr(4), cipi(4), cssr(3),
     * csrr(3), cyr(2), cyi(2)
      DATA zeror,zeroi,coner / 0.0d0, 0.0d0, 1.0d0 /
      DATA cipr(1),cipi(1),cipr(2),cipi(2),cipr(3),cipi(3),cipr(4),
     * cipi(4)/ 1.0d0,0.0d0, 0.0d0,1.0d0, -1.0d0,0.0d0, 0.0d0,-1.0d0/
      DATA hpi, aic  /
     1      1.57079632679489662d+00,     1.265512123484645396d+00/
C
      nz = 0
      nd = n
      nlast = 0
C-----------------------------------------------------------------------
C     COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
C     NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
C     EXP(ALIM)=EXP(ELIM)*TOL
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
C-----------------------------------------------------------------------
C     ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI
C-----------------------------------------------------------------------
      znr = zi
      zni = -zr
      zbr = zr
      zbi = zi
      cidi = -coner
      inu = int(sngl(fnu))
      ang = hpi*(fnu-dble(float(inu)))
      c2r = dcos(ang)
      c2i = dsin(ang)
      car = c2r
      sar = c2i
      in = inu + n - 1
      in = mod(in,4) + 1
      str = c2r*cipr(in) - c2i*cipi(in)
      c2i = c2r*cipi(in) + c2i*cipr(in)
      c2r = str
      IF (zi.GT.0.0d0) GO TO 10
      znr = -znr
      zbi = -zbi
      cidi = -cidi
      c2i = -c2i
   10 CONTINUE
C-----------------------------------------------------------------------
C     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
C-----------------------------------------------------------------------
      fn = dmax1(fnu,1.0d0)
      CALL zunhj(znr, zni, fn, 1, tol, phir, phii, argr, argi, zeta1r,
     * zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
      IF (kode.EQ.1) GO TO 20
      str = zbr + zeta2r
      sti = zbi + zeta2i
      rast = fn/xzabs(str,sti)
      str = str*rast*rast
      sti = -sti*rast*rast
      s1r = -zeta1r + str
      s1i = -zeta1i + sti
      GO TO 30
   20 CONTINUE
      s1r = -zeta1r + zeta2r
      s1i = -zeta1i + zeta2i
   30 CONTINUE
      rs1 = s1r
      IF (dabs(rs1).GT.elim) GO TO 150
   40 CONTINUE
      nn = min0(2,nd)
      DO 90 i=1,nn
        fn = fnu + dble(float(nd-i))
        CALL zunhj(znr, zni, fn, 0, tol, phir, phii, argr, argi,
     *   zeta1r, zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
        IF (kode.EQ.1) GO TO 50
        str = zbr + zeta2r
        sti = zbi + zeta2i
        rast = fn/xzabs(str,sti)
        str = str*rast*rast
        sti = -sti*rast*rast
        s1r = -zeta1r + str
        s1i = -zeta1i + sti + dabs(zi)
        GO TO 60
   50   CONTINUE
        s1r = -zeta1r + zeta2r
        s1i = -zeta1i + zeta2i
   60   CONTINUE
C-----------------------------------------------------------------------
C     TEST FOR UNDERFLOW AND OVERFLOW
C-----------------------------------------------------------------------
        rs1 = s1r
        IF (dabs(rs1).GT.elim) GO TO 120
        IF (i.EQ.1) iflag = 2
        IF (dabs(rs1).LT.alim) GO TO 70
C-----------------------------------------------------------------------
C     REFINE  TEST AND SCALE
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
        aphi = xzabs(phir,phii)
        aarg = xzabs(argr,argi)
        rs1 = rs1 + dlog(aphi) - 0.25d0*dlog(aarg) - aic
        IF (dabs(rs1).GT.elim) GO TO 120
        IF (i.EQ.1) iflag = 1
        IF (rs1.LT.0.0d0) GO TO 70
        IF (i.EQ.1) iflag = 3
   70   CONTINUE
C-----------------------------------------------------------------------
C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
C     EXPONENT EXTREMES
C-----------------------------------------------------------------------
        CALL zairy(argr, argi, 0, 2, air, aii, nai, idum)
        CALL zairy(argr, argi, 1, 2, dair, daii, ndai, idum)
        str = dair*bsumr - daii*bsumi
        sti = dair*bsumi + daii*bsumr
        str = str + (air*asumr-aii*asumi)
        sti = sti + (air*asumi+aii*asumr)
        s2r = phir*str - phii*sti
        s2i = phir*sti + phii*str
        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 80
        CALL zuchk(s2r, s2i, nw, bry(1), tol)
        IF (nw.NE.0) GO TO 120
   80   CONTINUE
        IF (zi.LE.0.0d0) s2i = -s2i
        str = s2r*c2r - s2i*c2i
        s2i = s2r*c2i + s2i*c2r
        s2r = str
        cyr(i) = s2r
        cyi(i) = s2i
        j = nd - i + 1
        yr(j) = s2r*csrr(iflag)
        yi(j) = s2i*csrr(iflag)
        str = -c2i*cidi
        c2i = c2r*cidi
        c2r = str
   90 CONTINUE
      IF (nd.LE.2) GO TO 110
      raz = 1.0d0/xzabs(zr,zi)
      str = zr*raz
      sti = -zi*raz
      rzr = (str+str)*raz
      rzi = (sti+sti)*raz
      bry(2) = 1.0d0/bry(1)
      bry(3) = d1mach(2)
      s1r = cyr(1)
      s1i = cyi(1)
      s2r = cyr(2)
      s2i = cyi(2)
      c1r = csrr(iflag)
      ascle = bry(iflag)
      k = nd - 2
      fn = dble(float(k))
      DO 100 i=3,nd
        c2r = s2r
        c2i = s2i
        s2r = s1r + (fnu+fn)*(rzr*c2r-rzi*c2i)
        s2i = s1i + (fnu+fn)*(rzr*c2i+rzi*c2r)
        s1r = c2r
        s1i = c2i
        c2r = s2r*c1r
        c2i = s2i*c1r
        yr(k) = c2r
        yi(k) = c2i
        k = k - 1
        fn = fn - 1.0d0
        IF (iflag.GE.3) GO TO 100
        str = dabs(c2r)
        sti = dabs(c2i)
        c2m = dmax1(str,sti)
        IF (c2m.LE.ascle) GO TO 100
        iflag = iflag + 1
        ascle = bry(iflag)
        s1r = s1r*c1r
        s1i = s1i*c1r
        s2r = c2r
        s2i = c2i
        s1r = s1r*cssr(iflag)
        s1i = s1i*cssr(iflag)
        s2r = s2r*cssr(iflag)
        s2i = s2i*cssr(iflag)
        c1r = csrr(iflag)
  100 CONTINUE
  110 CONTINUE
      RETURN
  120 CONTINUE
      IF (rs1.GT.0.0d0) GO TO 140
C-----------------------------------------------------------------------
C     SET UNDERFLOW AND UPDATE PARAMETERS
C-----------------------------------------------------------------------
      yr(nd) = zeror
      yi(nd) = zeroi
      nz = nz + 1
      nd = nd - 1
      IF (nd.EQ.0) GO TO 110
      CALL zuoik(zr, zi, fnu, kode, 1, nd, yr, yi, nuf, tol, elim, alim)
      IF (nuf.LT.0) GO TO 140
      nd = nd - nuf
      nz = nz + nuf
      IF (nd.EQ.0) GO TO 110
      fn = fnu + dble(float(nd-1))
      IF (fn.LT.fnul) GO TO 130
C      FN = CIDI
C      J = NUF + 1
C      K = MOD(J,4) + 1
C      S1R = CIPR(K)
C      S1I = CIPI(K)
C      IF (FN.LT.0.0D0) S1I = -S1I
C      STR = C2R*S1R - C2I*S1I
C      C2I = C2R*S1I + C2I*S1R
C      C2R = STR
      in = inu + nd - 1
      in = mod(in,4) + 1
      c2r = car*cipr(in) - sar*cipi(in)
      c2i = car*cipi(in) + sar*cipr(in)
      IF (zi.LE.0.0d0) c2i = -c2i
      GO TO 40
  130 CONTINUE
      nlast = nd
      RETURN
  140 CONTINUE
      nz = -1
      RETURN
  150 CONTINUE
      IF (rs1.GT.0.0d0) GO TO 140
      nz = n
      DO 160 i=1,n
        yr(i) = zeror
        yi(i) = zeroi
  160 CONTINUE
      RETURN
      END