cuni1.f

Go to the documentation of this file.

      SUBROUTINE cuni1(Z, FNU, KODE, N, Y, NZ, NLAST, FNUL, TOL, ELIM,
     * ALIM)
C***BEGIN PROLOGUE  CUNI1
C***REFER TO  CBESI,CBESK
C
C     CUNI1 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  CUCHK,CUNIK,CUOIK,R1MACH
C***END PROLOGUE  CUNI1
      COMPLEX CFN, CONE, CRSC, CSCL, CSR, CSS, CWRK, CZERO, C1, C2,
     * phi, rz, sum, s1, s2, y, z, zeta1, zeta2, cy
      REAL ALIM, APHI, ASCLE, BRY, C2I, C2M, C2R, ELIM, FN, FNU, FNUL,
     * rs1, tol, yy, r1mach
      INTEGER I, IFLAG, INIT, K, KODE, M, N, ND, NLAST, NN, NUF, NW, NZ
      dimension bry(3), y(n), cwrk(16), css(3), csr(3), cy(2)
      DATA czero, cone / (0.0e0,0.0e0), (1.0e0,0.0e0) /
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 = cmplx(1.0e0/tol,0.0e0)
      crsc = cmplx(tol,0.0e0)
      css(1) = cscl
      css(2) = cone
      css(3) = crsc
      csr(1) = crsc
      csr(2) = cone
      csr(3) = cscl
      bry(1) = 1.0e+3*r1mach(1)/tol
C-----------------------------------------------------------------------
C     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
C-----------------------------------------------------------------------
      fn = amax1(fnu,1.0e0)
      init = 0
      CALL cunik(z, fn, 1, 1, tol, init, phi, zeta1, zeta2, sum, cwrk)
      IF (kode.EQ.1) GO TO 10
      cfn = cmplx(fn,0.0e0)
      s1 = -zeta1 + cfn*(cfn/(z+zeta2))
      GO TO 20
   10 CONTINUE
      s1 = -zeta1 + zeta2
   20 CONTINUE
      rs1 = real(s1)
      IF (abs(rs1).GT.elim) GO TO 130
   30 CONTINUE
      nn = min0(2,nd)
      DO 80 i=1,nn
        fn = fnu + float(nd-i)
        init = 0
        CALL cunik(z, fn, 1, 0, tol, init, phi, zeta1, zeta2, sum, cwrk)
        IF (kode.EQ.1) GO TO 40
        cfn = cmplx(fn,0.0e0)
        yy = aimag(z)
        s1 = -zeta1 + cfn*(cfn/(z+zeta2)) + cmplx(0.0e0,yy)
        GO TO 50
   40   CONTINUE
        s1 = -zeta1 + zeta2
   50   CONTINUE
C-----------------------------------------------------------------------
C     TEST FOR UNDERFLOW AND OVERFLOW
C-----------------------------------------------------------------------
        rs1 = real(s1)
        IF (abs(rs1).GT.elim) GO TO 110
        IF (i.EQ.1) iflag = 2
        IF (abs(rs1).LT.alim) GO TO 60
C-----------------------------------------------------------------------
C     REFINE  TEST AND SCALE
C-----------------------------------------------------------------------
        aphi = cabs(phi)
        rs1 = rs1 + alog(aphi)
        IF (abs(rs1).GT.elim) GO TO 110
        IF (i.EQ.1) iflag = 1
        IF (rs1.LT.0.0e0) GO TO 60
        IF (i.EQ.1) iflag = 3
   60   CONTINUE
C-----------------------------------------------------------------------
C     SCALE S1 IF CABS(S1).LT.ASCLE
C-----------------------------------------------------------------------
        s2 = phi*sum
        c2r = real(s1)
        c2i = aimag(s1)
        c2m = exp(c2r)*real(css(iflag))
        s1 = cmplx(c2m,0.0e0)*cmplx(cos(c2i),sin(c2i))
        s2 = s2*s1
        IF (iflag.NE.1) GO TO 70
        CALL cuchk(s2, nw, bry(1), tol)
        IF (nw.NE.0) GO TO 110
   70   CONTINUE
        m = nd - i + 1
        cy(i) = s2
        y(m) = s2*csr(iflag)
   80 CONTINUE
      IF (nd.LE.2) GO TO 100
      rz = cmplx(2.0e0,0.0e0)/z
      bry(2) = 1.0e0/bry(1)
      bry(3) = r1mach(2)
      s1 = cy(1)
      s2 = cy(2)
      c1 = csr(iflag)
      ascle = bry(iflag)
      k = nd - 2
      fn = float(k)
      DO 90 i=3,nd
        c2 = s2
        s2 = s1 + cmplx(fnu+fn,0.0e0)*rz*s2
        s1 = c2
        c2 = s2*c1
        y(k) = c2
        k = k - 1
        fn = fn - 1.0e0
        IF (iflag.GE.3) GO TO 90
        c2r = real(c2)
        c2i = aimag(c2)
        c2r = abs(c2r)
        c2i = abs(c2i)
        c2m = amax1(c2r,c2i)
        IF (c2m.LE.ascle) GO TO 90
        iflag = iflag + 1
        ascle = bry(iflag)
        s1 = s1*c1
        s2 = c2
        s1 = s1*css(iflag)
        s2 = s2*css(iflag)
        c1 = csr(iflag)
   90 CONTINUE
  100 CONTINUE
      RETURN
C-----------------------------------------------------------------------
C     SET UNDERFLOW AND UPDATE PARAMETERS
C-----------------------------------------------------------------------
  110 CONTINUE
      IF (rs1.GT.0.0e0) GO TO 120
      y(nd) = czero
      nz = nz + 1
      nd = nd - 1
      IF (nd.EQ.0) GO TO 100
      CALL cuoik(z, fnu, kode, 1, nd, y, 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 + float(nd-1)
      IF (fn.GE.fnul) GO TO 30
      nlast = nd
      RETURN
  120 CONTINUE
      nz = -1
      RETURN
  130 CONTINUE
      IF (rs1.GT.0.0e0) GO TO 120
      nz = n
      DO 140 i=1,n
        y(i) = czero
  140 CONTINUE
      RETURN
      END