*--------------------------------------------------------------------*
      COMPLEX*16 FUNCTION HANKP(I,N,X)
*--------------------------------------------------------------------*
*------------------------------------------------*
*     THIS FUNCTION COMPUTES THE DERIVATIVE OF   *
*     THE HANKEL FUNCTION OF ORDER N AND SUPER   *
*     1 OR SUPER 2.                              *
*------------------------------------------------*
      COMPLEX*16 X,HANK
      HANKP=-HANK(I,N+1,X)+(FLOAT(N)/(X))*HANK(I,N,X)
      RETURN
      END
*--------------------------------------------------------------------*
      COMPLEX*16 FUNCTION BESSJP(N,X)
*--------------------------------------------------------------------*
*------------------------------------------------*
*     THIS FUNCTION COMPUTES THE DERIVATIVE OF   *
*     THE BESSEL FUNCTION OF THE FIRST KIND.     *
*------------------------------------------------*
      COMPLEX*16 X,BESSJ
      BESSJP=-BESSJ(N+1,X)+(FLOAT(N)/(X))*BESSJ(N,X)
      RETURN
      END
*--------------------------------------------------------------------*
      COMPLEX*16 FUNCTION HANK(I,N,X)
*--------------------------------------------------------------------*
*------------------------------------------------*
*   THIS FUNCTION RETURNS THE HANKEL FUNCTION    *
*   SUPER 1 AND 2 OF ORDER N AND ARGUMENT X      *
*   WRITTEN BY: LOUIS WILSON                     *
*------------------------------------------------*
      COMPLEX*16 J,X,BESSJ,BESSY
      PARAMETER(PI=3.1415926536)
      J=dcmplx(0.0,1.0)
      IF(I.EQ.1)THEN
         IF(N.GE.0)THEN
            HANK=BESSJ(N,X)+J*BESSY(N,X)
         ELSE
            N=-N
            HANK=BESSJ(N,X)+J*BESSY(N,X)
            HANK=CDEXP(N*PI*J)*HANK
         ENDIF
      ELSE
         IF(N.GE.0)THEN
            HANK=BESSJ(N,X)-J*BESSY(N,X)
c      write(*,*) 'hank,bessj,bessy',hank,bessj(n,x),bessy(n,x)	    
         ELSE
            N=-N
            HANK=BESSJ(N,X)-J*BESSY(N,X)
            HANK=CDEXP(-N*PI*J)*HANK
         ENDIF
      ENDIF
      RETURN
      END
*--------------------------------------------------------------------*
      COMPLEX*16 FUNCTION BESSJ(N,X)
*--------------------------------------------------------------------*
      COMPLEX*16 CBS(116),X,BESSJ0
      IF (N.EQ.0) THEN
        BESSJ=BESSJ0(X)
      ELSE
        CALL CBJNS(X,N+1,CBS)
        BESSJ=CBS(N+1)
      ENDIF
      RETURN
      END
*--------------------------------------------------------------------*
      COMPLEX*16 FUNCTION BESSJ0(X)
*--------------------------------------------------------------------*
*-------------------------------------------------*
*     FUNCTION TO COMPUTE THE BESSEL FUNCTION     *
*     J0 AS A FUNCTION OF X. USES POLYNOMIAL      *
*     APPROXIMATION FROM ABRAMOWITZ AND ATEGUN    *
*     WRITTEN BY : LOUIS WILSON                   *
*     COMMENTS--WORKS FINE.                       *
*     POLYNOMIAL COEFFICIENTS AS GIVEN IN THE     *
*     ABOVE REFERENCE.                            *
*-------------------------------------------------*
      COMPLEX*16 X,X3,FUNC,THETA
      DATA C1,C2,C3,C4,C5,C6/2.2499997,1.2656208,.3163866,
     $     .0444479,.0039444,.0002100/
      DATA F1,F2,F3,F4,F5,F6,F7/.79788456,.00000077,.00552740,
     $     .00009512,.00137237,.00072085,.00014476/
      DATA T1,T2,T3,T4,T5,T6,T7/.78539816,.04166397,.00003954,
     $     .00262573,.00054125,.00029333,.00013558/
*-------------------------------------------------------*
*   COMPUTATION OF BESSEL FUNCTION USING THE POLYNOMIAL *
*   APPROXIMATION, WHICH IS USUALLY GOOD ENOUGH.        *
*-------------------------------------------------------*
      IF(CDABS(X).LT.3)THEN
         X3=X/3.
         BESSJ0=1.-C1*X3**2+C2*X3**4-C3*X3**6+C4*X3**8-C5*X3**10
     $          +C6*X3**12
      ELSE
         X3=3./X
         FUNC=F1-F2*X3-F3*X3**2-F4*X3**3+F5*X3**4-F6*X3**5+F7*X3**6
         THETA=X-T1-T2*X3-T3*X3**2+T4*X3**3-T5*X3**4-T6*X3**5+T7*X3**6
         BESSJ0=(1./DSQRT(CDABS(X)))*FUNC*CDCOS(THETA)
      ENDIF
      RETURN
      END
*--------------------------------------------------------------------*
      COMPLEX*16 FUNCTION BESSY(N,X)
*--------------------------------------------------------------------*
      COMPLEX*16 X,X2,BY,BYP,BYM,BESSY0,BESSY1
      IF(N.EQ.0)THEN
         BY=BESSY0(X)
      ELSEIF(N.EQ.1) THEN
         BY=BESSY1(X)
      ELSE
         X2=2./X
         BY=BESSY1(X)
         BYM=BESSY0(X)
         DO 10 I=1,N-1
            BYP=I*X2*BY-BYM
            BYM=BY
            BY=BYP
10       CONTINUE
      ENDIF
      BESSY=BY
      RETURN
      END
*--------------------------------------------------------------------*
      COMPLEX*16 FUNCTION BESSY0(X)
*--------------------------------------------------------------------*
*     FUNCTION TO COMPUTE THE BESSEL FUNCTION    *
*     Y0 AS A FUNCTION OF X. USES POLYNOMIAL     *
*     APPROXIMATION FROM ABRAMOWITZ AND STEGUN   *
*     WRITTEN BY : LOUIS WILSON                  *
*     COMMENTS--WORKS FINE.                      *
*     POLYNOMIAL COEFFICIENTS AS GIVEN IN THE    *
*     ABOVE REFERENCE.                           *
*------------------------------------------------*
      COMPLEX*16 X,X3,FUNC,THETA,BJ0
      PARAMETER(PI=3.1415926536)
      DATA C1,C2,C3,C4,C5,C6/2.2499997,1.2656208,.3163866,
     $     .0444479,.0039444,.0002100/
      DATA D1,D2,D3,D4,D5,D6,D7/.36746691,.60559366,.74350384,
     $     .25300117,.04261214,.00427916,.00024846/
      DATA F1,F2,F3,F4,F5,F6,F7/.79788456,.00000077,.00552740,
     $     .00009512,.00137237,.00072085,.00014476/
      DATA T1,T2,T3,T4,T5,T6,T7/.78539816,.04166397,.00003954,
     $     .00262573,.00054125,.00029333,.00013558/
*-------------------------------------------------------*
*   COMPUTATION OF BESSEL FUNCTION USING THE POLYNOMIAL *
*   APPROXIMATION, WHICH IS USUALLY GOOD ENOUGH.        *
*-------------------------------------------------------*
      IF(CDABS(X).LT.3.)THEN
         X3=X/3.
         BJ0=1.-C1*X3**2+C2*X3**4-C3*X3**6+C4*X3**8-C5*X3**10
     $          +C6*X3**12
         BESSY0=(2/PI)*CDLOG(X/2)*BJ0+D1+D2*X3**2-D3*X3**4+D4*X3**6
     $          -D5*X3**8+D6*X3**10-D7*X3**12
      ELSE
         X3=3./X
         FUNC=F1-F2*X3-F3*X3**2-F4*X3**3+F5*X3**4-F6*X3**5+F7*X3**6
         THETA=X-T1-T2*X3-T3*X3**2+T4*X3**3-T5*X3**4-T6*X3**5+T7*X3**6
         BESSY0=(1./DSQRT(CDABS(X)))*FUNC*CDSIN(THETA)
      ENDIF
      RETURN
      END
*--------------------------------------------------------------------*
      COMPLEX*16 FUNCTION BESSY1(X)
*--------------------------------------------------------------------*
*-------------------------------------------------*
*     FUNCTION TO COMPUTE THE BESSEL FUNCTION     *
*     Y1 AS A FUNCTION OF X. USES POLYNOMIAL      *
*     APPROXIMATION FROM ABRAMOWITZ AND STEGUN    *
*     WRITTEN BY : LOUIS WILSON                   *
*     COMMENTS-- WORKS FINE.                      *
*     POLYNOMIAL COEFFICIENTS AS GIVEN IN THE     *
*     ABOVE REFERENCE.                            *
*-------------------------------------------------*
      COMPLEX*16 X,X3,FUNC,THETA,BJ1
      PARAMETER(PI=3.1415926536)
      DATA C1,C2,C3,C4,C5,C6/.56249985,.21093573,.03954289,
     $     .00443319,.00031716,.00001109/
      DATA D1,D2,D3,D4,D5,D6,D7/.6366198,.2212091,2.1682709,
     $     1.3164827,.3123951,.0400976,.0027873/
      DATA F1,F2,F3,F4,F5,F6,F7/.79788456,.00000156,.01659667,
     $     .00017105,.00249511,.00113653,.00020033/
      DATA T1,T2,T3,T4,T5,T6,T7/2.35619449,.12499612,.00005650,
     $     .00637879,.00074348,.00079824,.00029166/
*-------------------------------------------------------*
*   COMPUTATION OF BESSEL FUNCTION USING THE POLYNOMIAL *
*   APPROXIMATION, WHICH IS USUALLY GOOD ENOUGH.        *
*-------------------------------------------------------*
      IF(CDABS(X).LT.3.)THEN
         X3=X/3
         BJ1=X*(.5-C1*X3**2+C2*X3**4-C3*X3**6+C4*X3**8-C5*X3**10
     $          +C6*X3**12)
         BESSY1=(2/PI)*CDLOG(X/2)*BJ1+(-D1+D2*X3**2+D3*X3**4
     $          -D4*X3**6+D5*X3**8-D6*X3**10+D7*X3**12)/X
      ELSE
         X3=3/X
         FUNC=F1+F2*X3+F3*X3**2+F4*X3**3-F5*X3**4+F6*X3**5-F7*X3**6
         THETA=X-T1+T2*X3+T3*X3**2-T4*X3**3+T5*X3**4+T6*X3**5-T7*X3**6
         BESSY1=(1./DSQRT(CDABS(X)))*FUNC*CDSIN(THETA)
      ENDIF
      RETURN
      END
C-----------------------------------------------------------------------
C  IMSL Name:  CBJNS/DCBJNS (Single/Double precision version)
C
C  Computer:   VAX/SINGLE
C
C  Revised:    April 30, 1986
C
C  Purpose:    Evaluate a sequence of Bessel functions of the first
C              kind with integer order and complex arguments.
C
C  Usage:      CALL CBJNS (Z, N, CBS)
C
C  Arguments:
C     Z      - Complex argument for which the sequence of Bessel
C              functions is to be evaluated.  (Input)
C     N      - Number of elements in the sequence.  (Input)
C     CBS    - Vector of length N containing the values of the
C              function through the series.  (Output)
C
C  Remarks:
C  1. CBS(I) contains the value of the Bessel function of order I-1
C     at z for I = 1 to N.
C
C  2. The algorithm is based on a code which uses backward recurrsion
C     with strict error control.
C
C  GAMS:       C10a2
C
C  Chapter:    SFUN/LIBRARY Bessel Functions
C
C  Copyright:  1986 by IMSL, Inc.  All Rights Reserved.
C
C  Warranty:   IMSL warrants only that IMSL testing has been applied
C              to this code.  No other warranty, expressed or implied,
C              is applicable.
C
C---------------------------------------------------------------------
C
*--------------------------------------------------------------------*
      SUBROUTINE CBJNS (Z, N, CBS)
*--------------------------------------------------------------------*
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    N
      COMPLEX*16    Z, CBS(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    IMRECR, IPOS, JMAG, K, L, M, MLAST, MRECUR, NBMZ,
     $           NCALC, NEND, NEWNN, NN, NSTART
      REAL*8       SMALLZ, TEMP, TEMPB, TEMPC, TEST, TMPAR4, TOVER
      COMPLEX*16    CTEMP, CTEMPA, CTEMPB, P, PLAST, POLD, PSAVE, SUM,
     $           ZINV, ZTEMP
      LOGICAL    FLAG
C                                  SPECIFICATIONS FOR SAVE VARIABLES
      INTEGER    MAGZ
      REAL*8       EXPARG, RSIGN, RTEN, SIG, TEMPA
      LOGICAL    INITFL
      SAVE       EXPARG, INITFL, MAGZ, RSIGN, RTEN, SIG, TEMPA
C                                  SPECIFICATIONS FOR INTRINSICS
C     INTRINSIC  ABS,AIMAG,ALOG10,AMAX1,CABS,CMPLX,COS,EXP,INT,MIN0,
C    *REAL,SIN,SQRT
c      INTRINSIC  ABS, AIMAG, ALOG10, AMAX1, CABS, CMPLX, COS, EXP,
      INTRINSIC  ABS, AIMAG, ALOG10, AMAX1, CABS, COS, EXP,
     $           INT, MIN0, REAL, SIN, SQRT
      INTEGER    INT, MIN0
      REAL       ABS, AIMAG, ALOG10, AMAX1, CABS, COS, EXP, REAL, SIN,
     $           SQRT
c      COMPLEX    CMPLX
C                                  SPECIFICATIONS FOR SUBROUTINES
      EXTERNAL   CSET
C                                  SPECIFICATIONS FOR FUNCTIONS
C
      DATA INITFL/.FALSE./
C
C                                  Check N
      IF (N .LE. 0) THEN
         WRITE (*,*) 'The number of function values to be '//
     $               'calculated, N = %(I1), must be at least 1.'
         GO TO 9000
      END IF
C                                  Check Z
      IF (CDABS(Z) .GT. 1.0E+05) THEN
         WRITE (*,*) 'The complex absolute value of the input '//
     $               'argument, Z = %(C1), must not be greater than '//
     $               '10,000.'
         ZTEMP = dcmplx(1.7010200E+38,0.0E0)
         CALL CSET (N, ZTEMP, CBS, 1)
         GO TO 9000
      END IF
C                                  Initializations
      TEMPA = CDABS(Z)
      MAGZ = INT(TEMPA)
      IF (.NOT.INITFL) THEN
         SIG = -ALOG10(5.9618401E-08)
         EXPARG = ALOG10(1.7010200E+38)
         RTEN = INT(EXPARG)
         RSIGN = 1.0E0
         INITFL = .TRUE.
      END IF
      NCALC = N
      FLAG = .FALSE.
C                                  Use 2-term ascending series for
C                                  small Z
      TMPAR4 = TEMPA**4
      SMALLZ = .1E0**SIG
      IF (TMPAR4 .GE. SMALLZ) THEN
C                                  Initialize the calculation of the P*S
         NBMZ = N - MAGZ
         NN = MAGZ + 1
         ZINV = 1.0E0/Z
         PLAST = dcmplx(1.0E0,0.0E0)
         P = RSIGN*2*NN*ZINV
         TEST = 2.0E0*1.0E1**SIG
         M = 0
         IF (NBMZ .GE. 3) THEN
C                                  Calculate P*S until NN=N-1. Check for
C                                  possible overflow
            TOVER = 1.0E1**(RTEN-SIG)
            NSTART = MAGZ + 2
            NEND = N - 1
            DO 40  NN=NSTART, NEND
               POLD = PLAST
               PLAST = P
               P = RSIGN*(2*NN*PLAST*ZINV-POLD)
               IF (CDABS(P) .GT. TOVER) THEN
                  FLAG = .TRUE.
                  NEWNN = NN
                  NSTART = NEWNN + 1
C                                  To avoid overflow, normalize P*S by
C                                  dividing by TOVER. Calculate P*S
C                                  until unnormalized P would overflow
                  P = P/TOVER
                  PLAST = PLAST/TOVER
                  PSAVE = P
                  CTEMP = PLAST
                  TEST = 1.0E1**(SIG+SIG)
   10             NEWNN = NEWNN + 1
                  POLD = PLAST
                  PLAST = P
                  P = RSIGN*(2*NEWNN*PLAST*ZINV-POLD)
                  TMPAR4 = CDABS(P)
                  IF (TMPAR4**2 .LE. TEST) GO TO 10
C                                  Calculate backward test, and find
C                                  NCALC, the highest NEWNN such that
C                                  the test is passed.
                  TEMPB = CDABS(PLAST)/CDABS(POLD)
                  TEMPC = NEWNN/TEMPA
                  IF (TEMPB+1.0E0/TEMPB .GT. 2.0E0*TEMPC) TEMPB =
     $                TEMPC + SQRT(TEMPC**2-1.0E0)
                  TEST = .5E0*(1.0E0-1.0E0/TEMPB**2)/1.0E1**SIG
                  TEST = (CDABS(PLAST)*CDABS(POLD)*TEST)**2
                  P = PLAST*TOVER
                  NEWNN = NEWNN - 1
                  NEND = MIN0(N,NEWNN)
                  DO 20  NCALC=NSTART, NEND
                     POLD = CTEMP
                     CTEMP = PSAVE
                     PSAVE = RSIGN*(2*NEWNN*CTEMP*ZINV-POLD)
                     IF (CDABS(PSAVE)*CDABS(CTEMP) .GT. TEST) GO TO 30
   20             CONTINUE
                  NCALC = NEND + 1
   30             NCALC = NCALC - 1
                  GO TO 60
               END IF
   40       CONTINUE
            NN = NEND
C                                  Calculate special significance test
C                                  for NBMZ .GT. 2
            TEST = AMAX1(TEST,SQRT(2.0E0*CDABS(P))*SQRT(1.0E1**SIG*
     $             CDABS(PLAST)))
         END IF
C                                  Calculate P*S until significance
C                                  test is passed.
   50    NN = NN + 1
         POLD = PLAST
         PLAST = P
         P = RSIGN*(2*NN*PLAST*ZINV-POLD)
         IF (CDABS(P) .LT. TEST) GO TO 50
         IF (M .NE. 1) THEN
C                                  Calculate strict variant of
C                                  significance test and calculate
C                                  P*S until this test is passed.
            M = 1
            TEMPB = CDABS(P)/CDABS(PLAST)
            TEMPC = (NN+1)/TEMPA
            IF (TEMPB+1.0E0/TEMPB .GT. 2.0E0*TEMPC) TEMPB = TEMPC +
     $          SQRT(TEMPC**2-1.0E0)
            TEST = TEST/SQRT(TEMPB-1.0E0/TEMPB)
            IF (CDABS(P) .LT. TEST) GO TO 50
         END IF
C                                  The coefficient of B(NN) in the
C                                  normalization sum is
C                                  M*SQRT(-1)**JMAG, where M=-2,0,
C                                  or 2, and IMAG is 0 or 1. Calculate
C                                  recursion rules for M and JMAG, and
C                                  initialize them.
   60    IF (FLAG) NN = NEWNN
         NN = NN + 1
         TEMPB = DIMAG(Z)
         IPOS = 0
         IF (TEMPB.LT.0.0E0 .OR. TEMPB.GT.0.0E0) IPOS =
     $       INT(1.1E0*TEMPB/ABS(TEMPB))
         MRECUR = 4*((2+IPOS)/2) - 3 - 2*(IPOS)
         K = 2 + IPOS
         L = NN - 4*(NN/4)
         MLAST = 2 + 8*((K*L)/4) - 4*((K*L)/2)
         IF (IPOS.EQ.0 .AND. (L.EQ.1.OR.L.EQ.3)) MLAST = 0
         L = L + 3 - 4*((L+3)/4)
         M = 2 + 8*((K*L)/4) - 4*((K*L)/2)
         IF (IPOS.EQ.0 .AND. (L.EQ.1.OR.L.EQ.3)) M = 0
         IMRECR = IPOS*IPOS
         JMAG = IMRECR*(L-2*(L/2))
C                                  Initialize the backward recursion and
C                                  the normalization sum.
         CTEMPB = dcmplx(0.0E0,0.0E0)
         CTEMPA = 1.0/P
         IF (JMAG .EQ. 0) THEN
            SUM = M*CTEMPA
         ELSE
            SUM = M*dcmplx(0.0,1.0)*CTEMPA
         END IF
         NEND = NN - N
         IF (NEND .GT. 0) THEN
C                                  Recur backward via difference
C                                  equation calculation (but not
C                                  storing) CBS(NN) until NN = N.
            DO 70  L=1, NEND
               NN = NN - 1
               CTEMP = CTEMPB
               CTEMPB = CTEMPA
               P = 2*NN*ZINV
               CTEMPA = P*CTEMPB - RSIGN*CTEMP
               JMAG = (1-JMAG)*IMRECR
               K = MLAST
               MLAST = M
               M = K*MRECUR
               IF (JMAG .EQ. 0) THEN
                  SUM = SUM + M*CTEMPA
               ELSE
                  SUM = SUM + M*dcmplx(0.0E0,1.0E0)*CTEMPA
               END IF
   70       CONTINUE
         END IF
         IF (NEND .GE. 0) THEN
C                                  Store CBS(NN)
            CBS(NN) = CTEMPA
            IF (NN .LE. 1) THEN
C                                  N=1.  Since 2*CTEMPA was added to
C                                  SUM, we must subtract CTEMPA
               SUM = SUM - CTEMPA
               GO TO 90
            END IF
C                                  Calculate and store CBS(NN-1)
            NN = NN - 1
            P = 2*NN*ZINV
            CBS(NN) = P*CTEMPA - RSIGN*CTEMPB
            IF (NN .EQ. 1) THEN
               SUM = SUM + CBS(1)
               GO TO 90
            END IF
            JMAG = (1-JMAG)*IMRECR
            K = MLAST
            MLAST = M
            M = K*MRECUR
            IF (JMAG .EQ. 0) THEN
               SUM = SUM + M*CBS(NN)
            ELSE
               SUM = SUM + M*dcmplx(0.0E0,1.0E0)*CBS(NN)
            END IF
         ELSE
C                                  NN .LT. N, so store CBS(NN), and se
C                                  higher orders zero
            CBS(NN) = CTEMPA
            NEND = -NEND
            CALL CSET (NEND, dcmplx(0.0E0,0.0E0), CBS(NN+1), 1)
         END IF
         NEND = NN - 2
         IF (NEND .NE. 0) THEN
C                                  Calculate via difference equation and
C                                  store CBS(NN), UNTIL NN=2
            DO 80  L=1, NEND
               NN = NN - 1
               P = 2*NN*ZINV
               CBS(NN) = P*CBS(NN+1) - RSIGN*CBS(NN+2)
               JMAG = (1-JMAG)*IMRECR
               K = MLAST
               MLAST = M
               M = K*MRECUR
               IF (JMAG .EQ. 0) THEN
                  SUM = SUM + M*CBS(NN)
               ELSE
                  SUM = SUM + M*dcmplx(0.0E0,1.0E0)*CBS(NN)
               END IF
   80       CONTINUE
         END IF
C                                  Calculate and store CBS(1)
         CBS(1) = 2.0E0*(CBS(2)*ZINV) - RSIGN*CBS(3)
         SUM = SUM + CBS(1)
C                                  Calculate normalization factor,
C                                  CTEMPA
   90    CTEMP = (-IPOS)*dcmplx(0.0E0,1.0E0)*Z
         TEMP = EXP(REAL(CTEMP))
         CTEMPB = dcmplx(DCOS(DIMAG(CTEMP)),DSIN(DIMAG(CTEMP)))
         CTEMPA = TEMP*(CTEMPB/SUM)
C                                  Normalize
         DO 100  NN=1, N
            CBS(NN) = CBS(NN)*CTEMPA
  100    CONTINUE
      ELSE
C                                  Two-term ascending series for small Z
         CTEMPA = dcmplx(1.0E0,0.0E0)
         CTEMP = 0.25E0*Z*Z
         CBS(1) = CTEMPA - RSIGN*CTEMP
         IF (N .NE. 1) THEN
            DO 110  NN=2, N
               CTEMPA = (CTEMPA*Z)/(2*NN-2)
               CTEMPB = CTEMPA*CTEMP
               TEMP = RSIGN*REAL(CTEMPB)
               CTEMPB = dcmplx(TEMP,DIMAG(CTEMPB))
               CBS(NN) = CTEMPA - CTEMPB/NN
  110       CONTINUE
         END IF
      END IF
      IF (N .NE. NCALC) THEN
         WRITE (*,*) 'The function values CBS(I), (I = '//
     $               '1,%(I1)), are computed to machine precision, '//
     $               'but precision is lost for CBS(I), (I = '//
     $               '%(I2),%(I3)).  The input argument is Z = '//
     $               '%(C1).'
      END IF
C
 9000 RETURN
      END
C
*--------------------------------------------------------------------*
      SUBROUTINE CSET (N, CA, CX, INCX)
*--------------------------------------------------------------------*
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER    N, INCX
      COMPLEX*16    CA, CX(*)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER    I, NINCX
C
      IF (N .GT. 0) THEN
         IF (INCX .NE. 1) THEN
C                                  CODE FOR INCREMENT NOT EQUAL TO 1
            NINCX = N*INCX
            DO 10  I=1, NINCX, INCX
               CX(I) = CA
 10         CONTINUE
         ELSE
C                                  CODE FOR INCREMENT EQUAL TO 1
C                                    CLEAN-UP LOOP
            DO 20  I=1, N
               CX(I) = CA
 20         CONTINUE
         END IF
      END IF
      RETURN
      END
