zuni1.f

Go to the documentation of this file.

      SUBROUTINE zuni1(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL,
     * tol, elim, alim)
C***BEGIN PROLOGUE  ZUNI1
C***REFER TO  ZBESI,ZBESK
C
C     ZUNI1 COMPUTES I(FNU,Z)  BY MEANS OF THE UNIFORM ASYMPTOTIC
C     EXPANSION FOR I(FNU,Z) IN -PI/3.LE.ARG Z.LE.PI/3.
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  ZUCHK,ZUNIK,ZUOIK,D1MACH,XZABS
C***END PROLOGUE  ZUNI1
C     COMPLEX CFN,CONE,CRSC,CSCL,CSR,CSS,CWRK,CZERO,C1,C2,PHI,RZ,SUM,S1,
C    *S2,Y,Z,ZETA1,ZETA2
      DOUBLE PRECISION alim, aphi, ascle, bry, coner, crsc,
     * cscl, csrr, cssr, cwrki, cwrkr, c1r, c2i, c2m, c2r, elim, fn,
     * fnu, fnul, phii, phir, rast, rs1, rzi, rzr, sti, str, sumi,
     * sumr, s1i, s1r, s2i, s2r, tol, yi, yr, zeroi, zeror, zeta1i,
     * zeta1r, zeta2i, zeta2r, zi, zr, cyr, cyi, d1mach, xzabs
      INTEGER i, iflag, init, k, kode, m, n, nd, nlast, nn, nuf, nw, nz
      dimension bry(3), yr(n), yi(n), cwrkr(16), cwrki(16), cssr(3),
     * csrr(3), cyr(2), cyi(2)
      DATA zeror,zeroi,coner / 0.0d0, 0.0d0, 1.0d0 /
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     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
C-----------------------------------------------------------------------
      fn = dmax1(fnu,1.0d0)
      init = 0
      CALL zunik(zr, zi, fn, 1, 1, tol, init, phir, phii, zeta1r,
     * zeta1i, zeta2r, zeta2i, sumr, sumi, cwrkr, cwrki)
      IF (kode.EQ.1) go to 10
      str = zr + zeta2r
      sti = zi + zeta2i
      rast = fn/xzabs(str,sti)
      str = str*rast*rast
      sti = -sti*rast*rast
      s1r = -zeta1r + str
      s1i = -zeta1i + sti
      go to 20
   10 CONTINUE
      s1r = -zeta1r + zeta2r
      s1i = -zeta1i + zeta2i
   20 CONTINUE
      rs1 = s1r
      IF (dabs(rs1).GT.elim) go to 130
   30 CONTINUE
      nn = min0(2,nd)
      DO 80 i=1,nn
        fn = fnu + dble(float(nd-i))
        init = 0
        CALL zunik(zr, zi, fn, 1, 0, tol, init, phir, phii, zeta1r,
     *   zeta1i, zeta2r, zeta2i, sumr, sumi, cwrkr, cwrki)
        IF (kode.EQ.1) go to 40
        str = zr + zeta2r
        sti = zi + zeta2i
        rast = fn/xzabs(str,sti)
        str = str*rast*rast
        sti = -sti*rast*rast
        s1r = -zeta1r + str
        s1i = -zeta1i + sti + zi
        go to 50
   40   CONTINUE
        s1r = -zeta1r + zeta2r
        s1i = -zeta1i + zeta2i
   50   CONTINUE
C-----------------------------------------------------------------------
C     TEST FOR UNDERFLOW AND OVERFLOW
C-----------------------------------------------------------------------
        rs1 = s1r
        IF (dabs(rs1).GT.elim) go to 110
        IF (i.EQ.1) iflag = 2
        IF (dabs(rs1).LT.alim) go to 60
C-----------------------------------------------------------------------
C     REFINE  TEST AND SCALE
C-----------------------------------------------------------------------
        aphi = xzabs(phir,phii)
        rs1 = rs1 + dlog(aphi)
        IF (dabs(rs1).GT.elim) go to 110
        IF (i.EQ.1) iflag = 1
        IF (rs1.LT.0.0d0) go to 60
        IF (i.EQ.1) iflag = 3
   60   CONTINUE
C-----------------------------------------------------------------------
C     SCALE S1 IF CABS(S1).LT.ASCLE
C-----------------------------------------------------------------------
        s2r = phir*sumr - phii*sumi
        s2i = phir*sumi + phii*sumr
        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 70
        CALL zuchk(s2r, s2i, nw, bry(1), tol)
        IF (nw.NE.0) go to 110
   70   CONTINUE
        cyr(i) = s2r
        cyi(i) = s2i
        m = nd - i + 1
        yr(m) = s2r*csrr(iflag)
        yi(m) = s2i*csrr(iflag)
   80 CONTINUE
      IF (nd.LE.2) go to 100
      rast = 1.0d0/xzabs(zr,zi)
      str = zr*rast
      sti = -zi*rast
      rzr = (str+str)*rast
      rzi = (sti+sti)*rast
      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 90 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 90
        str = dabs(c2r)
        sti = dabs(c2i)
        c2m = dmax1(str,sti)
        IF (c2m.LE.ascle) go to 90
        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)
   90 CONTINUE
  100 CONTINUE
      RETURN
C-----------------------------------------------------------------------
C     SET UNDERFLOW AND UPDATE PARAMETERS
C-----------------------------------------------------------------------
  110 CONTINUE
      IF (rs1.GT.0.0d0) go to 120
      yr(nd) = zeror
      yi(nd) = zeroi
      nz = nz + 1
      nd = nd - 1
      IF (nd.EQ.0) go to 100
      CALL zuoik(zr, zi, fnu, kode, 1, nd, yr, yi, nuf, tol, elim, alim)
      IF (nuf.LT.0) go to 120
      nd = nd - nuf
      nz = nz + nuf
      IF (nd.EQ.0) go to 100
      fn = fnu + dble(float(nd-1))
      IF (fn.GE.fnul) go to 30
      nlast = nd
      RETURN
  120 CONTINUE
      nz = -1
      RETURN
  130 CONTINUE
      IF (rs1.GT.0.0d0) go to 120
      nz = n
      DO 140 i=1,n
        yr(i) = zeror
        yi(i) = zeroi
  140 CONTINUE
      RETURN
      END