Make your own free website on Tripod.com
Program for the Design of Recursive Digital Filters

Return to Main Software Page

C
C-----------------------------------------------------------------------
C MAIN PROGRAM:         D O R E D I
C
C               DESIGN AND OPTIMIZATION OF RECURSIVE DIGITAL FILTERS
C
C     BUTTERWORTH, CHEBYSHEV, AND ELLIPTIC APPROXIMATION;
C     OPTIMIZATION OF THE COEFFICIENT WORDLENGTH AND THE OUTPUT NOISE;
C     ANALYSIS OF THE DESIGNED FILTER
C
C VERSION:      V005-3 4-15-81
C
C AUTHOR:       DR.-ING. GUENTER F. DEHNER
C
C               INSTITUT FUER NACHRICHTENTECHNIK
C               UNIVERSITAET ERLANGEN-NUERNBERG
C
C               CAUERSTRASSE 7
C               D-8520 ERLANGEN, GERMANY
C
C INPUT:        MANUAL, SECTION 2
C-----------------------------------------------------------------------
C
      DOUBLE PRECISION DPI, DOMI
C
      COMPLEX ZP, ZPS
C
      COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /CONST2/ DPI, DOMI
      COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
C
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
     *    SC0(16)
      COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
      COMMON /CGRID/ GR(64), NGR(12)
      COMMON /CRECO/ JRCO, JECO, JDCO, JJDCO(5), JMAXV, JTRB2, LREF
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /COPTST/ LOPTS, ISTOR
      COMMON /CRENO/ LCNO, ICNO(16,5)
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /CREST2/ KSEQ(16,2)
      COMMON /CNOISE/ RI, RIN, RE, REN, FAC
C
      COMMON /SCRAT/ ADUM(32)
      COMMON /CPOL/ ZP(16,2), ZPS(16,2)
      COMMON /COPST2/ PN(100,2), TF(100,2), TFA(100,2)
      COMMON /CFFUNC/ PHI(5), BF2(5), BF1(5), BF0(5), IBB(5), ICOR
      COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
      COMMON /CPOW/ PNU, PNC, AND, ITCORP
      COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
      CALL DEFIN0
C
      NDDEL = -1
      CALL HEAD(KA2)
C
      IPRUN = 0
      IPCON = 0
      IPLOT = 0
  10  CALL DORINP
C
      IF (NDOUT.NE.(-1) .OR. NDDEL.NE.(-1)) GO TO 20
      CALL HEAD(KA4)
      NDDEL = 0
  20  IWLRR = IWLR
      GO TO (30, 40, 30, 80), IPRUN
C
  30  NGR(8) = 0
      CALL SYNTHE
      NINP = 1
      IF (IWLRR.EQ.100) IWLR = IWL
      IF (IPRUN.LT.0) GO TO 50
      IF (IPRUN.NE.3) GO TO 90
C
  40  IF (NINP.EQ.0) GO TO 100
      IF (NINP.GE.3) NGR(8) = 0
C
      CALL SEQSTA
C
      GO TO 60
  50  LOPTS = 0
      ISCAL = 0
      LSEQ = 0
C
  60  CALL SEQOPT
C
      IF (IPRUN.EQ.(-3)) GO TO 90
      IF (IPRUN.EQ.3) GO TO 70
      IF (IWLR.LT.100) CALL SEQSTA
      GO TO 90
C
C JSTRU .EQ. 1       SECOND COEFFICIENT OPTIMIZATION
C
  70  IF (JSTRU.EQ.2) GO TO 90
      IPRUN = -3
      GO TO 30
  80  IF (NINP.EQ.0) GO TO 100
      IF (NINP.GE.3) NGR(8) = 0
C
      CALL ANALYS
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C HERE INSERT THE CALLS TO YOUR OWN DESIGN AND ANALYSIS PROCEDURES
C DEFINE NEW LABELS AND EXTEND THE STATEMENT '100+1'
C IPRUN HAS TO BE EQUIVALENT TO THE DEFINITIONS IN 'DORINP'
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
  90  IF (LSPOUT.EQ.(-1)) CALL OUTSPE
      IF (NINP.NE.2) GO TO 10
      ENDFILE KA4
      GO TO 10
C
 100  CALL ERROR(26)
C
C THE PROGRAM WILL COME TO THE END IN THE SUBROUTINE DORINP
C OR IN THE SUBROUTINE ERROR
C
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DEFIN0
C DEFINITION OF THE MACHINE PARAMETERS
C-----------------------------------------------------------------------
C
      SUBROUTINE DEFIN0
C
      DOUBLE PRECISION DPI, DOMI, D1MACH
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /CONST2/ DPI, DOMI
      COMMON /CIOFOR/ IOFO(3), ION
C
      DIMENSION JOFO(3)
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      DATA JOFO /2H(4,2H0A,2H2)/, JON /40/
C
C THIS STATEMENT IS WRITTEN FOR THE PDP11/45 AND IS
C MACHINE DEPENDENT (SEE MANUAL 4.2.2.)
C E.G.:   CDC CYBER 172
C DATA IOFO /10H(8A10)    ,2*10H          /,ION/8/
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
      KA1 = I1MACH(1)
      KA2 = I1MACH(2)
      KA3 = 6
      KA4 = 1
      KA5 = 9
      LINE = 65
C
      PI = 4.*ATAN(1.0)
      FLMA = 2.**(I1MACH(13)-2)
      FLMI = 2.0*R1MACH(4)
      FLER = 1.0E02*FLMI
C
      DPI = 4.00*DATAN(1.0D00)
      DOMI = 2.0D0*D1MACH(4)
C
      MAXDEG = 32
      IWLMI = 2
      IWLMA = I1MACH(11)
      MBL = 16
C
      IOFO(1) = JOFO(1)
      IOFO(2) = JOFO(2)
      IOFO(3) = JOFO(3)
      ION = JON
C
      RETURN
C
      END
C
C ======================================================================
C
C DOREDI - SUBROUTINES:  ROOT ELEMENTS FOR OVERLAYS
C
C ======================================================================
C
C
C-----------------------------------------------------------------------
C SUBROUTINE:   SYNTHE
C SUBROUTINE TO HANDLE THE NECESSARY COMMONS FOR SYNNOR AND SYNOPT
C FOR OVERLAYS
C-----------------------------------------------------------------------
C
      SUBROUTINE SYNTHE
C
      DOUBLE PRECISION DK, DKS, DCAP02, DCAP04
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
      COMMON /TOLNOR/ VSNN, NDEGN, NBN
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
      COMMON /RESIN2/ DK, DKS, DCAP02, DCAP04
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
      CALL SYNNOR
      IF (LOPTW.EQ.0) RETURN
      CALL SYNOPT
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   SYNNOR
C SYNTHESIS OF RECURSIVE DIGITAL FILTERS
C CALCULATION OF THE ZEROS AND THE POLES
C-----------------------------------------------------------------------
C
      SUBROUTINE SYNNOR
C
      COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
C
      MOUT = NOUT
      IF (IPRUN.EQ.(-3)) NOUT = 0
      CALL DESIA
      NOUT = MOUT
      IF (LOPTW.NE.0) RETURN
      CALL DESIB
      IF (NDOUT.EQ.(-1)) CALL OUT012
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   SYNOPT
C SYNTHESIS OF RECURSIVE DIGITAL FILTERS WITH OPTIMIZED COEFFICIENTS
C CALCULATION OF THE POLES AND THE OPTIMIZATION PARAMETER
C-----------------------------------------------------------------------
C
      SUBROUTINE SYNOPT
C
      COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
      COMMON /COPTCP/ IS1, IS2, IABO, IWLGS, IWLGN, IWLGP, ADEPSS,
     *    ADEPSN, ADEPSP, ACXS, ACXN, ACXP
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
     *    SC0(16)
      COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
      FLAB = 10.*FLMI
      ACXN = ACXMI
      ACXP = ACXMA
      ADEPSS = FLMA
      ADEPSP = FLMA
      ADEPSN = -FLMA
      ITER = 0
      IWLGS = IWLMA
      IWLGP = IWLMA
      IWLGN = IWLMA
      IWLL = IWLMA
      IDEPSL = 0
      IS1 = 0
      IS2 = 0
      IABO = 0
      ITERM2 = ITERM
      IF (IPRUN.EQ.3 .AND. JSTRU.EQ.1) ITERM = ITERM1
      MOUT = NOUT
      IF (NOUT.LE.3) NOUT = NOUT - 2
C
  10  ITER = ITER + 1
      IF (IPRUN.EQ.(-3)) CALL BLNUMZ
      IF (NOUT.GE.2) CALL OUT023
      CALL DESIB
      IF (ITER.GT.1) GO TO 20
      CALL GRID01
      CALL GRID04
      IF (NOUT.EQ.5) CALL OUT042
      ACXS = ACX
  20  IF (LOPTW.EQ.0) GO TO 60
      CALL COEFW
      CALL COPY02
      CALL OPTPAR
      IF (NOUT.GE.2) CALL OUT045
      IF (ITER.GE.ITERM) GO TO 30
      IF (ACX-ACXMI.LT.FLAB) GO TO 40
      IF (ACXMA-ACX.LT.FLAB) GO TO 40
      IF (IABO.LT.3) GO TO 10
      IABO = 1
      GO TO 50
  30  IABO = 2
      GO TO 50
  40  IABO = 3
  50  NOUT = MOUT
      IF (NOUT.GE.1) CALL OUT022
      ACX = ACXS
      IF (IPRUN.EQ.3 .AND. JSTRU.EQ.1) GO TO 80
      IF (IPRUN.EQ.(-3)) CALL BLNUMZ
      CALL DESIB
      IF (NOUT.EQ.1) CALL OUT011
      IF (NOUT.EQ.2) CALL OUT038
  60  IF (NOUT/2.EQ.1) CALL OUT039
      IWL = IWLGS
      IF (IPRUN.EQ.3 .AND. JSTRUS.EQ.1) GO TO 70
      IF (NDOUT.EQ.(-1)) CALL OUT012
  70  LOPTW = 1
      MOUT = NOUT
      IF (NOUT.NE.0) NOUT = NOUT + 2
      CALL COEFW
      IF (IPRUN.NE.(-3)) CALL COPY02
      NOUT = MOUT
  80  ITERM = ITERM2
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   SEQSTA
C START NOISE ANALYSIS AND OPTIMIZATION
C-----------------------------------------------------------------------
C
      SUBROUTINE SEQSTA
C
      COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /COPTST/ LOPTS, ISTOR
C
      IF (LSEQ.NE.(-1)) GO TO 10
      IF (NOUT.GE.3) CALL OUT011
      CALL COPY01
      CALL ALLO02
      GO TO 20
C
  10  IF (NINP.GE.3 .AND. LOPTS.EQ.1) CALL FIXPAR
  20  L = LSTAB
      LSTAB = 0
      CALL STABLE
      IF (ISTAB.NE.(-1)) GO TO 50
      LSTAB = L
      IF (ISTRU.GE.30) GO TO 30
      CALL TSTR01
      GO TO 40
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C HERE INSERT TESTS FOR YOUR OWN STRUCTURE IMPLEMENTATIONS
C
  30  CALL ERROR(8)
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
  40  IF (IWLR.EQ.100) RETURN
      IF (IWLR.GT.IWLMA) RETURN
      L = IWL
      IWL = IWLR
      CALL RECO
      CALL ROUND(1, 5)
      CALL STABLE
      IWL = L
      RETURN
  50  CALL ERROR(30)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   SEQOPT
C OPTIMIZATION OF THE PAIRING AND ORDERING OF THE SECOND-ORDER BLOCKS
C OPTIMIZATION PROCEDURE WITH LIMITED STORAGE
C-----------------------------------------------------------------------
C
      SUBROUTINE SEQOPT
C
      COMPLEX ZP, ZPS
      COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /COPTST/ LOPTS, ISTOR
      COMMON /CGRID/ GR(64), NGR(12)
      COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
     *    SC0(16)
      COMMON /CPOL/ ZP(16,2), ZPS(16,2)
      COMMON /CNOISE/ RI, RIN, RE, REN, FAC
      COMMON /COPST2/ PN(100,2), TF(100,2), TFA(100,2)
      COMMON /CFFUNC/ PHI(5), BF2(5), BF1(5), BF0(5), IBB(5), ICOR
      COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
      COMMON /CPOW/ PNU, PNC, AND, ITCORP
C
      IF (ISCAL.EQ.0) GO TO 10
      CALL DESCAL
  10  CALL COPY01
      IF (NOUT.GE.3) CALL OUT011
      CALL POLLOC
      IF (NGR(8).EQ.0) CALL GRID02
      IF (NOUT.GE.5) CALL OUT042
      IF (NOUT.GE.3) CALL OUT025
      IF (LOPTS.NE.0) CALL OPTLST
      CALL ANANOI
      IF (NOUT.EQ.0) RETURN
      CALL OUT027
      IF (NOUT.EQ.1) RETURN
      IF (JSTRUS.EQ.2) CALL OUT010
      CALL OUT011
      IF (JSTRUS.EQ.2) CALL DENORM
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   ANALYS
C ANALYSIS OF THE DESIGNED RECURSIVE DIGITAL FILTER
C-----------------------------------------------------------------------
C
      SUBROUTINE ANALYS
C
      COMPLEX ZP, ZPS
      COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
      COMMON /COPTST/ LOPTS, ISTOR
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /CPOL/ ZP(16,2), ZPS(16,2)
      COMMON /CNOISE/ RI, RIN, RE, REN, FAC
C
      CALL ANAL01
      IF (LOPTS.EQ.0) GO TO 10
      CALL ANAL02
  10  IF (IPLOT.EQ.0) GO TO 20
      CALL PLOT00
C
  20  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   ANAL01
C FIRST PART OF ANALYSIS
C-----------------------------------------------------------------------
C
      SUBROUTINE ANAL01
C
      COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
      COMMON /COPTST/ LOPTS, ISTOR
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /CGRID/ GR(64), NGR(12)
      COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
      IF (LSEQ.NE.(-1)) GO TO 10
      IF (NOUT.GE.3) CALL OUT011
      CALL COPY01
      CALL ALLO02
C
  10  L = LSTAB
      LSTAB = 0
      CALL STABLE
      IF (ISTAB.NE.(-1)) GO TO 50
      LSTAB = L
      IF (LOPTW.EQ.0) GO TO 30
      IF (NGR(1).NE.0) GO TO 20
      IF (NINP.GE.3) CALL GRID03
      CALL GRID01
      IF (NINP.LE.2) CALL GRID04
  20  IF (NOUT.EQ.5) CALL OUT042
      CALL COEFW
      CALL COPY02
      LOPTW = 0
  30  IF (LOPTS.NE.(-1) .AND. IPLOT.EQ.0) GO TO 60
      IF (IWLR.EQ.100) GO TO 60
      IWL = IWLR
      IF (IWL.GT.IWLMA) GO TO 40
      CALL RECO
      CALL ROUND(1, 5)
      CALL STABLE
      GO TO 60
  40  CALL ERROR(16)
      GO TO 60
  50  CALL ERROR(30)
  60  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   ANAL02
C SECOND PART OF THE ANALYSIS
C-----------------------------------------------------------------------
C
      SUBROUTINE ANAL02
C
      COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
      COMMON /COPTST/ LOPTS, ISTOR
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /CNOISE/ RI, RIN, RE, REN, FAC
      COMMON /COPST2/ PN(100,2), TF(100,2), TFA(100,2)
      COMMON /CFFUNC/ PHI(5), BF2(5), BF1(5), BF0(5), IBB(5), ICOR
      COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
      COMMON /CPOW/ PNU, PNC, AND, ITCORP
C
      IF (ISTRU.GE.30) GO TO 10
      CALL TSTR01
      GO TO 20
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C HERE INSERT TESTS FOR YOUR OWN STRUCTURE IMPLEMENTATIONS
C
  10  CALL ERROR(8)
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
  20  IF (ISCAL.NE.0) CALL DESCAL
      CALL COPY01
      IF (NOUT.GE.3) CALL OUT011
      LOPTS = 0
      CALL POLLOC
      CALL GRID02
      IF (NOUT.GE.3) CALL OUT025
      CALL ANANOI
      CALL OUT027
      IF (JSTRUS.EQ.2) CALL OUT010
      CALL OUT011
      IF (JSTRUS.EQ.2) CALL DENORM
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   PLOT00
C PRINTER PLOTS
C-----------------------------------------------------------------------
C
      SUBROUTINE PLOT00
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CSTATV/ X1(16), X2(16)
      COMMON /CLINE/ Y, X, IZ(111)
C
      DIMENSION A(6)
      EQUIVALENCE (A(1),IZ(1))
C
      DATA ICHI, ICHP, ICHC, ICHM, ICHB /1HI,1H+,1H.,1H-,1H /
C
C NUMBER OF PRINTED LINES
C
      IF (IPLOT.GT.4) GO TO 230
      ML = IPAG*LINE - 6
      ADOM = (OMUP-OMLO)/FLOAT(ML-1)
C
C SCALING
C
      IF (LNOR.NE.(-1)) GO TO 60
C
      RRMAX = -FLMA
      RRMIN = FLMA
C
      IF (IPLOT.EQ.4) Y = RESP(1.,1)
      DO 30 L=1,ML
        IF (IPLOT.EQ.4) GO TO 20
        X = ADOM*FLOAT(L-1) + OMLO
        IF (IPLOT.EQ.3) GO TO 10
        Y = AMAGO(X)
        IF (IPLOT.EQ.1) GO TO 20
        IF (Y.LT.FLMI) Y = FLMI
        Y = -20.*ALOG10(Y)
        IF (Y.GT.RMAX) Y = RMAX
        GO TO 20
  10    Y = PHASE(X)
  20    RRMAX = AMAX1(RRMAX,Y)
        RRMIN = AMIN1(RRMIN,Y)
        IF (IPLOT.NE.4) GO TO 30
        Y = RESP(0.,0)
  30  CONTINUE
C
C SCALING OF THE PLOTS
C
      RR = ABS(RRMAX-RRMIN)
      RN = 0.005*RR
      RRMAX = RRMAX - RN
      RRMIN = RRMIN + RN
      RR = ABS(RRMAX-RRMIN)
      IF (RR.LE.0.) RR = ABS(RRMAX)
      IQ = INT(ALOG10(RR))
      RN = 10.**FLOAT(IQ)
      IF (RN.LT.RR) RN = RN*10.
      AN = RR/RN
      IF (AN.GT.0.0) AAN = 0.1
      IF (AN.GT.0.1) AAN = 0.2
      IF (AN.GT.0.2) AAN = 0.5
      IF (AN.GT.0.5) AAN = 1.0
      RR = AAN*RN
C
  40  RN = RR*0.2
      RMA = AINT(RRMAX/RN)*RN
      IF (RMA.LT.RRMAX) RMA = RMA + RN
      RMI = RMA - RR
      IF (RMI.LE.RRMIN) GO TO 50
      RMI = AINT(RRMIN/RN)*RN
      IF (RMI.GT.RRMIN) RMI = RMI - RN
      RMA = RMI + RR
      IF (RMA.GE.RRMAX) GO TO 50
      RR = RR*2.
      IF (AAN.EQ.0.2) RR = RR*1.25
      GO TO 40
C
  50  RRMAX = RMA
      RRMIN = RMI
      GO TO 110
C
  60  GO TO (70, 70, 80, 90), IPLOT
  70  RRMAX = RMAX
      RRMIN = 0.
      GO TO 100
  80  RRMAX = PI
      RRMIN = -RRMAX
      GO TO 100
  90  RRMAX = RMAX
      RRMIN = RMIN
 100  RR = RRMAX - RRMIN
C
 110  CALL OUT060
      NULL = INT(-RRMIN*100./RR+1.5)
      IF (NULL.LT.1 .OR. NULL.GT.111) NULL = 0
      Q = RR/5.
      DO 120 I=1,6
        A(I) = Q*FLOAT(I-1) + RRMIN
 120  CONTINUE
      CALL OUT061
      DO 130 I=1,111
        IZ(I) = ICHM
 130  CONTINUE
      IF (IPLOT.EQ.4) Y = RESP(1.,1)
      X = 1.
      DO 220 L=1,ML
        IF (IPLOT.EQ.4) GO TO 150
        X = ADOM*FLOAT(L-1) + OMLO
        IF (IPLOT.EQ.3) GO TO 140
        Y = AMAGO(X)
        IF (IPLOT.NE.2) GO TO 150
        IF (Y.LT.FLMI) Y = FLMI
        Y = -20.*ALOG10(Y)
        GO TO 150
 140    Y = PHASE(X)
 150    DO 160 I=1,101,10
          IZ(I) = ICHC
 160    CONTINUE
        IF (NULL.NE.0) IZ(NULL) = ICHC
        Y1 = (Y-RRMIN)/RR
        IY = INT(100.*Y1+1.5)
        IF (IY.GT.0) GO TO 170
        IZ(1) = ICHI
        GO TO 190
 170    IF (IY.LE.111) GO TO 180
        IZ(111) = ICHI
        GO TO 190
 180    IZ(IY) = ICHP
 190    CALL OUT062
        IF (IPLOT.NE.4) GO TO 200
        Y = RESP(0.,0)
        X = FLOAT(L)
 200    DO 210 I=1,111
          IZ(I) = ICHB
 210    CONTINUE
 220  CONTINUE
      GO TO 240
 230  CALL ERROR(8)
 240  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DORINP
C INPUT DATA COMPILER FOR THE DESIGN PROGRAM 'DOREDI'
C-----------------------------------------------------------------------
C
      SUBROUTINE DORINP
C
      COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
      COMMON /CARD/ LCODE, IZ(40)
      COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
      COMMON /INTDAT/ JINP(16)
      COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /COPTST/ LOPTS, ISTOR
C
      DIMENSION IPT1(5), IPT2(5)
      DATA IPT1(1), IPT1(2), IPT1(3), IPT1(4), IPT1(5) /1HS,1HS,1HS,1HA,
     *    1HF/
      DATA IPT2(1), IPT2(2), IPT2(3), IPT2(4), IPT2(5) /1HY,1HE,1HS,1HN,
     *    1HI/
      DATA IPMAX /5/
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C HERE INSERT YOUR OWN PROGRAM DEFINITIONS
C CHANGE DIMENSION AND DATA STATEMENT IN FRONT
C 'FI' HAS TO BE LAST DEFINITION
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
      LCODE = -1
      IF (IPLOT.NE.0) GO TO 50
      IF (IPCON.NE.0) GO TO 10
      CALL DEFIN1
      CALL DEFIN2
      CALL DEFIN3
      CALL DEFIN4
      CALL DEFIN5
      CALL DEFIN6
      GO TO 60
C
  10  IF (IPCON.EQ.99) STOP
      IPRUN = IPCON
      IF (IPCON.EQ.4) GO TO 40
      CALL DEFIN2
      IF (IPCON.EQ.2) GO TO 20
      CALL DEFIN3
      CALL DEFIN4
      IF (IPCON.EQ.1) GO TO 60
  20  CALL DEFIN5
  30  IF (IPRUN.EQ.2 .OR. IPRUN.EQ.3) LOPTS = -1
      IF (IPRUN.EQ.3) LOPTW = -1
      GO TO 60
  40  CALL DEFIN6
      GO TO 60
C
  50  IPLOT = 0
  60  CALL INP001
      IF (LCODE.EQ.0) GO TO 110
  70  DO 80 IP=1,IPMAX
        IF (ICODE.EQ.IPT1(IP) .AND. JCODE.EQ.IPT2(IP)) GO TO 90
  80  CONTINUE
      GO TO 190
  90  IF (IPRUN.EQ.0) GO TO 100
      IF (IP.EQ.IPMAX) IP = 99
      IPCON = IP
      IF (IPRUN.EQ.1 .OR. IPRUN.EQ.3) CALL INP031
      RETURN
C
 100  IPRUN = IP
      GO TO 30
C
 110  IF (IPRUN.EQ.0) GO TO 190
      CALL INP002
      IF (ICODE.EQ.0) GO TO 210
      GO TO (120, 130, 140, 150, 160, 170), ICODE
 120  CALL INP010
      GO TO 180
 130  IF (IPRUN.GE.4) GO TO 200
      CALL INP020
      GO TO 60
 140  IF (IPRUN.EQ.2) GO TO 200
      CALL INP030
      GO TO 60
 150  CALL INP040
      GO TO 180
 160  IF (IPRUN.LE.1) GO TO 200
      CALL INP050
      GO TO 180
 170  IF (IPRUN.NE.4) GO TO 200
      CALL INP060
      IF (IPLOT.EQ.0) GO TO 60
      RETURN
C
 180  IF (LCODE.NE.0) GO TO 70
      GO TO 60
C
 190  I = 3
      GO TO 220
 200  I = 4
      GO TO 220
 210  I = 2
 220  CALL ERROR(I)
      GO TO 60
      END
C
C ======================================================================
C
C DOREDI - SUBROUTINES: PART I
C
C ======================================================================
C
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DESIA
C FILTER DESIGN -- FIRST SECTION
C-----------------------------------------------------------------------
C
      SUBROUTINE DESIA
C
      COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
      COMMON /TOLNOR/ VSNN, NDEGN, NBN
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
      COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
C
      SFA = 0.
      IF (NOUT.GE.3) CALL OUT030
      CALL DESI00
      IF (NOUT.GE.4) CALL OUT031
      CALL DESI01
      IF (NOUT.GE.3) CALL OUT033
      GO TO (10, 20, 20, 30), IAPRO
  10  CALL DESI11
      GO TO 40
  20  CALL DESI12
      GO TO 40
  30  CALL DESI14
C
  40  VSNN = VSN
      NDEGN = NDEG
      NBN = NJ
      IF (NOUT.GE.5) CALL OUT034
      IF (NOUT.GE.3) CALL OUT035
      IF (NOUT.LT.4) GO TO 50
      IZ = 0
      IP = 0
      CALL OUT036
  50  CALL TRANZE
      IF (NOUT.LT.4) GO TO 60
      IZ = 1
      IP = 0
      CALL OUT036
  60  CALL TRBIZE
      CALL BLNUMZ
      CALL ROMEG
      IF (NOUT.GE.3) CALL OUT032
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DESIB
C FILTER DESIGN  --  SECOND SECTION
C-----------------------------------------------------------------------
C
      SUBROUTINE DESIB
C
      COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /TOLNOR/ VSNN, NDEGN, NBN
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
      COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
      VSN = VSNN
      NDEG = NDEGN
      NJ = NBN
      GO TO (10, 20, 20, 30), IAPRO
  10  CALL DESI21
      GO TO 40
  20  CALL DESI22
      GO TO 40
  30  CALL DESI24
C
  40  IF (NOUT.GE.3) CALL OUT037
      IF (NOUT.LT.4) GO TO 50
      IZ = 0
      IP = 1
      CALL OUT036
  50  CALL TRANPO
      IF (LSOUT.EQ.(-1)) GO TO 60
      IF (NOUT.LT.4) GO TO 70
  60  IZ = 1
      IP = 1
      CALL OUT036
      IF (LSOUT.EQ.(-1)) CALL OUT016
  70  CALL TRBIPO
      IF (NOUT.GE.3) CALL OUT038
      CALL BLDENZ
      IF (NOUT.GE.2) CALL OUT011
      IF (NOUT.GE.4) CALL OUT039
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DESI00
C TRANSFORM TOLERANCE SCHEME
C-----------------------------------------------------------------------
C
      SUBROUTINE DESI00
C
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
      IF (ITYP.GE.3) NDEG = (NDEG+1)/2
      IF (NDEG.NE.0) ADEG = FLOAT(NDEG)/(1.+EDEG)
      IF (LVSN.LT.0) GO TO 20
      IF (LVSN.EQ.0) GO TO 10
      CALL PARCHA
      CALL VSNMAX
  10  CALL TRANSN
  20  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DESI01
C DESIGN OF BUTTERWORTH, CHEBYSHEV (PASSBAND OR STOPBAND), AND
C ELLIPTIC FILTERS
C-----------------------------------------------------------------------
C
      SUBROUTINE DESI01
C
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
      IF (LVSN.NE.0) GO TO 20
      CALL PARCHA
      CALL DEGREE
      Q = ADEG*(1.+EDEG) + 0.5
      N = INT(Q)
      M = INT(ADEG)
      IF (FLOAT(M).LT.ADEG) M = M + 1
      N = MAX0(M,N)
      IF (NDEG.EQ.0) GO TO 10
      IF (NDEG.GE.N) GO TO 20
      CALL ERROR(15)
      IF (NDEGF.EQ.(-1)) STOP
C
  10  NDEG = N
C
  20  IF (NDEG.LE.MAXDEG) RETURN
      CALL ERROR(25)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   PARCHA
C COMPUTATION OF THE PARAMETERS OF THE CHARACTERISTIC FUNCTION
C-----------------------------------------------------------------------
C
      SUBROUTINE PARCHA
C
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
      GD1 = 0.
      GD2 = -1.
      IF (ADELP.GT.0.) GD1 = SQRT((2.-ADELP)*ADELP)/(1.-ADELP)
      IF (ADELS.GT.0.) GD2 = SQRT(1.-ADELS*ADELS)/ADELS
      ACAP12 = GD1/GD2
      IF (ACAP12.GT.0.) GO TO 60
      GO TO (10, 20, 20, 30), IAPRO
  10  ACAP12 = VSN**(-ADEG)
      GO TO 40
  20  Q = ARCOSH(VSN)*ADEG
      ACAP12 = 1./COSH(Q)
      GO TO 40
  30  CALL BOUND
  40  IF (GD2.EQ.(-1.)) GO TO 50
      GD1 = ACAP12*GD2
      ADELP = 1. - 1./SQRT(1.+GD1*GD1)
      GO TO 60
  50  GD2 = GD1/ACAP12
      ADELS = 1./SQRT(1.+GD2*GD2)
  60  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   VSNMAX
C COMPUTATION OF THE NORMALIZED CUTOFF FREQUENCY FOR A FILTER
C DEGREE AND THE RIPPLES GIVEN
C-----------------------------------------------------------------------
C
      SUBROUTINE VSNMAX
C
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
C
      AC = 1./ACAP12
      E = 1./ADEG
      GO TO (10, 20, 20, 30), IAPRO
  10  VSN = AC**E
      GO TO 40
  20  VSN = COSH(E*ARCOSH(AC))
      GO TO 40
  30  CALL BOUND
  40  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   BOUND
C CALCULATION OF A BOUND FOR VSN OR ACAP12 FOR ELLIPTIC FILTERS
C-----------------------------------------------------------------------
C
      SUBROUTINE BOUND
C
      DOUBLE PRECISION DPI, DOMI
      DOUBLE PRECISION DE, DKK, DDEG, DCAP12, DEG, DCAP14, DQ, DK1,
     *    DAB, DMAX, DDE
      DOUBLE PRECISION DK, DF, DD
      DOUBLE PRECISION DFF, DELLK
C
      COMMON /CONST2/ DPI, DOMI
      COMMON /TOLSN/ VSN, VD, A
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
C
      DIMENSION DK(3), DF(3)
C
      DATA DE /1.D00/
      DFF(DD) = (DELLK(DD)*DKK/DELLK(DSQRT(DE-DD*DD)))**II - DDEG
      IF (ACAP12.LE.0.) GO TO 10
      DCAP12 = DBLE(ACAP12)
      DEG = DBLE(1./ADEG)
      II = 1
      GO TO 20
  10  DCAP12 = DE/DBLE(VSN)
      DEG = DBLE(ADEG)
      II = -1
  20  DCAP14 = DSQRT(DE-DCAP12*DCAP12)
      DKK = DELLK(DCAP14)/DELLK(DCAP12)
      DQ = DEXP(-DPI*DKK*DEG)
      DK1 = 4.D00*DSQRT(DQ)
      IF (DK1.LT.DE) GO TO 30
      DQ = 2.D00*DQ
      DQ = (DE-DQ)/(DE+DQ)
      DQ = DQ*DQ
      DK1 = DSQRT(DE-DQ*DQ)
  30  DK(1) = DK1
      DK(2) = (DE+DK(1))/2.D00
      DDEG = DBLE(ADEG)
      DF(1) = DFF(DK(1))
      DF(2) = DFF(DK(2))
  40  DK(3) = DK(1) - DF(1)*(DK(1)-DK(2))/(DF(1)-DF(2))
      DF(3) = DFF(DK(3))
      IF (DABS(DF(3)).LT.1.D-6) GO TO 60
      DMAX = 0.D00
      DO 50 J=1,3
        DAB = DABS(DF(J))
        IF (DMAX.GT.DAB) GO TO 50
        JJ = J
        DMAX = DAB
  50  CONTINUE
      IF (JJ.EQ.3) GO TO 40
      DK(JJ) = DK(3)
      DF(JJ) = DF(3)
      GO TO 40
  60  IF (ACAP12.LE.0.) GO TO 70
      DDE = DE/DK(3)
      VSN = SNGL(DDE)
      RETURN
  70  ACAP12 = SNGL(DK(3))
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DEGREE
C COMPUTATION OF THE MINIMUM FILTER DEGREE (ADEG)
C-----------------------------------------------------------------------
C
      SUBROUTINE DEGREE
C
      DOUBLE PRECISION DE, DCAP02, DCAP04, DCAP12, DCAP14, DADEG
      DOUBLE PRECISION DELLK
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
C
      GO TO (10, 20, 20, 30), IAPRO
C
  10  ADEG = ALOG(1./ACAP12)/ALOG(VSN)
      RETURN
C
  20  ADEG = ARCOSH(1./ACAP12)/ARCOSH(VSN)
      RETURN
C
  30  DE = 1.D00
      DCAP02 = DE/DBLE(VSN)
      DCAP04 = DSQRT(DE-DCAP02*DCAP02)
      DCAP12 = DBLE(ACAP12)
      DCAP14 = DSQRT(DE-DCAP12*DCAP12)
      DADEG = (DELLK(DCAP02)*DELLK(DCAP14))/(DELLK(DCAP04)*DELLK(DCAP12)
     *    )
      ADEG = SNGL(DADEG)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DESI11
C BUTTERWORTH FILTER
C COMPUTATION OF THE ZEROS AND LOCATIONS OF THE EXTREMA
C-----------------------------------------------------------------------
C
      SUBROUTINE DESI11
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
      ADELTA = VSN**NDEG
C
      NH = NDEG/2
      NJ = (NDEG+1)/2
      FDEG = FLOAT(NDEG)
      FN = PI/2./FDEG
C
      DO 10 I=1,NJ
        NZERO(I) = 0
        III = I + I - 1
        Q = FN*FLOAT(III)
        PREN(I) = SIN(Q)
        PIMN(I) = COS(Q)
  10  CONTINUE
C
      FN = 2.*FN
      NZERO(1) = NDEG
      NZM(1) = 1
      SM(1,1) = 0.
      NZM(2) = 1
      SM(1,2) = 1.
      NZM(3) = 1
      SM(1,3) = VSN
      NZM(4) = 1
      SM(1,4) = FLMA
C
      UGC = GD2/ADELTA
      OGC = GD1
      SM(17,4) = 1.
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DESI12
C CHEBYSHEV FILTER (PASSBAND OR STOPBAND)
C COMPUTATION OF THE ZEROS AND THE LOCATIONS OF THE EXTREMA
C-----------------------------------------------------------------------
C
      SUBROUTINE DESI12
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
      ADELTA = COSH(FLOAT(NDEG)*ARCOSH(VSN))
C
      FA = 1.
      NH = NDEG/2
      NJ = (NDEG+1)/2
      FN = PI/(2.*FLOAT(NDEG))
C
      DO 10 I=1,NJ
        NZERO(I) = 0
        INJ = I + I - 1
        Q = FN*FLOAT(INJ)
        PREN(I) = SIN(Q)
        PIMN(I) = COS(Q)
  10  CONTINUE
C
      FN = 2.*FN
C
      IF (IAPRO.EQ.3) GO TO 40
C
      M = NJ + 1
      DO 20 I=1,NJ
        J = M - I
        SM(I,1) = PIMN(J)
  20  CONTINUE
      NZM(1) = NJ
      M = NH + 1
      DO 30 I=1,M
        MI = M - I
        Q = FLOAT(MI)*FN
        SM(I,2) = COS(Q)
  30  CONTINUE
      NZM(2) = M
      SM(1,3) = VSN
      NZM(3) = 1
      SM(1,4) = FLMA
      NZM(4) = 1
      NZERO(1) = NDEG
C
      UGC = GD2/ADELTA
      OGC = GD1
      GO TO 80
C
  40  SM(1,1) = 0.
      NZM(1) = 1
      SM(1,2) = 1.
      NZM(2) = 1
      DO 50 I=1,NJ
        INJ = NJ - I
        Q = FLOAT(INJ)*FN
        SM(INJ+1,3) = VSN/COS(Q)
  50  CONTINUE
      NZM(3) = NJ
      DO 60 I=1,NH
        NZERO(I) = 2
        Q = PIMN(I)
        FA = FA*Q*Q
        SM(I,4) = VSN/Q
  60  CONTINUE
      IF (NH.EQ.NJ) GO TO 70
      NZERO(NJ) = 1
      SM(NJ,4) = FLMA
  70  NZM(4) = NJ
C
      UGC = GD2
      OGC = GD1*ADELTA
  80  ACK = FA
      SM(17,4) = 1.
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DESI14
C ELLIPTIC FILTER
C COMPUTATION OF THE ZEROS AND THE LOCATIONS OF THE EXTREMA
C-----------------------------------------------------------------------
C
      SUBROUTINE DESI14
C
      DOUBLE PRECISION DPI, DOMI
      DOUBLE PRECISION DK, DKS, DCAP02, DCAP04, DE
      DOUBLE PRECISION DSK(16)
      DOUBLE PRECISION DCAP01, DQ, DN, DU, DM
      DOUBLE PRECISION DELLK, DSN2, DEL1, DEL2, DDE, DDELTA, DDELT
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CONST2/ DPI, DOMI
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
      COMMON /RESIN2/ DK, DKS, DCAP02, DCAP04
      EQUIVALENCE (PREN(1),DSK(1))
C
      DATA DE /1.D00/
C
      DCAP02 = DE/DBLE(VSN)
      DCAP01 = DSQRT(DCAP02)
      DCAP04 = DSQRT(DE-DCAP02*DCAP02)
      DK = DELLK(DCAP02)
      DKS = DELLK(DCAP04)
C
      DQ = DEXP(-DPI*DKS/DK)
C
      NH = NDEG/2
      M = NDEG + 1
      NJ = M/2
      MH = NH + 1
C
      DN = DK/DBLE(FLOAT(NDEG))
C
      DEL1 = DE
      IF (NH.EQ.0) GO TO 20
      DO 10 I=1,NH
        INH = M - I - I
        DU = DN*DBLE(FLOAT(INH))
        DM = DSN2(DU,DK,DQ)
        DEL1 = DEL1*DM*DCAP01
        DSK(I) = DM
        J = MH - I
        SM(J,1) = SNGL(DM)
        NZERO(I) = 2
        DDE = DE/(DCAP02*DM)
        SM(I,4) = SNGL(DDE)
  10  CONTINUE
      GO TO 30
  20  SM(1,1) = 0.
  30  KJ = NJ - 1
      MJ = NJ + 1
      DEL2 = DE
      IF (KJ.EQ.0) GO TO 50
      DO 40 I=1,KJ
        NDEGI = NDEG - I - I
        DU = DN*DBLE(FLOAT(NDEGI))
        DM = DSN2(DU,DK,DQ)
        J = MJ - I
        SM(J,2) = SNGL(DM)
        DDE = DE/(DCAP02*DM)
        SM(I+1,3) = SNGL(DDE)
        DEL2 = DEL2*DM*DCAP01
  40  CONTINUE
      GO TO 60
  50  SM(NDEG,2) = 1.
      SM(1,3) = VSN
C
  60  DDELT = DEL1*DEL1
      ADELTA = SNGL(DDELT)
      ACK = 1./ADELTA
      IF (NH.EQ.NJ) GO TO 80
      ACK = ACK*SNGL(DCAP01)
      DDELTA = DEL2*DEL2*DCAP01
      ADELTA = SNGL(DDELTA)
      DSK(NJ) = 0.D00
      NZERO(NJ) = 1
      SM(NJ,4) = FLMA
C
      IF (NH.EQ.0) GO TO 90
      DO 70 I=1,NH
        J = MJ - I
        SM(J,1) = SM(J-1,1)
        SM(I,2) = SM(I+1,2)
  70  CONTINUE
      SM(1,1) = 0.
      GO TO 90
C
  80  SM(MH,3) = FLMA
      SM(1,2) = 0.
C
  90  NZM(1) = NJ
      NZM(4) = NJ
      NZM(2) = MH
      NZM(3) = MH
      SM(MH,2) = 1.
      SM(1,3) = VSN
      UGC = GD2*ADELTA
      OGC = GD1/ADELTA
      SM(17,4) = 1.
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DESI21
C BUTTERWORTH FILTER
C COMPUTATION OF THE POLES
C-----------------------------------------------------------------------
C
      SUBROUTINE DESI21
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
C COMPUTATION OF CONSTANT C AND REDUCED TOLERANCE SCHEME
C
      IF (ACX.LT.999.) GO TO 20
      IF ((OGC-UGC).LT.FLMI) GO TO 10
      AC = (2.*ADELP/(ADELTA*ADELS))**0.33333
      ACX = ALOG10(AC/UGC)/ALOG10(OGC/UGC)
      IF (ACX.GE.0. .AND. ACX.LE.1.) GO TO 30
  10  ACX = 0.5
  20  AC = UGC*(OGC/UGC)**ACX
  30  RDELP = 1. - SQRT(1./(1.+AC*AC))
      Q = AC*ADELTA
      RDELS = SQRT(1./(1.+Q*Q))
C
C COMPUTATION OF FACTOR SFA AND POLES
C
      SFA = 1./AC
      Q = AC**(-1./FLOAT(NDEG))
C
      DO 40 I=1,NJ
        SPR(I) = -Q*PREN(I)
        SPI(I) = Q*PIMN(I)
  40  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DESI22
C CHEBYSHEV FILTER (PASSBAND OR STOPBAND)
C COMPUTATION OF THE POLES
C-----------------------------------------------------------------------
C
      SUBROUTINE DESI22
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
      IF (ACX.LT.999.) GO TO 20
      IF ((OGC-UGC).LT.FLMI) GO TO 10
      IF (IAPRO.EQ.2) Q = 1./ADELTA
      IF (IAPRO.EQ.3) Q = ADELTA*ADELTA
      AC = (2.*ADELP*Q/ADELS)**0.33333
      ACX = ALOG10(AC/UGC)/ALOG10(OGC/UGC)
      IF (ACX.GE.0. .AND. ACX.LE.1.) GO TO 30
  10  ACX = 0.5
  20  AC = UGC*(OGC/UGC)**ACX
C
C COMPUTATION OF THE REDUCED TOLERANCE SCHEME
C
  30  Q = AC
      IF (IAPRO.EQ.3) Q = Q/ADELTA
      Q = 1. + Q*Q
      RDELP = 1. - SQRT(1./Q)
      Q = AC
      IF (IAPRO.EQ.2) Q = Q*ADELTA
      Q = 1. + Q*Q
      RDELS = SQRT(1./Q)
C
C COMPUTATION OF THE FACTOR SFA AND THE POLES
C
      IF (IAPRO.EQ.3) GO TO 40
      SFA = 2./(AC*2.**NDEG)
      Q = -1./AC
      GO TO 50
C
  40  SFA = ACK
      Q = AC
  50  Q = ARSINH(Q)/FLOAT(NDEG)
      QR = SINH(Q)
      QI = COSH(Q)
      IF (IAPRO.EQ.3) GO TO 70
      DO 60 I=1,NJ
        SPR(I) = QR*PREN(I)
        SPI(I) = QI*PIMN(I)
  60  CONTINUE
      RETURN
C
  70  DO 80 I=1,NH
        Q = PIMN(I)*QI
        QA = PREN(I)*QR
        QQ = Q*Q
        QQA = QA*QA
        SFA = SFA/(QQ+QQA)
        SPR(I) = -VSN/(QQ/QA+QA)
        SPI(I) = VSN/(Q+QQA/Q)
  80  CONTINUE
      IF (NH.EQ.NJ) RETURN
      SPI(NJ) = 0.
      Q = VSN/QR
      SFA = SFA*Q
      SPR(NJ) = -Q
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DESI24
C ELLIPTIC FILTER
C COMPUTATION OF THE REDUCED TOLERANCE SCHEME, THE FACTOR SFA, AND
C THE POLES
C-----------------------------------------------------------------------
C
      SUBROUTINE DESI24
C
      DOUBLE PRECISION DPI, DOMI
      DOUBLE PRECISION DK, DKS, DCAP02, DCAP04, DE
      DOUBLE PRECISION DSK(16)
      DOUBLE PRECISION DU, DR, DQ, DUD, DUC, DRC, DRD, DM
      DOUBLE PRECISION DELLK, DSN2
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CONST2/ DPI, DOMI
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
      COMMON /RESIN2/ DK, DKS, DCAP02, DCAP04
      EQUIVALENCE (PREN(1),DSK(1))
C
      DATA DE /1.D00/
C
C IF ACX NOT DEFINED, COMPUTE A SYMMETRICAL USAGE OF THE TOLERANCE
C SCHEME
C
      IF (ACX.LT.999.) GO TO 20
      IF ((OGC-UGC).LT.FLMI) GO TO 10
      AC = (2.*ADELP/(ADELTA*ADELS))**0.333333333
      ACX = ALOG10(AC/UGC)/ALOG10(OGC/UGC)
      IF (ACX.GE.0. .AND. ACX.LE.1.) GO TO 30
  10  ACX = 0.5
C
  20  AC = UGC*(OGC/UGC)**ACX
C
C COMPUTATION OF THE REDUCED TOLERANCE SCHEME
C
  30  Q = AC*ADELTA
      DU = DE/DBLE(Q)
      RDELP = 1. - SQRT(1./(1.+Q*Q))
      Q = 1. + AC*AC/(ADELTA*ADELTA)
      RDELS = SQRT(1./Q)
C
C COMPUTATION OF THE FACTOR SFA AND THE POLES
C
      Q = AC*ACK
      IF (NH.EQ.NJ) Q = SQRT(1.+Q*Q)
      SFA = 1./Q
C
      DR = DBLE(ADELTA)
      DR = DR*DR
      DQ = DBLE(Q)
      CALL DELI1(DQ, DU, DR)
      DU = DQ
      DQ = DSQRT(DE-DR*DR)
      DQ = DELLK(DR)
      DU = DK*DU/(DQ*DBLE(FLOAT(NDEG)))
      DQ = DEXP(-DPI*DK/DKS)
      DU = -DSN2(DU,DKS,DQ)
      DQ = DU*DU
      DUD = DE - DCAP04*DCAP04*DQ
      DUD = DSQRT(DUD)
      DUC = DSQRT(DE-DQ)
      DO 40 I=1,NJ
        DR = DSK(I)
        DRC = DR*DR
        DRD = DE - DCAP02*DCAP02*DRC
        DRC = DSQRT(DE-DRC)
        DM = DE - DQ*DRD
        DRD = DSQRT(DRD)
        DRD = DRD*DU*DUC*DRC/DM
        SPR(I) = SNGL(DRD)
        DR = DR*DUD/DM
        SPI(I) = SNGL(DR)
  40  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   TRANZE
C REACTANCE TRANSFORMATION OF THE ZEROS AND THE LOCATIONS OF THE
C EXTREMA
C-----------------------------------------------------------------------
C
      SUBROUTINE TRANZE
C
      DOUBLE PRECISION DR, DQI
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
C
      DIMENSION MSM(4)
      FA = 1.
      IF (ITYP.EQ.1) GO TO 190
      IF (ITYP.EQ.3) GO TO 60
C
      ME = NZM(4)
      DO 10 I=1,ME
        Q = SM(I,4)
        IF (Q.LT.FLMA) FA = FA*Q
  10  CONTINUE
C
      FA = FA*FA
C
C LOWPASS - HIGHPASS
C
      DO 50 J=1,4
        ME = NZM(J)
        DO 40 I=1,ME
          QI = SM(I,J)
          IF (ABS(QI).LT.FLMI) GO TO 20
          QI = 1./QI
          GO TO 30
  20      QI = FLMA
  30      SM(I,J) = QI
  40    CONTINUE
  50  CONTINUE
      GO TO 90
  60  DO 80 J=1,2
        ME = NZM(J)
        MA = ME + 1
        ME = ME/2
        DO 70 I=1,ME
          QI = SM(I,J)
          II = MA - I
          SM(I,J) = SM(II,J)
          SM(II,J) = QI
  70    CONTINUE
  80  CONTINUE
C
  90  IF (ITYP.EQ.2) GO TO 190
C
C LOWPASS - BANDPASS TRANSFORMATION
C
      QA = 2.*A
      NN = NDEG + 1
      IF (ITYP.EQ.4) GO TO 110
C
      MSM(1) = 1
      IF (NZM(1).NE.1) MSM(1) = NDEG
      MSM(2) = 2
      IF (NZM(2).NE.1) MSM(2) = NN
      DO 100 J=3,4
        MSM(J) = 2*NZM(J)
 100  CONTINUE
      GO TO 130
C
 110  DO 120 J=1,2
        MSM(J) = 2*NZM(J)
 120  CONTINUE
      MSM(3) = 2
      IF (NZM(3).NE.1) MSM(3) = NN
      MSM(4) = 1
      IF (NZM(4).NE.1) MSM(4) = NDEG
C
 130  S = 1.
      DO 180 J=1,4
        ME = NZM(J)
        MA = MSM(J)
        NZM(J) = MA
        IF (J.EQ.3) S = -1.
        DO 170 I=1,ME
          QR = SM(I,J)
          NU = NZERO(I)
          IF (ABS(QR).LT.FLMA) GO TO 150
          IF (J.NE.4) GO TO 140
          FA = FA*(VD/A)**NU
 140      QI = QR
          GO TO 160
C
 150      QR = QR/QA
          DR = DBLE(QR)
          DQI = DSQRT(DR*DR+1.D00)
          QI = SNGL(DQI)
 160      SM(I,J) = QI - S*QR
          II = MA - I + 1
          IF (ABS(QR).LT.FLMI) NU = 2*NU
          IF (J.EQ.4) NZERO(II) = NU
          SM(II,J) = QI + S*QR
 170    CONTINUE
 180  CONTINUE
C
 190  DO 220 J=1,4
        ME = NZM(J)
        DO 210 I=1,ME
          Q = SM(I,J)
          IF (Q.LT.FLMA) GO TO 200
          IF (J.NE.4 .OR. ITYP.GE.3) GO TO 210
          NU = NZERO(I)
          FA = FA*VD**NU
          GO TO 210
 200      SM(I,J) = Q*VD
 210    CONTINUE
 220  CONTINUE
      SM(17,4) = SM(17,4)*FA
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   TRANPO
C REACTANCE TRANSFORMATION OF THE POLES
C-----------------------------------------------------------------------
C
      SUBROUTINE TRANPO
C
      DOUBLE PRECISION DR, DI, DQ
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
      IF (ITYP.EQ.1) GO TO 90
      IF (ITYP.EQ.3) GO TO 40
C
      DO 30 I=1,NJ
        QR = SPR(I)
        QI = SPI(I)
        QH = QR*QR + QI*QI
        IF (ABS(QI).GT.FLMI) GO TO 10
        SFA = -SFA/QR
        GO TO 20
  10    SFA = SFA/QH
  20    QI = QI/QH
        IF (ABS(QI).LT.FLMI) QI = 0.
        SPI(I) = QI
        SPR(I) = QR/QH
  30  CONTINUE
      IF (ITYP.EQ.2) GO TO 90
C
  40  QA = 2.*A
      NN = NJ
      NJ = NDEG
      NDEG = 2*NDEG
C
      ME = NJ
      DO 80 I=1,NN
        QR = SPR(I)/QA
        QI = SPI(I)/QA
        IF (ABS(QI).GE.FLMA) GO TO 70
        DR = DBLE(QR)
        DI = DBLE(QI)
        DQ = DI*DI
        DI = DI*DR*2.D00
        DR = DR*DR - DQ - 1.D00
        CALL DSQRTC(DR, DI, DR, DI)
        QZ = SNGL(DR)
        QN = SNGL(DI)
        IF (ABS(QN).GT.FLMI) GO TO 60
        JJ = NJ + ME
        DO 50 II=ME,NJ
          J = JJ - II
          SPR(J+1) = SPR(J)
          SPI(J+1) = SPI(J)
  50    CONTINUE
        NJ = NJ + 1
        ME = ME + 1
  60    SPR(I) = QR + QZ
        SPI(I) = QI + QN
        SPR(ME) = QR - QZ
        SPI(ME) = QN - QI
        ME = ME - 1
        GO TO 80
C
  70    SPR(I) = QR
        SPI(I) = FLMA
        NJ = NJ + 1
        SPR(NJ) = QR
        SPI(NJ) = 0.
  80  CONTINUE
C
  90  DO 100 I=1,NJ
        SPR(I) = SPR(I)*VD
        SPI(I) = SPI(I)*VD
 100  CONTINUE
C
      SFA = SFA*SM(17,4)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   TRBIZE
C BILINEAR TRANSFORMATION OF THE ZEROS AND THE LOCATIONS OF THE EXTREMA
C-----------------------------------------------------------------------
C
      SUBROUTINE TRBIZE
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
C
      FA = 1.
      DO 50 J=1,4
        ME = NZM(J)
        DO 40 I=1,ME
          QI = SM(I,J)
          ZM(I,J) = 2.*ATAN(QI)
          IF (J.NE.4) GO TO 40
          IF (QI.GE.FLMA) GO TO 10
          IF (QI.LT.FLMI) GO TO 20
          QQI = QI*QI
          Q = 1. + QQI
          ZZR(I) = (1.-QQI)/Q
          ZZI(I) = 2.*QI/Q
          NU = NZERO(I)/2
          FA = FA*Q**NU
          GO TO 40
C
  10      ZZR(I) = -1.
          GO TO 30
C
  20      ZZR(I) = 1.
  30      ZZI(I) = 0.
  40    CONTINUE
  50  CONTINUE
C
      SM(17,1) = FA
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   TRBIPO
C BILINEAR TRANSFORMATION OF THE POLES
C-----------------------------------------------------------------------
C
      SUBROUTINE TRBIPO
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
      ZFA = SFA*SM(17,1)
C
      DO 20 I=1,NJ
        QR = SPR(I)
        Q = 1. - QR
        QI = SPI(I)
        IF (ABS(QI).LT.FLMI) GO TO 10
        QQR = QR*QR
        QQI = QI*QI
        ZFA = ZFA/(Q-QR+QQR+QQI)
        Q = 1./(Q*Q+QQI)
        ZPR(I) = (1.-QQR-QQI)*Q
        ZPI(I) = 2.*QI*Q
        GO TO 20
C
  10    ZPR(I) = (1.+QR)/Q
        ZPI(I) = 0.
        ZFA = ZFA/Q
  20  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   BLNUMZ
C BUILD NUMERATOR BLOCKS OF SECOND ORDER
C-----------------------------------------------------------------------
C
      SUBROUTINE BLNUMZ
C
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
      DIMENSION NZE(16)
      COMMON /SCRAT/ ADUM(32)
      EQUIVALENCE (ADUM(1),NZE(1))
C
      N = 0
      ME = NZM(4)
      DO 10 I=1,ME
        NZE(I) = NZERO(I)
  10  CONTINUE
C
      DO 70 I=1,ME
        QR = ZZR(I)
        NZ = NZE(I)
C
  20    IF (NZ.EQ.0) GO TO 70
        N = N + 1
        B2(N) = 1.
        IF (NZ.EQ.1) GO TO 30
        B1(N) = -2.*QR
        B0(N) = 1.
        NZ = NZ - 2
        IF (NZ.GT.0) GO TO 20
        GO TO 70
C
  30    IF (I.EQ.ME) GO TO 60
        MA = I + 1
        DO 40 II=MA,ME
          IF (ZZI(II).EQ.0.) GO TO 50
  40    CONTINUE
        GO TO 60
C
  50    QRR = ZZR(II)
        B1(N) = -QR - QRR
        B0(N) = QR*QRR
        NZE(II) = NZE(II) - 1
        GO TO 70
C
  60    B1(N) = -QR
        B0(N) = 0.
C
  70  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   BLDENZ
C BUILD DENOMINATOR BLOCKS OF SECOND ORDER  --  Z-DOMAIN
C-----------------------------------------------------------------------
C
      SUBROUTINE BLDENZ
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
      NB = (NDEG+1)/2
      N = 0
      FACT = ZFA
C
      DO 40 I=1,NB
        N = N + 1
        QR = ZPR(N)
        QI = ZPI(N)
        IF (ABS(QI).LT.FLMI) GO TO 10
        C1(I) = -2.*QR
        C0(I) = QR*QR + QI*QI
        GO TO 40
C
  10    IF (N.GE.NJ) GO TO 20
        IF (ABS(ZPI(N+1)).LT.FLMI) GO TO 30
  20    C1(I) = -QR
        C0(I) = 0.
        GO TO 40
C
  30    N = N + 1
        QI = ZPR(N)
        C1(I) = -QR - QI
        C0(I) = QR*QI
  40  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   ROMEG
C REALIZED FREQUENCIES OMEGA
C-----------------------------------------------------------------------
C
      SUBROUTINE ROMEG
C
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
C
      N2 = NZM(2)
      N3 = NZM(3)
      GO TO (10, 20, 30, 40), ITYP
  10  ROM(1) = ZM(N2,2)
      ROM(2) = ZM(1,3)
      GO TO 50
  20  ROM(1) = ZM(1,3)
      ROM(2) = ZM(N2,2)
      GO TO 50
  30  ROM(1) = ZM(N3,3)
      ROM(2) = ZM(1,2)
      ROM(3) = ZM(N2,2)
      ROM(4) = ZM(1,3)
      GO TO 50
  40  N2 = N2/2
      ROM(1) = ZM(N2,2)
      ROM(4) = ZM(N2+1,2)
      ROM(3) = ZM(1,3)
      ROM(2) = ZM(N3,3)
  50  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   TRANSN
C COMPUTATION OF THE PARAMETERS OF THE NORMALIZED LOWPASS
C-----------------------------------------------------------------------
C
      SUBROUTINE TRANSN
C
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
      TAN2(AA) = SIN(AA/2.)/COS(AA/2.)
C
      V1 = TAN2(OM(1))
      V2 = TAN2(OM(2))
      IF (ITYP.LE.2) GO TO 210
      V3 = TAN2(OM(3))
      V4 = TAN2(OM(4))
      IF (ITYP.EQ.3) GO TO 10
      Q = V1
      V1 = -V4
      V4 = -Q
      Q = V2
      V2 = -V3
      V3 = -Q
      IF (LSYM.NE.0) LSYM = 5 - LSYM
      IF (LVSN.EQ.0) GO TO 10
      J = 5 - LVSN
      LVSN = LSYM
      LSYM = J
C
  10  J = LSYM + 1
      JJ = LVSN/2 + 1
      GO TO (20, 80, 120, 160, 180), J
  20  J = NORMA + 1
      GO TO (30, 30, 40, 70), J
  30  VDQ1 = V2*V3
      VSN1 = VDQ1/V1 - V1
      Q = V4 - VDQ1/V4
      IF (Q.LT.VSN1) VSN1 = Q
      A1 = 1./(V3-V2)
      VSN1 = VSN1*A1
      GO TO (40, 50, 40), J
  40  VDQ = V1*V4
      A = V2/(VDQ-V2*V2)
      Q = V3/(V3*V3-VDQ)
      IF (Q.LT.A) A = Q
      VSN = A*(V4-V1)
      IF (NORMA.EQ.2) GO TO 200
      IF (VSN.GE.VSN1) GO TO 200
  50  VDQ = VDQ1
  60  VSN = VSN1
      A = A1
      GO TO 200
C
  70  VDQ = SQRT(V1*V2*V3*V4)
      A1 = V3/(V3*V3-VDQ)
      VSN1 = (V4-VDQ/V4)*A1
      A = V2/(VDQ-V2*V2)
      VSN = (VDQ/V1-V1)*A
      IF (VSN.GE.VSN1) GO TO 200
      GO TO 60
C
  80  GO TO (90, 100, 110), JJ
  90  VDQ = V2*V3
      VSN = V4 - VDQ/V4
      GO TO 190
 100  V3 = V4*(V4+VSN*V2)/(V2+VSN*V4)
 110  VDQ = V2*V3
      A = 1./(V3-V2)
      GO TO 200
 120  GO TO (130, 150, 140), JJ
 130  VDQ = V1*V4
      A = V3/(V3*V3-VDQ)
      GO TO 170
 140  V4 = V3*(V1+VSN*V3)/(V3+VSN*V1)
 150  VDQ = V1*V4
      A = VSN/(V4-V1)
      GO TO 200
 160  VDQ = V1*V4
      A = V2/(VDQ-V2*V2)
 170  VSN = (V4-V1)*A
      GO TO 200
 180  VDQ = V2*V3
      VSN = VDQ/V1 - V1
 190  A = 1./(V3-V2)
      VSN = VSN*A
C
 200  VD = SQRT(VDQ)
      A = A*VD
      IF (ITYP.LE.3) GO TO 270
      A = A/VSN
      GO TO 270
C
 210  J = LVSN*2 + ITYP
      GO TO (220, 220, 230, 240, 250, 260), J
 220  VSN = V2/V1
      GO TO (250, 240), J
 230  VD = V2/VSN
      GO TO 270
 240  VD = V2
      GO TO 270
 250  VD = V1
      GO TO 270
 260  VD = V1*VSN
C
 270  RETURN
      END
C
C ======================================================================
C
C DOREDI - SUBROUTINES: PART II
C
C ======================================================================
C
C
C-----------------------------------------------------------------------
C SUBROUTINE:   GRID01
C BUILD A START GRID FOR THE SEARCH OF THE GLOBAL EXTREMES OF THE
C TRANSFER FUNCTION WITH ROUNDED COEFFICIENTS
C-----------------------------------------------------------------------
C
      SUBROUTINE GRID01
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
      COMMON /CGRID/ GR(64), NGR(12)
C
      M = NDEG
      IF (NDEG.GT.22) M = 22
      IF (ITYP.GE.3) M = M/2
      M1 = M/2
      FM = PI/FLOAT(2*M)
      M11 = 2*M1 + 1
      FM1 = PI/FLOAT(M11)
      M = M + 1
      PH2 = 1.5*PI
C
      NGR(8) = 1
C
C BANDPASS
C
      N = 0
      GO TO (70, 10, 20, 40), ITYP
C
C HIGHPASS
C
  10  IA = 1
      IB = NZM(2)
      GO TO 30
C
C BANDPASS
C
  20  IA = NZM(3)
      IB = 1
  30  N = N + 1
      GR(N) = 0.
      GO TO 50
C
C BANDSTOP
C
  40  IA = 1
      IB = NZM(2)/2 + 1
  50  J = 1
      PH1 = 0.
      JA = 3
      JB = 2
      GO TO 350
  60  NGR(1) = N + 1
  70  K = 2
      L = 1
      MAK = 1
      MAL = 1
      MEK = NZM(K)
      MEL = NZM(L)
      J = 1
      GO TO (80, 90, 100, 110), ITYP
C
C LOWPASS
C
  80  ID = 1
      NGR(1) = 1
      PH = 0.
      IF (MEK.GT.MEL) GO TO 370
      IQ = MEK
      MEK = MEL
      MEL = IQ
      K = 1
      L = 2
      GO TO 370
C
C HIGHPASS
C
  90  ID = -1
      GO TO 130
C
C BANDPASS
C
 100  MEK = MEK/2
      MEL = MEK
      GO TO 120
C
C BANDSTOP
C
 110  MAK = MEK/2 + 1
      MAL = MEL/2 + 1
 120  ID = 1
 130  PH = PH2
      GO TO 370
C
 140  NGR(2) = N
      GO TO (150, 160, 160, 160), ITYP
C
C LOWPASS
C
 150  IA = NZM(2)
      IB = 1
      JA = 2
      JB = 3
      PH1 = PH2
      J = 2
      GO TO 350
 160  NGR(3) = N
      IF (ITYP.LE.2) GO TO 230
      K = 2
      L = 1
      MEK = NZM(K)
      MEL = NZM(L)
      IF (ITYP.EQ.3) GO TO 180
      IF (MEK.GT.MEL) GO TO 170
      K = 1
      L = 2
 170  MEK = MEK/2
      MEL = MEL/2
      MAK = 1
      MAL = 1
      GO TO 190
 180  MAK = MEK/2 + 1
      MAL = MAK
 190  ID = 1
      PH = 0.
      J = 2
      GO TO 370
C
 200  NGR(4) = N
      IF (ITYP.EQ.3) GO TO 210
C
C BANDSTOP
C
      IA = NZM(2)/2
      IB = NZM(3)
      GO TO 220
C
C BANDPASS
C
 210  IA = NZM(2)
      IB = 1
 220  JA = 2
      JB = 3
      PH1 = PH2
      J = 3
C
      GO TO 350
 230  IF (ABS(GR(N)-PI).LT.FLMI) GO TO 240
      N = N + 1
      GR(N) = PI
 240  NGR(9) = N
      NGR(5) = N
C
C BANDSTOP
C
      K = 3
      L = 4
      MEK = NZM(K)
      MEL = NZM(L)
      MAK = 1
      MAL = 1
      J = 3
      ID = 1
      GO TO (280, 250, 270, 250), ITYP
C
C HIGHPASS
C
 250  IF (MEK.GT.MEL) GO TO 260
      IQ = MEK
      MEK = MEL
      MEL = IQ
      K = 4
      L = 3
 260  PH = 0.
      ID = -1
      IF (ITYP.NE.4) GO TO 330
C
C BANDSTOP
C
      MAK = MEK/2 + 1
      MAL = MEL/2 + 1
      GO TO 280
C
C BANDPASS
C
 270  MEK = MEK/2
      MEL = MEL/2
 280  PH = PH2
      GO TO 380
C
 290  NGR(6) = N
      J = 4
      GO TO (500, 500, 300, 320), ITYP
C
C BANDPASS
C
 300  IF (MEK.GT.MEL) GO TO 310
      IQ = MEK
      MEK = MEL
      MEL = IQ
      K = 4
      L = 3
 310  MAK = MEK + 1
      MAL = MEL + 1
      MEK = NZM(K)
      MEL = NZM(L)
      GO TO 330
C
C BANDSTOP
C
 320  MEK = MEK/2
      MEL = (MEL+1)/2
      K = 4
      L = 3
      MAK = 1
      MAL = 1
 330  PH = 0.
      GO TO 380
C
 340  NGR(7) = N
      GO TO 500
C
 350  QA = ZM(IA,JA)
      QD = ZM(IB,JB) - QA
      IF (PH1.NE.0.) QA = QA + QD
      DO 360 I=1,M1
        N = N + 1
        Q = FLOAT(I)*FM1 + PH1
        GR(N) = QA + SIN(Q)*QD
 360  CONTINUE
      GO TO (60, 160, 230), J
C
 370  GO TO (460, 390, 460, 390), IAPRO
 380  GO TO (460, 460, 390, 390), IAPRO
C
C MULTIPLE EXTREMA
C
 390  IF (ID.LT.0) GO TO 400
      IK = MAK
      IL = MAL
      GO TO 410
C
 400  IK = MEK
      IL = MEL
 410  DO 450 I=MAK,MEK
        N = N + 1
        GR(N) = ZM(IK,K)
        IF (IL.GT.MEL .OR. IL.LT.MAL) GO TO 440
        IF (NDEG.LE.22) GO TO 420
        IF (IL.NE.MEL .AND. IL.NE.MAL) GO TO 430
 420    N = N + 1
        GR(N) = ZM(IL,L)
 430    IL = IL + ID
 440    IK = IK + ID
 450  CONTINUE
      GO TO 490
C
C SINGLE EXTREMA
C
 460  QA = ZM(MEK,K)
      QD = ZM(MEL,L)
      IF (QD.GT.QA) GO TO 470
      Q = QA
      QA = QD
      QD = Q
 470  QD = QD - QA
      IF (PH.NE.0.) QA = QA + QD
      DO 480 I=1,M
        N = N + 1
        II = I - 1
        Q = FLOAT(II)*FM + PH
        GR(N) = QA + SIN(Q)*QD
 480  CONTINUE
 490  GO TO (140, 200, 290, 340), J
C
 500  NGR(10) = NGR(3) + 1
      NGR(11) = NGR(5) + 1
      NGR(12) = NGR(6) + 1
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   GRID03
C START VALUES FOR THE CALCULATION OF A GRID
C-----------------------------------------------------------------------
C
      SUBROUTINE GRID03
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
C
      IF (ITYP.LE.0 .OR. ITYP.GT.4) GO TO 70
      IF (NDEG.LE.0) GO TO 70
      IAPRO = 1
      DO 10 I=1,4
        NZM(I) = (ITYP+1)/2
  10  CONTINUE
      ME = 2
      IF (ITYP.GE.3) ME = 4
      DO 20 I=1,ME
        IF (OM(I).LE.0.) GO TO 80
  20  CONTINUE
C
      GO TO (30, 40, 50, 60), ITYP
  30  ZM(1,1) = 0.
      ZM(1,2) = OM(1)
      ZM(1,3) = OM(2)
      ZM(1,4) = PI
      GO TO 100
  40  ZM(1,1) = PI
      ZM(1,2) = OM(2)
      ZM(1,3) = OM(1)
      ZM(1,4) = 0.
      GO TO 100
  50  NZM(1) = 1
      ZM(1,1) = (OM(2)+OM(3))/2.
      ZM(1,2) = OM(2)
      ZM(2,2) = OM(3)
      ZM(1,3) = OM(4)
      ZM(2,3) = OM(1)
      ZM(1,4) = PI
      ZM(2,4) = 0.
      GO TO 100
  60  ZM(1,1) = 0.
      ZM(2,1) = PI
      ZM(1,2) = OM(1)
      ZM(2,2) = OM(4)
      ZM(1,3) = OM(3)
      ZM(2,3) = OM(2)
      NZM(4) = 1
      ZM(1,4) = (OM(2)+OM(3))/2.
      GO TO 100
C
  70  IE = 20
      GO TO 90
  80  IE = 23
  90  CALL ERROR(IE)
 100  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   GRID04
C CHANGE GRID FOR UNSYMMETRICAL TOLERANCE SCHEME
C-----------------------------------------------------------------------
C
      SUBROUTINE GRID04
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /CGRID/ GR(64), NGR(12)
C
      IF (ITYP.LE.2) GO TO 130
      DO 10 I=1,4
        IF (OM(I).EQ.0.) GO TO 130
  10  CONTINUE
C
      I = 1
      MA = NGR(1)
      ME = NGR(2)
      IF (ITYP.EQ.3) OMM = OM(2)
      IF (ITYP.EQ.4) OMM = OM(4)
      N = 1
      GO TO 50
C
  20  MA = NGR(10)
      ME = NGR(4)
      IF (ITYP.EQ.3) OMM = OM(3)
      IF (ITYP.EQ.4) OMM = OM(1)
      N = 4
      GO TO 90
C
  30  I = 2
      MA = NGR(11)
      ME = NGR(6)
      IF (ITYP.EQ.3) OMM = OM(4)
      IF (ITYP.EQ.4) OMM = OM(2)
      N = 11
      GO TO 50
C
  40  MA = NGR(12)
      ME = NGR(7)
      IF (ITYP.EQ.3) OMM = OM(1)
      IF (ITYP.EQ.4) OMM = OM(3)
      N = 7
      GO TO 90
C
  50  DO 60 II=MA,ME
        J = ME + MA - II
        IF ((GR(J)-OMM).LE.FLER) GO TO 70
  60  CONTINUE
      GO TO 80
  70  GR(J) = OMM
      NGR(N) = J
  80  IF (I.EQ.1) GO TO 20
      GO TO 40
C
  90  DO 100 II=MA,ME
        IF ((GR(II)-OMM).GE.(-FLER)) GO TO 110
 100  CONTINUE
      GO TO 120
 110  GR(II) = OMM
      NGR(N) = II
 120  IF (I.EQ.1) GO TO 30
C
 130  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   COEFW
C SEARCH OF THE MINIMUM COEFFICIENT WORDLENGTH
C COMPUTATION OF THE OPTIMIZATION CRITERION
C-----------------------------------------------------------------------
C
      SUBROUTINE COEFW
C
      COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
      COMMON /CEPSIL/ EPS(2), PMAX
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMPLEX ZP, ZPS
      COMMON /CPOL/ ZP(16,2), ZPS(16,2)
C
      IWLRR = IWLR
      IWLR = 100
      IWLG = IWLMA
      IWLD = 0
      IWLM = IWL
      II1 = 1
      II2 = 0
      II3 = 0
      CALL COPY01
      IF (IPRUN.GT.0) GO TO 20
      IF (LSEQ.NE.(-1)) GO TO 10
      CALL ALLO02
      CALL COPY01
      IF (NOUT.GE.2) CALL OUT011
  10  CALL ALLO01
      CALL COPY01
  20  IF (NOUT.EQ.4) CALL OUT043
      CALL RECO
  30  CALL COPY02
      IF (IPRUN.GT.0) GO TO 70
      CALL ROUND(4, 5)
      IF (ISCAL.NE.3) GO TO 40
      CALL POLLOC
      CALL COPY03
  40  AMAX = 1.
      FACT = 1.
      DO 60 J=1,NB
        IB = J
  50    CALL SCAL01
        CALL SCAL02
        IF (JSTRUD.EQ.1) GO TO 60
        IF (IB.NE.NB) GO TO 60
        IB = IB + 1
        GO TO 50
  60  CONTINUE
      CALL RECO
      CALL ROUND(1, 3)
      GO TO 80
  70  CALL ROUND(1, 5)
  80  CALL STABLE
      IF (NOUT.GE.5) CALL OUT046
      IF (ISTAB.NE.(-1)) GO TO 250
      CALL EPSILO
      IDEPS = 0
      DO 90 K=1,2
        IF (EPS(K).GT.1.) IDEPS = IDEPS + 1
  90  CONTINUE
      IF (NOUT.GE.5) CALL OUT043
      IF (NOUT.GE.4) CALL OUT044
      ADEPS = EPS(1) - EPS(2)
      IF (LOPTW.NE.(-1)) GO TO 260
      IF (IWL.NE.IWLL) GO TO 100
      ADEPSL = ADEPS
      IDEPSL = IDEPS
 100  IF (IDEPS.EQ.0) GO TO 110
      IF (IWL.LE.IWLD .OR. IWL.GT.IWLG) GO TO 120
      IWLD = IWL
      ADEPSD = ADEPS
      IDEPSD = IDEPS
      GO TO 120
 110  IF (IWL.GE.IWLG) GO TO 120
      IWLG = IWL
      ADEPSG = ADEPS
      IF (IWLD.LE.IWLG) GO TO 120
      IWLD = 0
      II3 = 0
 120  GO TO (130, 140, 180), II1
 130  II1 = 2
      IF (IDEPS.EQ.2) II1 = 3
      GO TO 120
 140  IF (IDEPS.NE.0) GO TO 170
 150  II2 = II2 + 1
 160  IWL = IWL - 1
      GO TO 220
 170  II3 = II3 + 1
      IF (II3.LT.2) GO TO 160
      IF (II2.GE.2) GO TO 270
      IWL = IWLM + 1
      II1 = 3
      GO TO 220
 180  IF (IDEPS.EQ.0) GO TO 210
 190  II3 = II3 + 1
 200  IWL = IWL + 1
      GO TO 220
 210  II2 = II2 + 1
      IF (II2.LT.2) GO TO 200
      IF (II3.GE.2) GO TO 270
      IWL = IWLM - 1
      II1 = 2
 220  IF (IWL.GE.IWLMI .AND. IWL.LE.IWLMA) GO TO 30
 230  IF (II1.EQ.3) GO TO 240
      IF (IDEPS.EQ.0) GO TO 270
      II1 = 3
      II2 = 2
      GO TO 190
 240  II1 = 2
      IF (IDEPS.NE.0) GO TO 270
      II3 = 2
      GO TO 150
 250  IF (LOPTW.NE.(-1)) GO TO 270
      IF (II1.NE.1) GO TO 230
      IWL = IWL + 1
      IF (IWL.LE.IWLMA) GO TO 30
 260  ADEPSD = ADEPS
      IWLG = IWL
      IWLD = IWL
 270  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   RECO
C TRANSFORM THE COEFFICIENTS FOR THE CHOSEN REALIZATION
C-----------------------------------------------------------------------
C
      SUBROUTINE RECO
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CRECO/ JRCO, JECO, JDCO, JJDCO(5), JMAXV, JTRB2, LREF
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
C
      DIMENSION ACO(16,5)
      EQUIVALENCE (ACO(1,1),B2(1))
      EQUIVALENCE (IRCO(1),IRB2), (IRCO(2),IRB1), (IRCO(3),IRB0)
      EQUIVALENCE (IRCO(4),IRC1), (IRCO(5),IRC0)
C
      ALOG2 = ALOG(2.)
C
C REPRESENTATION OF THE COEFFICIENTS WITH DIFFERENCES
C
      IF (JDCO.EQ.0) GO TO 70
      DO 60 J=1,5
        ACOM = 0.
        DO 10 I=1,NB
          ACOM = AMAX1(ABS(ACO(I,J)),ACOM)
  10    CONTINUE
        IF (ACOM.GE.FLMI) GO TO 20
        ID = -100
        GO TO 30
  20    ID = ALOG(ACOM)/ALOG2
        Q = 2.**ID
        IF (Q.GE.ACOM) GO TO 30
        ID = ID + 1
        Q = Q*2.
  30    JJ = JJDCO(J)
        IF (JJ.EQ.0) GO TO 60
        BCOM = 0.5*Q
        DO 50 I=1,NB
          AC = ACO(I,J)
          ACA = ABS(AC)
          IF (ACA.GE.BCOM) GO TO 40
          IF (JJ.NE.1) GO TO 50
  40      IDCO(I,J) = ID
  50    CONTINUE
  60  CONTINUE
C
C RANGE OF COEFFICIENT NUMBERS
C
  70  IF (JRCO.EQ.0) GO TO 180
      DO 120 J=1,5
        ACOM = 0.
        BCOM = 0.
        DO 80 I=1,NB
          ACS = ACOEFS(ACO(I,J),0,IDCO(I,J),0)
          ACOM = AMIN1(ACOM,ACS)
          BCOM = AMAX1(BCOM,ACS)
  80    CONTINUE
        ACOM = AMAX1(ABS(ACOM),BCOM)
        IF (ACOM.GE.FLMI) GO TO 90
        IQ = -100
        GO TO 110
  90    IQ = ALOG(ACOM)/ALOG2
        Q = 2.**IQ
        IF (Q.LT.ACOM) GO TO 100
        IF (JMAXV.LE.1) GO TO 110
        IF (JMAXV.EQ.2 .AND. Q.GT.BCOM) GO TO 110
        IF (JMAXV.EQ.3 .AND. Q.GT.ACOM) GO TO 110
 100    IQ = IQ + 1
 110    IF (JRCO.GT.0) IQ = MAX0(IQ,0)
        IRCO(J) = IQ
 120  CONTINUE
      JJ = IABS(JRCO)
      GO TO (180, 130, 130, 150, 160), JJ
 130  IQ = MAX0(IRB2,IRB1,IRB0)
      DO 140 J=1,3
        IRCO(J) = IQ
 140  CONTINUE
      IF (JJ.EQ.2) GO TO 180
      IQ = MAX0(IRC1,IRC0)
      IRC1 = IQ
      IRC0 = IQ
      GO TO 180
 150  IQ = MAX0(IRB1,IRC1)
      IRB1 = IQ
      IRC1 = IQ
      GO TO 180
 160  IQ = MAX0(IRB2,IRB1,IRB0,IRC1,IRC0)
      DO 170 J=1,5
        IRCO(J) = IQ
 170  CONTINUE
C
C COMPUTE "PSEUDO" FLOATING-POINT EXPONENTS
C
 180  IF (JECO.EQ.0) GO TO 230
      IECOM = 0
      DO 220 J=1,5
        DO 210 I=1,NB
          QA = ACOEFS(ACO(I,J),0,IDCO(I,J),IRCO(J))
          QA = ABS(QA)
          IF (QA.GE.FLMI) GO TO 190
          IQ = JECO
          GO TO 200
 190      IQ = -ALOG(QA)/ALOG2
          IF (IQ.GT.JECO) IQ = JECO
 200      IECO(I,J) = IQ
          IECOM = MAX0(IECOM,IQ)
 210    CONTINUE
 220  CONTINUE
 230  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   ROUND
C ROUNDING OF THE CHANGED COEFFICIENTS
C-----------------------------------------------------------------------
C
      SUBROUTINE ROUND(IRA, IRB)
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /CRECO/ JRCO, JECO, JDCO, JJDCO(5), JMAXV, JTRB2, LREF
      COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
C
      DIMENSION ACO(16,5)
      EQUIVALENCE (ACO(1,1),B2(1))
C
      ALSB = 2.**(1-IWL)
      ELSB = 1. - ALSB
      DO 120 I=1,NB
        DO 110 J=IRA,IRB
          AC = ACO(I,J)
          IF (J.NE.2 .AND. J.NE.3) GO TO 10
          IF (ACO(I,1).EQ.AC) GO TO 110
  10      ACS = ACOEFS(AC,IECO(I,J),IDCO(I,J),IRCO(J))
          ACSA = ABS(ACS)
          IF (ACSA.LT.ALSB .AND. J.EQ.1) GO TO 110
          IF (ACSA.GT.ELSB) GO TO 30
          R = 0.5
          IF (J.NE.1 .OR. JTRB2.LE.0) GO TO 20
          R = 0.
          IF (IDCO(I,J).GT.(-100)) R = 1.
  20      ACSA = AINT(ACSA/ALSB+R)
          ACSA = ACSA*ALSB
          GO TO 50
  30      IF (JMAXV.GE.1) GO TO 40
          IF (JMAXV.EQ.(-2) .AND. AC.LT.0.) GO TO 40
          ACSA = ELSB
          GO TO 50
  40      ACSA = 1.
  50      BC = SIGN(ACSA,ACS)
          BC = ACOEF(BC,IECO(I,J),IDCO(I,J),IRCO(J))
          ACO(I,J) = BC
          IF (J.NE.1) GO TO 110
          F = BC/AC
          JJ = IABS(JTRB2)
          GO TO (100, 80, 60), JJ
  60      DO 70 K=2,3
            ACO(I,K) = ACO(I,K)*F
  70      CONTINUE
  80      II = I + 1
          IF (II.GT.NB) GO TO 100
          DO 90 K=1,3
            ACO(II,K) = ACO(II,K)/F
  90      CONTINUE
          GO TO 110
 100      FACT = FACT/F
 110    CONTINUE
 120  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   STABLE
C CHECK STABILITY OF THE ROUNDED SYSTEM
C-----------------------------------------------------------------------
C
      SUBROUTINE STABLE
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
     *    SC0(16)
      COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
C
      DO 120 I=1,NB
        IT = 0
        CC0 = C0(I)
        CC1 = C1(I)
        IF (CC0.GE.1.) GO TO 20
        AC1 = ABS(CC1)
  10    IF (AC1.GE.1.+CC0) GO TO 30
        IF (IT.EQ.0) GO TO 120
        GO TO 110
  20    IT = IT + 1
        GO TO 40
  30    IT = IT + 2
  40    IF (LSTAB.NE.(-1)) GO TO 130
        ALSB = 2.**(1-IWL)
        Q0 = ACOEF(ALSB,IECO(I,5),-100,IRCO(5))
        IF (IT.GE.2) GO TO 50
        CC0 = CC0 - Q0
        GO TO 10
  50    IF (IT.GE.3) GO TO 80
        AC0 = CC0
  60    AC0 = AC0 + Q0
        IF (AC0.GE.1.) GO TO 70
        IF (AC1.GE.1.+AC0) GO TO 60
        Q3 = AC1 - ABS(SC1(I))
        Q2 = AC0 - SC0(I)
        A = Q3*Q3 + Q2*Q2
        A = SQRT(A)
        GO TO 80
  70    IT = 3
  80    Q1 = ACOEF(ALSB,IECO(I,4),-100,IRCO(4))
        BC1 = AC1
  90    BC1 = BC1 - Q1
        IF (BC1.GE.1.+AC0) GO TO 90
        IF (IT.EQ.3) GO TO 100
        Q3 = BC1 - ABS(SC1(I))
        Q2 = CC0 - SC0(I)
        B = Q3*Q3 + Q2*Q2
        B = SQRT(B)
        IF (A.GE.B) GO TO 100
        BC1 = AC1
        CC0 = AC0
 100    C1(I) = SIGN(BC1,CC1)
 110    C0(I) = CC0
 120  CONTINUE
      ISTAB = -1
      RETURN
 130  ISTAB = 0
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   EPSILO
C COMPUTATION OF THE GLOBAL EXTREMA OF THE TRANSFER FUNCTION
C WITH ROUNDED COEFFICIENTS
C COMPUTATION OF THE ERROR VALUES EPSILON
C-----------------------------------------------------------------------
C
      SUBROUTINE EPSILO
C
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /CGRID/ GR(64), NGR(12)
      COMMON /CEPSIL/ EPS(2), PMAX
C
C MAXIMUM IN THE PASSBAND(S)
C
      ME = 3
      IF (ITYP.EQ.3) ME = 5
      CALL SMAX(PMAX, 8, ME, 1)
      IF (ITYP.NE.4) GO TO 10
      CALL SMAX(QMAX, 10, 5, 1)
      PMAX = AMAX1(PMAX,QMAX)
C
C MINIMUM IN THE PASSBAND(S)
C
  10  ME = 2
      IF (ITYP.EQ.3) ME = 4
      CALL SMAX(PMIN, 1, ME, -1)
      IF (ITYP.NE.4) GO TO 20
      CALL SMAX(QMIN, 10, 4, -1)
      PMIN = AMIN1(PMIN,QMIN)
C
C MAXIMUM IN THE STOPBAND(S)
C
  20  ME = 6
      IF (ITYP.EQ.4) ME = 7
      CALL SMAX(STMAX, 11, ME, 1)
      IF (ITYP.NE.3) GO TO 30
      CALL SMAX(QMAX, 12, 7, 1)
      STMAX = AMAX1(STMAX,QMAX)
C
C COMPUTATION OF EPSILON
C
  30  EPS(1) = (1.-PMIN/PMAX)/ADELP
      EPS(2) = STMAX/PMAX/ADELS
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OPTPAR
C OPTIMIZATION OF THE COEFFICIENT WORDLENGTH BY VARIATION OF THE
C DESIGN PARAMETER ACX
C-----------------------------------------------------------------------
C
      SUBROUTINE OPTPAR
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /COPTCP/ IS1, IS2, IABO, IWLGS, IWLGN, IWLGP, ADEPSS,
     *    ADEPSN, ADEPSP, ACXS, ACXN, ACXP
C
      IF (LOPTW.NE.(-1)) GO TO 40
      IF (IWLD.GT.IWLL .AND. IDEPSL.LT.2) GO TO 20
      IF (IDEPSD.NE.2) GO TO 10
      ADEPSD = ADEPSG
      IWLD = IWLG
  10  IWLL = IWLD
      GO TO 30
  20  ADEPSD = ADEPSL
      IWLD = IWLL
  30  IWL = IWLD
  40  IABO = IABO + 1
      IF (ADEPSD.LT.0.) GO TO 80
      IF (IS1.NE.(-1)) GO TO 50
      IF (IS2.NE.(-1)) GO TO 70
      GO TO 60
  50  IF (IS2.NE.(-1)) GO TO 60
      IS1 = -1
  60  IS2 = -((IS2+2)/2)
      ACXA = (ACX+ACXN)/2.
      GO TO 120
  70  ACXA = (ACX+ACXMI)/2.
      GO TO 120
C
  80  IF (IS1.EQ.(-1)) GO TO 90
      IF (IS2.NE.(-1)) GO TO 110
      GO TO 100
  90  IF (IS2.NE.(-1)) GO TO 100
      IS1 = 0
 100  IS2 = -((IS2+2)/2)
      ACXA = (ACX+ACXP)/2.
      GO TO 120
 110  ACXA = (ACX+ACXMA)/2.
C
 120  IF (IWLG.LT.IWLGS) GO TO 130
      IF (ADEPSS-ABS(ADEPSD).LE.FLMI .OR. IWLG.GT.IWLGS) GO TO 140
 130  IWLGS = IWLG
      ADEPSS = ABS(ADEPSD)
      ACXS = ACX
 140  IF (ADEPSD.LT.0.) GO TO 160
      IF (IWLG.LT.IWLGP) GO TO 150
      IF (ADEPSP-ADEPSD.LE.FLMI .OR. IWLG.GT.IWLGP) GO TO 190
 150  IWLGP = IWLG
      ADEPSP = ADEPSD
      ACXP = ACX
      GO TO 180
 160  IF (IWLG.LT.IWLGN) GO TO 170
      IF (ADEPSN-ADEPSD.GE.FLMI .OR. IWLG.GT.IWLGN) GO TO 190
 170  IWLGN = IWLG
      ADEPSN = ADEPSD
      ACXN = ACX
 180  IABO = 0
 190  ACX = ACXA
      RETURN
      END
C
C ======================================================================
C
C DOREDI - SUBROUTINES: PART III
C
C ======================================================================
C
C
C-----------------------------------------------------------------------
C SUBROUTINE:   GRID02
C BUILD A START GRID OUT OF THE POLE LOCATIONS FOR THE SEARCH OF THE
C GLOBAL EXTREMES OF THE TRANSFER FUNCTION
C-----------------------------------------------------------------------
C
      SUBROUTINE GRID02
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /CGRID/ GR(64), NGR(12)
C
      NGR(1) = 0
      J = 1
      GR(1) = 0.
      DO 10 I=1,NB
        Q = -C1(I)/2.
        QN = C0(I)
        IF (ABS(QN).LT.FLMI) GO TO 10
        QN = QN - Q*Q
        IF (QN.LE.0.) GO TO 10
        J = J + 1
        GR(J) = ATAN2(SQRT(QN),Q)
  10  CONTINUE
      J = J + 1
      GR(J) = PI
      NGR(8) = 1
      NGR(9) = J
C
C SORTING OF THE LOCATIONS
C
      J = J - 1
  20  JJ = 0
      DO 30 I=1,J
        Q1 = GR(I)
        Q2 = GR(I+1)
        IF (Q1.LT.Q2) GO TO 30
        JJ = 1
        GR(I) = Q2
        GR(I+1) = Q1
  30  CONTINUE
      IF (JJ.NE.0) GO TO 20
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   ANANOI
C NOISE ANALYSIS OF A GIVEN CASCADE
C-----------------------------------------------------------------------
C
      SUBROUTINE ANANOI
C
      COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /COPTST/ LOPTS, ISTOR
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /CPOW/ PNU, PNC, AND, ITCORP
      COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
      COMMON /CFFUNC/ PHI(5), BF2(5), BF1(5), BF0(5), IBB(5), ICOR
      COMMON /CNOISE/ RI, RIN, RE, REN, FAC
      COMMON /COPST2/ PN(100,2), TF(100,2), TFA(100,2)
C
      IQ = IWLR
      IF (IQ-1.GT.IWLMA) IQ = IWLMA + 1
      IQ = IQ - 2
      ALSBI = 2.**IQ
      IB = NB
      IF (LOPTS.NE.0) GO TO 20
      DO 10 I=1,NB
        JSEQN(I) = I
        JSEQD(I) = I
  10  CONTINUE
      GO TO 40
  20  DO 30 I=1,NB
        JSEQN(I) = ISEQ(I,1)
        JSEQD(I) = ISEQ(I,2)
  30  CONTINUE
  40  SAND = 0.
      SPNU = 0.
      SPNC = 0.
      CALL ALLOND
      CALL SMAX(Q, 8, 9, 1)
      FAC1 = FACT/Q
      IF (ISCAL.NE.0 .OR. JSTRUD.EQ.1) FACT = 1.
      AMAX = 1.
      ITCORP = -1
C
      DO 80 I=1,NB
        IB = I
        IF (ISTRU.GE.30) GO TO 60
  50    IF (ISCAL.EQ.0) SCA = 1.
        IF (ISCAL.NE.0) CALL SCAL01
        CALL ALNS01
        CALL ALNS02
        GO TO 70
  60    CONTINUE
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C HERE INSERT THE CALLS TO YOUR OWN STRUCTURE IMPLEMENTATIONS
C SEE SUBROUTINE OPTLST
C
        CALL ERROR(8)
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
  70    Q = FAC1/FACT
        CALL POWER
        AND = AND*Q
        PNC = Q
        Q = Q*Q
        PNU = PNU*Q
        CALL PCORP
        PNC = PNC*Q
        SAND = SAND + AND
        SPNU = SPNU + PNU
        SPNC = SPNC + PNC
        IF (NOUT.GE.3) CALL OUT026
        IF (ISCAL.NE.0) CALL SCAL02
        IF (JSTRUD.EQ.1 .OR. IB.NE.1) FAC1 = FAC1/SCA
        IF (JSTRUD.EQ.1) GO TO 80
        IF (IB.NE.NB) GO TO 80
        IB = IB + 1
        GO TO 50
  80  CONTINUE
      CALL TCORP
      RIN = SPNU + SAND*SAND + PNC
      FAC = FACT/FAC1
      Q = FAC*FAC
      RI = RIN*Q
      IF (NOUT.LT.3) GO TO 90
      IB = 0
      AND = SAND
      PNU = SPNU
      CALL OUT026
  90  CALL SPOW(RE)
      REN = RE/Q
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OPTLST
C OPTIMIZATION OF THE PAIRING AND ORDERING BY THE PROCEDURE WITH
C LIMITED STORAGE
C DYNAMIC PROGRAMMING FOR SYSTEMS OF AN ORDER LESS THAN OR EQUAL TO TEN
C DIRECT PROCEDURE FOR THE PARAMETER ISTOR EQUAL TO ONE
C-----------------------------------------------------------------------
C
      SUBROUTINE OPTLST
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /COPTST/ LOPTS, ISTOR
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
      COMMON /COPST1/ LSEQN(16,100,2), LSEQD(16,100,2)
      COMMON /COPST2/ PN(100,2), TF(100,2), TFA(100,2)
      COMMON /CPOW/ PNU, PNC, AND, ITCORP
      COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
      COMMON /CFFUNC/ PHI(5), BF2(5), BF1(5), BF0(5), IBB(5), ICOR
      COMMON /CNOISE/ RI, RIN, RE, REN, FAC
C
      DIMENSION ISEQN(16), ISEQD(16)
      EQUIVALENCE (ISEQN(1),ISEQ(1,1)), (ISEQD(1),ISEQ(1,2))
C
      ITCORP = 0
      ALSBI = 2.**(IWLMA-1)
      II = 1
      PN(1,2) = 0.
      TF(1,2) = 1.
      TFA(1,2) = 1.
      CALL SMAX(Q, 8, 9, 1)
      FAC1 = FACT/Q
      JS = 1
      IS = 2
      IN = 1
      PNTM = FLMA
      MA = NB
      IF (LOPTS.NE.(-1)) MA = 1
C
      DO 220 IK=1,NB
        IIK = IK
        MK = IK - 1
        IQ = IS
        IS = JS
        JS = IQ
        JN = IN
        IN = 0
C
C EXTENSION OF THE ASSIGNED ORDERING
C
        DO 210 JK=1,JN
          IF (IK.EQ.1) GO TO 10
          CALL CODE4(JSEQD, LSEQD(1,JK,JS), MK)
          IF (LOPTS.EQ.(-1)) CALL CODE4(JSEQN, LSEQN(1,JK,JS), MK)
  10      DO 200 I=1,NB
            IF (IK.EQ.1) GO TO 30
            DO 20 II=1,MK
              IF (JSEQD(II).EQ.I) GO TO 200
  20        CONTINUE
  30        JSEQD(IK) = I
            DO 190 J=1,MA
              IF (LOPTS.NE.(-1)) GO TO 60
              IF (IK.EQ.1) GO TO 50
              DO 40 II=1,MK
                IF (JSEQN(II).EQ.J) GO TO 190
  40          CONTINUE
  50          JSEQN(IK) = J
C
  60          FACT = TF(JK,JS)
              AMAX = TFA(JK,JS)
              PNT = PN(JK,JS)
              IB = IK
              CALL ALLOND
  70          IF (ISTRU.GE.30) GO TO 80
              CALL SCAL01
              CALL ALNS01
              CALL ALNS02
              GO TO 90
  80          CONTINUE
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C HERE INSERT THE CALLS TO YOUR OWN STRUCTURE IMPLEMENTATIONS
C SCALING OF THE IB-TH BLOCK
C ALLOCATIONS OF THE NOISE SOURCES IN THE IB-TH BLOCK
C ALLOCATION OF THE INPUT TRANSFER FUNCTIONS
C  FACT : PRODUCT OF THE IB-1 SCALING FACTORS IN FRONT
C  AMAX : SCALING OF THE (IB-1)-TH BLOCK
C  JSTRU  = 1 : SECOND COEFFICIENT OPTIMIZATION
C  JSTRU  = 2 : NO SECOND OPTIMIZATION
C  JSTRUS = 1 : NO SEPARATE SCALING FACTOR
C  JSTRUS = 2 : SEPARATE SCALING FACTOR
C  JSTRUD = 1 : NUMERATOR FIRST
C  JSTRUD = 2 : DENOMINATOR FIRST
C  BN2(I),BN1(I),BN0(I) : NOISE TRANSFER FUNCTION IN THE IB-TH BLOCK
C  BF2(I),BF1(I),BF0(I) : TRANSFER FUNCTIONS IN THE IB-TH BLOCK
C                         FROM THE INPUT TO THE NONLINEARITIES
C  IBB(I) .NE. IB : BLOCK IB ONLY NUMERATOR
C  PHI(I) : VARIANCE OF THE INPUT TRANSFER FUNCTIONS
C
C REMOVE THE FOLLOWING STATEMENT
C
              CALL ERROR(8)
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
  90          Q = FAC1/FACT
              CALL POWER
              CALL PCORP
              PNT = (PNU+PNC+AND*AND)*Q*Q + PNT
              IF (JSTRUD.EQ.1) GO TO 100
              IF (IB.NE.NB) GO TO 100
              IB = IB + 1
              GO TO 70
 100          IF (IK.EQ.NB) GO TO 170
C
C STORE INTERMEDIATE RESULT
C CHECK IF COMBINATION IS IN STORAGE
C
              IF (IK.EQ.1) GO TO 130
              IF (IN.EQ.0) GO TO 150
              DO 120 IL=1,IN
                CALL CODE5(ISEQD, LSEQD(1,IL,IS), IIK, JSEQD, II)
                IF (II.EQ.2) GO TO 120
                IF (LOPTS.NE.(-1)) GO TO 110
                CALL CODE5(ISEQN, LSEQN(1,IL,IS), IIK, JSEQN, II)
                IF (II.EQ.2) GO TO 120
 110            CONTINUE
                IF (PN(IL,IS).LE.PNT) GO TO 190
                CALL CODE3(JSEQD, LSEQD(1,IL,IS), IIK)
                IF (LOPTS.EQ.(-1)) CALL CODE3(JSEQN, LSEQN(1,IL,IS),
     *              IIK)
                PN(IL,IS) = PNT
                TF(IL,IS) = FACT
                TFA(IL,IS) = AMAX
                GO TO 190
 120          CONTINUE
 130          IF (IN.LT.ISTOR) GO TO 150
C
C SEARCH IN THE STORAGE FOR THE ALLOCATION WITH GREATEST NOISE POWER
C
              PNM = PNT
              IQ = 0
              DO 140 II=1,ISTOR
                IF (PN(II,IS).LT.PNM) GO TO 140
                IQ = II
                PNM = PN(II,IS)
 140          CONTINUE
              IF (IQ.EQ.0) GO TO 190
              GO TO 160
 150          IN = IN + 1
              IQ = IN
 160          PN(IQ,IS) = PNT
              TF(IQ,IS) = FACT
              TFA(IQ,IS) = AMAX
              CALL CODE3(JSEQD, LSEQD(1,IQ,IS), IIK)
              IF (LOPTS.EQ.(-1)) CALL CODE3(JSEQN, LSEQN(1,IQ,IS), IIK)
              GO TO 190
C
C CHECK RESULT
C
 170          IF (PNT.GE.PNTM) GO TO 190
              PNTM = PNT
              FAC = FACT
              DO 180 II=1,NB
                ISEQD(II) = JSEQD(II)
                ISEQN(II) = JSEQN(II)
 180          CONTINUE
 190        CONTINUE
 200      CONTINUE
          IF (IK.EQ.1) GO TO 220
 210    CONTINUE
 220  CONTINUE
      FAC = FAC/FAC1
      RIN = PNTM
      RI = RIN*FAC*FAC
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   ALLOND
C ALLOCATION OF THE NUMERATOR AND DENOMINATOR BLOCKS
C-----------------------------------------------------------------------
C
      SUBROUTINE ALLOND
C
      COMPLEX ZP, ZPS
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
     *    SC0(16)
      COMMON /CPOL/ ZP(16,2), ZPS(16,2)
      COMMON /COPTST/ LOPTS, ISTOR
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
      DO 40 I=1,IB
        J = JSEQD(I)
        C1(I) = SC1(J)
        C0(I) = SC0(J)
        ZP(I,1) = ZPS(J,1)
        ZP(I,2) = ZPS(J,2)
        IF (LOPTS.NE.1) GO TO 20
        IF (JSTRUD.EQ.1) GO TO 30
        IF (I.NE.IB) GO TO 10
        J = JSEQD(1)
        GO TO 30
  10    J = JSEQD(I+1)
        GO TO 30
  20    J = JSEQN(I)
  30    B2(I) = SB2(J)
        B1(I) = SB1(J)
        B0(I) = SB0(J)
        JSEQN(I) = J
  40  CONTINUE
      IF (IB.GE.NB) RETURN
      II = IB
      JJ = IB
      DO 80 I=1,NB
        DO 50 J=1,IB
          IF (JSEQD(J).EQ.I) GO TO 60
  50    CONTINUE
        II = II + 1
        C1(II) = SC1(I)
        C0(II) = SC0(I)
        ZP(II,1) = ZPS(I,1)
        ZP(II,2) = ZPS(I,2)
  60    DO 70 J=1,IB
          IF (JSEQN(J).EQ.I) GO TO 80
  70    CONTINUE
        JJ = JJ + 1
        B2(JJ) = SB2(I)
        B1(JJ) = SB1(I)
        B0(JJ) = SB0(I)
  80  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   SCAL01
C SCALING OF THE IB-TH BLOCK
C-----------------------------------------------------------------------
C
      SUBROUTINE SCAL01
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
      NNB = NB
      BMAX = 0.
      J = 1
      IF (IB.GT.NB) GO TO 30
      NB = IB
      IF (JSTRUD.EQ.1) GO TO 30
      J = 2
      BB1 = B1(IB)
      BB0 = B0(IB)
      B1(IB) = 0.
      B0(IB) = 0.
      GO TO 30
  10  IF (JSTRUS.EQ.1) GO TO 20
      J = 3
      B1(IB) = BB1
      B0(IB) = BB0
      GO TO 30
  20  IF (ISTRU.LT.25) GO TO 80
      J = 3
      B1(IB) = 0.
      B0(IB) = 1.
  30  GO TO (40, 50, 60), ISCAL
  40  CALL SMIMP(QMAX)
      GO TO 70
  50  CALL SMAX(QMAX, 8, 9, 1)
      GO TO 70
  60  CALL SPOW(QMAX)
      QMAX = SQRT(QMAX)
  70  BMAX = AMAX1(BMAX,QMAX)
      GO TO (90, 10, 80), J
  80  B1(IB) = BB1
      B0(IB) = BB0
  90  SCA = SCALM/BMAX
      IF (JSTRUS.EQ.1) GO TO 100
      IF (JSTRUD.EQ.2) GO TO 100
      IF (AMAX*SCA.GT.SCALM) SCA = SCALM/AMAX
C
C SCALING WITH A POWER OF TWO
C
 100  IF (LPOT2.NE.(-1)) GO TO 110
      I = ALOG(SCA)/ALOG(2.)
      Q = 2.**I
      IF (Q.GT.SCA) Q = Q/2.
      SCA = Q
      GO TO 120
 110  IF (IWLR.GE.100) GO TO 120
      BB0 = 2.*ALSBI
      Q = SCA*BB0
      SCA = AINT(Q)/BB0
 120  AMAX = BMAX*SCA
      FACT = FACT*SCA
      NB = NNB
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   SCAL02
C CALCULATE SCALING FACTOR INTO NUMERATOR
C-----------------------------------------------------------------------
C
      SUBROUTINE SCAL02
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
C
      I = IB
      IF (JSTRUD.EQ.1) GO TO 10
      IF (I.EQ.1) GO TO 20
      I = I - 1
  10  B2(I) = SCA
      B1(I) = B1(I)*SCA
      B0(I) = B0(I)*SCA
      FACT = FACT/SCA
  20  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   ALNS01
C ALLOCATION OF THE NOISE SOURCES FOR THE MODIFIED STRUCTURES OF THE
C FIRST AND SECOND CANONIC FORMS
C-----------------------------------------------------------------------
C
      SUBROUTINE ALNS01
C
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
      COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
C
      DO 10 I=1,5
        AQC(I) = 1.
  10  CONTINUE
C
      BB2 = SCA
      IF (IB.GT.NB) GO TO 120
      CALL QUAN(C1(IB), AQC1)
      CALL QUAN(C0(IB), AQC0)
      BB1 = B1(IB)
      BB0 = B0(IB)
C
      J = ISTRU/10
      IF (J.EQ.2) GO TO 130
      BB2 = BB2*B2(IB)
      CALL QUAN(BB2, AQB2)
      IF (JSTRU.EQ.2) GO TO 20
      BB1 = BB1*SCA
      BB0 = BB0*SCA
  20  CALL QUAN(BB1, AQB1)
      CALL QUAN(BB0, AQB0)
      J = ISTRU - 10
      IIS = 1
      GO TO (30, 30, 50, 50, 60, 70, 90, 100), J
  30  AQC(5) = AQC0
      AQC(4) = AQC1
      AQC(3) = AQB0
      AQC(2) = AQB1
  40  AQC(1) = AQB2
      GO TO 270
  50  AQC(4) = AMIN1(AQB1,AQC1)
      AQC(5) = AMIN1(AQB0,AQC0)
      GO TO 40
  60  AQC(3) = AQB0
      GO TO 80
  70  AQC(1) = AQB2
  80  AQC(2) = AQB1
      AQC(4) = AQC1
      AQC(5) = AQC0
      GO TO 270
  90  AQC(5) = AMIN1(AQB0,AQC0)
      GO TO 110
 100  AQC(5) = AQC0
 110  AQC(4) = AMIN1(AQB1,AQC1)
      GO TO 40
C
 120  AQC1 = 1.
      AQC0 = 1.
      AQB1 = 1.
      AQB0 = 1.
      GO TO 140
 130  J = ISTRU - 20
      CALL QUAN(BB1, AQB1)
      CALL QUAN(BB0, AQB0)
 140  IIS = 2
      IF (IB.NE.1) BB2 = BB2*B2(IB-1)
      CALL QUAN(BB2, AQB2)
      GO TO (170, 150, 210, 200, 190, 160, 210, 220), J
 150  AQC(3) = AQB0
 160  AQC(2) = AQB1
 170  AQC(1) = AQB2
 180  AQC(4) = AQC1
      AQC(5) = AQC0
      GO TO 230
 190  IF (IB.EQ.1) GO TO 170
      GO TO 180
 200  AQC(2) = AMIN1(AQB1,AQB0)
 210  AQC(4) = AMIN1(AQB2,AQC1,AQC0)
      IF (IB.LE.NB) GO TO 230
      AQC(1) = AQC(4)
      AQC(4) = 1.
      GO TO 270
 220  AQC(2) = AQB1
      GO TO 210
C
 230  IF (IB.EQ.1) GO TO 270
      BB1 = B1(IB-1)
      BB0 = B0(IB-1)
      IF (JSTRUS.EQ.2) GO TO 240
      BB1 = BB1*SCA
      BB0 = BB0*SCA
 240  CALL QUAN(BB1, AQB1)
      CALL QUAN(BB0, AQB0)
      GO TO (250, 270, 260, 270, 250, 270, 260, 270), J
 250  AQC(3) = AQB0
      AQC(2) = AQB1
      GO TO 270
 260  AQC(4) = AMIN1(AQC(4),AQB1,AQB0)
C
 270  DO 460 I=1,5
        BN2(I) = 0.
        BN1(I) = 0.
        BN0(I) = 0.
        IF (AQC(I).EQ.1.) GO TO 460
        IF (IB.GT.NB) GO TO 420
        IF (IIS.EQ.2) GO TO 340
        GO TO (280, 290, 300, 310, 320, 330, 300, 310), J
 280    GO TO (420, 430, 440, 430, 440), I
 290    GO TO (400, 430, 440, 430, 440), I
 300    GO TO (420, 460, 460, 430, 440), I
 310    GO TO (400, 460, 460, 430, 440), I
 320    GO TO (460, 430, 410, 430, 440), I
 330    GO TO (400, 430, 460, 430, 440), I
C
 340    GO TO (390, 350, 360, 370, 390, 380, 360, 370), J
 350    GO TO (390, 450, 450, 390, 390), I
 360    GO TO (460, 460, 460, 390, 460), I
 370    GO TO (460, 450, 460, 390, 460), I
 380    GO TO (390, 450, 460, 390, 390), I
C
 390    BN2(I) = B2(IB)
        BN1(I) = B1(IB)
        BN0(I) = B0(IB)
        GO TO 460
 400    BN2(I) = 1.
        BN1(I) = B1(IB)/B2(IB)
        BN0(I) = B0(IB)/B2(IB)
        GO TO 460
 410    BN0(I) = 1.
 420    BN2(I) = 1.
        GO TO 460
 430    BN1(I) = 1.
        GO TO 460
 440    BN0(I) = 1.
        GO TO 460
 450    BN2(I) = B2(IB)
        BN1(I) = B2(IB)*C1(IB)
        BN0(I) = B2(IB)*C0(IB)
 460  CONTINUE
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ALNS02
C ALLOCATION OF THE TRANSFER FUNCTIONS FROM THE INPUT
C TO THE NONLINEARITIES
C-----------------------------------------------------------------------
C
      SUBROUTINE ALNS02
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /CRENO/ LCNO, ICNO(16,5)
      COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
      COMMON /CFFUNC/ PHI(5), BF2(5), BF1(5), BF0(5), IBB(5), ICOR
C
      NNB = NB
      IF (IB.GT.NNB) GO TO 10
      BB2 = B2(IB)
      BB1 = B1(IB)
      BB0 = B0(IB)
  10  ICOR = 0
      IS = ISTRU - 10
      IIS = 1
      IF (IS.LE.8) GO TO 20
      IIS = 2
      IS = IS - 10
C
  20  DO 450 N=1,5
        BF2(N) = 0.
        BF1(N) = 0.
        BF0(N) = 0.
        IBB(N) = 0
        AQ = AQC(N)
        IF (AQ.EQ.1.) GO TO 400
        IF (LCNO.EQ.(-1)) AQ = 0.
        K = ICNO(IB,N)
        GO TO (400, 30, 400, 400, 40, 50, 400), K
  30    AQ = SQRT(3.)*AQ
        GO TO 60
  40    AQ = -SQRT(3.)
        GO TO 60
  50    AQ = SQRT(3.)*(AQ-1.)
  60    IF (AQ.EQ.0.) GO TO 440
        IBB(N) = IB
        IF (IB.GT.NNB) GO TO 380
        IF (IIS.EQ.2) GO TO 120
        GO TO (70, 70, 80, 90, 100, 110, 90, 90), IS
  70    GO TO (200, 210, 220, 240, 250), N
  80    GO TO (260, 400, 400, 280, 290), N
  90    GO TO (200, 400, 400, 280, 290), N
 100    GO TO (400, 210, 220, 240, 250), N
 110    GO TO (200, 210, 400, 240, 250), N
C
 120    GO TO (130, 140, 150, 160, 170, 180, 150, 190), IS
 130    GO TO (300, 310, 320, 340, 350), N
 140    GO TO (330, 360, 370, 340, 350), N
 150    GO TO (330, 400, 400, 380, 400), N
 160    GO TO (330, 390, 400, 380, 400), N
 170    GO TO (330, 310, 320, 340, 350), N
 180    GO TO (330, 360, 400, 340, 350), N
 190    GO TO (330, 360, 400, 380, 400), N
C
 200    Q = BB2
        GO TO 230
 210    Q = BB1
        GO TO 230
 220    Q = BB0
 230    BF2(N) = Q
        BF1(N) = Q*C1(IB)
        BF0(N) = Q*C0(IB)
        GO TO 410
 240    Q = -C1(IB)
        GO TO 270
 250    Q = -C0(IB)
        GO TO 270
 260    Q = 1.
 270    BF2(N) = Q*BB2
        BF1(N) = Q*BB1
        BF0(N) = Q*BB0
        GO TO 410
 280    BF2(N) = BB1 - BB2*C1(IB)
        BF1(N) = BB0 - BB2*C0(IB)
        GO TO 410
 290    BF2(N) = BB0 - BB2*C0(IB)
        BF1(N) = BB0*C1(IB) - BB1*C0(IB)
        GO TO 410
 300    IF (IB.EQ.1) GO TO 330
        IQ = IB - 1
        IBB(N) = IQ
        BF2(N) = B2(IQ)
        GO TO 410
 310    IQ = IB - 1
        IBB(N) = IQ
        BF1(N) = B1(IQ)
        GO TO 410
 320    IQ = IB - 1
        IBB(N) = IQ
        BF0(N) = B0(IQ)
        GO TO 410
 330    Q = 1.
        GO TO 230
 340    BF1(N) = -C1(IB)
        GO TO 410
 350    BF0(N) = -C0(IB)
        GO TO 410
 360    BF1(N) = BB1/BB2
        GO TO 410
 370    BF0(N) = BB0/BB2
        GO TO 410
 380    BF2(N) = 1.
        GO TO 410
 390    BF1(N) = BB1/BB2
        GO TO 370
C
 400    AQ = 0.
        GO TO 440
C
 410    NB = IBB(N)
        IF (IB.GT.NNB) GO TO 430
        IF (NB.EQ.IB) GO TO 420
        BB22 = B2(NB)
        BB12 = B1(NB)
        BB02 = B0(NB)
 420    B2(NB) = BF2(N)
        B1(NB) = BF1(N)
        B0(NB) = BF0(N)
 430    CALL SPOW(Q)
        PHI(N) = SQRT(Q)
        ICOR = 1
        IF (NB.EQ.IB) GO TO 440
        B2(NB) = BB22
        B1(NB) = BB12
        B0(NB) = BB02
C
 440    BF2(N) = BF2(N)*AQ
        BF1(N) = BF1(N)*AQ
        BF0(N) = BF0(N)*AQ
 450  CONTINUE
C
      IF (IB.GT.NNB) GO TO 460
      B2(IB) = BB2
      B1(IB) = BB1
      B0(IB) = BB0
 460  NB = NNB
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   TSTR01
C TEST ROUTINE FOR THE STRUCTURES NOS. 11-28
C-----------------------------------------------------------------------
C
      SUBROUTINE TSTR01
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /COPTST/ LOPTS, ISTOR
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
      IF (LOPTS.NE.1) GO TO 10
      IF (ISTRU.GT.20 .AND. JSTRUS.EQ.2) GO TO 30
  10  L = ISTRU
      IF (L.GT.20) L = L - 10
      IF (L.LT.15) GO TO 60
      DO 20 L=1,NB
        Q = ABS(B0(L))
        IF (Q.LT.FLMI) Q = ABS(B1(L))
        IF ((B2(L)-Q).GT.FLMI) GO TO 40
  20  CONTINUE
      GO TO 60
  30  L = 31
      GO TO 50
  40  L = 32
  50  CALL ERROR(L)
  60  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   PCORP
C COMPUTATION OF THE CORRELATED NOISE OF THE BLOCK IB BY A
C LINEARIZED MODEL
C-----------------------------------------------------------------------
C
      SUBROUTINE PCORP
C
      COMPLEX ZBL1, ZQQ, ZOM, Z, ZA, ZAA
      COMPLEX ZPHI(300)
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /CPOW/ PNU, PNC, AND, ITCORP
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
      COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
      COMMON /CFFUNC/ PHI(5), BF2(5), BF1(5), BF0(5), IBB(5), ICOR
      COMMON /COPST2/ PN(100,2), TF(100,2), TFA(100,2)
C
      EQUIVALENCE (ZPHI(1),PN(1,1))
C
      ZBL1(ZQQ,A2,A1,A0) = (ZQQ*A2+A1)*ZQQ + A0
C
      FAC = PNC
      IF (ITCORP.GE.0) GO TO 20
      DO 10 I=1,300
        ZPHI(I) = CMPLX(0.,0.)
  10  CONTINUE
      ITCORP = 1
C
  20  PNC = 0.
      IF (ICOR.EQ.0) GO TO 80
      ADOM = PI/299.
      DO 70 LL=1,300
        OM = ADOM*FLOAT(LL-1)
        ZOM = CMPLX(COS(OM),SIN(OM))
        Z = CMPLX(FACT,0.)
C
        DO 30 I=1,NB
          IF (I.NE.IB) Z = Z*ZBL1(ZOM,B2(I),B1(I),B0(I))
          Z = Z/ZBL1(ZOM,1.,C1(I),C0(I))
  30    CONTINUE
C
        IF (IB.LE.NB) Z = Z/ZBL1(ZOM,1.,C1(IB),C0(IB))
        ZA = CMPLX(0.,0.)
        DO 60 N=1,5
          IF (IBB(N).EQ.0) GO TO 60
          IF (IB.GT.NB) GO TO 50
          ZAA = ZBL1(ZOM,BF2(N),BF1(N),BF0(N))*ZBL1(ZOM,BN2(N),BN1(N),
     *        BN0(N))
          ZAA = ZAA/PHI(N)
          IF (IBB(N).EQ.IB) GO TO 40
          ZAA = ZAA*ZBL1(ZOM,1.,C1(IB),C0(IB))/ZBL1(ZOM,B2(IB),B1(IB),
     *        B0(IB))
  40      ZA = ZA + ZAA
          GO TO 60
  50      ZA = ZA + BF2(N)
  60    CONTINUE
C
        Z = Z*ZA
        IF (ITCORP.EQ.1) ZPHI(LL) = ZPHI(LL) + Z*FAC
        Q = CABS(Z)
        PNC = PNC + Q*Q
  70  CONTINUE
C
      PNC = PNC/PI/150.
  80  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   TCORP
C COMPUTATION OF THE TOTAL CORRELATED NOISE BY A LINEARIZED MODEL
C-----------------------------------------------------------------------
C
      SUBROUTINE TCORP
C
      COMPLEX ZPHI(300)
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CPOW/ PNU, PNC, AND, ITCORP
      COMMON /COPST2/ PN(100,2), TF(100,2), TFA(100,2)
      COMMON /CRENO/ LCNO, ICNO(16,5)
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
      EQUIVALENCE (ZPHI(1),PN(1,1))
C
      PNC = 0.
C
      DO 20 I=1,NB
        DO 10 J=1,5
          K = ICNO(I,J)
          IF (K.EQ.2) GO TO 30
          IF (K.EQ.5) GO TO 30
          IF (K.EQ.6) GO TO 30
  10    CONTINUE
  20  CONTINUE
      GO TO 50
C
  30  DO 40 I=1,300
        Q = CABS(ZPHI(I))
        PNC = PNC + Q*Q
  40  CONTINUE
C
      PNC = PNC/PI/150.
  50  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DESCAL
C REMOVE SCALING OF SECOND-ORDER BLOCKS
C-----------------------------------------------------------------------
C
      SUBROUTINE DESCAL
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
      DO 10 I=1,NB
        Q = B2(I)
        B2(I) = 1.
        B1(I) = B1(I)/Q
        B0(I) = B0(I)/Q
        FACT = FACT*Q
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DENORM
C TRANSFORM STRUCTURE TO STRUCTURE WITHOUT SEPARATE SCALING FACTOR
C-----------------------------------------------------------------------
C
      SUBROUTINE DENORM
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
      DO 10 I=1,NB
        Q = B2(I)
        B1(I) = B1(I)*Q
        B0(I) = B0(I)*Q
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   ALLO01
C ALLOCATION OF NUMERATORS AND DENOMINATORS ACCORDING TO THE
C ALLOCATION FIELD  ISEQ(I,J)
C-----------------------------------------------------------------------
C
      SUBROUTINE ALLO01
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
     *    SC0(16)
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
C
      NB = NBS
      FACT = SFACT
      DO 10 I=1,NB
        J = ISEQ(I,1)
        B2(I) = SB2(J)
        B1(I) = SB1(J)
        B0(I) = SB0(J)
        J = ISEQ(I,2)
        C1(I) = SC1(J)
        C0(I) = SC0(J)
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   ALLO02
C ALLOCATION OF NUMERATORS AND DENOMINATORS ACCORDING TO THE
C ALLOCATION FIELD  KSEQ(I,J), FOR A GIVEN START ORDERING
C-----------------------------------------------------------------------
C
      SUBROUTINE ALLO02
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
     *    SC0(16)
      COMMON /CREST2/ KSEQ(16,2)
C
      NB = NBS
      FACT = SFACT
      DO 10 I=1,NB
        J = KSEQ(I,1)
        B2(I) = SB2(J)
        B1(I) = SB1(J)
        B0(I) = SB0(J)
        J = KSEQ(I,2)
        C1(I) = SC1(J)
        C0(I) = SC0(J)
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   FIXPAR
C FIXED PAIRING OF POLES AND ZEROS
C USED TOGETHER WITH A FREE INPUT OF THE COEFFICIENTS
C-----------------------------------------------------------------------
C
      SUBROUTINE FIXPAR
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
     *    SC0(16)
      COMMON /SCRAT/ ADUM(32)
C
      DIMENSION IPOL(16,2)
      EQUIVALENCE (ADUM(1),IPOL(1,1))
C
      DO 10 I=1,NB
        IPOL(I,1) = 0
        IPOL(I,2) = 0
  10  CONTINUE
C
      DO 30 I=1,NB
        PM = 0.
        DO 20 J=1,NB
          IF (IPOL(J,1).NE.0) GO TO 20
          P2O = POW2O(1.0,0.,0.,C1(J),C0(J))
          IF (P2O.LE.PM) GO TO 20
          PM = P2O
          IPM = J
  20    CONTINUE
        IPOL(IPM,1) = I
  30  CONTINUE
C
      DO 70 I=1,NB
        DO 40 K=1,NB
          IF (IPOL(K,1).NE.I) GO TO 40
          KK = K
          GO TO 50
  40    CONTINUE
C
  50    CC0 = C0(KK)
        CC1 = C1(KK)
        DO 60 J=1,NB
          IF (IPOL(J,2).NE.0) GO TO 60
          P2O = POW2O(1.0,BB1,BB0,CC1,CC0)
          PM = P2O
          IZM = J
  60    CONTINUE
        IPOL(IZM,2) = KK
  70  CONTINUE
C
      CALL COPY01
      DO 80 I=1,NB
        IPM = IPOL(I,2)
        B2(IPM) = SB2(I)
        B1(IPM) = SB1(I)
        B0(IPM) = SB0(I)
  80  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   POLLOC
C CALCULATE THE POLES IN THE CLOSED UPPER Z-DOMAIN
C-----------------------------------------------------------------------
C
      SUBROUTINE POLLOC
C
      COMPLEX ZP, ZPS
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /CPOL/ ZP(16,2), ZPS(16,2)
C
      DO 20 I=1,NB
        Q = -C1(I)/2.
        QN = C0(I)
        QN = Q*Q - QN
        IF (QN.LT.0.) GO TO 10
        QN = SQRT(QN)
        ZPS(I,1) = CMPLX(Q+QN,0.)
        ZPS(I,2) = CMPLX(Q-QN,0.)
        GO TO 20
  10    QN = SQRT(-QN)
        ZPS(I,1) = CMPLX(Q,QN)
        ZPS(I,2) = CMPLX(Q,-QN)
  20  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION:     POW2O
C NOISE POWER OF A SECOND-ORDER BLOCK
C-----------------------------------------------------------------------
C
      FUNCTION POW2O(B2, B1, B0, C1, C0)
C
      BB0 = B2*B2 + B1*B1 + B0*B0
      BB1 = 2.*B1*(B2+B0)
      BB2 = 2.*B2*B0
      E1 = 1. + C0
      Q = BB0*E1 - BB1*C1 + BB2*(C1*C1-C0*E1)
      Q = Q/((1.-C0*C0)*E1-(C1-C1*C0)*C1)
      POW2O = Q
      RETURN
      END
C
C ======================================================================
C
C DOREDI - SUBROUTINES: PART IV
C
C ======================================================================
C
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DEFIN1
C DEFAULT VALUES FOR CLASS 01 INPUT DATA
C-----------------------------------------------------------------------
C
      SUBROUTINE DEFIN1
C
      COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /CGRID/ GR(64), NGR(12)
C
      NINP = 0
      NOUT = 3
      NDOUT = 0
      LSPOUT = 0
      NSPOUT = 0
      DO 10 I=1,12
        NGR(I) = 0
  10  CONTINUE
      DO 20 I=1,4
        ROM(I) = 0.
  20  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   INP010
C COMMAND CARDS FOR INPUT AND OUTPUT DIRECTIVES
C PROCESSING OF THE COMMAND CARDS CLASS 01
C-----------------------------------------------------------------------
C
      SUBROUTINE INP010
C
      COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
      COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
      COMMON /CGRID/ GR(64), NGR(12)
C
      J = JCODE/10
      IF (J.EQ.0) GO TO 110
      IF (J.GT.4) GO TO 110
      GO TO (10, 80, 90, 100), J
C
C *0110  B,NINP
C
  10  IF (IINP.GT.4) GO TO 120
      GO TO (20, 30, 40, 50), IINP
  20  IF (NINP.EQ.0) GO TO 130
      GO TO 70
  30  CALL INP012
      GO TO 70
  40  NINP = IINP
      CALL INP013
      GO TO 60
  50  NINP = IINP
      CALL INP014
  60  NGR(1) = 0
      RETURN
  70  NINP = IINP
      RETURN
C
C *0120  B,NOUT
C
  80  IF (IINP.LT.0 .OR. IINP.GT.5) GO TO 120
      NOUT = IINP
      RETURN
C
C *0130  NDOUT
C
  90  NDOUT = LINP
      RETURN
C
C *0140  LSPOUT,ISPOUT
C
 100  LSPOUT = LINP
      NSPOUT = IINP
      RETURN
C
 110  I = 1
      GO TO 140
 120  I = 4
      GO TO 140
 130  I = 12
 140  CALL ERROR(I)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   INP012
C INPUT OF THE FILTER DESCRIPTION AND DATA FROM DISK: (CHANNEL KA4)
C-----------------------------------------------------------------------
C
      SUBROUTINE INP012
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
      COMMON /CRECO/ JRCO, JECO, JDCO, JJDCO(5), JMAXV, JTRB2, LREF
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
C
9999  FORMAT (10I5)
9998  FORMAT (4E15.7)
C
      DO 10 I=1,3
        READ (KA4,9999)
  10  CONTINUE
      READ (KA4,9999) ITYP, IAPRO, NDEG
      READ (KA4,9998) SF, OM, ADELP, ADELS
      READ (KA4,9998) AC, ROM, RDELP, RDELS
      READ (KA4,9999) NZM, M
      DO 20 I=1,M
        READ (KA4,9998) (ZM(I,J),J=1,4)
  20  CONTINUE
      READ (KA4,9999)
      READ (KA4,9999) NB
      READ (KA4,9998) FACT
      DO 30 I=1,NB
        READ (KA4,9998) B2(I), B1(I), B0(I), C1(I), C0(I)
  30  CONTINUE
      READ (KA4,9999) IRCO
      DO 40 I=1,NB
        READ (KA4,9999) (IECO(I,J),J=1,5)
  40  CONTINUE
      DO 50 I=1,NB
        READ (KA4,9999) (IDCO(I,J),J=1,5)
  50  CONTINUE
      READ (KA4,9999) IWL, IECOM, JMAXV, JTRB2
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   INP013
C INPUT OF SECOND-ORDER BLOCKS
C-----------------------------------------------------------------------
C
      SUBROUTINE INP013
C
      COMMON /CARD/ LCODE, IZ(40)
      COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
      CALL INP001
      IF (LCODE.NE.0) RETURN
      CALL INP002
      IF (ICODE.NE.0 .OR. JCODE.NE.10) CALL ERROR(2)
      NB = IINP
      FACT = AINP(1)
C
      MA = 2*NB
      DO 20 I=1,MA
        CALL INP001
        IF (LCODE.NE.0) RETURN
        CALL INP002
        IF (ICODE.NE.0 .OR. JCODE.NE.10) CALL ERROR(2)
        IF (I.GT.NB) GO TO 10
        B2(I) = AINP(1)
        B1(I) = AINP(2)
        B0(I) = AINP(3)
        GO TO 20
  10    J = I - NB
        C1(J) = AINP(2)
        C0(J) = AINP(3)
  20  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   INP014
C INPUT OF POLES AND ZEROS
C-----------------------------------------------------------------------
C
      SUBROUTINE INP014
C
      COMMON /CARD/ LCODE, IZ(40)
      COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
      CALL INP001
      IF (LCODE.NE.0) RETURN
      CALL INP002
      IF (ICODE.NE.0 .OR. JCODE.NE.10) CALL ERROR(2)
      NB = IINP
      FACT = AINP(1)
      NN = 0
      JJ = NB
C
  10  CALL INP001
      IF (LCODE.NE.0) RETURN
      CALL INP002
      IF (ICODE.NE.0 .OR. JCODE.NE.10) CALL ERROR(2)
      QR = AINP(1)
      QI = AINP(2)
      IF (QI.EQ.0.) GO TO 20
      A0 = QR*QR + QI*QI
      A1 = -2.*QR
      GO TO 40
  20  IF (JJ.NE.NB) GO TO 30
      JJ = NB - 1
      R = QR
      GO TO 10
  30  JJ = NB
      A0 = R*QR
      A1 = -R - QR
  40  NN = NN + 1
      IF (K.EQ.2) GO TO 50
      B2(NN) = 1.
      B1(NN) = A1
      B0(NN) = A0
      GO TO 60
  50  C1(NN) = A1
      C0(NN) = A0
  60  IF (NN.LT.JJ) GO TO 10
      IF (JJ.EQ.NB) GO TO 70
      A1 = -R
      A0 = 0.
      JJ = JJ + 1
      GO TO 40
  70  IF (K.EQ.2) RETURN
      K = 2
      NN = 0
      GO TO 10
C
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DEFIN2
C DEFAULT VALUES FOR THE COMMAND CARDS OF CLASS 02
C-----------------------------------------------------------------------
C
      SUBROUTINE DEFIN2
C
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST2/ KSEQ(16,2)
      COMMON /COPTST/ LOPTS, ISTOR
C
      LOPTW = 0
      LSTAB = -1
      IWLR = 100
      ACXMI = 0.
      ACXMA = 1.
      LOPTS = 0
      ISTOR = 100
      LSEQ = 0
      ISCAL = 2
      SCALM = 1.
      ITERM = 10
      ITERM1 = 4
      DO 20 J=1,2
        DO 10 I=1,MBL
          KSEQ(I,J) = I
          ISEQ(I,J) = I
  10    CONTINUE
  20  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   INP020
C COMMAND CARDS FOR CONTROLLING THE OPTIMIZATION
C PROCESSING OF THE COMMAND CARDS OF CLASS 02
C-----------------------------------------------------------------------
C
      SUBROUTINE INP020
C
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
      COMMON /INTDAT/ JINP(16)
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST2/ KSEQ(16,2)
      COMMON /COPTST/ LOPTS, ISTOR
C
      IF (JCODE.LT.10) GO TO 110
      J = JCODE/10
      IF (J.GT.4) GO TO 110
      GO TO (10, 50, 60, 70), J
C
C *0210  LWLF,ITERM,ACX,ACXMI,ACXMA
C
  10  LOPTW = -1
      IF (LINP.EQ.(-1)) LOPTW = 1
      IF (IINP.EQ.0) GO TO 20
      IF (IINP.LT.0) GO TO 120
      ITERM = IINP
  20  IF (NPAR.EQ.0) GO TO 30
      Q = AINP(1)
      IF (Q.LT.0. .OR. Q.GT.1.) GO TO 120
      ACX = Q
  30  Q = AINP(2)
      IF (Q.EQ.0.) GO TO 40
      IF (Q.LT.0. .OR. Q.GE.1.) GO TO 120
      ACXMI = Q
  40  Q = AINP(3)
      IF (Q.EQ.0.) RETURN
      IF (Q.LE.0. .OR. Q.GT.1.) GO TO 120
      ACXMA = Q
      RETURN
C
C *0220  B,ITERM1
C
  50  IF (IINP.LT.0) GO TO 120
      ITERM1 = IINP
      RETURN
C
C *0230  PAIRF,ISTOR
C
  60  LOPTS = -1
      IF (LINP.EQ.(-1)) LOPTS = 1
      IF (IINP.GT.100) GO TO 120
      IF (IINP.LT.0) GO TO 120
      IF (IINP.NE.0) ISTOR = IINP
      RETURN
C
C *0240  LSEQ,ISCAL,SCALM
C
  70  IF (LINP.NE.(-1)) GO TO 100
      LSEQ = -1
      DO 90 J=1,2
        CALL INP001
        IF (ICODE.NE.0) GO TO 130
        IF (JCODE.NE.40) GO TO 130
        CALL INP003
        DO 80 I=1,MBL
          KSEQ(I,J) = JINP(I)
  80    CONTINUE
  90  CONTINUE
C
 100  IF (IINP.NE.0) ISCAL = IINP
      IF (ISCAL.LT.0) ISCAL = 0
      IF (AINP(1).NE.0.) SCALM = AINP(1)
      RETURN
C
 110  I = 4
      GO TO 140
 120  I = 5
      GO TO 140
 130  I = 2
 140  CALL ERROR(I)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DEFIN3
C SET DEFAULT VALUES
C-----------------------------------------------------------------------
C
      SUBROUTINE DEFIN3
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
      ITYP = 1
      IAPRO = 4
      NDEG = 0
      NDEGF = 0
      MDEG = 0
      NORMA = 0
      LSOUT = 0
C
      DO 10 I=1,4
        OM(I) = 0.
  10  CONTINUE
      SF = 0.
      A = 0.
      ADELP = 0.
      ADELS = 1./FLMA
      VSN = 0.
      VD = 0.
      ACX = 1000.
      EDEG = ACX
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   INP030
C INPUT ROUTINE FOR THE FILTER TYPE AND FOR THE TOLERANCE SCHEME
C PROCESSING OF THE COMMAND CARDS OF CLASS 03
C-----------------------------------------------------------------------
C
      SUBROUTINE INP030
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
      SOM(AA) = 2.*ATAN(AA)
C
      J = JCODE/10
      IF (J.GT.7) GO TO 240
      GO TO (10, 20, 30, 40, 170, 220, 230), J
C
C *0310  B,FILTER TYPE,SAMPLING FREQUENCY
C
  10  IF (IINP.NE.0) ITYP = IINP
      IF (AINP(1).NE.0.) SF = AINP(1)
      RETURN
C
C *0320  B,IAPRO,,ACX
C
  20  IF (IINP.NE.0) IAPRO = IINP
      IF (NPAR.NE.0) ACX = AINP(1)
      RETURN
C
C *0330  NDEGF,EDEG
C
  30  IF (LINP.EQ.(-1)) NDEGF = -1
      IF (IINP.NE.0) NDEG = IINP
      IF (NPAR.NE.0) EDEG = AINP(1)
      RETURN
C
C CUTOFF FREQUENCIES OF THE TOLERANCE SCHEME
C
  40  K = JCODE - J*10
      IF (K.GT.7) GO TO 250
      GO TO (50, 50, 100, 110, 120, 140, 160), K
C
C *0341  B,B,FR(1),FR(2),    ::  FREQUENCY DOMAIN
C *0342  B,B,FR(3),FR(4),
C
  50  IF (AINP(3).EQ.0.) GO TO 70
      IF (SF.NE.0.) GO TO 260
  60  SF = AINP(3)
  70  IF (SF.EQ.0.) GO TO 270
      Q = 2.*PI/SF
      GO TO (80, 90), K
  80  OM(1) = AINP(1)*Q
      OM(2) = AINP(2)*Q
      RETURN
C
  90  OM(3) = AINP(1)*Q
      OM(4) = AINP(2)*Q
      RETURN
C
C *0343  B,B,OM(1),OM(2)        ::  Z-DOMAIN
C
 100  OM(1) = AINP(1)
      OM(2) = AINP(2)
      RETURN
C
C *0344  B,B,OM(3),OM(4)       ::  Z-DOMAIN SEC. CARD
C
 110  OM(3) = AINP(1)
      OM(4) = AINP(2)
      RETURN
C
C *0345  B,B,S(1),S(2)         ::  S-DOMAIN
C
 120  DO 130 I=1,2
        OM(I) = SOM(AINP(I))
 130  CONTINUE
      RETURN
C
C *0346  B,B,S(3),S(4)         :: S-DOMAIN  SEC. CARD
C
 140  DO 150 I=1,2
        OM(I+2) = SOM(AINP(I))
 150  CONTINUE
      RETURN
C
C *0347  B,B,VSN,VD,A           ::  NORMALIZED PARAMETER
C
 160  VSN = AINP(1)
      VD = AINP(2)
      A = AINP(3)
      RETURN
C
C DELTA P, DELTA S
C
 170  K = JCODE - J*10
      IF (K.GT.3) GO TO 250
      GO TO (180, 190, 210), K
C
C *0351  B,B,DELP,DELS
C
 180  ADELP = AINP(1)
      ADELS = AINP(2)
      RETURN
C
C *0352  B,B,AP,AS (IN DB)
C
 190  ADELP = 1. - 10.**(-0.05*AINP(1))
 200  ADELS = 10.**(-0.05*AINP(2))
      RETURN
C
C *0353  B,B,P,AS    (AS IN DB)
C
 210  Q = AINP(1)
      ADELP = 1. - SQRT(1.-Q*Q)
      GO TO 200
C
C *0360  B,NORMA
C
 220  IF (IINP.LT.0 .OR. IINP.GT.3) GO TO 250
      NORMA = IINP
      RETURN
C
C *0370  LSOUT
C
 230  LSOUT = LINP
      RETURN
C
C OUTPUT OF ERROR MESSAGES
C
 240  I = 4
      GO TO 280
 250  I = 5
      GO TO 280
 260  I = 7
      GO TO 280
 270  I = 6
 280  CALL ERROR(I)
      IF (I.EQ.7) GO TO 60
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   INP031
C CHECK ROUTINE FOR THE FILTER INPUT DATA
C CHECK DATA CARDS OF CLASS 03
C-----------------------------------------------------------------------
C
      SUBROUTINE INP031
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
      IF (ITYP.GT.4) GO TO 90
      IF (IAPRO.GT.4) GO TO 100
C
      IF (EDEG.GT.999.) EDEG = 0.2
C
C CHECK OM(1) ... OM(4)
C
      LVSN = 0
      LSYM = 0
      J = 0
      Q = 0.
      ME = 2
      IF (ITYP.GE.3) ME = 4
      DO 20 I=1,ME
        QQ = OM(I)
        IF (QQ.EQ.0.) GO TO 10
        IF (QQ.GT.PI) GO TO 120
        IF (QQ.LT.Q) GO TO 110
        Q = QQ
        GO TO 20
  10    LSYM = LVSN
        LVSN = I
        J = J + 1
  20  CONTINUE
      IF (J.EQ.0) GO TO 80
      IF (ITYP.GT.2) GO TO 30
      IF (J.EQ.1) GO TO 70
      GO TO 50
  30  GO TO (40, 60, 50, 50), J
  40  LSYM = LVSN
      LVSN = 0
      J = 0
      GO TO 80
C
  50  IF (VSN.EQ.0.) GO TO 120
      IF (VD.EQ.0.) GO TO 120
      IF (ITYP.GT.2 .AND. A.EQ.0.) GO TO 110
      LVSN = -1
      GO TO 80
C
  60  J = 1
      IQ = LSYM + LVSN
      IF (IQ.LE.3 .OR. IQ.GE.7) GO TO 120
  70  IF (NDEG.EQ.0) GO TO 120
C
C CHECK FOR PASSBAND AND STOPBAND RIPPLE
C
  80  IF (ADELP.EQ.0.) J = J + 1
      IF (ADELS.LE.(1./FLMA)) J = J + 1
      IF (J.EQ.0) RETURN
      IF (J.NE.1 .OR. NDEG.EQ.0) GO TO 130
      RETURN
C
C OUTPUT OF ERROR MESSAGE
C
  90  I = 20
      GO TO 140
 100  I = 21
      GO TO 140
 110  I = 22
      GO TO 140
 120  I = 23
      GO TO 140
 130  I = 24
 140  CALL ERROR(I)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DEFIN4
C SET DEFAULT VALUES FOR COMMAND CARDS OF CLASS 04
C-----------------------------------------------------------------------
C
      SUBROUTINE DEFIN4
C
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /CRECO/ JRCO, JECO, JDCO, JJDCO(5), JMAXV, JTRB2, LREF
      COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
C
      LSTAB = 0
      LREF = 0
      IWL = 16
      JRCO = 2
      JECO = 0
      JDCO = 0
      JMAXV = -3
      JTRB2 = 3
      DO 20 J=1,5
        IRCO(J) = 0
        DO 10 I=1,MBL
          IECO(I,J) = 0
          IDCO(I,J) = -100
  10    CONTINUE
  20  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   INP040
C INPUT ROUTINE FOR THE SPECIFICATION OF THE COEFFICIENT REALIZATION
C-----------------------------------------------------------------------
C
      SUBROUTINE INP040
C
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /CARD/ LCODE, IZ(40)
      COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
      COMMON /INTDAT/ JINP(16)
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /CRECO/ JRCO, JECO, JDCO, JJDCO(5), JMAXV, JTRB2, LREF
      COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
C
      J = JCODE/10
      IF (J.EQ.0) GO TO 170
      IF (J.GT.6) GO TO 170
      GO TO (10, 20, 50, 90, 150, 160), J
C
C *0410  LSTAB,IWL
C
  10  LSTAB = LINP
      IF (IINP.NE.0) IWL = IINP
      RETURN
C
C *0420  RCOINP,JRCO
C
  20  IF (LINP.NE.(-1)) GO TO 40
      CALL INP001
      IF (LCODE.EQ.1) RETURN
      CALL INP003
      IF (ICODE.NE.0 .OR. JCODE.NE.20) CALL ERROR(2)
      DO 30 I=1,5
        IRCO(I) = JINP(I)
  30  CONTINUE
      JRCO = 0
      IF (IINP.NE.0) GO TO 180
  40  JRCO = IINP
      RETURN
C
C *0430  ECOINP,JECO
C
  50  IF (LINP.NE.(-1)) GO TO 80
      JECO = 0
      IECOM = 0
      DO 70 J=1,5
        CALL INP001
        IF (LCODE.EQ.1) RETURN
        CALL INP003
        IF (ICODE.NE.0 .OR. JCODE.NE.30) CALL ERROR(2)
        DO 60 I=1,MBL
          IECOM = MAX0(IECOM,JINP(I))
          IECO(I,J) = JINP(I)
  60    CONTINUE
  70  CONTINUE
      RETURN
C
  80  JECO = IINP
      RETURN
C
C *0440  DCOINP,JDCO
C
  90  IF (LINP.EQ.(-1)) GO TO 120
      DO 110 J=1,5
        II = IINP/10
        M = 6 - J
        JJ = IINP - 10*II
        IINP = II
        JJDCO(M) = JJ
        IF (II.NE.0) GO TO 110
        DO 100 I=1,16
          IDCO(I,J) = -100
 100    CONTINUE
 110  CONTINUE
      JDCO = -1
      RETURN
C
 120  JDCO = 0
      DO 140 J=1,5
        CALL INP001
        IF (LCODE.EQ.1) RETURN
        CALL INP003
        IF (ICODE.NE.0 .OR. JCODE.NE.40) CALL ERROR(2)
        DO 130 I=1,16
          IDCO(I,J) = JINP(I)
 130    CONTINUE
 140  CONTINUE
      RETURN
C
C *0450  B,JMAXV
C
 150  JMAXV = IINP
      RETURN
C
C *0460  B,JTRB2
C
 160  JTRB2 = IINP
      RETURN
C
C OUTPUT OF ERROR MESSAGES
C
 170  I = 1
      GO TO 190
 180  I = 9
 190  CALL ERROR(I)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DEFIN5
C DEFAULT VALUES FOR COMMAND CARDS OF CLASS 05
C-----------------------------------------------------------------------
C
      SUBROUTINE DEFIN5
C
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /CRENO/ LCNO, ICNO(16,5)
C
      LCNO = 0
      LPOT2 = 0
      ISTRU = 13
      JSTRU = 1
      JSTRUS = 1
      JSTRUD = 1
      IWLR = 100
      DO 20 J=1,5
        DO 10 I=1,MBL
          ICNO(I,J) = 3
  10    CONTINUE
  20  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   INP050
C INPUT DATA FOR THE STRUCTURE REALIZATION
C COMMAND CARDS OF CLASS 05
C-----------------------------------------------------------------------
C
      SUBROUTINE INP050
C
      COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
      COMMON /CARD/ LCODE, IZ(40)
      COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
      COMMON /INTDAT/ JINP(16)
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CREST1/ JSTRUS, JSTRUD
      COMMON /CRENO/ LCNO, ICNO(16,5)
C
      J = JCODE/10
      IF (J.EQ.0) GO TO 100
      IF (J.GT.3) GO TO 100
      GO TO (10, 20, 90), J
C
C *0510  LPOT2,ISTRU
C
  10  LPOT2 = LINP
      IF (IINP.EQ.0) RETURN
      IF (IINP.GE.30) GO TO 110
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C HERE INSERT THE PARAMETERS JSTRU,JSTRUS,JSTRUD FOR YOUR STRUCTURES
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
      ISTRU = IINP
      JSTRU = 1
      IF (ISTRU/2.EQ.(ISTRU+1)/2) JSTRU = 2
      JSTRUS = JSTRU
      JSTRUD = 1
      IF (ISTRU.GE.20) JSTRUD = 2
      RETURN
C
C *0520  CNOINP,JCNO
C
  20  IF (LINP.NE.(-1)) GO TO 50
      DO 40 J=1,5
        CALL INP001
        IF (LCODE.NE.0) RETURN
        CALL INP003
        IF (ICODE.NE.0 .OR. JCODE.NE.20) CALL ERROR(2)
        DO 30 I=1,MBL
          ICNO(I,J) = JINP(I)
  30    CONTINUE
  40  CONTINUE
      IF (IINP.NE.0) GO TO 130
      RETURN
C
  50  IF (IINP.EQ.0) GO TO 120
      DO 80 J=1,5
        II = IINP/10
        M = 6 - J
        JJ = IINP - 10*II
        IINP = II
        IF (JJ.NE.0) GO TO 60
        IF (J.EQ.1) GO TO 130
        JJ = ICNO(1,M+1)
  60    DO 70 I=1,16
          ICNO(I,M) = JJ
  70    CONTINUE
  80  CONTINUE
      RETURN
C
C *0530  LCNO,IWLR
C
  90  LCNO = LINP
      IF (IINP.LT.0) GO TO 110
      IF (IINP.NE.0) IWLR = IINP
      RETURN
C
 100  I = 1
      GO TO 140
 110  I = 5
      GO TO 140
 120  I = 9
      GO TO 140
 130  I = 4
 140  CALL ERROR(I)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DEFIN6
C DEFAULT VALUES FOR THE COMMAND CARDS OF CLASS 06
C-----------------------------------------------------------------------
C
      SUBROUTINE DEFIN6
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /COPTST/ LOPTS, ISTOR
      COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
C
      LOPTS = 0
      LOPTW = 0
      IPOT = 0
      LNOR = 0
      IPAG = 2
      OMLO = 0.
      OMUP = PI
      RMAX = 1.
      RMIN = -1.
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   INP060
C COMMAND CARDS FOR THE ANALYSIS
C PROCESSING OF THE COMMAND CARDS OF CLASS 06
C-----------------------------------------------------------------------
C
      SUBROUTINE INP060
C
      COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
      COMMON /COPTST/ LOPTS, ISTOR
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
      COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
C
      J = JCODE/10
      IF (J.EQ.0) GO TO 90
      IF (J.GT.6) GO TO 90
      GO TO (10, 20, 30, 60, 70, 80), J
C
C *0610  LWLM,IWL
C
  10  LOPTW = 1
      IF (LINP.EQ.(-1)) LOPTW = -1
      IF (IINP.LT.0) GO TO 100
      IF (IINP.NE.0) IWL = IINP
      RETURN
C
C *0620  LSEQ,LSCAL,SCALM
C
  20  LOPTS = -1
      JCODE = 40
      CALL INP020
      RETURN
C
C *0630  LNOR,IPAG,OMLO,OMUP,RMAX
C
  30  IPLOT = 1
  40  LNOR = LINP
      IF (IINP.GT.0) IPAG = IINP
      IF (IPLOT.EQ.4) RETURN
      IF (NPAR.NE.0) OMLO = AINP(1)
      Q = AINP(2)
      IF (Q.EQ.0.) GO TO 50
      IF (Q.LE.OMLO) GO TO 110
      OMUP = Q
  50  IF (AINP(3).NE.0.) RMAX = AINP(3)
      RETURN
C
C *0640  LNOR,IPAG,OMLO,OMUP,RMAX
C
  60  IPLOT = 2
      RMAX = 100.
      GO TO 40
C
C *0650  LNOR,IPAG,OMLO,OMUP
C
  70  IPLOT = 3
      GO TO 40
C
C *0660  LNOR,IPAG,RMIN,RMAX
C
  80  IPLOT = 4
      RMAX = 1.
      IF (NPAR.NE.0) RMIN = AINP(1)
      Q = AINP(2)
      IF (Q.EQ.0.) GO TO 40
      IF (Q.LE.RMIN) GO TO 110
      RMAX = Q
      GO TO 40
C
  90  I = 1
      GO TO 120
 100  I = 5
      GO TO 120
 110  I = 10
      IPLOT = 0
 120  CALL ERROR(I)
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   INP001
C READ COMMAND CARD
C
C IF 'DECODE' DOES NOT CONFORM TO YOUR COMPILER,
C SET THE ALTERNATIVE INPUT VERSION
C
C ELIMINATE THE STATEMENTS ENCLOSED BY THE COMMENT LINES C111
C
C AND REPLACE THE COMMENT MARKS 'C222' WITH '    '
C
C-----------------------------------------------------------------------
C
      SUBROUTINE INP001
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CARD/ LCODE, IZ(40)
      COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
C111
      COMMON /CIOFOR/ IOFO(3), ION
C111
C
      DIMENSION IZ1(2)
      DATA IPR, ISL /1H*,1H//
C
9999  FORMAT (1X, 4A1)
9998  FORMAT (1X//11H DATA INPUT)
C
      IF (LCODE.EQ.(-1)) WRITE (KA2,9998)
C
  10  CONTINUE
C111
      READ (KA1,IOFO) (IZ(I),I=1,ION)
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C DECODE IS A MACHINE-DEPENDENT FUNCTION
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
      DECODE(2,21,IZ) IZ1(1),IZ1(2)
C111
C
C222  READ (KA1,21) (IZ1(I),I=1,2)
  21  FORMAT (2A1)
      IF (IZ1(1).EQ.IPR) GO TO 30
      IF (IZ1(1).EQ.ISL .AND. IZ1(2).EQ.ISL) GO TO 20
      CALL ERROR(1)
      GO TO 10
  20  CONTINUE
C111
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C DECODE IS A MACHINE-DEPENDENT FUNCTION
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
      DECODE (4,22,IZ) ICODE,JCODE
C111
C
C222  READ (KA1,22) ICODE,JCODE
  22  FORMAT (2X, 2A1)
      WRITE (KA2,9999) IZ1(1), IZ1(2), ICODE, JCODE
      LCODE = 1
      RETURN
  30  LCODE = 0
      ICODE = 0
      RETURN
C
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   INP002
C DECODE COMMAND CARD
C
C IF 'DECODE' DOES NOT CONFORM TO YOUR COMPILER, SEE INP001
C
C-----------------------------------------------------------------------
C
      SUBROUTINE INP002
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CARD/ LCODE, IZ(40)
      COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
      COMMON /SCRAT/ ADUM(32)
      DIMENSION IZZ(15)
      EQUIVALENCE (ADUM(1),IZZ(1))
C
      DATA IBL, IT /1H ,1HT/
C
      NPAR = 0
C111
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C DECODE IS A MACHINE-DEPENDENT FUNCTION
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
      DECODE (30,11,IZ) IZZ
  11  FORMAT (15X, 15A1)
      DO 10 I=1,15
        IF (IZZ(I).NE.IBL) NPAR = 1
  10  CONTINUE
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C DECODE IS A MACHINE-DEPENDENT FUNCTION
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
      DECODE (60,12,IZ) ICODE,JCODE,JINP,IINP,AINP
C111
C
C222  READ (KA1,12) ICODE,JCODE,JINP,IINP,(AINP(I),I=1,3)
C222  NPAR=1
C
      WRITE (KA2,9999) ICODE, JCODE, JINP, IINP, AINP
C
  12  FORMAT (1X, 2I2, 2X, 1A1, I6, 1X, 3E15.0)
9999  FORMAT (2H *, 2I4, 2X, 1A1, I6, 1X, 3E15.6)
C
      LINP = 0
      IF (JINP.EQ.IT) LINP = -1
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   INP003
C INTEGER DATA INPUT
C
C
C IF 'DECODE' DOES NOT CONFORM TO YOUR COMPILER, SEE INP001
C-----------------------------------------------------------------------
C
      SUBROUTINE INP003
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CARD/ LCODE, IZ(40)
      COMMON /INTDAT/ JINP(16)
C
  12  FORMAT (1X, 2I2, 3X, 16I4)
9999  FORMAT (1H*, 2I3, 1X, 16I4)
C111
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C DECODE IS A MACHINE-DEPENDENT FUNCTION
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
      DECODE (72,12,IZ) ICODE,JCODE,(JINP(I),I=1,16)
C111
C
C222  READ (KA1,12) ICODE,JCODE,(JINP(I),I=1,16)
C
      WRITE (KA2,9999) ICODE, JCODE, JINP
      RETURN
      END
C
C ======================================================================
C
C DOREDI - SUBROUTINES: PART V
C
C ======================================================================
C
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT010
C LEADER TO OUT011 FOR STRUCTURES WITH SEPARATE SCALING
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT010
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
      WRITE (KA2,9999)
9999  FORMAT (1X//46H STRUCTURE WITH SEPARATE SCALING FACTORS B2(L)//
     *    45H FOR STRUCTURES 21 TO 28  B2(0) = GAIN FACTOR)
C
      DO 10 I=1,NB
        Q = B2(I)
        B1(I) = B1(I)/Q
        B0(I) = B0(I)/Q
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT011
C OUTPUT OF BLOCKS OF SECOND ORDER
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT011
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
      WRITE (KA2,9999) FACT
9999  FORMAT (//23H BLOCKS OF SECOND ORDER//24H CONSTANT GAIN FACTOR = ,
     *    E15.7//3H  L, 5X, 5HB2(L), 8X, 5HB1(L), 8X, 5HB0(L), 10X,
     *    5HC1(L), 8X, 5HC0(L)/)
      WRITE (KA2,9998) (I,B2(I),B1(I),B0(I),C1(I),C0(I),I=1,NB)
9998  FORMAT (1X, I3, 3F13.8, 2X, 2F13.8)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT012
C OUTPUT OF THE RESULTS TO THE DISK: (CHANNEL KA4)
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT012
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
      COMMON /CRECO/ JRCO, JECO, JDCO, JJDCO(5), JMAXV, JTRB2, LREF
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
C
9999  FORMAT (10I5)
9998  FORMAT (4E15.7)
C
      WRITE (KA4,9997)
9997  FORMAT (19H FILTER DESCRIPTION)
      WRITE (KA4,9999) ITYP, IAPRO, NDEG
      WRITE (KA4,9998) SF, OM, ADELP, ADELS
      WRITE (KA4,9998) AC, ROM, RDELP, RDELS
      M = 0
      DO 10 I=1,4
        M = MAX0(M,NZM(I))
  10  CONTINUE
      WRITE (KA4,9999) NZM, M
      DO 20 I=1,M
        WRITE (KA4,9998) (ZM(I,J),J=1,4)
  20  CONTINUE
      WRITE (KA4,9996)
9996  FORMAT (5H DATA)
      WRITE (KA4,9999) NB
      WRITE (KA4,9998) FACT
      DO 30 I=1,NB
        WRITE (KA4,9998) B2(I), B1(I), B0(I), C1(I), C0(I)
  30  CONTINUE
      WRITE (KA4,9999) IRCO
      DO 40 I=1,NB
        WRITE (KA4,9999) (IECO(I,J),J=1,5)
  40  CONTINUE
      DO 50 I=1,NB
        WRITE (KA4,9999) (IDCO(I,J),J=1,5)
  50  CONTINUE
      WRITE (KA4,9999) IWL, IECOM, JMAXV, JTRB2
      WRITE (KA4,9995)
9995  FORMAT (12H END OF DATA)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT016
C OUTPUT OF BLOCKS OF SECOND ORDER IN THE S-DOMAIN
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT016
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
      COMMON /SCRAT/ ADUM(32)
      DIMENSION NZE(16)
      EQUIVALENCE (ADUM(1),NZE(1))
C
      WRITE (KA2,9999) SFA
9999  FORMAT (//39H BLOCKS OF SECOND ORDER IN THE S-DOMAIN//9H CONSTANT,
     *    15H GAIN FACTOR = , E15.7//3H  L, 5X, 5HA2(L), 8X, 5HA1(L),
     *    8X, 5HA0(L), 10X, 5HD1(L), 8X, 5HD0(L)/)
C
      N = 0
      ME = NZM(4)
      DO 10 I=1,ME
        NZE(I) = NZERO(I)
  10  CONTINUE
      NB = (NDEG+1)/2
      II = 1
      NZ = NZE(II)
C
      DO 180 I=1,NB
  20    IF (NZ.GT.0) GO TO 30
        II = II + 1
        NZ = NZE(II)
        GO TO 20
  30    QI = SM(II,4)
        IF (NZ.EQ.1) GO TO 70
        IF (QI.GE.FLMA) GO TO 50
  40    B2 = 1.
        B1 = 0.
        B0 = QI*QI
        GO TO 60
  50    B2 = 0.
        B1 = 0.
        B0 = 1.
  60    NZ = NZ - 2
        GO TO 130
  70    IF (II.EQ.ME) GO TO 120
        MA = II + 1
        DO 80 J=MA,ME
          IF (SM(J,4).GE.FLMA) GO TO 90
          IF (SM(J,4).LE.FLMI) GO TO 100
  80    CONTINUE
        GO TO 120
  90    NZE(J) = NZE(J) - 1
        IF (QI.GE.FLMA) GO TO 50
        GO TO 110
 100    NZE(J) = NZE(J) - 1
        IF (QI.LE.FLMI) GO TO 40
 110    B2 = 0.
        B1 = 1.
        B0 = 0.
        NZ = NZ - 1
        GO TO 130
 120    IF (QI.GE.FLMA) GO TO 50
        GO TO 110
C
 130    N = N + 1
        QR = SPR(N)
        QI = SPI(N)
        IF (ABS(QI).LT.FLMI) GO TO 140
        C1 = -2.*QR
        C0 = QR*QR + QI*QI
        GO TO 170
 140    IF (N.GE.NJ) GO TO 150
        IF (ABS(SPI(N+1)).LT.FLMI) GO TO 160
 150    C1 = -QR
        C0 = 0.
        GO TO 170
 160    N = N + 1
        QI = SPR(N)
        C1 = -QR - QI
        C0 = QR*QI
C
 170    WRITE (KA2,9998) I, B2, B1, B0, C1, C0
9998    FORMAT (1X, I3, 3F13.8, 2X, 2F13.8)
 180  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT022
C OUTPUT OF THE TERMINATION KIND
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT022
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
      COMMON /COPTCP/ IS1, IS2, IABO, IWLGS, IWLGN, IWLGP, ADEPSS,
     *    ADEPSN, ADEPSP, ACXS, ACXN, ACXP
C
      WRITE (KA2,9999) ITER
9999  FORMAT (///38H OPTIMIZATION IS TERMINATED AFTER THE , I4,
     *    6H. STEP)
      GO TO (10, 20, 30), IABO
  10  WRITE (KA2,9998)
9998  FORMAT (25H THREE UNSUCCESSFUL STEPS)
      GO TO 40
  20  WRITE (KA2,9997)
9997  FORMAT (24H MAXIMUM NUMBER OF STEPS)
      GO TO 40
  30  WRITE (KA2,9996)
9996  FORMAT (54H OPTIMIZATION PARAMETER IS COMING TO THE END OF THE RE,
     *    4HGION)
  40  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT023
C OUTPUT OF THE ITERATION NUMBER
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT023
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
C
      WRITE (KA2,9999) ITER
9999  FORMAT (1X////1X, I3, 13H-TH ITERATION)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT025
C OUTPUT OF THE OPTIMIZATION OPTIONS
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT025
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /COPTST/ LOPTS, ISTOR
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /CRENO/ LCNO, ICNO(16,5)
C
      DIMENSION IES(4), NO(4)
      DATA IES(1), IES(2), IES(3), IES(4) /1H ,1HY,1HE,1HS/
      DATA NO(1), NO(2), NO(3), NO(4) /1H ,1H ,1HN,1HO/
C
      IF (LOPTS.EQ.0) GO TO 10
      WRITE (KA2,9999)
9999  FORMAT (1X//41H OPTIMIZATION OF THE PAIRING AND ORDERING//)
      WRITE (KA2,9998) ISTOR
9998  FORMAT (46H STORAGE FOR THE INTERMEDIATE RESULTS  ISTOR =, I4/)
      IF (LOPTS.EQ.(-1)) WRITE (KA2,9997) NO
      IF (LOPTS.NE.(-1)) WRITE (KA2,9997) IES
9997  FORMAT (/14H FIXED PAIRING, 32X, 4A1)
C
  10  WRITE (KA2,9996) ISTRU
9996  FORMAT (1X/20H REALIZED STRUCTURE , 19X, 7HISTRU =, I4/)
      WRITE (KA2,9995) ISCAL
9995  FORMAT (/16H SCALING  OPTION, 23X, 7HISCAL =, I4)
      IF (LPOT2.EQ.(-1)) WRITE (KA2,9994) IES
      IF (LPOT2.NE.(-1)) WRITE (KA2,9994) NO
9994  FORMAT (10X, 24HBY A FACTOR OF POWER TWO, 12X, 4A1)
      WRITE (KA2,9993) SCALM
9993  FORMAT (10X, 21HCHOSEN MAXIMUM OF THE/10X, 15HOVERFLOW POINTS,
     *    14X, 7HSCALM =, F6.3/)
      WRITE (KA2,9992)
9992  FORMAT (/33H REALIZATION OF THE COEFFICIENTS )
      IF (LCNO.EQ.(-1)) WRITE (KA2,9991) NO, IWLR
      IF (LCNO.NE.(-1)) WRITE (KA2,9991) IES, IWLR
9991  FORMAT (5X, 33HCONSIDERATION OF THE QUANTIZATION, 8X, 4A1/5X,
     *    10HWORDLENGTH, 24X, 7HIWLR  =, I4/24X, 19H(FOR IWLR=100 NO RO,
     *    7HUNDING)/)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT026
C OUTPUT OF THE NOISE VALUES OF THE SECOND-ORDER BLOCKS
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT026
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
     *    LPOT2, JSTRU
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
      COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
      COMMON /CRENO/ LCNO, ICNO(16,5)
      COMMON /CPOW/ PNU, PNC, AND, ITCORP
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
      DIMENSION NCHA(5), NQUAN(5)
C
      IF (IB.LE.0) GO TO 50
      IF (IB.EQ.1) WRITE (KA2,9999)
9999  FORMAT (1X//47H NUM  DEN   SCALING   UNCOR. NOISE   COR. NOISE,
     *    14H     DC-OFFSET, 12H CHA CO-QUAN)
      ALOG2 = ALOG(2.)
      IF (IB.GT.NB) GO TO 10
      II = IB
      IN = ISEQ(IB,1)
      ID = ISEQ(IB,2)
      GO TO 20
  10  IN = 0
      ID = 0
      II = NB
  20  DO 40 J=1,5
        Q = AQC(J)
        IF (Q.EQ.1.) GO TO 30
        Q = ALOG(Q)/ALOG2
        NQUAN(J) = -IFIX(Q-0.5)
        NCHA(J) = ICNO(II,J)
        GO TO 40
  30    NQUAN(J) = 0
        NCHA(J) = 0
  40  CONTINUE
      WRITE (KA2,9998) IN, ID, SCA, PNU, PNC, AND, (NCHA(J),NQUAN(J),J=
     *    1,5)
9998  FORMAT (/1X, 2(I3, 2X), F9.6, 1X, 3E13.5, 2I4, 4H BIT/(60X, 2I4,
     *    4H BIT))
      RETURN
C
  50  WRITE (KA2,9997) AND
9997  FORMAT (46X, 14(1H-)/46X, E14.5/)
      Q = AND*AND
      PNT = PNU + PNC + Q
      WRITE (KA2,9996) PNU, PNC, Q, PNT
9996  FORMAT (12H TOTAL NOISE, 9X, 3E13.5, 1H=, E11.4)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT027
C OUTPUT OF THE DIFFERENT NOISE VALUES
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT027
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CNOISE/ RI, RIN, RE, REN, FAC
C
      ANP = RI/12.
      ALANP = 10.*ALOG10(ANP)
      WRITE (KA2,9999) ANP, ALANP
9999  FORMAT (1H //22H ABSOLUTE OUTPUT NOISE, 7X, 3HANP, 3X, 1H=,
     *    E13.5, 12H  (Q*Q)    =, F6.2, 3H DB/)
      WRITE (KA2,9998) FAC
9998  FORMAT (/22H SCALING AT THE OUTPUT, 7X, 3HSCA, 3X, 1H=, E13.5/)
      RNP = RIN/12.
      ALRNP = 10.*ALOG10(RNP)
      WRITE (KA2,9997) RNP, ALRNP
9997  FORMAT (/22H RELATIVE OUTPUT NOISE, 7X, 3HRNP, 3X, 1H=, E13.5,
     *    12H  (Q*Q)    ,F6.2, 3H DB)
      WRITE (KA2,9995)
      ALRIN = 10.*ALOG10(RIN)
      WRITE (KA2,9996) RIN, ALRIN
9996  FORMAT (/19H INNER NOISE FIGURE, 10X, 3HRIN, 3X, 1H=, E13.5,
     *    12H (Q*Q)/12 = , F6.2, 3H DB)
      WRITE (KA2,9995)
9995  FORMAT (/23H RELATED TO MAX /H/ = 1/)
      WRITE (KA2,9994) REN
9994  FORMAT (/22H ENTRANCE NOISE FIGURE, 7X, 3HREN, 3X, 1H=, E13.5/)
      WRITE (KA2,9995)
      ADW = RIN/(1.+REN)
      ADW = ALOG(ADW)/(2.*ALOG(2.))
      IDW = IFIX(ADW)
      Q = FLOAT(IDW)
      IF (Q.LT.ADW) IDW = IDW + 1
      WRITE (KA2,9993) ADW, IDW
9993  FORMAT (/22H ADDITIONAL WORDLENGTH, 6X, 11HDELTA W = /, F7.3,
     *    4H/ = , I4/)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT030
C OUTPUT OF THE TOLERANCE SCHEME
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT030
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
C
      DIMENSION FE(4)
C
      AL(ADEL) = -20.*ALOG10(ADEL)
C
      ME = 2
      IF (ITYP.GT.2) ME = 4
      WRITE (KA2,9999)
9999  FORMAT (//17H TOLERANCE SCHEME)
      WRITE (KA2,9998)
9998  FORMAT (//13H FILTER-TYPE )
      GO TO (10, 20, 30, 40), ITYP
  10  WRITE (KA2,9997)
9997  FORMAT (25X, 7HLOWPASS)
      GO TO 50
  20  WRITE (KA2,9996)
9996  FORMAT (25X, 8HHIGHPASS)
      GO TO 50
  30  WRITE (KA2,9995)
9995  FORMAT (25X, 8HBANDPASS)
      GO TO 50
  40  WRITE (KA2,9994)
9994  FORMAT (25X, 8HBANDSTOP)
  50  WRITE (KA2,9993)
9993  FORMAT (/15H APPROXIMATION )
      GO TO (60, 70, 80, 90), IAPRO
  60  WRITE (KA2,9992)
9992  FORMAT (25X, 11HBUTTERWORTH)
      GO TO 100
  70  WRITE (KA2,9991)
9991  FORMAT (25X, 11HCHEBYSHEV I)
      GO TO 100
  80  WRITE (KA2,9990)
9990  FORMAT (25X, 12HCHEBYSHEV II)
      GO TO 100
  90  WRITE (KA2,9989)
9989  FORMAT (25X, 8HELLIPTIC)
C
 100  IF (SF.EQ.0.) GO TO 120
      WRITE (KA2,9988) SF
9988  FORMAT (/18H SAMPLING FREQ.   , F15.6)
      Q = SF*0.5/PI
      DO 110 I=1,ME
        FE(I) = OM(I)*Q
 110  CONTINUE
      WRITE (KA2,9987) (FE(I),I=1,ME)
9987  FORMAT (20H CUTOFF FREQUENCIES , 4F13.6)
C
 120  WRITE (KA2,9986) (OM(I),I=1,ME)
9986  FORMAT (/20H NORM. CUTOFF FREQ. , 4F13.6)
      DO 130 I=1,ME
        Q = OM(I)*0.5
        FE(I) = SIN(Q)/COS(Q)
 130  CONTINUE
      WRITE (KA2,9985) (FE(I),I=1,ME)
9985  FORMAT (/20H CUTOFF FREQ. S-DOM., 4F13.6)
      P = 1. - ADELP
      Q = AL(P)
      P = SQRT(1.-P*P)
      WRITE (KA2,9984) ADELP, Q, P
9984  FORMAT (//20H PASSBAND RIPPLE(S) , F15.6, F12.4, 9H DB   P =,
     *    F12.4)
      Q = FLMA
      IF (ADELS.GT.0.) Q = AL(ADELS)
      WRITE (KA2,9983) ADELS, Q
9983  FORMAT (20H STOPBAND RIPPLE(S) , F15.6, F12.4, 3H DB)
      IF (NDEG.EQ.0) RETURN
      WRITE (KA2,9982) NDEG
9982  FORMAT (//32H FILTER DEGREE ASSIGNED TO ORDER, I3)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT031
C OUTPUT OF THE NORMALIZED PARAMETER IN THE S-DOMAIN
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT031
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
C
      WRITE (KA2,9999)
9999  FORMAT (//37H NORMALIZED PARAMETER IN THE S-DOMAIN)
      WRITE (KA2,9998) VD, VSN
9998  FORMAT (/3H VD, 6X, 1H=, F13.6/4H VSN, 5X, 1H=, F13.6)
      IF (ITYP.LE.2) RETURN
C
      WRITE (KA2,9997) A
9997  FORMAT (2H A, 7X, 1H=, F13.6)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT032
C OUTPUT OF THE REALIZED CUTOFF FREQUENCIES
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT032
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
C
      ME = 2
      IF (ITYP.GE.3) ME = 4
      WRITE (KA2,9999) (ROM(I),I=1,ME)
9999  FORMAT (/9H REALIZED/20H NORM. CUTOFF FREQ. , 4F13.6)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT033
C OUTPUT OF THE FILTER DEGREE
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT033
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLSN/ VSN, VD, A
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
      Q = ADEG*EDEG
      WRITE (KA2,9999) ADEG, Q, NDEG
9999  FORMAT (//20H MIN. FILTER DEGREE , F12.4//18H DEGREE EXTENSION ,
     *    F14.4//20H CHOSEN FILTER DEG. , I7)
      IF (ITYP.LE.2) RETURN
      IQ = 2*NDEG
      WRITE (KA2,9998) IQ
9998  FORMAT (26H FOR THE REFERENCE LOWPASS//21H ACTUAL FILTER DEGREE,
     *    I6)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT034
C OUTPUT OF THE CHARACTERISTIC FUNCTION  /K(JV)/**2
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT034
C
      DOUBLE PRECISION DSK(16), DK, DKS, DCAP02, DCAP04
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
      COMMON /RESIN2/ DK, DKS, DCAP02, DCAP04
      EQUIVALENCE (PREN(1),DSK(1))
C
      WRITE (KA2,9999) ADELTA
9999  FORMAT (/11H CAP. DELTA, 10X, F15.6)
      IF (IAPRO.LT.4) RETURN
      WRITE (KA2,9998)
9998  FORMAT (//50H ZEROS OF THE CHARACTERISTIC FUNCTION  /K(J*V)/**2/)
      IP = 0
      IZ = 2
      ZIM = 0.
      DO 10 I=1,NJ
        ZRE = SNGL(DSK(I))
        CALL OUT002
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT035
C OUTPUT OF THE TOLERANCE
C OUTPUT OF THE BOUND PAIR OF THE DESIGN PARAMETER AC
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT035
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
      WRITE (KA2,9999) UGC, OGC
9999  FORMAT (//38H BOUND PAIR OF THE DESIGN PARAMETER C , F9.4,
     *    16H  .LE.  C  .LE. , F9.4)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT036
C OUTPUT OF THE ZEROS AND POLES OF THE NORMALIZED REFERENCE LOWPASS
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT036
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
      NZ = NZM(4)
      IF (IZ.NE.0) GO TO 10
      WRITE (KA2,9999)
9999  FORMAT (////50H POLES AND ZEROS OF THE NORMALIZED REFERENCE LOWPA,
     *    18HSS IN THE S-DOMAIN)
      IF (IP.NE.0) NZ = 0
      GO TO 20
  10  WRITE (KA2,9998)
9998  FORMAT (////50H POLES AND ZEROS OF THE REFERENCE FILTER IN THE S-,
     *    6HDOMAIN)
  20  CONTINUE
      PRE = SFA
      CALL OUT001
      ME = MAX0(NZ,NJ)
      DO 60 I=1,ME
        IF (IP.EQ.0) GO TO 30
        IP = 2
        Q = SPI(I)
        IF (ABS(Q).LT.FLMI) IP = 1
        PRE = SPR(I)
        PIM = Q
        IF (NZ.LT.I) GO TO 40
        IF (NJ.LT.I) IP = 0
  30    IZ = NZERO(I)
        ZRE = 0.
        ZIM = SM(I,4)
        GO TO 50
  40    IZ = 0
  50    CALL OUT002
  60  CONTINUE
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT037
C OUTPUT OF THE CHOSEN DESIGN PARAMETER
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT037
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
      COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
C
      WRITE (KA2,9999) ACX, AC
9999  FORMAT (/26H CHOSEN DESIGN PAR.   CX =, F10.7, 10H  (=)  C =,
     *    F15.7)
      QP = 100.*RDELP/ADELP
      QS = 100.*RDELS/ADELS
      WRITE (KA2,9998) RDELP, QP, RDELS, QS
9998  FORMAT (/40H UTILIZATION OF THE PASSBAND  DELTA P  =, F10.7,
     *    3H  =, F10.5, 8H PERCENT/20X, 20HSTOPBAND  DELTA S  =, F10.7,
     *    3H  =, F10.5, 8H PERCENT)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT038
C OUTPUT OF THE POLES AND ZEROS IN THE Z-DOMAIN
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT038
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
      COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
      WRITE (KA2,9999)
9999  FORMAT (////32H POLES AND ZEROS IN THE Z-DOMAIN)
      NZ = NZM(4)
      PRE = ZFA
      CALL OUT001
      ME = MAX0(NJ,NZ)
C
      DO 30 I=1,ME
        IP = 0
        IF (NJ.LT.I) GO TO 10
        IP = 2
        Q = ZPI(I)
        IF (ABS(Q).LT.FLMI) IP = 1
        PRE = ZPR(I)
        PIM = Q
  10    IZ = 0
        IF (NZ.LT.I) GO TO 20
        IZ = NZERO(I)
        ZRE = ZZR(I)
        ZIM = ZZI(I)
  20    CALL OUT002
  30  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT039
C OUTPUT OF THE EXTREMA OF THE MAGNITUDE OF THE TRANSFER FUNCTION
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT039
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
      COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
      COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
C
      WRITE (KA2,9999)
9999  FORMAT (///51H EXTREMES OF THE MAGNITUDE OF THE TRANSFER FUNCTION,
     *    1H , 20H(COEFS NOT ROUNDED) )
      FN = 180./PI
      WRITE (KA2,9998)
9998  FORMAT (/1X, 26HMAXIMA IN THE PASSBAND (P)/1X, 15HMAXIMA IN THE S,
     *    11HTOPBAND (S)/1X, 14HBAND  S-DOMAIN, 7X, 8HZ-DOMAIN, 10X,
     *    9HMAGNITUDE/18X, 6HIN RAD, 3X, 9HIN DEGREE/)
C
      II = 1
      JJ = 3
  10  I = NZM(II)
      J = NZM(JJ)
      K = MAX0(I,J)
      DO 40 L=1,K
        IF (I.LT.L) GO TO 20
        OMEG1 = FN*ZM(L,II)
        AMAG1 = AMAGO(ZM(L,II))
        IF (J.LT.L) GO TO 30
        OMEG2 = FN*ZM(L,JJ)
        AMAG2 = AMAGO(ZM(L,JJ))
        WRITE (KA2,9997) SM(L,II), ZM(L,II), OMEG1, AMAG1, SM(L,JJ),
     *      ZM(L,JJ), OMEG2, AMAG2
9997    FORMAT (2H P, 3(F10.4, 2X), F12.6/2H S, 3(F10.4, 2X), F12.6)
9996    FORMAT (2H P, 3(F10.4, 2X), F12.6)
        GO TO 40
  20    OMEG2 = FN*ZM(L,JJ)
        AMAG2 = AMAGO(ZM(L,JJ))
        WRITE (KA2,9995) SM(L,JJ), ZM(L,JJ), OMEG2, AMAG2
9995    FORMAT (2H S, 3(F10.4, 2X), F12.6)
        GO TO 40
  30    WRITE (KA2,9996) SM(L,II), ZM(L,II), OMEG1, AMAG1
  40  CONTINUE
      IF (II.EQ.2) RETURN
      WRITE (KA2,9994)
9994  FORMAT (/1X, 26HMINIMA IN THE PASSBAND (P)/1X, 15HMINIMA IN THE S,
     *    11HTOPBAND (S)/1X, 14HBAND  S-DOMAIN, 7X, 8HZ-DOMAIN, 10X,
     *    9HMAGNITUDE/18X, 6HIN RAD, 3X, 9HIN DEGREE/)
      II = 2
      JJ = 4
      GO TO 10
C
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT042
C OUTPUT OF THE START GRID FOR THE SEARCH OF THE GLOBAL EXTREMA
C OF THE TRANSFER FUNCTION
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT042
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CGRID/ GR(64), NGR(12)
C
      WRITE (KA2,9999)
9999  FORMAT (///36H GRID FOR THE SEARCH OF THE EXTREMES)
      J = 1
      IF (NGR(1).NE.0) GO TO 10
      MA = NGR(8)
      ME = NGR(9)
      J = 4
  10  GO TO (20, 30, 40, 50, 80), J
  20  WRITE (KA2,9998)
9998  FORMAT (/10H PASS BAND/)
      MA = 1
      ME = NGR(3)
      GO TO 60
  30  MA = ME + 1
      ME = NGR(5)
      IF (ME.LT.MA) GO TO 70
      GO TO 60
  40  WRITE (KA2,9997)
9997  FORMAT (/10H STOP BAND/)
      MA = MAX0(NGR(3),NGR(5)) + 1
      ME = NGR(6)
      GO TO 60
  50  MA = ME + 1
      ME = NGR(7)
      IF (ME.LT.MA) GO TO 80
  60  WRITE (KA2,9996) (GR(I),I=MA,ME)
9996  FORMAT (1X, 7F10.4)
  70  J = J + 1
      GO TO 10
  80  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT043
C HEADLINE FOR THE SEARCH OF THE MINIMUM WORDLENGTH
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT043
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
C
      WRITE (KA2,9999)
9999  FORMAT (//29H SEARCH OF MINIMUM WORDLENGTH//4H IWL, 6X, 5HEPS-P,
     *    7X, 5HEPS-S, 7X, 4HPMAX)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT044
C OUTPUT OF ERROR EPSILON AND MAXIMUM MAGNITUDE
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT044
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
      COMMON /CEPSIL/ EPS(2), PMAX
C
9999  FORMAT (1X, I4, 3F12.6)
C
      WRITE (KA2,9999) IWL, EPS, PMAX
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT045
C OUTPUT OF MINIMUM WORDLENGTH
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT045
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
C
      WRITE (KA2,9999) IWLG, ADEPSG, IWLD, ADEPSD
9999  FORMAT (/19H MINIMUM WORDLENGTH, I4, 13H  DELTA EPS =,
     *    F12.6/19H TEST    WORDLENGTH, I4, 13H  DELTA EPS =, F12.6)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT046
C OUTPUT OF THE LAYOUT OF THE ROUNDED COEFFICIENTS
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT046
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
     *    ISTAB, IDEPSL, IDEPSD
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
      COMMON /SCRAT/ ADUM(32)
C
      DIMENSION ACO(16,5)
      DIMENSION ITEXT(5), JTEXT(5)
      DIMENSION IOCT(12)
      EQUIVALENCE (B2(1),ACO(1,1))
      EQUIVALENCE (ADUM(1),IOCT(1))
C
      DATA ITEXT(1), ITEXT(2), ITEXT(3), ITEXT(4), ITEXT(5) /1HB,1HB,
     *    1HB,1HC,1HC/
      DATA JTEXT(1), JTEXT(2), JTEXT(3), JTEXT(4), JTEXT(5) /1H2,1H1,
     *    1H0,1H1,1H0/
C
      WRITE (KA2,9999) IWL
9999  FORMAT (///35H LAYOUT OF THE ROUNDED COEFFICIENTS//11H WORDLENGTH,
     *    2H  , I4///11X, 4HCOEF, 10X, 5HCOEFS, 10X, 2HIR, 3X, 2HIE,
     *    10X, 2HID, 5X, 5HOCTAL/)
      DO 20 J=1,5
        WRITE (KA2,9998)
9998    FORMAT (/)
        DO 10 I=1,NB
          AC = ACO(I,J)
          ACS = ACOEFS(AC,IECO(I,J),IDCO(I,J),IRCO(J))
          ID = SIGN(2.,AC)
          Q = ACS
          CALL OCTAL(Q, IOCT, 12)
          WRITE (KA2,9997) ITEXT(J), JTEXT(J), I, AC, ACS, IRCO(J),
     *        IECO(I,J), ID, IDCO(I,J), (IOCT(II),II=1,12)
9997  FORMAT (1X, 2A1, 1H(, I2, 2H) , F13.8, 1H=, F13.8, 5H*2**(, I3,
     *        2H -, I3, 1H), 2X, I3, 2H**, I4, 1X, 12I1)
  10    CONTINUE
  20  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT060
C NEW PAGE
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT060
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      WRITE (KA2,9999)
9999  FORMAT (1H1)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT061
C HEAD OF PRINTER PLOT
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT061
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
      COMMON /CLINE/ Y, X, IZ(111)
      DIMENSION A(6)
      EQUIVALENCE (A(1),IZ(1))
      GO TO (10, 20, 30, 40), IPLOT
  10  WRITE (KA2,9999)
9999  FORMAT (1H /30X, 35HMAGNITUDE OF THE FREQUENCY RESPONSE/)
      GO TO 50
  20  WRITE (KA2,9998)
9998  FORMAT (1H /28X, 37HATTENUATION OF THE FREQUENCY RESPONSE/)
      GO TO 50
  30  WRITE (KA2,9997)
9997  FORMAT (1H /31X, 14HPHASE RESPONSE/)
      GO TO 50
  40  WRITE (KA2,9996)
9996  FORMAT (1H /31X, 16HIMPULSE RESPONSE/)
  50  WRITE (KA2,9995) A(1), A(3), A(5)
9995  FORMAT (10H MAGNITUDE, 2X, 1HX, F8.3, 2(11X, F9.3)/)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT062
C OUTPUT OF THE PLOTTER LINE
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT062
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
      COMMON /CLINE/ Y, X, IZ(111)
      DIMENSION IZDUM(56)
      DATA ICHP, ICHB /1H+,1H /
      DO 10 I=1,110,2
        J = (I+1)/2
        IZDUM(J) = ICHB
        IF (IZ(I).NE.ICHB) IZDUM(J) = IZ(I)
        IF (IZ(I+1).NE.ICHB) IZDUM(J) = IZ(I+1)
        IF ((IZ(I).EQ.ICHP) .OR. (IZ(I+1).EQ.ICHP)) IZDUM(J) = ICHP
  10  CONTINUE
      IZDUM(56) = IZ(111)
      IF (IPLOT.EQ.4) GO TO 20
      WRITE (KA2,9999) Y, X, IZDUM
9999  FORMAT (1X, E8.2, F6.3, 1X, 56A1)
      GO TO 30
  20  IX = INT(X)
      WRITE (KA2,9998) Y, IX, IZDUM
9998  FORMAT (1X, E8.2, I6, 1X, 56A1)
  30  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   HEAD
C OUTPUT OF A HEAD LINE
C-----------------------------------------------------------------------
C
      SUBROUTINE HEAD(KAN)
C
      WRITE (KAN,9999)
9999  FORMAT (41H     ===  DOREDI  ===  VERSION  V005 === )
C
C HERE ON LINE FOR THE OUTPUT OF DATE AND TIME;
C REMOVE THE FOLLOWING STATEMENT
      WRITE (KAN,9998)
9998  FORMAT (1H )
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT001
C HEADLINE OF THE POLES AND ZERO OUTPUT
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT001
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
C
      IF (PRE.NE.0.) WRITE (KA2,9999) PRE
9999  FORMAT (/24H CONSTANT GAIN FACTOR = , E15.7)
      WRITE (KA2,9998)
9998  FORMAT (//3X, 4HNUM., 11X, 5HPOLES, 15X, 4HNUM., 11X, 5HZEROS/)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUT002
C OUTPUT OF POLE AND ZERO
C-----------------------------------------------------------------------
C
      SUBROUTINE OUT002
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
      COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
C
      IF (IP.EQ.0) GO TO 20
      IF (IZ.EQ.0) GO TO 10
      WRITE (KA2,9999) IP, PRE, PIM, IZ, ZRE, ZIM
9999  FORMAT (2(1X, I5, 1X, F11.6, 5H +-J*, F11.6, 1X))
      RETURN
  10  WRITE (KA2,9999) IP, PRE, PIM
      RETURN
  20  IF (IZ.EQ.0) RETURN
      WRITE (KA2,9998) IZ, ZRE, ZIM
9998  FORMAT (36X, I5, 1X, F11.6, 5H +-J*, F11.6, 1X)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OCTAL
C OCTAL CONVERSION
C-----------------------------------------------------------------------
C
      SUBROUTINE OCTAL(AC, IOCT, N)
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
C
      DIMENSION IOCT(1)
C
      BC = AC
      IF (BC.LT.(-1.0)) BC = -1.0
      IF (BC.LT.0.) BC = -BC
      IF (BC.GE.1.) BC = 0.
C
      IOCT(1) = 0
      DO 10 J=2,N
        Q = 8.*BC
        IQ = INT(Q)
        IOCT(J) = IQ
        BC = Q - FLOAT(IQ)
  10  CONTINUE
C
      IF (AC.GE.FLMI) RETURN
      IA = 1
      IOCT(1) = 1
      IQ = N - 1
      DO 20 J=1,IQ
        JJ = N - J + 1
        IF (IA.EQ.1 .AND. IOCT(JJ).EQ.0) GO TO 20
        IOCT(JJ) = 7 - IOCT(JJ) + IA
        IA = 0
  20  CONTINUE
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   OUTSPE
C     DUMMY SUBROUTINE FOR 'SPECIAL OUTPUT'
C-----------------------------------------------------------------------
C
      SUBROUTINE OUTSPE
C
      COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, ISPOUT
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
C
C     DATA TRANSFER BY NAMED COMMON BLOCKS:
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
C     ISPMAX: MAXIMUM NUMBER OF IMPLEMENTED OPTIONS
      ISPMAX=0
C
      IF(ISPOUT.LE.0) RETURN
      IF(ISPOUT.GT.ISPMAX) GO TO 900
C
C     EXAMPLE FOR ISPMAX=3
C
C     GO TO (100,200,300),ISPOUT
C 100 CONTINUE
C     .
C     OPTION NO.  1
C     .
C     RETURN
C
C 200 CONTINUE
C     .
C     OPTION NO.  2
C     .
C     RETURN
C
C 300 CONTINUE
C     .
C     OPTION NO.  3
C     .
C     RETURN
C
  900 CALL ERROR (8)
      RETURN
      END
C
C ======================================================================
C
C DOREDI - SUBROUTINES: PART VI
C
C ======================================================================
C
C
C-----------------------------------------------------------------------
C SUBROUTINE:   SPOW
C CALCULATION OF THE L-2 NORM
C-----------------------------------------------------------------------
C
      SUBROUTINE SPOW(BMAX)
C
      COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
      COMMON /CPOW/ PNU, PNC, AND, ITCORP
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /CRENO/ LCNO, ICNO(16,5)
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
      BB2 = BN2(1)
      BB1 = BN1(1)
      BB0 = BN0(1)
      BN2(1) = B2(1)
      BN1(1) = B1(1)
      BN0(1) = B0(1)
      LQ = LCNO
      LCNO = 2
      IQ = IB
      IB = 1
      Q = PNU
      CALL POWER
      BMAX = PNU*FACT*FACT
      PNU = Q
      IB = IQ
      LCNO = LQ
      BN2(1) = BB2
      BN1(1) = BB1
      BN0(1) = BB0
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   POWER
C CALCULATION OF THE NOISE POWER GENERATED IN THE IB-TH BLOCK
C-----------------------------------------------------------------------
C
      SUBROUTINE POWER
C
      COMPLEX ZA, ZB, ZP, ZPS, ZQ, ZQI, Z0, Z1, ZBL, ZBL1, ZZ, ZZ1
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CPOL/ ZP(16,2), ZPS(16,2)
      COMMON /CRENO/ LCNO, ICNO(16,5)
      COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
      COMMON /CPOW/ PNU, PNC, AND, ITCORP
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
      DIMENSION EU(5)
C
      DATA Z0,Z1 /(0.,0.),(1.,0.)/
C
      ZBL(ZZ,A1,A0) = (ZZ+A1)*ZZ + A0
      ZBL1(ZZ1,A21,A11,A01) = (A21*ZZ1+A11)*ZZ1 + A01
C
      ZA = Z0
      PNU = 0.
      IF (LCNO.NE.2) GO TO 20
      EU(1) = 1.
      DO 10 I=2,5
        EU(I) = 0.
  10  CONTINUE
      GO TO 160
  20  DO 130 N=1,5
        AQ = AQC(N)
        IF (AQ.EQ.1.) GO TO 120
        IF (LCNO.EQ.(-1)) AQ = 0.
        KK = IB
        IF (KK.GT.NB) KK = NB
        K = ICNO(KK,N)
        QS = 1.
        QC1 = 0.
        GO TO (30, 40, 50, 60, 70, 80, 90), K
  30    QC2 = -1.
        Q = AQ/2.
        GO TO 100
  40    QC2 = 2.*(1.-3./PI)
        GO TO 110
  50    QC2 = 2.
        GO TO 110
  60    QC2 = -1.
        Q = (AQ-1.)/2.
        GO TO 100
  70    QS = 4. - 6./PI
        GO TO 30
  80    QC1 = -6.*(1.-2./PI)
        GO TO 40
  90    QS = 4.
        GO TO 30
 100    ZA = ZA + Q*ZBL1(Z1,BN2(N),BN1(N),BN0(N))
 110    EU(N) = QS + (QC2*AQ+QC1)*AQ
        GO TO 130
 120    EU(N) = 0.
 130  CONTINUE
      IF (IB.GT.NB) GO TO 340
      IF (CABS(ZA).EQ.0.) GO TO 150
      ZA = ZA/ZBL(Z1,C1(IB),C0(IB))
      IF (IB.EQ.NB) GO TO 150
      MA = IB + 1
      DO 140 I=MA,NB
        ZA = ZA*ZBL1(Z1,B2(I),B1(I),B0(I))/ZBL(Z1,C1(I),C0(I))
 140  CONTINUE
 150  AND = REAL(ZA)
C
 160  MARK = 0
      R1 = 0.
      DO 170 N=1,5
        IF (EU(N).EQ.0.) GO TO 170
        R1 = R1 + BN2(N)*BN0(N)*EU(N)
 170  CONTINUE
      R2 = 1.
      DO 190 J=IB,NB
        IF (J.NE.IB) R2 = R2*B2(J)*B0(J)
        Q = C0(J)
        IF (ABS(Q).GT.FLMI) GO TO 180
        IF (MARK.NE.0) GO TO 360
        MARK = J
        Q = C1(J)
 180    R2 = R2/Q
 190  CONTINUE
      IF (MARK.NE.0) GO TO 200
      PNU = PNU + R1*R2
      GO TO 280
C
 200  DO 270 I=IB,NB
        IF (I.EQ.MARK) GO TO 210
        Q1 = C0(I)
        Q2 = C1(I)
        GO TO 220
 210    Q1 = C1(I)
        Q2 = 1.
 220    IF (I.NE.IB) GO TO 240
        R = 0.
        DO 230 N=1,5
          IF (EU(N).EQ.0.) GO TO 230
          QQ1 = BN2(N)*(BN1(N)*Q1-BN0(N)*Q2)/Q1
          QQ2 = BN0(N)*(BN1(N)-BN2(N)*C1(IB))
          R = R + (QQ1+QQ2)*EU(N)
 230    CONTINUE
        R = R*R2
        GO TO 260
 240    R = R1*(B2(I)*(B1(I)*Q1-B0(I)*Q2)/Q1+B0(I)*(B1(I)-C1(I)*B2(I)))
        DO 250 J=IB,NB
          IF (J.NE.I .AND. J.NE.IB) R = R*B2(J)*B0(J)
          Q = C0(J)
          IF (J.EQ.MARK) Q = C1(J)
          R = R/Q
 250    CONTINUE
 260    PNU = PNU + R
 270  CONTINUE
C
 280  DO 330 J=IB,NB
        DO 320 K=1,2
          ZQ = ZP(J,K)
          IF (CABS(ZQ).LT.FLMI) GO TO 320
          KK = 3 - K
          ZQI = ZP(J,KK)
          ZB = Z1/(ZQ-ZQI)
          ZQI = Z1/ZQ
          ZB = ZB*ZQI
          ZA = Z0
          DO 290 N=1,5
            IF (EU(N).EQ.0.) GO TO 290
            ZA = ZA + EU(N)*ZBL1(ZQ,BN2(N),BN1(N),BN0(N))*ZBL1(ZQI,
     *          BN2(N),BN1(N),BN0(N))
 290      CONTINUE
          ZB = ZB*ZA
          DO 310 JJ=IB,NB
            IF (JJ.EQ.IB) GO TO 300
            ZB = ZB*ZBL1(ZQ,B2(JJ),B1(JJ),B0(JJ))*ZBL1(ZQI,B2(JJ),B1(JJ)
     *          ,B0(JJ))
C MAY 15, 1980 Q CHANGED TO ZQ IN FOLLOWING STATEMENT 300
 300        IF  (JJ.NE.J) ZB = ZB/ZBL(ZQ,C1(JJ),C0(JJ))
            ZB = ZB/ZBL(ZQI,C1(JJ),C0(JJ))
 310      CONTINUE
          PNU = PNU + REAL(ZB)
 320    CONTINUE
 330  CONTINUE
      RETURN
C
 340  AND = REAL(ZA)
      DO 350 N=1,5
        PNU = PNU + EU(N)
 350  CONTINUE
      RETURN
C
 360  CALL ERROR(33)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   SMAX
C COMPUTE THE MAXIMUM (IC=1) OR MINIMUM (IC=-1) OF THE TRANSFER FUNCTION
C-----------------------------------------------------------------------
C
      SUBROUTINE SMAX(BMAX, MMA, MME, IC)
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CGRID/ GR(64), NGR(12)
C
      DIMENSION O(5), B(5)
C
      FLM = 1./FLMA
      FLER1 = FLER*0.05
      MA = NGR(MMA)
      ME = NGR(MME)
      BMAX = 0.
      MAX = 0
      DO 30 I=MA,ME
        QR = AMAGO(GR(I))
        IF (IC.GT.0) GO TO 10
        IF (QR.LT.FLM) GO TO 220
        QR = 1./QR
  10    IF (QR.LT.BMAX) GO TO 20
        BMAX1 = BMAX
        MAX1 = MAX
        BMAX = QR
        MAX = I
        GO TO 30
  20    IF (QR.LT.BMAX1) GO TO 30
        BMAX1 = QR
        MAX1 = I
  30  CONTINUE
      JJ = 1
      IF ((ME-MA).GT.1) GO TO 40
      O(1) = GR(MA)
      O(3) = (GR(MA)+GR(ME))/2.
      O(5) = GR(ME)
      JJ = 3
      GO TO 60
  40  IF (IABS(MAX-MAX1).LT.2) JJ = 3
  50  IF (MAX.EQ.MA) MAX = MA + 1
      IF (MAX.EQ.ME) MAX = MAX - 1
      O(1) = GR(MAX-1)
      O(3) = GR(MAX)
      O(5) = GR(MAX+1)
  60  DO 80 I=1,5,2
        Q = AMAGO(O(I))
        IF (IC.GT.0) GO TO 70
        IF (Q.LT.FLM) GO TO 220
        Q = 1./Q
  70    B(I) = Q
  80  CONTINUE
      GO TO 150
C
  90  EMAX = 0.
      DO 100 I=1,5
        Q = B(I)
        IF (Q.LT.EMAX) GO TO 100
        MAX = I
        EMAX = Q
 100  CONTINUE
      IF (EMAX.EQ.BMAX) GO TO 110
      IF (EMAX/BMAX-1..LT.FLER) GO TO 180
 110  IF (MAX.EQ.3) GO TO 130
      IF (MAX.LT.3) GO TO 120
      O(1) = O(3)
      B(1) = B(3)
      O(3) = O(4)
      B(3) = B(4)
      GO TO 140
 120  O(5) = O(3)
      B(5) = B(3)
      O(3) = O(2)
      B(3) = B(2)
      GO TO 140
 130  O(1) = O(2)
      B(1) = B(2)
      O(5) = O(4)
      B(5) = B(4)
 140  BMAX = EMAX
 150  O(2) = (O(3)+O(1))/2.
      Q = AMAGO(O(2))
      IF (IC.GT.0) GO TO 160
      IF (Q.LT.FLM) GO TO 220
      Q = 1./Q
 160  B(2) = Q
      O(4) = (O(3)+O(5))/2.
      Q = AMAGO(O(4))
      IF (IC.GT.0) GO TO 170
      IF (Q.LT.FLM) GO TO 220
      Q = 1./Q
 170  B(4) = Q
      IF (O(5)-O(1).GE.FLER1) GO TO 90
 180  GO TO (190, 200, 210), JJ
 190  EMAX = BMAX
      BMAX = BMAX1
      BMAX1 = EMAX
      MAX = MAX1
      JJ = 2
      GO TO 50
 200  BMAX = AMAX1(BMAX,BMAX1)
 210  IF (IC.LT.0) BMAX = 1./BMAX
      RETURN
 220  BMAX = 0.
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   SMIMP
C CALCULATION OF THE SUM OF THE IMPULSE RESPONSE MAGNITUDES
C ABSOLUTE SCALING CRITERION
C-----------------------------------------------------------------------
C
      SUBROUTINE SMIMP(SUM)
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
      DIMENSION X1(16), X2(16)
      COMMON /SCRAT/ ADUM(32)
      EQUIVALENCE (ADUM(1),X1(1)), (ADUM(17),X2(1))
C
      DO 10 I=1,NB
        X1(I) = 0.
        X2(I) = 0.
  10  CONTINUE
C
      BMAX = 0.
      SUM = 0.
      IN = 0
      JN = 0
      Y = FACT
      IM = 0
      GO TO 30
  20  Y = 0.
  30  DO 40 I=1,NB
        U = Y
        Y = U*B2(I) + X1(I)
        X1(I) = U*B1(I) - Y*C1(I) + X2(I)
        X2(I) = U*B0(I) - Y*C0(I)
  40  CONTINUE
      JM = SIGN(1.,Y)
      Y = ABS(Y)
      SUM = SUM + Y
      BMAX = AMAX1(BMAX,Y)
      IF (IM.NE.0) GO TO 50
      EMAX = BMAX
      IM = JM
      GO TO 20
  50  IF (IM.EQ.JM) GO TO 70
      IF (JN.GT.3) RETURN
      IN = 0
      Q = EMAX/BMAX
      EMAX = Y
      IM = JM
      IF (Q.LT.FLER) GO TO 60
      JN = 0
      GO TO 20
  60  JN = JN + 1
      GO TO 20
  70  IF (IN.GT.3) RETURN
      EMAX = AMAX1(EMAX,Y)
      IF (Y/EMAX.LT.FLER) IN = IN + 1
      GO TO 20
      END
C
C-----------------------------------------------------------------------
C FUNCTION:     AMAGO
C COMPUTATION OF THE MAGNITUDE OF THE TRANSFER FUNCTION
C PARAMETER OMEGA IN RADIANS
C-----------------------------------------------------------------------
C
      FUNCTION AMAGO(OMEG)
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
      COMPLEX C, CH
C
      C = CMPLX(COS(OMEG),SIN(OMEG))
      CH = CMPLX(FACT,0.)
C
      DO 10 I=1,NB
        CH = CH*((B2(I)*C+B1(I))*C+B0(I))/((C+C1(I))*C+C0(I))
  10  CONTINUE
      AMAGO = CABS(CH)
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION:     PHASE
C PHASE OF OMEGA
C-----------------------------------------------------------------------
C
      FUNCTION PHASE(OM)
C
      COMPLEX CQ, CF
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      CF = CMPLX(FACT,0.)
      CQ = CMPLX(COS(OM),SIN(OM))
      DO 10 I=1,NB
        CF = CF*((B2(I)*CQ+B1(I))*CQ+B0(I))/((CQ+C1(I))*CQ+C0(I))
  10  CONTINUE
      Y = AIMAG(CF)
      X = REAL(CF)
      PHASE = ATAN2(Y,X)
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION:     RESP
C RESPONSE OF AN INPUT SIGNAL
C-----------------------------------------------------------------------
C
      FUNCTION RESP(X, IS)
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /CSTATV/ X1(16), X2(16)
C
      IF (IS.EQ.0) GO TO 20
C CLEAR STATE VARIABLES
      DO 10 I=1,NB
        X1(I) = 0.
        X2(I) = 0.
  10  CONTINUE
C
  20  Y = X*FACT
      DO 30 I=1,NB
        U = Y
        Y = U*B2(I) + X1(I)
        X1(I) = U*B1(I) - Y*C1(I) + X2(I)
        X2(I) = U*B0(I) - Y*C0(I)
  30  CONTINUE
      RESP = Y
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   CODE3
C STORE PAIRING AND ORDERING
C-----------------------------------------------------------------------
C
      SUBROUTINE CODE3(L, J, N)
C
      DIMENSION J(16), L(16)
      DO 10 I=1,N
        J(I) = L(I)
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   CODE4
C RESTORE PAIRING AND ORDERING
C-----------------------------------------------------------------------
C
      SUBROUTINE CODE4(L, J, N)
C
      DIMENSION J(16), L(16)
      DO 10 I=1,N
        L(I) = J(I)
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   CODE5
C SEARCH FOR A GIVEN CODE ARRAY IN LL
C RESTORE PAIRING AND ORDERING TO L
C-----------------------------------------------------------------------
C
      SUBROUTINE CODE5(L, J, N, LL, M)
C
      DIMENSION L(16), J(16), LL(16)
      M = 2
      DO 30 I=1,N
        IQ = J(I)
        DO 10 II=1,N
          IF (IQ.EQ.LL(II)) GO TO 20
  10    CONTINUE
        RETURN
  20    L(I) = IQ
  30  CONTINUE
      M = 1
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   QUAN
C CHECK THE ACTUAL COEFFICIENT QUANTIZATION
C-----------------------------------------------------------------------
C
      SUBROUTINE QUAN(ACO, ACQ)
C
      COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
      Q = ABS(ACO)
      IF (Q.NE.AINT(Q)) GO TO 10
      ACQ = 1.
      RETURN
  10  Q = Q*ALSBI
      ACQ = 0.5/ALSBI
      GO TO 30
  20  Q = Q*0.5
      ACQ = ACQ*2.
  30  IF (Q.EQ.AINT(Q)) GO TO 20
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION:     ACOEF
C RETRANSFORM ONE COEFFICIENT
C-----------------------------------------------------------------------
C
      FUNCTION ACOEF(AS, IE, ID, IR)
C
      Q = AS*2.**(IR-IE)
      ACOEF = Q
      IF (ID.LE.(-100)) RETURN
      AD = 2.**ID
      ACOEF = Q - SIGN(AD,AS)
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION:     ACOEFS
C TRANSFORM ONE COEFFICIENT
C-----------------------------------------------------------------------
C
      FUNCTION ACOEFS(AC, IE, ID, IR)
C
      AD = 0.
      IF (ID.GT.(-100)) AD = 2.**ID
      ACOEFS = (AC-SIGN(AD,AC))*2.**(IE-IR)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   COPY01
C COPY COEFFICIENT FIELD /FILT/ INTO THE FIELD /SFILT/
C-----------------------------------------------------------------------
C
      SUBROUTINE COPY01
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
     *    SC0(16)
C
      DIMENSION ACO(16,5), SACO(16,5)
      EQUIVALENCE (ACO(1,1),B2(1)), (SACO(1,1),SB2(1))
C
      NBS = NB
      SFACT = FACT
      DO 20 J=1,5
        DO 10 I=1,NB
          SACO(I,J) = ACO(I,J)
  10    CONTINUE
  20  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   COPY02
C COPY COEFFICIENT FIELD /SFILT/ INTO THE FIELD /FILT/
C-----------------------------------------------------------------------
C
      SUBROUTINE COPY02
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
     *    SC0(16)
C
      DIMENSION ACO(16,5), SACO(16,5)
      EQUIVALENCE (ACO(1,1),B2(1)), (SACO(1,1),SB2(1))
C
      NB = NBS
      FACT = SFACT
      DO 20 J=1,5
        DO 10 I=1,NB
          ACO(I,J) = SACO(I,J)
  10    CONTINUE
  20  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   COPY03
C COPY POLE LOCATIONS
C-----------------------------------------------------------------------
C
      SUBROUTINE COPY03
C
      COMPLEX ZP, ZPS
C
      COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
      COMMON /CPOL/ ZP(16,2), ZPS(16,2)
C
      DO 20 I=1,NB
        DO 10 J=1,2
          ZP(I,J) = ZPS(I,J)
  10    CONTINUE
  20  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION:     DELLK
C CALCULATE COMPLETE ELLIPTIC INTEGRAL OF FIRST KIND
C-----------------------------------------------------------------------
C
      DOUBLE PRECISION FUNCTION DELLK(DK)
C
      DOUBLE PRECISION DPI, DOMI
      DOUBLE PRECISION DE, DGEO, DK, DRI, DARI, DTEST
C
      COMMON /CONST/ PI, FLMA, FLMI, FLER
      COMMON /CONST2/ DPI, DOMI
C
      DATA DE /1.D00/
C
      DGEO = DE - DK*DK
      IF (DGEO) 10, 10, 20
  10  DELLK = DBLE(FLMA)
      RETURN
C
  20  DGEO = DSQRT(DGEO)
      DRI = DE
  30  DARI = DRI
      DTEST = DARI*DOMI
      DRI = DGEO + DRI
C
      IF (DARI-DGEO-DTEST) 50, 50, 40
  40  DGEO = DSQRT(DARI*DGEO)
      DRI = 0.5D00*DRI
      GO TO 30
  50  DELLK = DPI/DRI
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION:     DSN2
C CALCULATION OF THE JACOBI'S ELLIPTIC FUNCTION SN(U,K)
C
C EXTERNAL CALCULATION OF THE PARAMETER NECESSARY
C DK = K($K)
C DQ = EXP(-PI*K'/K) ... (JACOBI'S NOME)
C-----------------------------------------------------------------------
C
      DOUBLE PRECISION FUNCTION DSN2(DU, DK, DQ)
C
      DOUBLE PRECISION DPI, DOMI
      DOUBLE PRECISION DE, DZ, DPI2, DQ, DM, DU, DK, DC, DQQ, DH, DQ1,
     *    DQ2
C
      COMMON /CONST2/ DPI, DOMI
C
      DATA DE, DZ /1.D00,2.D00/
C
      DPI2 = DPI/DZ
      IF (DABS(DQ).GE.DE) GO TO 30
C
      DM = DPI2*DU/DK
      DC = DZ*DM
      DC = DCOS(DC)
C
      DM = DSIN(DM)*DK/DPI2
      DQQ = DQ*DQ
      DQ1 = DQ
      DQ2 = DQQ
C
      DO 10 I=1,100
        DH = (DE-DQ1)/(DE-DQ2)
        DH = DH*DH
        DH = DH*(DE-DZ*DQ2*DC+DQ2*DQ2)
        DH = DH/(DE-DZ*DQ1*DC+DQ1*DQ1)
        DM = DM*DH
C
        DH = DABS(DE-DH)
        IF (DH.LT.DOMI) GO TO 20
C
        DQ1 = DQ1*DQQ
        DQ2 = DQ2*DQQ
  10  CONTINUE
C
      GO TO 30
C
  20  DSN2 = DM
      RETURN
C
  30  DSN2 = 0.D00
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DELI1
C ELLIPTIC FUNCTION
C-----------------------------------------------------------------------
C
      SUBROUTINE DELI1(RES, X, CK)
C
      DOUBLE PRECISION RES, X, CK, ANGLE, GEO, ARI, PIM, SQGEO, AARI,
     *    TEST, DPI
      DOUBLE PRECISION DOMI
C
      COMMON /CONST2/ DPI, DOMI
      IF (X) 20, 10, 20
  10  RES = 0.D0
      RETURN
C
  20  IF (CK) 40, 30, 40
  30  RES = DLOG(DABS(X)+DSQRT(1.D0+X*X))
      GO TO 130
C
  40  ANGLE = DABS(1.D0/X)
      GEO = DABS(CK)
      ARI = 1.D0
      PIM = 0.D0
  50  SQGEO = ARI*GEO
      AARI = ARI
      ARI = GEO + ARI
      ANGLE = -SQGEO/ANGLE + ANGLE
      SQGEO = DSQRT(SQGEO)
      IF (ANGLE) 70, 60, 70
C
C REPLACE 0 BY A SMALL VALUE, TEST
C
  60  ANGLE = SQGEO*DOMI
  70  TEST = AARI*DOMI*1.D+05
      IF (DABS(AARI-GEO)-TEST) 100, 100, 80
  80  GEO = SQGEO + SQGEO
      PIM = PIM + PIM
      IF (ANGLE) 90, 50, 50
  90  PIM = PIM + DPI
      GO TO 50
 100  IF (ANGLE) 110, 120, 120
 110  PIM = PIM + DPI
 120  RES = (DATAN(ARI/ANGLE)+PIM)/ARI
 130  IF (X) 140, 150, 150
 140  RES = -RES
 150  RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION:     SINH
C-----------------------------------------------------------------------
C
      FUNCTION SINH(X)
C
      SINH = (EXP(X)-EXP(-X))/2.
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION:     COSH
C-----------------------------------------------------------------------
C
      FUNCTION COSH(X)
C
      COSH = (EXP(X)+EXP(-X))/2.
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION:     ARSINH
C-----------------------------------------------------------------------
C
      FUNCTION ARSINH(X)
C
      ARSINH = ALOG(X+SQRT(X*X+1.))
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION:     ARCOSH
C-----------------------------------------------------------------------
C
      FUNCTION ARCOSH(X)
C
      IF (X.LT.1.) GO TO 10
      ARCOSH = ALOG(X+SQRT(X*X-1.))
      RETURN
  10  ARCOSH = 0.
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   DSQRTC
C COMPUTATION OF THE COMPLEX SQUARE ROOT IN DOUBLE PRECISION
C     DU + J*DV = SQRT ( DX + J*DY )
C-----------------------------------------------------------------------
C
      SUBROUTINE DSQRTC(DX, DY, DU, DV)
C
      DOUBLE PRECISION DPI, DOMI, D1MACH
      DOUBLE PRECISION DX, DU, DY, DV, DQ, DP
      COMMON /CONST2/ DPI, DOMI
C
      DQ = DX
      DP = DY
C
      DV = 0.5D00*DQ
      DU = DQ*DQ + DP*DP
      DU = DSQRT(DU)
      DU = 0.5D00*DU
      DV = DU - DV
      DU = DV + DQ
      IF (DABS(DU).LE.3.D0*D1MACH(3)) DU = 0.D0
      DU = DSQRT(DU)
      IF (DABS(DV).LE.3.D0*D1MACH(3)) DV = 0.D0
      DV = DSQRT(DV)
      IF (DP.LT.(-DOMI)) DU = -DU
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:   ERROR
C OUTPUT OF AN ERROR MESSAGE
C-------------------------------------------------------------------
C
      SUBROUTINE ERROR(IERCOD)
C
      COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
C
      WRITE (KA2,9999) IERCOD
9999  FORMAT (22H *** ERROR *** NUMBER , I3)
      IF (IERCOD.GE.20) STOP
      RETURN
      END
Return to Main Software Page