C=======================================================================
      SUBROUTINE BESLPN (H,M,GN,N,FPK,AP,DT,IERR)
C ----------------------------------------------------------------------
C      BESSEL LOW PASS FILTER COEFFICIENT
C ----------------------------------------------------------------------
C **** ARGUMENTS*****
C ( O )  H      : FILTER COEFFICIENTS
C ( O )  M      : ORDER OF FILTER  (M=(N+1)/2)
C ( O )  GN     : GAIN FACTOR
C ( I )  N      : ORDER OF BESSEL FILTER FUNCTION
C ( I )  FPK    : PASS BAND FREQUENCY  (DIMENSIONAL)  ( HZ )
C ( I )  AP     : MAX. ATTENUATION IN PASS BAND
C ( I )  DT     : SAMPLING INTERVAL
C ---------------------------------------------------------------------
C      A. KATSUMATA ( 1989/5/20 )
C ---------------------------------------------------------------------
      DIMENSION H(*)
      DIMENSION COEFBS(0:50),WCOEF(0:50)
      COMPLEX   ZROOT(50)
      DATA  PI/3.141593/,HP/1.570796/
C **** ARGUMENTS  *****
      FP = FPK * DT
      WP = ABS(FP)*PI
      IF( WP.EQ.0. .OR. WP.GE.HP )  GO TO  100
C****  DETERMINE  C  ******
      TP = TAN(WP)
      PA = ABS(AP)
      IF( PA.EQ.0. )  PA = 0.5
C----- GET BESSEL FILTER POLYNOMIAL 
      CALL GETBPL(N,COEFBS(0),MAXDBS,IERROR)
      IF(IERROR.NE.0) GOTO 100
      CALL DETBTX(PA,COEFBS(0),MAXDBS,TX,IERROR)
      IF(IERROR.NE.0) GOTO 100
      C=TX/TP
      W1=1.0
      IF(COEFBS(MAXDBS).LT.0.99999.OR.COEFBS(MAXDBS).GT.1.00001)
     +             W1=1.0/COEFBS(MAXDBS)
C
      DO 20 J=0,MAXDBS-1
         WCOEF(MAXDBS-J)=COEFBS(J)*W1
   20 CONTINUE
C------ GET ROOTS OF POLYNAMIAL
      CALL BAIHIT(MAXDBS,WCOEF(1),ZROOT(1))
C ---------------------------------------------------------------------
      M  = MAXDBS/2
      G  = 1.
C ---------------------------------------------------------------------
      DO 1 J=1,M
        SR = REAL(ZROOT(2*J-1))
        SI = AIMAG(ZROOT(2*J-1))
        SA2= SR*SR+SI*SI
        A  = 1./(SA2+C*C-C*2.0*SR)
        G  = G*A*SA2
          H(J*4-3) = 2.
          H(J*4-2) = 1.
          H(J*4-1) = (2.0*SA2-2.0*C*C)*A
          H( J*4 ) = (C*C+SA2+C*2.0*SR)*A
    1 CONTINUE
C****  FOR ODD N  *****
      IF( MOD(MAXDBS,2).NE.0 ) THEN
        M = M + 1
        SR=REAL(ZROOT(MAXDBS))
        A = 1.0/(C-SR)
        G = G*A*(-SR)
         H(M*4-3) = 1.
         H(M*4-2) = 0.
         H(M*4-1) = (-C-SR)*A
         H( M*4 ) = 0.
      ENDIF
      GN = G
      IERR=0
      RETURN
C****  ERROR  ******
  100 WRITE(6,101)  FP
  101  FORMAT(/1X,5('?'),'   (BESLPN)   INVALID INPUT   FP =',
     *        1PE14.6,3X,5('?')//)
      IERR=1
      RETURN
      END
C====================================================================
      SUBROUTINE BESHPN (H,M,GN,N,FPK,AP,DT,IERR)
C ----------------------------------------------------------------------
C      BESSEL HIGH PASS FILTER COEFFICIENT
C ----------------------------------------------------------------------
C **** ARGUMENTS*****
C ( O )  H      : FILTER COEFFICIENTS
C ( O )  M      : ORDER OF FILTER  (M=(N+1)/2)
C ( O )  GN     : GAIN FACTOR
C ( I )  N      : ORDER OF BESSLE FILTER FUNCTION
C ( I )  FPK    : PASS BAND FREQUENCY  (DIMENSIONAL)  ( HZ )
C ( I )  AP     : MAX. ATTENUATION IN PASS BAND
C ( I )  DT     : SAMPLING INTERVAL (SEC)
C ----------------------------------------------------------------------
C      A. KATSUMATA (1989/5/20)
C ----------------------------------------------------------------------
      DIMENSION H(*)
      DIMENSION COEFBS(0:50),WCOEF(0:50)
      COMPLEX   ZROOT(50)
      DATA  PI/3.141593/,HP/1.570796/
C
C------ NON-DIMENSIONAL FREQUENCY
      FP = FPK       * DT
      WP = ABS(FP)*PI
      IF( WP.EQ.0..OR. WP.GE.HP)  GO TO  100
C****  DETERMINE C  *****
      TP = TAN(WP)
      PA = ABS(AP)
      IF( PA.EQ.0. )  PA = 0.5
C----- GET BESSEL FILTER POLYNOMIAL 
      CALL GETBPL(N,COEFBS(0),MAXDBS,IERROR)
      IF(IERROR.NE.0) GOTO 100
      CALL DETBTX(PA,COEFBS(0),MAXDBS,TX,IERROR)
      IF(IERROR.NE.0) GOTO 100
      C  = TX*TP
C----- GET BESSEL FILTER POLYNOMIAL ---------------------
      CALL GETBPL(N,COEFBS(0),MAXDBS,IERROR)
      W1=1.0
      IF(COEFBS(MAXDBS).LT.0.99999.OR.COEFBS(MAXDBS).GT.1.00001)
     +             W1=1.0/COEFBS(MAXDBS)
C
      DO 20 J=0,MAXDBS-1
         WCOEF(MAXDBS-J)=COEFBS(J)*W1
   20 CONTINUE
C------ GET ROOTS OF POLYNAMIAL
      CALL BAIHIT(MAXDBS,WCOEF(1),ZROOT(1))
C ---------------------------------------------------------------------
      M  = MAXDBS/2
      G  = 1.
C ---------------------------------------------------------------------
      DO 1 J=1,M
        SR = REAL(ZROOT(2*J-1))
        SI = AIMAG(ZROOT(2*J-1))
        SA2= SR*SR+SI*SI
        A  = 1./(SA2+C*C-C*2.0*SR)
        G  = G*A*SA2
          H(J*4-3) = -2.
          H(J*4-2) = 1.
          H(J*4-1) =-(2.0*SA2-2.0*C*C)*A
          H( J*4 ) = (C*C+SA2+C*2.0*SR)*A
    1 CONTINUE
      IF( MOD(MAXDBS,2).NE.0 ) THEN
C****  FOR ODD N  *****
        M = M + 1
        SR=REAL(ZROOT(MAXDBS))
        A = 1.0/(C-SR)
        G = G*A*(-SR)
         H(M*4-3) = -1.
         H(M*4-2) = 0.
         H(M*4-1) =-(-C-SR)*A
         H( M*4 ) = 0.
      ENDIF
      GN = G
      IERR=0
      RETURN
C****  ERROR  *****
  100 WRITE(6,101)  FP
  101  FORMAT(/1X,5('?'),'   (BESHPN)   INVALID INPUT   FP =',
     *        1PE14.6,3X,5('?')//)
      IERR=1
      RETURN
      END
C=======================================================================
      SUBROUTINE BESPSN (H,M,GN,N,FLK,FHK,AP,DT,IERR)
C ----------------------------------------------------------------------
C      BESSEL BAND PASS FILTER COEFFICIENT
C ----------------------------------------------------------------------
C **** ARGUMENTS*****
C ( O )  H      : FILTER COEFFICIENTS
C ( O )  M      : ORDER OF FILTER  (M=(N+1)/2)
C ( O )  GN     : GAIN FACTOR
C ( I )  N      : ORDER OF BESSEL FILTER FUNCTION
C ( I )  FLK    : LOW  FREQUENCY CUT-OFF  (DIMENSIONAL) ( HZ )
C ( I )  FHK    : HIGH FREQUENCY CUT-OFF  (DIMENSIONAL) ( HZ )
C ( I )  AP     : MAX. ATTENUATION IN PASS BAND
C ( I )  DT     : SAMPLING INTERVAL
C ----------------------------------------------------------------------
C      A. KATSUMATA (1989/5/22)
C ----------------------------------------------------------------------
      DIMENSION H(*)
      DIMENSION COEFBS(0:50),WCOEF(0:50)
      COMPLEX  ZROOT(50)
     +        ,SIO,SIC,R02,CR1,CR2
C
      DATA  PI/3.141593/,HP/1.570796/
      FL = FLK       * DT
      FH = FHK       * DT
      WL = AMIN1(ABS(FL),ABS(FH))*PI
      WH = AMAX1(ABS(FL),ABS(FH))*PI
      PA = ABS( AP )
      IF( PA.EQ.0. )  PA = 0.5
C
      IF( WL.LE.0. .OR. WL.EQ.WH .OR. WH.GE.HP ) GOTO 100
C
C----- GET BESSEL FILTER POLYNOMIAL ---------------------
      CALL GETBPL(N,COEFBS(0),MAXDBS,IERROR)
      IF(IERROR.NE.0) GOTO 100
C
      CALL DETBTX(PA,COEFBS(0),MAXDBS,TX,IERROR)
      IF(IERROR.NE.0) GOTO 100
C****  DETERMINE  C *****
      C=TX/(TAN(WH)-TAN(WL))
      WA = TAN(WL)*TAN(WH)
C------------------------------------------------------------------
      W1=1.0
      IF(COEFBS(MAXDBS).LT.0.99999.OR.COEFBS(MAXDBS).GT.1.00001)
     +             W1=1.0/COEFBS(MAXDBS)
      DO 20 J=0,MAXDBS-1
         WCOEF(MAXDBS-J)=COEFBS(J)*W1
   20 CONTINUE
C------ GET ROOTS OF POLYNAMIAL
      CALL BAIHIT(MAXDBS,WCOEF(1),ZROOT(1))
C------------------------------------------------------------------
      K  = MAXDBS/2
      M  = K*2
      L  = 0
      G  = 1.
C ---------------------------------------------------------------------
      DO  2  J=1,K
         SIO = ZROOT(J*2-1)
         SIC = ZROOT(J*2)
         ZABS = SQRT(SIO*SIC)
         R02 = CMPLX(WA,0.0)*CMPLX(C,0.0)*CMPLX(C,0.0)
         CR1 = (SIO + CSQRT( SIO*SIO -4.0*R02)) / 2.0
         CR2 = (SIO - CSQRT( SIO*SIO -4.0*R02)) / 2.0
C
         A21 = C/ZABS+2.0*REAL(CR1)/ZABS+CABS(CR1)*CABS(CR1)/C/ZABS
         A11 = -2.0*C/ZABS+2.0*CABS(CR1)*CABS(CR1)/C/ZABS
         A01 = C/ZABS-2.0*REAL(CR1)/ZABS+CABS(CR1)*CABS(CR1)/C/ZABS
C
         A22 = C/ZABS+2.0*REAL(CR2)/ZABS+CABS(CR2)*CABS(CR2)/C/ZABS
         A12 = -2.0*C/ZABS+2.0*CABS(CR2)*CABS(CR2)/C/ZABS
         A02 = C/ZABS-2.0*REAL(CR2)/ZABS+CABS(CR2)*CABS(CR2)/C/ZABS
C         
         H(J*8-7) = 0.0
         H(J*8-6) =-1.0
         H(J*8-5) = A11/A01
         H(J*8-4) = A21/A01
         H(J*8-3) = 0.0
         H(J*8-2) =-1.0
         H(J*8-1) = A12/A02
         H(J*8  ) = A22/A02
         G = G/(A01*A02)
    2 CONTINUE
C****  FOR ODD N  ******
      IF( MAXDBS.NE.M ) THEN
         M = M + 1
         WB=C/REAL(ZROOT(MAXDBS))
         WW0=-WA*WB-WB+1.0
         WW1=-2.0*(WA-1.0)*WB
         WW2=-WA*WB-WB-1.0
         G  = G / WW0
          H(M*4-3) = 0.
          H(M*4-2) =-1.
          H(M*4-1) = WW1/WW0
          H(M*4  ) = WW2/WW0
      ENDIF
      GN = G
      IERR=0
      RETURN
C****  ERROR  *****
  100 WRITE(6,101)  FL,FH
  101 FORMAT(/1X,5('?'),'   (BESPSN)   INVALID INPUT   FL =',
     *        1PE14.6,3X,'FH =',E14.6,3X,5('?')//)
      IERR=1
      RETURN
      END
C**********************************************
C
      SUBROUTINE GETBPL(NN,COEFBS,MAXBS,IERR)
C
C      GET THE BESSEL FILTER POLYNOMIAL
C
C       NN     :  ORDER OF POLYNOMIAL
C       COEFBS : COEFFICIENT OF POLYNOMIAL
C       MAXBS  : MAXIMUM ORDER OF POLYNOMIAL
C       IERR   : .NE.0 ==> ERROR
C
C**********************************************
C
      DIMENSION COEFBS(0:50)
C
C       COEFG  : COEFFICIENT OF POLYNOMIAL  ( FOR REAL PART )
C       COEFU  : COEFFICIENT OF POLYNOMIAL  ( FOR IMAGINARY PART )
      DIMENSION COEFG(0:50),COEFU(0:50)
C
C==== CLEAR OUTPUT BUFFER ========================
      DO 10 I=0,50
         COEFBS(I)=0.0
   10 CONTINUE
      CALL GETBGU(NN,COEFG(0),MAXDG,COEFU(0),MAXDU,IERROR)
      IF(IERROR.NE.0) THEN
         IERR=IERROR
         RETURN
      ENDIF
C------ REAL PART OF POLYNOMIAL
      DO 20 I=0,MAXDG
         IF(MOD(I,4).EQ.0) THEN
            COEFBS(I)=COEFBS(I)+COEFG(I)
         ELSEIF(MOD(I,4).EQ.2) THEN
            COEFBS(I)=COEFBS(I)-COEFG(I)
         ENDIF
   20 CONTINUE
C------ IMAGINARY PART OF POLYNOMIAL
      DO 30 I=0,MAXDU
         IF(MOD(I,4).EQ.1) THEN
            COEFBS(I)=COEFBS(I)+COEFU(I)
         ELSEIF(MOD(I,4).EQ.3) THEN
            COEFBS(I)=COEFBS(I)-COEFU(I)
         ENDIF
   30 CONTINUE
C
      MAXBS=MAXDG
      IF(MAXBS.LT.MAXDU) MAXBS=MAXDU
C
      RETURN
      END
C**********************************************
C
      SUBROUTINE GETBGU(NN,COEFG,MAXDG,COEFU,MAXDU,IERR)
C
C      GET BESSEL FILTER POLYNOMIAL ( SEPARATLY ABOUT G,U )
C
C       NN  :  ORDER OF POLYNOMIAL
C       COEFG : COEFFICIENT OF POLYNOMIAL OF REAL PART      ( G )
C       MAXDG : MAXIMUM ORDER OF COEFG
C       COEFU : COEFFICIENT OF POLYNOMIAL OF IMAGINARY PART ( U )
C       MAXDU : MAXIMUM ORDER OF COEFU
C       IERR  : .NE.0 ==> ERROR
C
C**********************************************
C
      DIMENSION COEFG(0:50),COEFU(0:50)
      REAL*8    DCOEFG(0:50),DCOEFU(0:50)
C
C  MRKEXG   MARK FOR NEXT EXTENSION OF POLYNOMIAL ORDER  ( REAL PART )
C  MRKEXU                                                  ( IM. PART )
C  WCOEFG    COEFFICIENT OF EACH TERM   ( REAL PART )
C  WCOEFU                               ( IMAGINARY PART )
C  IDEGG    ORDER OF EACH TERM        ( FOR REAL PART )
C  IDEGU                               ( FOR IMAGINARY PART )
C
      DIMENSION MRKEXG(1000),MRKEXU(1000),IDEGG(1000),IDEGU(1000)
      REAL*8    WCOEFG(1000),WCOEFU(1000)
C
C====== CLEAR OUTPUT BUFFER ==========================
      DO 10 I=0,50
         COEFG(I)=0.0
         COEFU(I)=0.0
         DCOEFG(I)=0.0D0
         DCOEFU(I)=0.0D0
   10 CONTINUE
C
      IERR=0
C----- NN ILLEGAL ? -------------------------
      IF(NN.LE.0.OR.NN.GT.20) THEN
         IERR=1
         WRITE(6,'(/4X,40(''#''),''   GETBGU ERROR NN='',I8,3X,
     +        40(''#''))') NN
         RETURN
C----- NN=1 ----------------------------------
      ELSEIF(NN.EQ.1) THEN
         COEFG(0)=1.0
         COEFU(1)=1.0
         MAXDG=0
         MAXDU=1
         RETURN
      ENDIF
C
C====== CLEAR WORK BUFFER ============================
      DO 20 I=1,1000
         MRKEXG(I)=0
         MRKEXU(I)=0
         WCOEFG(I)=0.0D0
         WCOEFU(I)=0.0D0
         IDEGG(I)=0
         IDEGU(I)=0
   20 CONTINUE
      MRKEXG(1)=1
      WCOEFG(1)=1.0D0
      IDEGG(1)=0
      WCOEFU(1)=1.0D0
      IDEGU(1)=1
C..... NUMBER OF TERM
      NTERMG=1
      NTERMU=1
C
C==========================================
C  CALCULATE EACH TERM COEFFICIENT
C==========================================
      DO 30 I=2,NN
C
C===== REAL PART =============================
         NG=NTERMG
         DO 40 J=1,NTERMG
C------ NO EXTENSION TERM --------------------------
            IF(MRKEXG(J).EQ.0) THEN
               MRKEXG(J)=1
               WCOEFG(J)=WCOEFG(J)*DBLE(2*I-1)
C------ EXTENSION TERM
            ELSE
C     ( CREATE NEW TERM )
               NG=NG+1
               WCOEFG(NG)=-WCOEFG(J)/DBLE(2*I-3)
               IDEGG(NG)=IDEGG(J)+2
               MRKEXG(NG)=0
C
               WCOEFG(J)=WCOEFG(J)*DBLE(2*I-1)
               MRKEXG(J)=1
            ENDIF
   40    CONTINUE
         NTERMG=NG
C
C======= IMAGINARY PART =============================
         NU=NTERMU
         DO 50 J=1,NTERMU
C------ NO EXTENSION TERM --------------------------
            IF(MRKEXU(J).EQ.0) THEN
               MRKEXU(J)=1
               WCOEFU(J)=WCOEFU(J)*DBLE(2*I-1)
C------ EXTENSION TERM
            ELSE
C     ( CREATE NEW TERM )
               NU=NU+1
               WCOEFU(NU)=-WCOEFU(J)/DBLE(2*I-3)
               IDEGU(NU)=IDEGU(J)+2
               MRKEXU(NU)=0
C
               WCOEFU(J)=WCOEFU(J)*DBLE(2*I-1)
               MRKEXU(J)=1
            ENDIF
   50    CONTINUE
         NTERMU=NU
C
   30 CONTINUE
C
C================================================
C    GROUPING EACH ORDER
C================================================
C----- REAL PART --------------------------
      MAXDG=0
      DO 60 I=1,NTERMG
         IF(IDEGG(I).GT.MAXDG) MAXDG=IDEGG(I)
   60    DCOEFG(IDEGG(I))=DCOEFG(IDEGG(I))+WCOEFG(I)
      DO 65 I=0,MAXDG
   65    COEFG(I)=DCOEFG(I)
C----- IMAGINARY PART -------------------------
      MAXDU=0
      DO 70 I=1,NTERMU
         IF(IDEGU(I).GT.MAXDU) MAXDU=IDEGU(I)
   70    DCOEFU(IDEGU(I))=DCOEFU(IDEGU(I))+WCOEFU(I)
C--- REAL*8 ==> REAL*4
      DO 75 I=0,MAXDU
   75    COEFU(I)=DCOEFU(I)
C
      RETURN
      END
C***************************************************
C
      SUBROUTINE DETBTX(AP,COEFFB,MAXDBS,TX,IERR)
C
C      DETERMINE BESSEL FILTER POLYNOMIAL X VALUE FOR GIVEN AP
C
C ( I ) AP   : PASS BAND MAXIMUM ATTENATION
C ( I ) COEFFB : COEFFICIENTS OF POLYNOMIAL
C ( I ) MAXDBS : MAX. ORDER OF THE POLYNOMIAL
C ( O ) TX   : THE VALUE WHICH GIVES AP IN BESSEL POLYNOMIAL
C ( O ) IERR : ERROR FLAG
C
C***************************************************
C
      DIMENSION COEFBS(0:50),COEFFB(0:50)
C
C----- ATTENATION VALUE
      BP=SQRT(1.0+AP*AP)
C
C----- NORMALIZING OF POLYNOMIAL
         WW=COEFFB(0)
         DO 20 J=0,MAXDBS
            COEFBS(J)=COEFFB(J)/WW
   20    CONTINUE
C----- GET TX BY THE BISECTION METHOD
         TX1=0.0
         TX2=1.0
         DO 30 J=1,50
            CALL VALBPL(TX2,COEFBS(0),MAXDBS,VAL)
            IF(VAL-BP.GT.0.0) GOTO 40
            TX2=TX2*2.0
   30    CONTINUE
   40    CONTINUE
C
         DO 50 J=1,150
            TX=(TX1+TX2)/2.0
            CALL VALBPL(TX,COEFBS(0),MAXDBS,VAL)
            IF(VAL-BP.GT.0.0) THEN
               TX2=TX
            ELSE
               TX1=TX
            ENDIF
            IF(ABS(TX2-TX1)/ABS(TX2+TX1).LT.1.0E-5) GOTO 60
   50    CONTINUE
         WRITE(6,'('' GETBSN ERROR  NO CONVERGENCE OF X'')')
         GOTO 100
C
   60    CONTINUE
C
      IERR=0
      RETURN
C---- ERROR
  100 CONTINUE
      IERR=1
      RETURN
      END
C**************************************************
C
      SUBROUTINE VALBPL(X,COEFF,N,VAL)
C
C       GET OUTPUT VALUE FOR BESSEL FILTER POLYNOMIAL FUNCTION
C
C         X     : GIVEN VALUE FOR CALCURATION
C         COEFF : COEFFICIENTS OF POLYNOMIAL
C         N     : THE ORDER OF THE POLYNOMIAL
C         VAL   : OUTPUT
C                VAL=ABS( SIGMA (COEFF(J)*(IX)**J) )
C
C***************************************************
C
      DIMENSION COEFF(0:N)
C
      REAL*8 VR,VI,XN
C----- INITIALIZE
      VR=0.0D0
      VI=0.0D0
      XN=1.0D0
C
      DO 10 I=0,N
C.... REAL PART
         IF(MOD(I,4).EQ.0) THEN
            VR=VR+DBLE(COEFF(I))*XN
         ELSEIF(MOD(I,4).EQ.2) THEN
            VR=VR-DBLE(COEFF(I))*XN
C,.... IMAGINARY PART
         ELSEIF(MOD(I,4).EQ.1) THEN
            VI=VI+DBLE(COEFF(I))*XN
         ELSEIF(MOD(I,4).EQ.3) THEN
            VI=VI-DBLE(COEFF(I))*XN
         ELSE
            WRITE(6,'('' VALBPL  ERROR I='',I5)') I
         ENDIF
C...... XN=X**(I)
         XN=XN*DBLE(X)
   10 CONTINUE
C
      VAL=DSQRT(VR*VR+VI*VI)
C
      RETURN
      END
C************************************************
C
      SUBROUTINE BAIHIT(N,A,Z)
C
C      CALCULATE ROOTS OF POLYNOMIAL BY BAIRSTOW-HITCHICOCK METHOD
C        ( 1989/5/12 PROGRAMMED BY KATSUMATA )
C       ( AFTER 'SUUCHI-KENSAN-HOU' BY HAYATO TOGAWA
C          CORONA PUBLISHING CO., LTD.   1981 )
C
C          N : ORDER OF POLYNOMIAL
C          A(I) : COEFFIENTS OF POLYNOMIAL
C                   A(I) COEFFICIENT FOR N-I POWERS OF X
C                   COEFFIENT OF N POWERS OF X IS 1.0
C          Z(I) : COMPLEX ROOTS (I : 1 -> N )
C
C*************************************************
C
      DIMENSION A(N)
      COMPLEX Z(N)
C P,Q : FACTOR POLYNOMIAL
C            PBUR(I):  COEFFICIENT OF 1 POWER OF X
C            QBUF(I):  COEFFICIENT OF 0 POWER OF X
C            I: INDEX FOR POLYNOMIAL
      REAL*8    PBUF(20),QBUF(20)
      REAL*8    AWRK(50),B(50),BWRK(50),C(50)
     +         ,P,Q,BN,BNM1
     +         ,CNM1,CNM2,CNM3,A11,A12,A21,A22,B1,B2,DP,DQ,H
     +         ,RA,SA,RB,SB
C
C------- PRE-CHECK FOR POLYNIMIAL DIVIDE
      IF(N.EQ.1) THEN
         Z(1)=CMPLX(-A(1),0.0)
         RETURN
      ELSEIF(N.EQ.2) THEN
         NPQ=1
         PBUF(NPQ)=A(1)
         QBUF(NPQ)=A(2)
         GOTO 20
      ELSEIF(N.LT.1) THEN
         WRITE(6,'('' BAIHIT  ORDER ERROR N='',I3)') N
         RETURN
      ENDIF
C
C------ NUMBER OF POLYNOMIALS WITH 2 POWERS
      NPQ=0
C------ ORDER OF REMAINED POLYNOMIAL
      NA=N
C
      DO 40 I=1,N
         AWRK(I)=A(I)
   40 CONTINUE
C
   10 CONTINUE
C======================================================
C     GET P & Q ( POLYNOMIAL COEFFICIENT OF 2 POWERS )
C             BY ITERATION METHOD
C======================================================
   30 CONTINUE
C------ INITIAL SETTING OF P & Q
         IF(ABS(AWRK(NA)).GT.1.0E-2) THEN
            P=(ABS(AWRK(NA)))**(1.0/REAL(NA))
            Q=P
            RJUGDE=P
         ELSE
            P=1.0D0
            Q=1.0D0
            RJUGDE=P
         ENDIF
C
         DLIMIT=1.0E-12 
  100    CONTINUE
         DO 50 KKK=1,50
C------- DIVIDE BY (X**2+P*X+Q)
            CALL DDVPLN(AWRK(1),NA,B(1),P,Q,RA,SA)
C....... END?
            IF(ABS(RA)/RJUGDE.LT.DLIMIT.AND.ABS(SA)/RJUGDE.LT.DLIMIT)
     +                        THEN
               NA=NA-2
               NPQ=NPQ+1
               PBUF(NPQ)=P
               QBUF(NPQ)=Q
C
               IF(NA.GE.3) THEN
                  DO 80 I=1,NA
                     AWRK(I)=B(I)
   80             CONTINUE
                  GOTO 30
               ELSEIF(NA.EQ.2) THEN
                  NPQ=NPQ+1
                  PBUF(NPQ)=B(1)
                  QBUF(NPQ)=B(2)
                  GOTO 20
               ELSEIF(NA.EQ.1) THEN
                  Z(N)=DCMPLX(-B(1),0.0D0)
                  GOTO 20
               ENDIF
            ENDIF
C====== GET NEXT P & Q VALUE ============================
            NB=NA-2
            DO 70 I=1,NB
   70          BWRK(I)=B(I)
            BWRK(NB+1)=0.0
C
            CALL DDVPLN(BWRK(1),NB+1,C(1),P,Q,RB,SB)
C
            BN=B(NB)
            IF(NB.EQ.1) THEN
               BNM1=1.0
            ELSE
               BNM1=B(NB-1)
            ENDIF
C
            IF(NB.EQ.1) THEN
               CNM1=1.0
               CNM2=0.0
               CNM3=0.0
            ELSEIF(NB.EQ.2) THEN
               CNM1=C(NB-1)
               CNM2=1.0
               CNM3=0.0
            ELSEIF(NB.EQ.3) THEN
               CNM1=C(NB-1)
               CNM2=C(NB-2)
               CNM3=1.0
            ELSE
               CNM1=C(NB-1)
               CNM2=C(NB-2)
               CNM3=C(NB-3)
            ENDIF
C
            A11=BN-CNM2*Q-CNM1*P
            A12=BNM1-CNM3*Q-CNM2*P
            A21=-CNM1*Q
            A22=BN-CNM2*Q
            B1=RA
            B2=SA
C
            DP=(B1*A22-A12*B2)/(A11*A22-A12*A21)
            DQ=(A11*B2-A21*B1)/(A11*A22-A12*A21)
            P=P+DP
            Q=Q+DQ
C------- GOTO NEXT ITERATION
   50    CONTINUE
         DLIMIT=DLIMIT*2.0
         GOTO 100
C
C=======================================================
C     CALCULATE ROOTS
C=======================================================
   20 CONTINUE
      DO 60 I=1,NPQ
         H=PBUF(I)*PBUF(I)-4.0*QBUF(I)
         IF(H.GE.0.0) THEN
            Z(2*I-1)=DCMPLX((-PBUF(I)+DSQRT(H))/2.0,0.0D0)
            Z(2*I)  =DCMPLX((-PBUF(I)-DSQRT(H))/2.0,0.0D0)
         ELSE
            Z(2*I-1)=DCMPLX(-PBUF(I)/2.0, DSQRT(-H)/2.0D0)
            Z(2*I)  =DCMPLX(-PBUF(I)/2.0,-DSQRT(-H)/2.0D0)
         ENDIF
   60 CONTINUE
C
      RETURN
      END
C************************************************
C
      SUBROUTINE DDVPLN(A,N,B,P,Q,R,S)
C
C     DIVIDE POLYNOMIAL BY (X*X+PX+Q)  
C      (A*X) = (X*X+P*X+Q)*(B*X) + (RX+S)
C       ( A*X = X**N + A1*X**(N-1) + ...... + AN )
C
C************************************************
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(N),B(N)
C
      IF(N.EQ.1) THEN
         R=1.0
         S=A(1)
      ELSEIF(N.EQ.2) THEN
         R=A(1)-P
         S=A(2)-Q
      ELSE
         DO 10 I=1,N-2
            IF(I.EQ.1) THEN
               B(I)=A(1)-P
            ELSEIF(I.EQ.2) THEN
               B(I)=A(2)-P*B(1)-Q
            ELSE
               B(I)=A(I)-P*B(I-1)-Q*B(I-2)
            ENDIF
   10    CONTINUE
C------ AMARI -----------------------------
         IF(N.EQ.3) THEN
            R=A(2)-Q-P*B(1)
            S=A(3)-Q*B(1)
         ELSE
            R=A(N-1)-Q*B(N-3)-P*B(N-2)
            S=A(N)-Q*B(N-2)
         ENDIF
      ENDIF
C
      RETURN
      END
