cuni2.f

Go to the documentation of this file.

      SUBROUTINE cuni2(Z, FNU, KODE, N, Y, NZ, NLAST, FNUL, TOL, ELIM,
     * ALIM)
C***BEGIN PROLOGUE  CUNI2
C***REFER TO  CBESI,CBESK
C
C     CUNI2 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  CAIRY,CUCHK,CUNHJ,CUOIK,R1MACH
C***END PROLOGUE  CUNI2
      COMPLEX AI, ARG, ASUM, BSUM, CFN, CI, CID, CIP, CONE, CRSC, CSCL,
     * csr, css, cy, czero, c1, c2, dai, phi, rz, s1, s2, y, z, zb,
     * zeta1, zeta2, zn, zar
      REAL AARG, AIC, ALIM, ANG, APHI, ASCLE, AY, BRY, CAR, C2I, C2M,
     * c2r, elim, fn, fnu, fnul, hpi, rs1, sar, tol, yy, r1mach
      INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST,
     * nn, nuf, nw, nz, idum
      dimension bry(3), y(n), cip(4), css(3), csr(3), cy(2)
      DATA czero,cone,ci/(0.0e0,0.0e0),(1.0e0,0.0e0),(0.0e0,1.0e0)/
      DATA cip(1),cip(2),cip(3),cip(4)/
     1 (1.0e0,0.0e0), (0.0e0,1.0e0), (-1.0e0,0.0e0), (0.0e0,-1.0e0)/
      DATA hpi, aic  /
     1      1.57079632679489662e+00,     1.265512123484645396e+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 = 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
      yy = aimag(z)
C-----------------------------------------------------------------------
C     ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI
C-----------------------------------------------------------------------
      zn = -z*ci
      zb = z
      cid = -ci
      inu = int(fnu)
      ang = hpi*(fnu-float(inu))
      car = cos(ang)
      sar = sin(ang)
      c2 = cmplx(car,sar)
      zar = c2
      in = inu + n - 1
      in = mod(in,4)
      c2 = c2*cip(in+1)
      IF (yy.GT.0.0e0) GO TO 10
      zn = conjg(-zn)
      zb = conjg(zb)
      cid = -cid
      c2 = conjg(c2)
   10 CONTINUE
C-----------------------------------------------------------------------
C     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
C-----------------------------------------------------------------------
      fn = amax1(fnu,1.0e0)
      CALL cunhj(zn, fn, 1, tol, phi, arg, zeta1, zeta2, asum, bsum)
      IF (kode.EQ.1) GO TO 20
      cfn = cmplx(fnu,0.0e0)
      s1 = -zeta1 + cfn*(cfn/(zb+zeta2))
      GO TO 30
   20 CONTINUE
      s1 = -zeta1 + zeta2
   30 CONTINUE
      rs1 = real(s1)
      IF (abs(rs1).GT.elim) GO TO 150
   40 CONTINUE
      nn = min0(2,nd)
      DO 90 i=1,nn
        fn = fnu + float(nd-i)
        CALL cunhj(zn, fn, 0, tol, phi, arg, zeta1, zeta2, asum, bsum)
        IF (kode.EQ.1) GO TO 50
        cfn = cmplx(fn,0.0e0)
        ay = abs(yy)
        s1 = -zeta1 + cfn*(cfn/(zb+zeta2)) + cmplx(0.0e0,ay)
        GO TO 60
   50   CONTINUE
        s1 = -zeta1 + zeta2
   60   CONTINUE
C-----------------------------------------------------------------------
C     TEST FOR UNDERFLOW AND OVERFLOW
C-----------------------------------------------------------------------
        rs1 = real(s1)
        IF (abs(rs1).GT.elim) GO TO 120
        IF (i.EQ.1) iflag = 2
        IF (abs(rs1).LT.alim) GO TO 70
C-----------------------------------------------------------------------
C     REFINE  TEST AND SCALE
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
        aphi = cabs(phi)
        aarg = cabs(arg)
        rs1 = rs1 + alog(aphi) - 0.25e0*alog(aarg) - aic
        IF (abs(rs1).GT.elim) GO TO 120
        IF (i.EQ.1) iflag = 1
        IF (rs1.LT.0.0e0) 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 cairy(arg, 0, 2, ai, nai, idum)
        CALL cairy(arg, 1, 2, dai, ndai, idum)
        s2 = phi*(ai*asum+dai*bsum)
        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 80
        CALL cuchk(s2, nw, bry(1), tol)
        IF (nw.NE.0) GO TO 120
   80   CONTINUE
        IF (yy.LE.0.0e0) s2 = conjg(s2)
        j = nd - i + 1
        s2 = s2*c2
        cy(i) = s2
        y(j) = s2*csr(iflag)
        c2 = c2*cid
   90 CONTINUE
      IF (nd.LE.2) GO TO 110
      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 100 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 100
        c2r = real(c2)
        c2i = aimag(c2)
        c2r = abs(c2r)
        c2i = abs(c2i)
        c2m = amax1(c2r,c2i)
        IF (c2m.LE.ascle) GO TO 100
        iflag = iflag + 1
        ascle = bry(iflag)
        s1 = s1*c1
        s2 = c2
        s1 = s1*css(iflag)
        s2 = s2*css(iflag)
        c1 = csr(iflag)
  100 CONTINUE
  110 CONTINUE
      RETURN
  120 CONTINUE
      IF (rs1.GT.0.0e0) GO TO 140
C-----------------------------------------------------------------------
C     SET UNDERFLOW AND UPDATE PARAMETERS
C-----------------------------------------------------------------------
      y(nd) = czero
      nz = nz + 1
      nd = nd - 1
      IF (nd.EQ.0) GO TO 110
      CALL cuoik(z, fnu, kode, 1, nd, y, 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 + float(nd-1)
      IF (fn.LT.fnul) GO TO 130
C      FN = AIMAG(CID)
C      J = NUF + 1
C      K = MOD(J,4) + 1
C      S1 = CIP(K)
C      IF (FN.LT.0.0E0) S1 = CONJG(S1)
C      C2 = C2*S1
      in = inu + nd - 1
      in = mod(in,4) + 1
      c2 = zar*cip(in)
      IF (yy.LE.0.0e0)c2=conjg(c2)
      GO TO 40
  130 CONTINUE
      nlast = nd
      RETURN
  140 CONTINUE
      nz = -1
      RETURN
  150 CONTINUE
      IF (rs1.GT.0.0e0) GO TO 140
      nz = n
      DO 160 i=1,n
        y(i) = czero
  160 CONTINUE
      RETURN
      END