C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C                        POPFNC.FOR
C
C  THE POPULATION FUNCTION A.K.A. MAXIMUM LIKELIHOOD FUNCTION.
C  THE IDEA BEHIND THIS ROUTINE IS TO ALLOCATE THE WORKSPACE VARIABLES
C  AND FILTER THEM DOWN TO THE ROUTINES WHICH NEED THEM!!!
C
C***********************************************************************
C  INPUT TO THE MLE ROUTINE
C        NX------LENGTH OF DATA TO ANALYZE                     (INTEGER)
C        X-------THE TIME SERIES DATA SEQUENCE                    (REAL)
C        NCPHI---FOR AN ARMA(P,Q) THE P ASSOCIATED WITH THIS   (INTEGER)
C        NCTHET--FOR AN ARMA(P,Q) THE Q ASSOCIATED WITH THIS   (INTEGER)
C        ETALIM--THE MIN AND MAX VALUES OF THE ETA'S              (REAL)
C                MATRIX WHERE FIRST ROW IS THE MINIMUM AND
C                SECOND ROW IS THE MAXIMUM. I.E.
C                ETALIM(1,J)---MINIMUM JTH COEFFICIENT VALUE
C                ETALIM(2,J)---MAXIMUM JTH COEFFICIENT VALUE
C                1 <= J <= (P+Q+1)
C                WHEN J=(P+Q+1) ==> STANDARD DEVIATION OF THE NOISE
C        C1------DISTRIBUTIONAL CONSTANT NUMBER 1                 (REAL)
C        C2------DISTRIBUTIONAL CONSTANT NUMBER 2                 (REAL)
C        NDIV----GRID DENSITY FOR SEARCH (NUMBER OF DIVISIONS) (INTEGER)
C        NREP----NUMBER OF GRID REDUCTIONS (NUMBER OF REPETITIONS) (INT)
C        NSUP----SUPPRESS THE SIMULATED ANNEALING AND GENETIC ALGORITHM
C                SUPSAG =1 FOR THE SUPPRESSION                 (INTEGER)
C        RSPC----REAL WORK SPACE VARIABLE WHOSE MIN DIMENSION IS  (REAL)
C                24*(NP)+(NP+1)*NPOP+2*NX+10*(MAXL+1)+10*NP*NP)
C                WHERE NP=P+Q+1, MAXL=100, NPOP=100
C        CSPC----COMPLEX WORK SPACE VARIABLE WHOSE MIN DIMENSION IS
C                4*(NP)                                        (COMPLEX)
C        ISPC----INTEGER WORK SPACE WHOSE MIN DIMENSION IS     (INTEGER)
C                4*(NP)+2*NPOP
C  OUTPUT OF THE MLE ROUTINE
C        Z-------MOVING AVERAGE SEQUENCE                          (REAL)
C        NZ------LENGTH OF Z SEQUENCE                          (INTEGER)
C        XMLE----THE MLE FUNCTION                                 (REAL)
C        ETA-----THE COEFFICIENTS OF THE MLE PROCESS          (REAL VEC)
C        COVVEC--THE COVARIANCE VECTOR ASSOCIATED WITH ETA    (REAL VEC)
C        STDEV---THE STANDARD DEVIATION OF THE COEFFICIENTS OF ETA(RVEC)
C        NNZC----THE NUMBER OF NON-ZERO COEFFICIENTS           (INTEGER)
C***********************************************************************
      SUBROUTINE  NGMLE(NX,X,NCPHI,NCTHET,ETALIM,C1,C2,NDIV,NREP,NSUP,
     1                 RSPC,CSPC,ISPC,Z,NZ,XMLE,ETA,COVVEC,STDEV,NNZC)

        INTEGER NX,NZ,NZT,NNZC,MAXL,I,IQ,NCPHI,NCTHET,ICOEFF,IRW,IW,ICW
        INTEGER NPI,NPO,NTI,NTO,MP1,NP2,NDIV,NREP,NSUP,NPOP,NPARM,NP,J
        REAL ACURCY
        PARAMETER(NPOP=100,MAXL=100)
        PARAMETER(ACURCY=.00000010)
        REAL X(NX),ETALIM(2,NCPHI+NCTHET+1),C1,C2,Z(NX),XMLE
        REAL ETA(NCPHI+NCTHET+1),STDEV(NCPHI+NCTHET+1)
        REAL COVVEC((NCPHI+NCTHET+1)*(NCPHI+NCTHET+2)/2)
        INTEGER I1,I2,I3,I4,I5,I6,ISPC(4*(NCPHI+NCTHET+1)+2*NPOP)
        REAL RSPC(24*(NCPHI+NCTHET+1)+(NCPHI+NCTHET+1+1)*NPOP+
     1     2*NX+10*(MAXL+1)+10*(NCPHI+NCTHET+1)*(NCPHI+NCTHET+1))
        COMPLEX CSPC(NCPHI+NCTHET+1)
        REAL SIGMIN,SIGMAX,SEQUEN
C  THESE ARE THE VARIABLE NAMES ALLOCATED TO THE WORKING ARRAYS
C        INTEGER IRNGP(NPARM),IRNGPB(NPARM),IPOP(NPOP),IDU(NPOP)
C        INTEGER IPLIM(NPARM),IP(NPARM)
C        REAL SMA(NPARM),GA(NPARM)
C        REAL RPARM(NPARM,NPOP),POP(NPOP),PL(2,NPARM)
C        REAL PINT(NPARM),PR1(NPARM),PR2(NPARM),PA(NPARM)
C        REAL DU(NPARM),SA1(NPARM),SA2(NPARM)
C        REAL DELX(NPARM),TCOEFF(NPARM)
C        REAL COEFF(NPARM)
C        REAL X(NX),Z(NX),ZS(NPARM+MAXL+MAXX),Y(NPARM+MAXL+MAXX)
C        REAL COEFR(MAXL),COEFPI(MAXL),THET(MAXL),COEFPO(MAXL)
C        REAL PI(NPARM),PO(NPARM),TI(NPARM),TO(NPARM)
C        REAL COEFCC(NPARM),PW(NPARM)
C        REAL WORK(NPARM*NPARM+2*NPARM)
C        COMPLEX ZPHI(NPARM),ZTHETA(NPARM)
C        COMPLEX ROOTIN(NPARM),ROOTOT(NPARM)
        
C  ASSIGN SOME CONSTANTS
        NZ=NX
        NPARM=NCPHI+NCTHET+1
        ICOEFF=NCPHI+NCTHET
        IQ=SQRT(NX*1.0)
        NP=NPARM
        NP2=NP*NP
        MP1=MAXL+1
        SIGMIN=ETALIM(1,NPARM)
        SIGMAX=ETALIM(2,NPARM)
        I1=1+14*NP+NPOP+NP*NPOP
        I2=I1+2*NX+2*NP+2*MAXL
        I3=I2+4*MAXL
        I4=I3+8*NP+NP2
C  CHANGE OF VARIABLE TO COMPLEX        I5=I4+8*NP
        I5=I4
        I6=I5+4*MP1+8*NP2

C  INITIALIZE THE WORK SPACE VARIABLES
        IRW=24*NP+(NP+1)*NPOP+2*NX+10*(MAXL+1)+10*NP*NP
        IW=4*NP+2*NPOP
        ICW=4*NP
        DO 5  I=1,IRW
               RSPC(I)=0.0
5       CONTINUE
        DO 10 I=1,IW
               ISPC(I)=0
10      CONTINUE
        DO 15 I=1,ICW
               CSPC(I)=CMPLX(0,0)
15      CONTINUE

C  THIS IS THE ALLOCATION OF THE WORKING ARRAYS FOR THE SEARCH ROUTINE
C  PERFORM A SEQUENTIAL OR STOCHASTIC SEARCH
C       CALL REPSQS(NPOP,ICOEFF,NDIV,NREP,NSUP,IRNGP,
C     1            IRNGPB,IPOP,IDU,IPLIM,IP,
C     2            COEFF,PL,SMA,GA,RPARM,POP,ETALIM,PINT,
C     3            PR1,PR2,PA,DU,SA1,SA2,DELX,TCOEFF,
C     4   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
C     5   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,
C     6   COEFR,COEFPI,THET,COEFPO,PI,PO,TI,TO,
C     7   COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT,XMLE)

        CALL REPSQS(NPOP,ICOEFF,NDIV,NREP,NSUP,
     1  ISPC(1     ),ISPC(1+NP   ),ISPC(1+2*NP),ISPC(1+2*NP+NPOP),
     2  ISPC(1+2*NP+2*NPOP),ISPC(1+3*NP+2*NPOP),RSPC(1)     ,
     3  RSPC(1+NP) ,RSPC(1+3*NP),RSPC(1+4*NP),RSPC(1+5*NP)     ,
     4  RSPC(1+5*NP+NP*NPOP),ETALIM,RSPC(1+5*NP+NPOP+NP*NPOP),
     5  RSPC(1+6*NP+NPOP+NP*NPOP),RSPC(1+7*NP+NPOP+NP*NPOP),
     6  RSPC(1+8*NP+NPOP+NP*NPOP),RSPC(1+9*NP+NPOP+NP*NPOP),
     7  RSPC(1+10*NP+NPOP+NP*NPOP),RSPC(1+11*NP+NPOP+NP*NPOP),
     8  RSPC(1+12*NP+NPOP+NP*NPOP),RSPC(1+13*NP+NPOP+NP*NPOP),
     9  NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,ACURCY,SIGMIN,SIGMAX,
     1  X       ,Z          ,RSPC(I1),RSPC(I1+NX+NP+MAXL),
     1  RSPC(I2)   ,RSPC(I2+MAXL),RSPC(I2+2*MAXL),RSPC(I2+3*MAXL),
     2  RSPC(I3) ,RSPC(I3+NP),RSPC(I3+2*NP),RSPC(I3+3*NP),
     3  RSPC(I3+4*NP),RSPC(I3+5*NP),RSPC(I3+6*NP),
     4  CSPC(1),CSPC(1+NP),CSPC(1+2*NP),CSPC(1+3*NP),XMLE)

C  AGAIN THE WORKING ARRAY ALLOCATION
C  CALCULATE THE INPUT SEQUENCE USING THE OUTPUT SEQUENCE AND THE
C  ARMA COEFFICIENTS
C                CALL PUCOEF(COEFF,NCPHI,NCTHET,
C     1                   PI,NPI,PO,NPO,
C     2                   TI,NTI,TO,NTO,NPARM,
C     3                   COEFCC,PW,
C     4                   ZPHI,ZTHETA,ROOTIN,ROOTOT,WORK)
C        CALL PSIOPI(TI,NTI,TO,NTO,PI,NPI,PO,NPO,NPARM,
C     1                         X,NX,Z,NZ,MAXL,ACURCY,
C     2                  COEFR,COEFPI,THET,COEFPO,ZS,Y)
C        COEFF(NPARM)=SEQUEN(SIGMAX,SIGMIN,Z,NZ,IQ)

                CALL PUCOEF(RSPC(1),NCPHI,NCTHET,
     1  RSPC(I3) ,NPI,RSPC(I3+NP),NPO,RSPC(I3+2*NP),NTI,RSPC(I3+3*NP),
     2  NTO,NPARM,RSPC(I3+4*NP),RSPC(I3+5*NP),CSPC(1),CSPC(1+NP),
     3  CSPC(1+2*NP),CSPC(1+3*NP),RSPC(I3+6*NP))
 
                CALL PSIOPI(RSPC(I3+2*NP),NTI,RSPC(I3+3*NP),NTO,
     1                      RSPC(I3) ,NPI,RSPC(I3+NP),NPO,NPARM,
     2                      X,NX,Z,NZ,MAXL,ACURCY,
     3      RSPC(I2),RSPC(I2+MAXL),RSPC(I2+2*MAXL),RSPC(I2+3*MAXL),
     4      RSPC(I1),RSPC(I1+NX+NP+MAXL))
        RSPC(NPARM)=SEQUEN(SIGMAX,SIGMIN,Z,NZ,IQ)

C PLACE THE TRUNCATED SERIES IN THE ARRAY
        NZT=0
        DO 20 I=IQ+1,NX-IQ
                NZT=NZT+1
                Z(NZT)=Z(I)
20        CONTINUE
        NZ=NZT

C  THE FINAL WORK SPACE ALLOCATIONS
C  FIND THE COVARIANCE MATRIX AND PARAMETER ESTIMATES
C        CALL COVAR(C1,C2,COEFF,NCPHI,NCTHET,
C     1       MAXL,ACURCY,NX,IQ,NPARM,ETA,COVMAT,STDEV,NNZC,
C     2       NPI,NPO,NTI,NTO,PI,PO,TI,TO,
C     3       RSPC(I5),            RSPC(I5+  MP1),
C     4       RSPC(I5+2*MP1),      RSPC(I5+3*MP1),
C     6       RSPC(I5+4*MP1),      RSPC(I5+4*MP1+NP2),
C     7       RSPC(I5+4*MP1+2*NP2),RSPC(I5+4*MP1+3*NP2),
C     8       RSPC(I5+4*MP1+4*NP2),RSPC(I5+4*MP1+5*NP2),
C     9       RSPC(I5+4*MP1+6*NP2),RSPC(I5+4*MP1+7*NP2))

        CALL COVAR(C1,C2,RSPC(1),NCPHI,NCTHET,
     1       MAXL,ACURCY,NX,IQ,NPARM,ETA,RSPC(I6),STDEV,NNZC,
     2       NPI,NPO,NTI,NTO,
     3       RSPC(I3) ,RSPC(I3+NP),RSPC(I3+2*NP),RSPC(I3+3*NP),
     4       RSPC(I5),            RSPC(I5+  MP1),
     5       RSPC(I5+2*MP1),      RSPC(I5+3*MP1),
     6       RSPC(I5+4*MP1),      RSPC(I5+4*MP1+NP2),
     7       RSPC(I5+4*MP1+2*NP2),RSPC(I5+4*MP1+3*NP2),
     8       RSPC(I5+4*MP1+4*NP2),RSPC(I5+4*MP1+5*NP2),
     9       RSPC(I5+4*MP1+6*NP2),RSPC(I5+4*MP1+7*NP2))
 
C  PUT THE MATRIX DATA INTO A VECTOR THEN EXIT THE SUBROUTINE
        NZT=0
        DO 40 I=1,NNZC
              DO 30 J=I,NNZC
              NZT=NZT+1
              COVVEC(NZT)=RSPC(I6+NPARM*(J-1)+I-1)
30            CONTINUE
40      CONTINUE

        RETURN
        END

C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C  THESE ROUTINES CALCULATE THE LOG-LIKELIHOOD FUNCTION AS DEFINED
C  IN DR. LII'S PAPER
C
C***********************************************************************
C        POPFNC--POPULATION FUNCTION WHICH IS MAXIMIZED BY THE 
C                GENETIC/SIMULATED ANNEALING/SEQUENTIAL SEARCH ROUTINE
C        PSIOPI--IN ARMA TERMINOLOGY THIS IS AN AR(INFINITY)COEFFICIENT
C                DEFINITION.  AS IT IS USED HERE IT CALCULATES PI WHICH
C                IS AN MA(INFINITY) REPRESENTATION.
C***********************************************************************
C
C        CALCULATE THE POPULATION FUNCTION TO MAXIMIZE
C***********************************************************************
C        INPUT
C                COEFF--THE ARMA COEFFICIENTS
C        OUTPUT
C                POPFNC-LOG/LIKELIHOOD FUNCTION
C        OTHER
C                OTHER VARIABLES ARE JUST WORKSPACE
C***********************************************************************

      FUNCTION POPFNC(COEFF,
     1                NP,NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     2                ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,
     3                COEFR,COEFPI,THET,COEFPO,PI,PO,TI,TO,
     4                COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT)

        INTEGER NX,MAXL,IQ,NCPHI,NCTHET,NP,NZ,NPI,NPO,NTI,NTO,I
        REAL ACURCY,SIGMIN,SIGMAX,SIGMA,SEQUEN,XTL,XTP,SUM,TOT
        REAL FSIG,XLFSMN,POPFNC
        REAL COEFF(NP)
        REAL X(NX),Z(NX),ZS(NP+MAXL+NX),Y(NP+MAXL+NX)
        REAL COEFR(MAXL),COEFPI(MAXL),THET(MAXL),COEFPO(MAXL)
        REAL PI(NP),PO(NP),TI(NP),TO(NP)
        REAL COEFCC(NP),PW(NP),WORK(NP*NP+2*NP)
        COMPLEX ZPHI(NP),ZTHETA(NP)
        COMPLEX ROOTIN(NP),ROOTOT(NP)

C  CALCULATE THE INPUT SEQUENCE USING THE OUTPUT SEQUENCE AND THE
C  ARMA COEFFICIENTS
                CALL PUCOEF(COEFF,NCPHI,NCTHET,
     1                   PI,NPI,PO,NPO,
     2                   TI,NTI,TO,NTO,NP,
     3                   COEFCC,PW,
     4                   ZPHI,ZTHETA,ROOTIN,ROOTOT,WORK)

        CALL PSIOPI(TI,NTI,TO,NTO,PI,NPI,PO,NPO,NP,
     1                         X,NX,Z,NZ,MAXL,ACURCY,
     2                  COEFR,COEFPI,THET,COEFPO,ZS,Y)


        SIGMA=SEQUEN(SIGMAX,SIGMIN,Z,NZ,IQ)

C  CALCULATE THE LOG-LIKELIHOOD FUNCTION
        XTL=ALOG(ABS(TI(NTI+1)))
        XTP=ALOG(ABS(PI(NPI+1)))

        SUM=0
        TOT=NZ-2*IQ
        DO 10 I=IQ+1,NZ-IQ
        SUM=SUM+ALOG(FSIG(Z(I),SIGMA))
10      CONTINUE
        XLFSMN=SUM/TOT

        POPFNC=XLFSMN+XTP-XTL

        RETURN
        END

C***********************************************************************
C        PSIOPI--IN ARMA TERMINOLOGY THIS IS AN AR(INFINITY)COEFFICIENT
C                DEFINITION.  AS IT IS USED HERE IT CALCULATES PI WHICH
C                IS AN MA(INFINITY) REPRESENTATION.
C***********************************************************************
C        CALCULATE THE OUTPUT USING THE PSI COEFFICIENTS DEVELOPED 
C        FROM SEGREGATED POLYNOMIALS
C***********************************************************************
C        PSI=THETAIN*THETAOUT/PHIIN/PHIOUT
C        X=PSI*Z
C  INPUT
C        PI------PHI POLYNOMIAL COEFFICIENTS INSIDE UNIT CIRCLE
C        NPI-----DEGREE OF PI POLYNOMIAL
C        PO------PHI POLYNOMIAL COEFFICIENTS OUTSIDE UNIT CIRCLE
C        NPO-----DEGREE OF PO POLYNOMIAL
C        TI------THETA POLYNOMIAL COEFFICIENTS INSIDE UNIT CIRCLE
C        NTI-----DEGREE OF TI POLYNOMIAL
C        TO------THETA POLYNOMIAL COEFFICIENTS OUTSIDE UNIT CIRCLE
C        NTO-----DEGREE OF TO POLYNOMIAL
C        NP------DIMENSION OF THE ABOVE POLYNOMIALS
C        Z-------INPUT STATIONARY SEQUENCE
C        NZ------LENGTH OF INPUT STATIONARY SEQUENCE
C        NX------LENGTH OF OUTPUT STATIONARY SEQUENCE
C        MAXL----MAXIMUM NUMBER OF COEFFICIENTS IN THE INVERTED
C                POLYNOMIALS
C        ACURCY--ACCURACY OF THE INVERTED POLYNOMIALS
C  OUTPUT
C        X-------OUTPUT STATIONARY SEQUENCE
C  OTHER
C        JUST WORKSPACE VARIABLES
C******************************************************************

      SUBROUTINE PSIOPI(PI,NPI,PO,NPO,TI,NTI,TO,NTO,NP,
     1                  Z,NZ,X,NX,MAXL,ACURCY,
     2                  COEFR,COEFPI,THET,COEFPO,ZS,Y)

        INTEGER NPI,NPO,NTI,NTO,NP,NZ,NX,MAXL,NCOFPO,NCOFPI
        INTEGER LNGTHX,LCOFPI,LINIT,I,J
        REAL ACURCY,SUM
        REAL PI(NP),PO(NP),TI(NP),TO(NP),X(NX),Z(NZ)
        REAL COEFR(MAXL),COEFPI(MAXL),THET(MAXL),COEFPO(MAXL)
        REAL ZS(NP+MAXL+NZ),Y(NP+MAXL+NZ)

C  MULTIPLY THE POLYNOMIALS 
        CALL PMULT(TI,NTI,TO,NTO,THET,NP)
        CALL POLYIV(1,PO,NPO,NP,
     1               COEFPO,MAXL,MAXL,
     2               ACURCY,NCOFPO)
        IF(NCOFPO.GT.(NTI+NTO)) THEN
                NCOFPO=NCOFPO-NTI-NTO
        ELSE
                NCOFPO=0
        ENDIF
        CALL PMULT(COEFPO,NCOFPO,THET,NTI+NTO,COEFR,100)
        CALL POLYIV(2,PI,NPI,NP,
     1               COEFPI,MAXL,MAXL,
     2               ACURCY,NCOFPI)

C  DETERMINE THE LENGTHS OF THE COEFFICIENTS
        LNGTHX=NTI+NTO+NCOFPO+1
        LCOFPI=NCOFPI+1
        LINIT=LNGTHX+NZ
        IF(LINIT.LT.(LCOFPI+NZ)) LINIT=LCOFPI

C INITIALIZE THE INTERMEDIATE AND WORKING ARRAYS
        DO 10 I=1,LINIT
                Y(I)=0
                ZS(I)=0
10        CONTINUE
        DO 20 I=1,NZ
                ZS(I+LNGTHX)=Z(I)
20        CONTINUE

C  CALCULATE THE INTERMEDIATE VARIABLE Y
                DO 40 I=NZ,1,-1
                        SUM=0.0
                        DO 30 J=1,LNGTHX
                        SUM=SUM+COEFR(J)*ZS(I+LNGTHX-J+1)
30                      CONTINUE
                Y(I)=SUM
40              CONTINUE
        
C CALCULATE THE FINAL VALUE
                DO 60 I=1,NX,1
                        SUM=0.0
                        DO 50 J=1,LCOFPI
                        SUM=SUM+COEFPI(J)*Y(I+J-1+NPI)
50                      CONTINUE
                X(I)=SUM
60              CONTINUE

          RETURN
        END

C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012

C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C  THESE ROUTINES PERFORM THE NECESSARY FUNCTIONS WITH RESPECT TO 
C  POLYNOMIAL MANIPULATION
C***********************************************************************
C        PUCOEF--PUTS COEFFICIENTS INTO THE AR AND MA POLYNOMIALS
C        OROOTS--BREAKS THE AR AND MA POLYNOMIALS INTO INSIDE AND
C                OUTSIDE UNIT CIRCLE
C        PMULT---MULTIPLIES POLYNOMIALS TOGETHER TO CREATE A THIRD
C        PMULTI--PERFORMS POLYNOMIAL MULTIPLICATION IN PLACE
C        POLYIV--CREATES AN INVERTED POLYNOMIAL SEQUECE
C***********************************************************************
C
C        PUTS THE COEFFICIENTS INTO THEIR AR AND MA POLYNOMIALS
C***********************************************************************
C
C        INPUT
C                  COEFF-------THE COEFFICIENTS PASSED TO THIS ROUTINE
C                  NCPHI-------NUMBER OF COEFFICIENTS IN AR
C                  NCTHET------NUMBER OF COEFFICIENTS IN MA
C                  N-----------DIMENSION OF POLYNOMIALS
C        OUTPUT
C                  PHI---------INSIDE AND OUTSIDE AR POLYNOMIAL
C                  NIP,NOP-----DEGREE OF POLYNOMIALS PHI INSIDE AND 
C                            OUTSIDE OF UNIT CIRCLE
C                  THETA-------INSIDE AND OUTSIDE MA POLYNOMIAL
C                  NIT,NOT-----DEGREE OF POLYNOMIALS THETA INSIDE AND 
C                            OUTSIDE OF UNIT CIRCLE
C***********************************************************************

      SUBROUTINE PUCOEF(COEFF,NCPHI,NCTHET,
     1                   PHIIN,NIP,PHIOUT,NOP,
     2                   THETAI,NIT,THETAO,NOT,N,
     3                   COEFCC,PW,
     4                   ZPHI,ZTHETA,ROOTIN,ROOTOT,WORK)
        INTEGER NCPHI,NCTHET,NIP,NOP,NIT,NOT,N,I,IERR
        REAL COEFF(N),PHIIN(N),PHIOUT(N),THETAI(N),THETAO(N)
        REAL COEFCC(N),PW(N),WORK(N*N+2*N)
        COMPLEX ZPHI(N),ZTHETA(N)
        COMPLEX ROOTIN(N),ROOTOT(N)

C  INITIALIZE THE VECTORS
        DO 5 I=1,N
                PHIIN(I)=0.0
                PHIOUT(I)=0.0
                THETAI(I)=0.0
                THETAO(I)=0.0
5        CONTINUE

C  FIND THE ROOTS OF PHI THEN PUT INTO POLYNOMIAL
        IF(NCPHI.GT.0) THEN

C  FILL THE COMPLEX COEFFICIENTS FOR THE NUMERICAL RECIPE ROUTINE
        COEFCC(NCPHI+1)=1.0
        DO 10 I=1,NCPHI
                COEFCC(NCPHI+1-I)=COEFF(I)
10        CONTINUE

        CALL RPQR79(NCPHI,COEFCC,ZPHI,IERR,WORK)
        CALL OROOTS(ZPHI,NCPHI,PHIIN,NIP,PHIOUT,NOP,N,
     1                  PW,ROOTIN,ROOTOT)
        ELSE
                NIP=0
                NOP=0
                PHIIN(1)=1.0
                PHIOUT(1)=1.0
        ENDIF
                

C  FIND THE ROOTS OF THETAT THEN PUT INTO POLYNOMIAL
        IF(NCTHET.GT.0) THEN

C  FILL THE COMPLEX COEFFICIENTS FOR THE NUMERICAL RECIPE ROUTINE
        COEFCC(NCTHET+1)=1.0
        DO 20 I=1,NCTHET
                COEFCC(NCTHET+1-I)=COEFF(NCPHI+I)
20        CONTINUE

        CALL RPQR79(NCTHET,COEFCC,ZTHETA,IERR,WORK)        
        CALL OROOTS(ZTHETA,NCTHET,THETAI,NIT,THETAO,NOT,N,
     1                  PW,ROOTIN,ROOTOT)

        ELSE
                NIT=0
                NOT=0
                THETAI(1)=1.0
                THETAO(1)=1.0
        ENDIF

        RETURN
        END

C***********************************************************************
C       THIS ROUTINE CREATES POLYNOMIALS INSIDE AND OUTSIDE UNIT CIRCLE
C***********************************************************************
C        INPUT
C                ROOTS-------THE ROOTS TO PLACE INTO POLYNOMIALS
C                M-----------NUMBER OF ROOTS
C        OUTPUT
C                  POLYIN------POLYNOMIAL WITH ROOTS INSIDE UNIT CIRCLE
C                  NI----------DEGREE OF POLYNOMIAL INSIDE UNIT CIRCLE
C                  POLYOT-----POLYNOMIAL WITH ROOTS OUTSIDE UNIT CIRCLE
C                  NO----------DEGREE OF POLYNOMIAL OUTSIDE UNIT CIRCLE
C                  N-----------DIMENSION OF POLYNOMIALS
C***********************************************************************

      SUBROUTINE OROOTS(ROOTS,M,POLYIN,NI,POLYOT,NO,N,
     1                  PW,ROOTIN,ROOTOT)
        INTEGER M,NI,NO,N,I,ICOUNT
        COMPLEX ROOTS(N),ROOTIN(N),ROOTOT(N)
        REAL TOL,DENOM
        REAL POLYIN(N),POLYOT(N),PW(N),POLY(3)
        PARAMETER(TOL=1.0E-6)

C  FIND THE MAGNITUDE OF EACH ROOT
        NI=0
        NO=0
        DO 10 I=1,M
        IF(CABS(ROOTS(I)).GT.1.0) THEN
                NO=NO+1
                ROOTOT(NO)=ROOTS(I)
        ELSE
                NI=NI+1
                ROOTIN(NI)=ROOTS(I)
        ENDIF
        POLYIN(I)=0
        POLYOT(I)=0
10      CONTINUE
        POLYIN(1)=1.0
        POLYOT(1)=1.0

C  CREATE THE POLYNOMIAL OUTSIDE THE UNIT CIRCLE
        ICOUNT=1
20      IF(ICOUNT.LE.NO) THEN
                IF(ABS(AIMAG(ROOTOT(ICOUNT))).LT.TOL) THEN
                        POLY(1)=1.0
                        POLY(2)=-1.0/ROOTOT(ICOUNT)
                        POLY(3)=0.0
                        CALL PMULTI(POLY,1,POLYOT,ICOUNT-1,N,PW)
                        ICOUNT=ICOUNT+1
                ELSE
                        DENOM=1.0/CABS(ROOTOT(ICOUNT))**2
                        POLY(1)=1.0
                        POLY(2)=-2.0*REAL(ROOTOT(ICOUNT))*DENOM
                        POLY(3)=DENOM
                        CALL PMULTI(POLY,2,POLYOT,ICOUNT-1,N,PW)
                        ICOUNT=ICOUNT+2
                ENDIF
        ENDIF
        IF(ICOUNT.LE.NO) GOTO 20

C  CREATE THE POLYNOMIAL INSIDE THE UNIT CIRCLE
        ICOUNT=1
30      IF(ICOUNT.LE.NI) THEN
                IF(ABS(AIMAG(ROOTIN(ICOUNT))).LT.TOL) THEN
                        POLY(1)=1.0
                        POLY(2)=-1.0/ROOTIN(ICOUNT)
                        POLY(3)=0.0
                        CALL PMULTI(POLY,1,POLYIN,ICOUNT-1,N,PW)
                        ICOUNT=ICOUNT+1
                ELSE
                        DENOM=CABS(ROOTIN(ICOUNT))**2
                        POLY(1)=1.0
                        POLY(2)=-2.0*REAL(ROOTIN(ICOUNT))*DENOM
                        POLY(3)=DENOM
                        CALL PMULTI(POLY,2,POLYIN,ICOUNT-1,N,PW)    
                        ICOUNT=ICOUNT+2
                ENDIF
        ENDIF
        IF(ICOUNT.LE.NI) GOTO 30

        RETURN
        END

C***********************************************************************
C          THIS SUBROUTINE MULTIPLIES TWO POLYNOMIALS INTO A THIRD
C***********************************************************************
C        INPUT
C                  POLY1-------FIRST POLYNOMIAL IN THE MULTIPLICATION
C                  M1----------ORDER OF FIRST POLYNOMIAL
C                  POLY2-------SECOND POLYNOMIAL IN THE MULTIPLICATION
C                  M2----------ORDER OF SECOND POLYNOMIAL
C                  N-----------DIMENSION OF POLYNOMIALS
C        OUTPUT
C                  POLY3-------PRODUCT OF THE ABOVE
C***********************************************************************

      SUBROUTINE PMULT(POLY1,M1,POLY2,M2,POLY3,N)

        INTEGER M1,M2,N,LENGTH,I,J
        REAL POLY1(N),POLY2(N),POLY3(N)

C  LENGTH OF NEW POLYNOMIAL
        LENGTH=M1+M2+1

C  INITIALIZE THE RESULTANT POLYNOMIAL
        DO 10 I=1,LENGTH
                POLY3(I)=0
10      CONTINUE

C  PERFORM THE MULTIPLICATION
        DO 30 J=1,M2+1
                DO 20 I=1,M1+1
                        POLY3(I+J-1)=POLY3(I+J-1)+POLY1(I)*POLY2(J)
20              CONTINUE
30      CONTINUE

        RETURN
        END

C***********************************************************************
C          THIS SUBROUTINE MULTIPLIES TWO POLYNOMIALS IN PLACE
C***********************************************************************
C        INPUT
C                  POLY1-------FIRST POLYNOMIAL IN THE MULTIPLICATION
C                  M1----------ORDER OF FIRST POLYNOMIAL
C                  POLY2-------SECOND POLYNOMIAL IN THE MULTIPLICATION
C                  M2----------ORDER OF SECOND POLYNOMIAL
C                  N-----------DIMENSION OF POLYNOMIALS
C        OUTPUT
C                  POLY2-------PRODUCT OF THE ABOVE
C***********************************************************************

      SUBROUTINE PMULTI(POLY1,M1,POLY2,M2,N,POLY3)

        INTEGER M1,M2,N,I,J,LENGTH
        REAL POLY1(N),POLY2(N),POLY3(N)

C  LENGTH OF NEW POLYNOMIAL
        LENGTH=M1+M2+1

C  INITIALIZE THE RESULTANT POLYNOMIAL
        DO 10 I=1,LENGTH
                POLY3(I)=0
10      CONTINUE

C  PERFORM THE MULTIPLICATION
        DO 30 J=1,M2+1
                DO 20 I=1,M1+1
                        POLY3(I+J-1)=POLY3(I+J-1)+POLY1(I)*POLY2(J)
20              CONTINUE
30      CONTINUE

        DO 40 I=1,LENGTH
                POLY2(I)=POLY3(I)
40        CONTINUE

        RETURN
        END

C***********************************************************************
C          THIS SUBROUTINE TAKES THE INVERSE OF A POLYNOMIAL
C***********************************************************************
C        INPUT
C                  ITYPE-----TYPE OF POLYNOMIAL
C                       1 FOR OUTSIDE UNIT CIRCLE
C                       2 FOR INSIDE UNIT CIRCLE
C                  POLY------INPUT POLYNOMIAL    
C                  NP--------OF DEGREE HAVING NP+1 TERMS
C                  N---------DIMENSION OF POLY
C        OUTPUT
C                  COEFF-----OUTPUT POLYNOMIAL
C                  NCMAX-----NUMBER OF TERMS
C                  NC--------DIMENSION OF NC
C                  ACURCY--SMALLEST TERM OF COEF BIGGER THAN
C                  NDEG------DEGREE OF RESULTANT POLYNOMIAL
C***********************************************************************

      SUBROUTINE POLYIV(ITYPE,POLY,NP,N,
     1                     COEFF,NCMX,NC,
     2                     ACURCY,NDEG)

        INTEGER ITYPE,N,NC,NP,NCMX,NDEG,I,NCMAX,LCOEFF,NTIMES,J
        REAL POLY(N),COEFF(NC),ACURCY,DIV,ACCDIV

C  INITIALIZE THE RESULTANT POLYNOMIAL
        DO 10 I=1,NC
        COEFF(I)=0
10      CONTINUE
        COEFF(1)=1

C  LIMIT THE MAXIMUM TO THE ALLOWABLE
        NCMAX=NCMX
        IF(NCMAX.GT.(NC-NP)) NCMAX=NC-NP

C  DON'T DO ANYTHING IF THIS POLYNOMIAL IS OF ORDER ZERO
        IF(NP.EQ.0) THEN
                COEFF(1)=POLY(1)
                LCOEFF=1
        ELSE

C  WHAT TYPE OF INVERSE IS DESIRED
        IF(ITYPE.EQ.1) THEN

C  DO THE INVERSE OF A POLYNOIMIAL ON THE OUTSIDE OF UNIT CIRCLE
        I=1
        NTIMES=0
15      IF((NTIMES.LT.NP).AND.(I.LT.NCMAX)) THEN
                IF(ABS(COEFF(I)).LT.ACURCY) THEN
                        NTIMES=NTIMES+1
                ELSE
                        NTIMES=0
                ENDIF
                DO 20 J=2,NP+1
                COEFF(I+J-1)=COEFF(I+J-1)-COEFF(I)*POLY(J)
20              CONTINUE
        I=I+1
        ENDIF
        IF((NTIMES.LT.NP).AND.(I.LT.NCMAX)) GOTO 15

C  BOOKKEEPING PROCEDURES
        LCOEFF=I

        DO 30 J=1,NP
                COEFF(I+J)=0
30      CONTINUE

        ELSE

C  DO THE INVERSE OF A POLYNOMIAL ON THE INSIDE OF A UNIT CIRCLE
C  DIVIDE THE ORIGNAL POLYNOMIAL FIRST
        I=1
        DIV=POLY(NP+1)
        ACCDIV=ACURCY*ABS(DIV)
        DO 40 J=1,NP+1
        POLY(J)=POLY(J)/DIV
40      CONTINUE

        NTIMES=0
45      IF((NTIMES.LT.NP).AND.(I.LT.NCMAX)) THEN
                IF(ABS(COEFF(I)).LT.ACCDIV) THEN
                        NTIMES=NTIMES+1
                ELSE
                        NTIMES=0
                ENDIF
                DO 50 J=NP,1,-1
                COEFF(I+NP-J+1)=COEFF(I+NP-J+1)-COEFF(I)*POLY(J)
50              CONTINUE
        I=I+1
        ENDIF
        IF((NTIMES.LT.NP).AND.(I.LT.NCMAX)) GOTO 45

C  NOW DO THE BOOKKEEPING
        LCOEFF=I

        DO 60 J=1,NP
                COEFF(I+J)=0
60      CONTINUE

C  RETURN THE ORIGINAL POLYNOMIAL TO THE BEGINING
        DO 70 J=1,NP+1
        POLY(J)=POLY(J)*DIV
70      CONTINUE

C  REDUCE THE POLYNOMIAL BY THE FIRST COEFFICENT OF THE POLYNOMIAL
        DO 80 J=1,LCOEFF
        COEFF(J)=COEFF(J)/DIV
80      CONTINUE

        ENDIF

        ENDIF
        NDEG=LCOEFF-1

        RETURN
        END

C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012
