C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C                        SQSRCH.LIB
C
C  THESE ROUTINES PERFORM A SEQUENTIAL SEARCH ON MULTIPLE VARIABLES
C  AND SAVES ANY IMPORTANT REAL OR INTEGER VARIABLES FOR LATTER PRINTING
C  IF REQUIRED.  A SIMULATED ANNEALING OR GENETIC ALGORITHM SEARCH WILL
C  BE MADE IF NOT SUPPRESSED.
C
C***********************************************************************
C        REPSQS---REPTITIVELY CALL THE SEQUENTIAL SEARCH ROUTINE
C        SQSRCH---SEQUENTIAL SEARCH ROUTINETO FIND THE MAXIMUM
C        DEFPAR---DEFINE THE PARAMETERS OVER THE SEARCH GRID
C        SMANGA---SIMULATED ANNEALING AND GENETIC ALGORITHM
C***********************************************************************
C
C        FIND THE MAXIMUM MLE WITH A REPITITION TO A SEQUENTIAL SEARCH
C***********************************************************************
C        INPUT
C                NDIV---NUMBER OF INTERVALS TO SEARCH BETWEEN LIMITS
C                NREP---NUMBER OF REDUCED PARAMTER LIMIT REPETITIONS
C                ETOL---STOPPING TOLERANCE TESTED ON MAX LIKELIHOOD
C***********************************************************************

      SUBROUTINE REPSQS(NPOP,NPARM,NDIV,NREP,NSUP,
     1                  IRNGP,IRNGPB,IPOP,IDU,IPLIM,IP,
     2                  COEFF,PL,SMA,GA,RPARM,POP,PLIM,PINT,
     3                  PR1,PR2,PA,DU,SA1,SA2,DELX,TCOEFF,
     4   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     5   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,
     6   COEFR,COEFPI,THET,COEFPO,PI,PO,TI,TO,
     7   COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT,XMAX)


        INTEGER NDIV,NREP,NSUP
        INTEGER NPOP,NPARM
        INTEGER NUMCOF,I,ICOUNT
        INTEGER IRNGP(NPARM),IRNGPB(NPARM),IPOP(NPOP),IDU(NPOP)
        INTEGER IPLIM(NPARM),IP(NPARM)
        REAL XMAX,TEST0,TOTNUM,RANGE,RLOW,RHIGH
        REAL COEFF(NPARM),PL(2,NPARM),SMA(NPARM),GA(NPARM)
        REAL RPARM(NPARM,NPOP),POP(NPOP),PLIM(2,NPARM)
        REAL PINT(NPARM),PR1(NPARM),PR2(NPARM),PA(NPARM)
        REAL DU(NPARM),SA1(NPARM),SA2(NPARM)
        REAL DELX(NPARM),TCOEFF(NPARM)

        INTEGER NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO
        REAL ACURCY,SIGMIN,SIGMAX
        REAL X(NX),Z(NX),ZS(NPARM+MAXL+NX),Y(NPARM+MAXL+NX)
        REAL COEFR(MAXL),COEFPI(MAXL),THET(MAXL),COEFPO(MAXL)
        REAL PI(NPARM),PO(NPARM),TI(NPARM),TO(NPARM)
        REAL COEFCC(NPARM),PW(NPARM)
        REAL WORK(NPARM*NPARM+2*NPARM)
        COMPLEX ZPHI(NPARM),ZTHETA(NPARM)
        COMPLEX ROOTIN(NPARM),ROOTOT(NPARM)

C  CALCULATE IF A STOCHASTIC SEARCH IS REQUIRED DUE TO TIME 
C  CONSTRAINTS
        NUMCOF=0
        DO 5 I=1,NPARM
                TEST0=PLIM(2,I)-PLIM(1,I)
                IF(TEST0.NE.0) THEN
                        NUMCOF=NUMCOF+1
                ENDIF
5         CONTINUE
        TOTNUM=NDIV**NUMCOF
        IF(TOTNUM.LT.120000) NSUP=1
        
C  PERFORM THE INITIAL SEQUENTIAL SEARCH OR STOCHASTIC SEARCH
        IF(NSUP.EQ.1) THEN
           CALL SQSRCH(COEFF,PLIM,NDIV,NPARM,DELX,IPLIM,IP,TCOEFF,
     1   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     2   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     3   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT,XMAX)

        ELSE

           CALL SMANGA(NPARM,NPOP,IRNGP,IRNGPB,IPOP,IDU,
     1   COEFF,SMA,GA,RPARM,POP,PLIM,PINT,PR1,PR2,PA,DU,SA1,SA2,
     2   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     3   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     4   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT,XMAX)

        ENDIF

C  MODIFY THE PARAMETER LIMITS TO SEARCH OVER
        DO 10 I=1,NPARM
                RANGE=(PLIM(2,I)-PLIM(1,I))/NDIV
                RLOW=COEFF(I)-RANGE
                RHIGH=COEFF(I)+RANGE
                IF(RLOW.GT.PLIM(1,I)) THEN
                        PL(1,I)=RLOW
                ELSE
                        PL(1,I)=PLIM(1,I)
                ENDIF
                IF(RHIGH.LT.PLIM(2,I)) THEN
                        PL(2,I)=RHIGH
                ELSE
                        PL(2,I)=PLIM(2,I)
                ENDIF
10        CONTINUE   

C  NOW REPEAT THE SEQUENTIAL SEARCH AS OFTEN AS DESIRED
        ICOUNT=1
        IF(ICOUNT.GE.NREP) GOTO 30
15              ICOUNT=ICOUNT+1
                IF(NSUP.EQ.1) THEN
           CALL SQSRCH(COEFF,PL,NDIV,NPARM,DELX,IPLIM,IP,TCOEFF,
     1   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     2   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     3   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT,XMAX)

                ELSE

           CALL SMANGA(NPARM,NPOP,IRNGP,IRNGPB,IPOP,IDU,
     1   COEFF,SMA,GA,RPARM,POP,PL,PINT,PR1,PR2,PA,DU,SA1,SA2,
     2   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     3   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     4   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT,XMAX)

                ENDIF
                DO 20 I=1,NPARM
                        RANGE=(PL(2,I)-PL(1,I))/NDIV
                        RLOW=COEFF(I)-RANGE
                        RHIGH=COEFF(I)+RANGE
                        IF(RLOW.GT.PL(1,I)) PL(1,I)=RLOW
                        IF(RHIGH.LT.PL(2,I)) PL(2,I)=RHIGH
20                CONTINUE
        IF(ICOUNT.LT.NREP) GOTO 15        

30      RETURN
        END

C***********************************************************************
C        FIND THE MAXIMUM MLE WITH A SEQUENTIAL SEARCH
C***********************************************************************
C        INPUT
C                COEFF---COEFFICIENTS TO SEARCH THROUGH
C                PLIM----PARAMETER LIMITS
C                NDIV----NUMBER OF INTERVALS TO SEARCH BETWEEN LIMITS
C                ICOEFF--PARAMETERS TO SEARCH OVER
C***********************************************************************

      SUBROUTINE SQSRCH(COEFF,PLIM,NDIV,NPARM,DELX,IPLIM,IP,TCOEFF,
     1   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     2   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     3   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT,XMAX)

        INTEGER NPARM,NDIV,ISTOP,I
        INTEGER IPLIM(NPARM),IP(NPARM)
        REAL XMAX,POPFNC,XMLE
        REAL COEFF(NPARM),PLIM(2,NPARM),DELX(NPARM),TCOEFF(NPARM)
        INTEGER NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO
        REAL ACURCY,SIGMIN,SIGMAX
        REAL X(NX),Z(NX),ZS(NPARM+MAXL+NX),Y(NPARM+MAXL+NX)
        REAL COEFR(MAXL),COEFPI(MAXL),THET(MAXL),COEFPO(MAXL)
        REAL PI(NPARM),PO(NPARM),TI(NPARM),TO(NPARM)
        REAL COEFCC(NPARM),PW(NPARM)
        REAL WORK(NPARM*NPARM+2*NPARM)
        COMPLEX ZPHI(NPARM),ZTHETA(NPARM)
        COMPLEX ROOTIN(NPARM),ROOTOT(NPARM)

C  INITIALIZE THE PARAMETERS TO START THE SEARCH
        ISTOP=-1
        CALL DEFPAR(NPARM,PLIM,NDIV,COEFF,ISTOP,DELX,IPLIM,IP)
        DO 5 I=1,NPARM
              TCOEFF(I)=COEFF(I)
5       CONTINUE
        XMAX=POPFNC(COEFF,
     1                NPARM,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)

C  NOW SEARCH OVER THE GRID
        IF(ISTOP.NE.0) GOTO 20
10              CALL DEFPAR(NPARM,PLIM,NDIV,TCOEFF,ISTOP,DELX,IPLIM,IP)
                XMLE=POPFNC(TCOEFF,
     1                NPARM,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)

                IF(XMAX.LT.XMLE) THEN
                        XMAX=XMLE
                              DO 15 I=1,NPARM
                                     COEFF(I)=TCOEFF(I)
15                            CONTINUE
                ENDIF
        IF(ISTOP.EQ.0) GOTO 10

20      RETURN
        END

C***********************************************************************
C        DEFINES THE GRID TO PERFORM THE SEQUENTIAL SEARCH
C***********************************************************************
C        INPUT
C                COEFF---COEFFICIENTS TO SEARCH THROUGH
C                PLIM----PARAMETER LIMITS
C                NDIV----NUMBER OF INTERVALS TO SEARCH BETWEEN LIMITS
C                ICOEFF--WHICH PARAMETER TO SEARCH
C                ISTOP---FLAG WHICH INITIALIZES OR ENDS THE PARAMETER
C                        CALCULATIONS
C***********************************************************************

      SUBROUTINE DEFPAR(NPARM,PLIM,NDIV,COEFF,ISTOP,DELX,IPLIM,IP)

        INTEGER NPARM,NDIV,ISTOP,I,IP(NPARM),IPLIM(NPARM),IDSTOP
        REAL RANGE,COEFF(NPARM),PLIM(2,NPARM),DELX(NPARM)

C  INITIALIZE THE GRID PARAMETERS
        IF(ISTOP.EQ.-1) THEN
                ISTOP=0
                DO 10 I=1,NPARM
                        IP(I)=0
                        RANGE=PLIM(2,I)-PLIM(1,I)
                        DELX(I)=RANGE/NDIV
                        IF(DELX(I).GT.0) THEN
                                IPLIM(I)=NDIV
                        ELSE
                                IPLIM(I)=0
                        ENDIF
                COEFF(I)=PLIM(1,I)
10              CONTINUE

C  NOW STEP ACROSS THE GRID EVERTIME THIS ROUTINE IS CALLED
        ELSE
                I=1
                IDSTOP=0
20                      IP(I)=IP(I)+1
                        COEFF(I)=PLIM(1,I)+DELX(I)*IP(I)
                        IF(IP(I).GT.IPLIM(I)) THEN
                                IP(I)=0
                                COEFF(I)=PLIM(1,I)+DELX(I)*IP(I)
                                I=I+1
                                IF(I.GT.NPARM) THEN
                                        ISTOP=1
                                        IDSTOP=1
                                ENDIF
                        ELSE
                                IDSTOP=1
                        ENDIF
                 IF(IDSTOP.EQ.0) GOTO 20
        ENDIF

        RETURN
        END

C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012

C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C                        SCALE.LIB
C
C
C  THESE ROUTINES FIND THE OPTIMUM SCALING PARAMETER
C
C***********************************************************************
C        SEQUEN--PERFORMS A SEQUENCTIAL SEARCH
C        BISECT--MODIFIED BISECTION TECHNIQUE
C        TOMAX---THE FUNCTION TO MAXIMIZE, TRUNCATED LOG LIKELIHOOD
C        FSIG----THE SCALED PROBABILITY DENSITY FUNCTION
C***********************************************************************
C
C        THIS ROUTINE PERFORMS A SEQUENTIAL SEARCH USING 10 STEPS
C***********************************************************************
C        INPUT
C                BIG----ENDPOINT OF SEARCH
C                SMALL--BEGINING OF SEARCH
C                Z------SEQUENCE TO SEARCH
C                NZ-----DIMENSION OF SEQUENCE VECTOR
C                IQ-----TRUNCATION OF SEQUENCE AT BOTH ENDS
C***********************************************************************

      FUNCTION SEQUEN(BIG,SMALL,Z,NZ,IQ)
        INTEGER NZ,IQ,NS,I,ISTEP
        REAL Z(NZ),BIG,SMALL,WIDTH,STEP,SIGSTP
        REAL XMAX,TOMAX,TEST,ULIM,BLIM,SEQUEN,BISECT
        PARAMETER(NS=10)

C  FIND THE INTERVAL WHERE THE MAXIMUM OCCURS
        WIDTH=BIG-SMALL
        STEP=WIDTH/NS
        ISTEP=0
        TEST=TOMAX(Z,NZ,IQ,SMALL)
        DO 10 I=1,NS
                SIGSTP = SMALL + STEP * I
                XMAX=TOMAX(Z,NZ,IQ,SIGSTP)
                IF(XMAX.GT.TEST) THEN
                        TEST=XMAX
                        ISTEP=I
                ENDIF
10      CONTINUE

C  REFINE THE RESULTS WITH A BISECTION CALL
        IF(ISTEP.EQ.0) THEN
               ULIM=SMALL+STEP
               BLIM=SMALL
        ELSEIF(ISTEP.EQ.NS) THEN
               ULIM=BIG
               BLIM=BIG-STEP
        ELSE
               ULIM = SMALL + STEP * (ISTEP + 1)
               BLIM = SMALL + STEP * (ISTEP - 1)
        ENDIF
        SEQUEN = BISECT(BLIM,ULIM,Z,NZ,IQ)
                
        RETURN
        END

C***********************************************************************
C        THIS ROUTINE PERFORMS A MODIFIED BISECTION USING NO MORE
C        THAN 20 ITERATIONS
C***********************************************************************
C        INPUT
C                SH----ENDPOINT OF SEARCH
C                SL----BEGINING OF SEARCH
C                Z------SEQUENCE TO SEARCH
C                NZ-----DIMENSION OF SEQUENCE VECTOR
C                IQ-----TRUNCATION OF SEQUENCE AT BOTH ENDS
C***********************************************************************

      FUNCTION BISECT(SL,SH,Z,NZ,IQ)
        INTEGER NZ,IQ,NITER,ICOUNT
        REAL SL,SH,Z(NZ),TOL,S1,S2,FZ1,TOMAX,FZ2,DEL,HALF,FZH,BISECT
        PARAMETER(TOL=1.0E-6,NITER=20)

C  INIITIALIZE THE BISECTION PARAMETERS
        S1=SL
        S2=SH
        IF(S1.LE.0.0) THEN
                S1=TOL
        ENDIF
        FZ1=TOMAX(Z,NZ,IQ,S1)
        FZ2=TOMAX(Z,NZ,IQ,S2)
        ICOUNT=0
        DEL=2*TOL

C  FIND THE SCALE PARAMTER WHICH MAXIMUMIZES THE LOG-LIKLIHOOD SUM
10      IF ((ICOUNT.LT.NITER).AND.(DEL.GT.TOL)) THEN
                ICOUNT=ICOUNT+1
                HALF=(S2+S1)/2.0
                FZH=TOMAX(Z,NZ,IQ,HALF)
                IF(FZ1.GE.FZ2) THEN
                        FZ2=FZH
                        S2=HALF
                ELSE
                        FZ1=FZH
                        S1=HALF
                ENDIF
                DEL=ABS(FZ2-FZ1)
        ENDIF
        IF ((ICOUNT.LT.NITER).AND.(DEL.GT.TOL)) GOTO 10
        BISECT=HALF

        RETURN
        END

C***********************************************************************
C        THIS ROUTINE FINDS THE TRUNCATED LOG LIKELIHOOD
C***********************************************************************
C        INPUT
C                Z------SEQUENCE
C                NZ-----DIMENSION OF SEQUENCE VECTOR
C                IQ-----TRUNCATION OF SEQUENCE AT BOTH ENDS
C                SIG----SCALE PARAMETER
C***********************************************************************
 
      FUNCTION TOMAX(Z,NZ,IQ,SIG)
        INTEGER NZ,IQ,I
        REAL Z(NZ),SIG,SUM,FSIG,TOMAX

        SUM=0
        DO 10 I=IQ+1,NZ-IQ
                SUM=SUM+ALOG(FSIG(Z(I),SIG))
10      CONTINUE

        TOMAX=SUM

        RETURN
        END

C***********************************************************************
C        THIS ROUTINE PERFORMS THE SCALED DENSITY FUNCTION
C***********************************************************************
C        INPUT
C                X------SEQUENCE
C                SIG----SCALE PARAMETER
C***********************************************************************

      FUNCTION FSIG(X,SIG)
      REAL X,SIG,FSIG,DNSITY

        FSIG=1.0/SIG*DNSITY(X/SIG)

        RETURN
        END

C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012

C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C                        SIMGEN.LIB
C
C  NOTE:  SUBROUTINE SMANGA IS THE INTERFACE TO THE SIMULATED ANNEALING
C          AND GENETIC ALGORITHM ROUTINE.  THE SETTINGS HAVE BEEN
C          OPTIMIZED FOR THE NON-GAUSSIAN NGMLE ROUTINES
C
C***********************************************************************

      SUBROUTINE SMANGA(NPARM,NPOP,IRNGP,IRNGPB,IPOP,IDU,
     1   COEFF,SMA,GA,RPARM,POP,PLIM,PINT,PR1,PR2,PA,DU,SA1,SA2,
     2   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     3   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     4   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT,XMAX)

        INTEGER NPOP,NPARM,NGEN,NPROP,NMUT
        INTEGER NTAUSM,ISEED,NSMAGA,INIT,I
        INTEGER IRNGP(NPARM),IRNGPB(NPARM),IPOP(NPOP),IDU(NPOP)
        REAL XMAX,XSMA,POPFNC,XGA,PRBSIM,ZTSIM,DELSIM
        REAL COEFF(NPARM),SMA(NPARM),GA(NPARM)
        REAL RPARM(NPARM,NPOP),POP(NPOP),PLIM(2,NPARM)
        REAL PINT(NPARM),PR1(NPARM),PR2(NPARM),PA(NPARM)
        REAL DU(NPARM),SA1(NPARM),SA2(NPARM)

        INTEGER NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO
        REAL ACURCY,SIGMIN,SIGMAX
        REAL X(NX),Z(NX),ZS(NPARM+MAXL+NX),Y(NPARM+MAXL+NX)
        REAL COEFR(MAXL),COEFPI(MAXL),THET(MAXL),COEFPO(MAXL)
        REAL PI(NPARM),PO(NPARM),TI(NPARM),TO(NPARM)
        REAL COEFCC(NPARM),PW(NPARM)
        REAL WORK(NPARM*NPARM+2*NPARM)
        COMPLEX ZPHI(NPARM),ZTHETA(NPARM)
        COMPLEX ROOTIN(NPARM),ROOTOT(NPARM)

        SAVE INIT
        DATA INIT /0/

        IF(INIT.EQ.0) THEN
                INIT   = 1
                NGEN   = 600
                NPROP  = 20
                NMUT   = 2
                ZTSIM  = 1.0 
                DELSIM = 0.05
                NTAUSM = 100
                ISEED = -(4*17*19*23*29+1)
                DO 10 I=1,NPARM
                        IRNGPB(I)=15
10              CONTINUE
        ENDIF

        NSMAGA=1
              CALL SIMGEN(SMA,NPOP,NPARM,NGEN,NPROP,NMUT,
     1                    PLIM,IRNGP,IRNGPB,RPARM,POP,IPOP,
     2                    PRBSIM,ZTSIM,DELSIM,NTAUSM,NSMAGA,ISEED,
     3                    PINT,PR1,PR2,PA,DU,SA1,SA2,IDU,
     4   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     5   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     6   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT)

        XSMA=POPFNC(SMA,
     1      NPARM,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)

        NSMAGA=2
              CALL SIMGEN(GA,NPOP,NPARM,NGEN,NPROP,NMUT,
     1                    PLIM,IRNGP,IRNGPB,RPARM,POP,IPOP,
     2                    PRBSIM,ZTSIM,DELSIM,NTAUSM,NSMAGA,ISEED,
     3                    PINT,PR1,PR2,PA,DU,SA1,SA2,IDU,
     4   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     5   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     6   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT)

        XGA = POPFNC(GA,
     1      NPARM,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)

        IF(XSMA.GT.XGA) THEN
                XMAX=XSMA
                DO 30 I=1,NPARM
                        COEFF(I)=SMA(I)
30              CONTINUE
        ELSE
                XMAX=XGA
                DO 40 I=1,NPARM
                        COEFF(I)=GA(I)
40              CONTINUE
        ENDIF

        RETURN
        END

C**********************************************************************
C       SIMULATED ANNEALING & GENETIC ALGORITHM PROGRAM VERSION 3
C               WRITTEN BY MARK E.LEHR
C                          UCR STAT DEPT
C                UPDATED
C                       JULY 29, 1995
C       THESE ALGORITHMS WERE DESCRIBED IN
C               SCIENTIFIC AMERICA, JULY 92--GENETIC ALGORITHMS
C                       BY JOHN H. HOLLAND
C               SIMULATED ANEALING, THEORY AND APPLICATIONS
C                       BY P.J.M. VANHAARHOVEN &
C                          E.H.L. AARTS
C       THE FOLLOWING PROGRAM IS AN INTERPRETATION/IMPLEMENTATION
C       OF THE ALGORITHMS DESCIRBED IN THOSE ARTICLES/BOOKS
C
C       TO ALL USERS
C               PLEASE PROVIDE ANY COMMENTS ON THE
C               UTILIZATION, CODING, OR IMPLEMENTATION
C               OF THIS ALGORITHM TO YOUR'S TRULY
C
C       HISTORICAL NOTE:
C               THIS PROGRAM GREW FROM THE 7TH VERSION OF THE 
C               GENETIC ALGORITHM IN DEVELOPMENT SINCE MARCH 93.
C               THE FIRST VERSION OF THE SIMULATED ANNEALING
C               ALGORITHM IS CONTAINED HEREIN.
C        THIS VERSION
C                USED IN THE ARTICLE TO BE PUBLISHED ON NON-
C                GAUSSIAN ARMA ESTIMATION
C**********************************************************************

C        PROGRAM SIMGEN_VERSION_3

      SUBROUTINE SIMGEN(COEFF,NPOP,NPARM,NGEN,NPROP,NMUT,
     1                   PLIM,IRNGP,IRNGPB,RPARM,POP,IPOP,
     2                   PRBSIM,ZTSIM,DELSIM,NTAUSM,NSMAGA,ISEED,
     3                   PINT,PR1,PR2,PA,DU,SA1,SA2,IDU,
     4   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     5   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     6   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT)


        INTEGER NPOP,NPARM,NGEN,NPROP,NMUT,NTAUSM,NSMAGA,ISEED,I
        INTEGER IRNGP(NPARM),IRNGPB(NPARM),IPOP(NPOP),IDU(NPOP)
        REAL PRBSIM,ZTSIM,DELSIM
        REAL PLIM(2,NPARM),RPARM(NPARM,NPOP),POP(NPOP)
        REAL PINT(NPARM),PR1(NPARM),PR2(NPARM),PA(NPARM)
        REAL DU(NPARM),SA1(NPARM),SA2(NPARM),COEFF(NPARM)

        INTEGER NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO
        REAL ACURCY,SIGMIN,SIGMAX
        REAL X(NX),Z(NX),ZS(NPARM+MAXL+NX),Y(NPARM+MAXL+NX)
        REAL COEFR(MAXL),COEFPI(MAXL),THET(MAXL),COEFPO(MAXL)
        REAL PI(NPARM),PO(NPARM),TI(NPARM),TO(NPARM)
        REAL COEFCC(NPARM),PW(NPARM)
        REAL WORK(NPARM*NPARM+2*NPARM)
        COMPLEX ZPHI(NPARM),ZTHETA(NPARM)
        COMPLEX ROOTIN(NPARM),ROOTOT(NPARM)

C  INITIALIZE THE POPULATION UNLESS RESTART IS DESIRED        
        CALL INTPOP(NPOP,NPARM,
     1              PLIM,IRNGP,IRNGPB,RPARM,POP,IPOP,ISEED,PINT,
     2   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     3   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     4   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT)

C  DIFFERENT LOOPS FOR SIMULATED ANNEALING
C  OR THE GENETIC ALGORITHM
        
        IF(NSMAGA.EQ.1) THEN

C  INNER LOOP USED IN THE GENERATION PROCESS
        DO 10 I=1,NGEN,1
                CALL SIMANL(I,NPOP,NPARM,
     1                    PLIM,IRNGPB,RPARM,POP,
     2                    PRBSIM,ZTSIM,DELSIM,NTAUSM,SA1,SA2,
     3   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     4   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     5   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT)

                CALL SORTIT(POP,IPOP,NPOP,NPOP)
10      CONTINUE
        
C  START THE GENETIC ALGORITHM LOOP

        ELSE        

C  INNER LOOP USED IN THE GENERATION PROCESS
        DO 20 I=1,NGEN,1
                CALL PROPOP(NPOP,NPARM,NPROP,
     1                    PLIM,IRNGP,IRNGPB,RPARM,POP,IPOP,PR1,PR2,
     2   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     3   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     4   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT)

                CALL PARMUT(NPOP,NPARM,NPROP,NMUT,
     1                    PLIM,IRNGP,RPARM,POP,IPOP,PA,
     2   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     3   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     4   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT)

                CALL DUPMUT(NPOP,NPARM,NPROP,
     1                    PLIM,IRNGP,RPARM,POP,IPOP,DU,IDU,
     2   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     3   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     4   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT)

20      CONTINUE
        ENDIF

C  NOW JUST FINISH WRITING THE RESULTS TO VARIOUS FILES
        CALL SORTIT(POP,IPOP,NPOP,NPOP)

        DO 30 I=1,NPARM
                COEFF(I)=RPARM(I,IPOP(I))
30      CONTINUE
        
        RETURN        
        END

C**********************************************************************
C       THE INITILIALIZATION OF THE POPULATION IS ACCOMPLISHED
C       BY DEFINING THE PARAMETERS UNIFORMLY BETWEEN THE LIMITS
C       SPECIFIED IN THE INPUT DECK
C**********************************************************************

        SUBROUTINE INTPOP(NPOP,NPARM,
     1                    PLIM,IRNGP,IRNGPB,RPARM,POP,IPOP,ISEED,V,
     2   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     3   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     4   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT)

 
        INTEGER NPOP,NPARM,ISEED,INIT,IDUM,I,J,IXRAND
        INTEGER IRNGP(NPARM),IRNGPB(NPARM),IPOP(NPOP)
        REAL XUNIF,SUNIF,A,B,C,D,POPFNC
        REAL V(NPARM),PLIM(2,NPARM),RPARM(NPARM,NPOP),POP(NPOP)

        INTEGER NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO
        REAL ACURCY,SIGMIN,SIGMAX
        REAL X(NX),Z(NX),ZS(NPARM+MAXL+NX),Y(NPARM+MAXL+NX)
        REAL COEFR(MAXL),COEFPI(MAXL),THET(MAXL),COEFPO(MAXL)
        REAL PI(NPARM),PO(NPARM),TI(NPARM),TO(NPARM)
        REAL COEFCC(NPARM),PW(NPARM)
        REAL WORK(NPARM*NPARM+2*NPARM)
        COMPLEX ZPHI(NPARM),ZTHETA(NPARM)
        COMPLEX ROOTIN(NPARM),ROOTOT(NPARM)

        DATA INIT/0/

C  INITIALIZE THE RANDOM NUMBER GENERATORS

        IF(INIT.EQ.0) THEN
                INIT=1
                IDUM=-ISEED
                XUNIF=SUNIF(IDUM)
        ENDIF

C  PARAMETER DEFINITION USING A UNIFORM DISTRIBUTION
        DO 10 I=1,NPARM,1
        IRNGP(I)=2**IRNGPB(I)
        A=PLIM(1,I)
        B=PLIM(2,I)
        D=IRNGP(I)
        DO 5 J=1,NPOP
        C=IXRAND(0,IRNGP(I))
        RPARM(I,J)=A+(B-A)*C/D
5       CONTINUE
10      CONTINUE

C  DETERMINE THE POPULATION FUNCTION VALUES FOR ALL MEMBERS
C  OF THE POPULATION
        DO 20 I=1,NPOP,1
                DO 15 J=1,NPARM
                V(J)=RPARM(J,I)
15              CONTINUE
        POP(I)=POPFNC(V,
     1      NPARM,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)

20      CONTINUE

C  INITIALIZE THE POINTER INDEX OF THE POPULATION AND
C  PARAMTERS
        DO 25 I=1,NPOP,1
                IPOP(I)=I
25      CONTINUE

        RETURN
        END

C**********************************************************************
C       ONLY ALLOW THOSE MEMBERS OF THE POPULATION TO CHANGE
C       IF THEY MEET THE SIMULATED ANNEALING REQUIREMENTS
C**********************************************************************

        SUBROUTINE SIMANL(IGEN,NPOP,NPARM,
     1                    PLIM,IRNGPB,RPARM,POP,
     2                    PRBSIM,ZTSIM,DELSIM,NTAUSM,DELPRM,V3,
     3   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     4   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     5   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT)


        INTEGER NPOP,NPARM,NTAUSM,IXRAND,K,I,J,IGEN
        INTEGER IRNGPB(NPARM)
        REAL PRBSIM,ZTSIM,DELSIM,PROB,OLD,XNEW,ACCEPT,DELTA,POPFNC
        REAL PLIM(2,NPARM),RPARM(NPARM,NPOP),POP(NPOP)
        REAL DELPRM(NPARM),V3(NPARM)

        INTEGER NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO
        REAL ACURCY,SIGMIN,SIGMAX
        REAL X(NX),Z(NX),ZS(NPARM+MAXL+NX),Y(NPARM+MAXL+NX)
        REAL COEFR(MAXL),COEFPI(MAXL),THET(MAXL),COEFPO(MAXL)
        REAL PI(NPARM),PO(NPARM),TI(NPARM),TO(NPARM)
        REAL COEFCC(NPARM),PW(NPARM)
        REAL WORK(NPARM*NPARM+2*NPARM)
        COMPLEX ZPHI(NPARM),ZTHETA(NPARM)
        COMPLEX ROOTIN(NPARM),ROOTOT(NPARM)

C  CALCULATE THE PROBABILITY OF ACCEPTING A
C  DECREASING POPULATION FUNCTION FOR 
C  SIMULATED ANNNEALING

        PRBSIM=ZTSIM*EXP(-1.0*(IGEN-1)/NTAUSM)

C  MODIFY EACH MEMBER IN THE POPULATION ACCORDING
C  TO THE RULES OF SIMULATED ANNEALING

        DO 20 J=1,NPOP,1

C  CALCULATE PERCENTAGE CHANGE FOR EACH PARAMETER
        
        I=IXRAND(1,NPARM)
        DELTA=DELSIM/IRNGPB(I)
        DELPRM(I)=DELTA*IXRAND(-IRNGPB(I),IRNGPB(I))
        DO 10 K=1,NPARM
                V3(K)=RPARM(K,J)
10      CONTINUE

C  DETERIME THE CHANGE IN THE PARAMETER AND LIMIT

        V3(I)=DELPRM(I)*(PLIM(2,I)-PLIM(1,I))+RPARM(I,J)
        IF     (V3(I).GT.PLIM(2,I)) THEN
                V3(I)=PLIM(2,I)
        ELSEIF (V3(I).LT.PLIM(1,I)) THEN
                V3(I)=PLIM(1,I)
        ENDIF

C  COMPARE THE CHANGE IN THE POPULATION FUNCTION VALUE
        
        ACCEPT=0
        XNEW=POPFNC(V3,
     1      NPARM,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)

        OLD=POP(J)
        IF(XNEW.GE.OLD) ACCEPT=1

C  ALSO COMPARE TO SIMULATED ANNEALING CRITERIA

        PROB=IXRAND(1,32000)/32000.0
        IF(PROB.LT.PRBSIM) ACCEPT=1

C  NOW CHANGE THE POPULATION BY THE ABOVE
C  ACCEPTED CRITERIA

        IF(ACCEPT.EQ.1) THEN
        POP(J)=XNEW
        DO 15 I=1,NPARM,1
                RPARM(I,J)=V3(I)
15      CONTINUE
        ENDIF
20      CONTINUE

        RETURN
        END

C**********************************************************************
C       ONLY ALLOW THOSE MEMBERS OF THE POPULATION MEETING
C       THE FITNESS REQUIREMENTS TO PROPAGATE
C**********************************************************************

        SUBROUTINE PROPOP(NPOP,NPARM,NPROP,
     1                    PLIM,IRNGP,IRNGPB,RPARM,POP,IPOP,V3,V4,
     2   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     3   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     4   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT)


        INTEGER NPOP,NPARM,NPROP
        INTEGER IR1MOD,IR2MOD,IR3,IR4,J,IXRAND,IWPARM,INST,I,INSI
        INTEGER IRNGP(NPARM),IRNGPB(NPARM),IPOP(NPOP)
        REAL POPFNC,P1,P2,P3,P4,A1,A2,SPLIT,SPARM,B,C
        REAL PLIM(2,NPARM),RPARM(NPARM,NPOP),POP(NPOP)
        REAL V3(NPARM),V4(NPARM)

        INTEGER NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO
        REAL ACURCY,SIGMIN,SIGMAX
        REAL X(NX),Z(NX),ZS(NPARM+MAXL+NX),Y(NPARM+MAXL+NX)
        REAL COEFR(MAXL),COEFPI(MAXL),THET(MAXL),COEFPO(MAXL)
        REAL PI(NPARM),PO(NPARM),TI(NPARM),TO(NPARM)
        REAL COEFCC(NPARM),PW(NPARM)
        REAL WORK(NPARM*NPARM+2*NPARM)
        COMPLEX ZPHI(NPARM),ZTHETA(NPARM)
        COMPLEX ROOTIN(NPARM),ROOTOT(NPARM)

C  SORTIT TO FIND THOSE MEMBERS OF THE POPULATION WHICH ARE
C  ALLOWED TO PROPAGATE
        CALL SORTIT(POP,IPOP,NPOP,NPROP)
        
C  RANDOMIZE THE MATING AMONGST THOSE ALLOWED TO PROPAGATE
        CALL RANMAT(NPOP,NPROP,IPOP)

C  PROPAGATE BY TWO'S AND REMOVE THE LEAST SIGNIFICANT
C  MEMBERS OF THE POPULATION
C  THE NEW OFFSPRING ARE DETERMINED BY CROSSOVER
        INSI=NPOP-NPROP
        DO 30 I=1,NPROP,2
        INST=INSI+I

C  WHICH PARAMETER AND WHERE TO SPLIT
        IWPARM=IXRAND(1,NPARM)
        SPARM=IXRAND(1,IRNGPB(IWPARM))
        SPLIT=2**SPARM             

C  DETERIME THE PERCENTAGE OF SPLIT
        A1=RPARM(IWPARM,IPOP(I  ))
        A2=RPARM(IWPARM,IPOP(I+1))
        B=PLIM(1,IWPARM)
        C=PLIM(2,IWPARM)
        IF(C.EQ.B) THEN
                P1=0
                P2=0
        ELSE
                P1=(A1-B)/(C-B)
                P2=(A2-B)/(C-B)
        ENDIF

C  DETERMINE NEW PERCENTAGE BY CROSSOVER
        IR1MOD=MOD(IRNGP(IWPARM)*P1,SPLIT)
        IR2MOD=MOD(IRNGP(IWPARM)*P2,SPLIT)
        IR3=IRNGP(IWPARM)*P1-IR1MOD+IR2MOD
        IR4=IRNGP(IWPARM)*P2-IR2MOD+IR1MOD
        P3=REAL(IR3)/IRNGP(IWPARM)
        P4=REAL(IR4)/IRNGP(IWPARM)

C  PUT THE UNCHANGED PARAMETERS UPTO THE MODIFIED PARAMETER
C  INTO THE CHILDRENS GENES
        DO 10 J=1,IWPARM,1
                RPARM(J,IPOP(INST  ))=RPARM(J,IPOP(I  ))
                V3(J)=RPARM(J,IPOP(INST))
                RPARM(J,IPOP(INST+1))=RPARM(J,IPOP(I+1))
                V4(J)=RPARM(J,IPOP(INST+1))
10      CONTINUE

C  DO THE SAME FOR THE NEW PARAMETERS
        RPARM(IWPARM,IPOP(INST  ))=B+(C-B)*P3
        V3(IWPARM)=RPARM(IWPARM,IPOP(INST))
        RPARM(IWPARM,IPOP(INST+1))=B+(C-B)*P4
        V4(IWPARM)=RPARM(IWPARM,IPOP(INST+1))
                
C  NOW FINISH THE PROCESS FOR THOSE PARAMETERS ABOVE THE
C  MODIFIED PARAMTER
        DO 20 J=IWPARM+1,NPARM,1
                RPARM(J,IPOP(INST  ))=RPARM(J,IPOP(I  ))
                V3(J)=RPARM(J,IPOP(INST))
                RPARM(J,IPOP(INST+1))=RPARM(J,IPOP(I+1))
                V4(J)=RPARM(J,IPOP(INST+1))
20      CONTINUE

C  DEFINE THE POPULATION RANKING FOR THE NEW CHILDREN
        POP(IPOP(INST  ))=POPFNC(V3,
     1      NPARM,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)


        POP(IPOP(INST+1))=POPFNC(V4,
     1      NPARM,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)


30      CONTINUE
        
        RETURN
        END

C**********************************************************************
C       THIS ROUTINE ALLOWS RANDOM MATING OF THOSE WHICH ARE
C       FIT TO MATE AND IS SELF EXPLANATORY
C**********************************************************************

        SUBROUTINE RANMAT(NPOP,NPROP,IPOP)

        INTEGER NPOP,NPROP,IPOP(NPOP),I,IRANDM,IXRAND,TEMP

        DO 10 I=1,NPROP
               IRANDM=IXRAND(1,NPROP)
               TEMP=IPOP(I)
               IPOP(I)=IPOP(IRANDM)
               IPOP(IRANDM)=TEMP
10      CONTINUE

        RETURN
        END

C**********************************************************************
C       MUTATE THOSE CAPABLE OF PROPAGATION BY ONE PARAMETER
C       THIS NUMBER IS SPECIFIED IN THE INPUT DECK
C       BITWISE MUTATION IS UNNECESSARY
C**********************************************************************

        SUBROUTINE PARMUT(NPOP,NPARM,NPROP,NMUT,
     1                    PLIM,IRNGP,RPARM,POP,IPOP,V,
     2   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     3   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     4   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT)

        INTEGER NPOP,NPARM,NPROP,NMUT,I,J
        INTEGER IXRAND,INSPOP,INSPRP,INPOP,INPROP
        INTEGER IRNGP(NPARM),IPOP(NPOP)
        REAL POPFNC,A,B,C,D
        REAL V(NPARM),POP(NPOP),RPARM(NPARM,NPOP),PLIM(2,NPARM)

        INTEGER NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO
        REAL ACURCY,SIGMIN,SIGMAX
        REAL X(NX),Z(NX),ZS(NPARM+MAXL+NX),Y(NPARM+MAXL+NX)
        REAL COEFR(MAXL),COEFPI(MAXL),THET(MAXL),COEFPO(MAXL)
        REAL PI(NPARM),PO(NPARM),TI(NPARM),TO(NPARM)
        REAL COEFCC(NPARM),PW(NPARM)
        REAL WORK(NPARM*NPARM+2*NPARM)
        COMPLEX ZPHI(NPARM),ZTHETA(NPARM)
        COMPLEX ROOTIN(NPARM),ROOTOT(NPARM)

                
C  SORTIT AND RANOMLY CHOOSE THOSE OFFSPRING TO MUTATE
        CALL SORTIT(POP,IPOP,NPOP,NPROP)
        CALL RANMAT(NPOP,NPROP,IPOP)

C  AGAIN PUT THE MUTATIONS AT THE POPULATION BOTTOM UNTIL
C  SORTITING HAS OCCURED
        INPROP=NPROP-NMUT
        INPOP=NPOP-NMUT
        DO 10 J=1,NMUT                          
               INSPRP=INPROP+J
               INSPOP=INPOP+J
               DO  5 I=1,NPARM
                 RPARM(I,IPOP(INSPOP))=RPARM(I,IPOP(INSPRP))
5              CONTINUE
10      CONTINUE

C  DO THE SAME FOR THE MUTATED PARAMTER
        DO 15 J=1,NMUT
               INSPOP=INPOP+J
               I=IXRAND(1,NPARM)
               A=PLIM(1,I)
               B=PLIM(2,I)
               D=IRNGP(I)
               C=IXRAND(0,IRNGP(I))
               RPARM(I,IPOP(INSPOP))=A+(B-A)*C/D
15      CONTINUE

C  ESTABLISH THE POPULATION RANKING OF THE MUTATION
        DO 25 I=1,NMUT,1
                INSPOP=INPOP+I
                DO 20 J=1,NPARM
                       V(J)=RPARM(J,IPOP(INSPOP))
20              CONTINUE
        POP(IPOP(INSPOP))=POPFNC(V,
     1      NPARM,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)

25      CONTINUE

        RETURN
        END

C**********************************************************************
C       ANY POPULATION MEMBERS WHICH ARE DUPLICATES SHOULD
C       BE MUTATED OR EXECUTED ON THE SPOT DUE TO THE LACK
C       OF ANY ADDITIONAL INFORMATION PROVIDED
C**********************************************************************

        SUBROUTINE DUPMUT(NPOP,NPARM,NPROP,
     1                    PLIM,IRNGP,RPARM,POP,IPOP,V,ISAME,
     2   NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO,
     3   ACURCY,SIGMIN,SIGMAX,X,Z,ZS,Y,COEFR,COEFPI,THET,COEFPO,
     4   PI,PO,TI,TO,COEFCC,PW,WORK,ZPHI,ZTHETA,ROOTIN,ROOTOT)

        INTEGER NPOP,NPARM,NPROP,IXRAND,J,I,NSAME
        INTEGER ISAME(NPOP),IRNGP(NPARM),IPOP(NPOP)
        REAL A,B,C,D,POPFNC
        REAL V(NPARM),RPARM(NPARM,NPOP),POP(NPOP),PLIM(2,NPARM)

        INTEGER NX,NZ,MAXL,IQ,NCPHI,NCTHET,NPI,NPO,NTI,NTO
        REAL ACURCY,SIGMIN,SIGMAX
        REAL X(NX),Z(NX),ZS(NPARM+MAXL+NX),Y(NPARM+MAXL+NX)
        REAL COEFR(MAXL),COEFPI(MAXL),THET(MAXL),COEFPO(MAXL)
        REAL PI(NPARM),PO(NPARM),TI(NPARM),TO(NPARM)
        REAL COEFCC(NPARM),PW(NPARM)
        REAL WORK(NPARM*NPARM+2*NPARM)
        COMPLEX ZPHI(NPARM),ZTHETA(NPARM)
        COMPLEX ROOTIN(NPARM),ROOTOT(NPARM)

C  SORTIT AND CHOOSE THOSE DUPLICATES FOR ELIMINATION FROM
C  THOSE ALLOWED TO PROPAGATE
        CALL SORTIT(POP,IPOP,NPOP,NPROP)

C  FIND THE DUPLICATES
        NSAME=0
        DO 5 I=2,NPROP,1
              ISAME(I)=0
              IF(POP(IPOP(I-1)).EQ.POP(IPOP(I))) THEN
                ISAME(I)=IPOP(I)
                NSAME=NSAME+1
              ENDIF
5       CONTINUE

C  MUTATE THE DUPLICATES AND PLACE THE PARAMETERS AT THE
C  BOTTOM OF THE POPULATION STACK
        DO 15 I=1,NPARM,1
               A=PLIM(1,I)
               B=PLIM(2,I)
               D=IRNGP(I)
               DO 10 J=2,NPROP,1
                      IF(ISAME(J).NE.0) THEN
                           C=IXRAND(0,IRNGP(I))
                           RPARM(I,ISAME(J))=A+(B-A)*C/D
                      ENDIF
10             CONTINUE
15      CONTINUE

C  PERFORM POPULATION RANKING
        DO 25 I=1,NPROP,1
                IF(ISAME(I).NE.0) THEN
                      DO 20 J=1,NPARM
                             V(J)=RPARM(J,ISAME(I))
20                    CONTINUE

            POP(ISAME(I))=POPFNC(V,
     1          NPARM,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)
 
                ENDIF
25      CONTINUE

        RETURN
        END

C**********************************************************************
C       FIND THE MAXIMUM AND MINIMUM IN A REAL VALUED ARRAY
C**********************************************************************

        SUBROUTINE RMXMIN(RARRAY,IARRAY,IDIM,ISTART,IMAX,IMIN)

        INTEGER IDIM,ISTART,IMAX,IMIN,ISP1,IEND,I
        REAL MAX,MIN,RARRAY(IDIM),TEST
        INTEGER IARRAY(IDIM)

        MAX=RARRAY(IARRAY(ISTART))
        MIN=MAX
        IMAX=ISTART
        IMIN=IMAX
        ISP1=ISTART+1
        IEND=IDIM-ISTART+1

        DO 10 I=ISP1,IEND,1
                TEST=RARRAY(IARRAY(I))
                IF(TEST.GT.MAX) THEN
                        IMAX=I
                        MAX=TEST
                ELSEIF(TEST.LT.MIN) THEN
                        IMIN=I
                        MIN=TEST
                ENDIF
10      CONTINUE

        RETURN
        END

C**********************************************************************
C       SORTIT AN ARRAY FROM THE TOP AND BOTTOM SIMULTANEOUSLY
C       STOPPING AT A VALUE SPECIFIED
C**********************************************************************

        SUBROUTINE SORTIT(RARRAY,IARRAY,IDIM,NORDER)

        INTEGER IDIM,NORDER,ILOOP,IDIM2,I,IFRONT,IEND
        REAL RARRAY(IDIM)
        INTEGER IARRAY(IDIM)
        INTEGER IMAX,IMIN,ITEMP

        ILOOP=NORDER
        IDIM2=IDIM/2
        IF(NORDER.GT.IDIM2) THEN
                ILOOP=IDIM2
        ENDIF
        DO 10 I=1,ILOOP,1
                IFRONT=I
                IEND=IDIM-I+1
                CALL RMXMIN(RARRAY,IARRAY,IDIM,IFRONT,IMAX,IMIN)
                IF     ((IMIN.EQ.IFRONT).AND.(IMAX.EQ.IEND)) THEN
                        ITEMP=IARRAY(I)
                        IARRAY(I)=IARRAY(IMAX)
                        IARRAY(IMAX)=ITEMP
                ELSEIF (IMIN.EQ.IFRONT)  THEN
                        ITEMP=IARRAY(IEND)
                        IARRAY(IEND)=IARRAY(IMIN)
                        IARRAY(IMIN)=ITEMP
                        ITEMP=IARRAY(I)
                        IARRAY(I)=IARRAY(IMAX)
                        IARRAY(IMAX)=ITEMP
                ELSE                                          
                        ITEMP=IARRAY(I)
                        IARRAY(I)=IARRAY(IMAX)
                        IARRAY(IMAX)=ITEMP
                        ITEMP=IARRAY(IEND)
                        IARRAY(IEND)=IARRAY(IMIN)
                        IARRAY(IMIN)=ITEMP
                ENDIF
10      CONTINUE
        
        RETURN
        END

C**********************************************************************
C        DEVELOP AN INTEGER RANDOM NUMBER
C**********************************************************************

      FUNCTION IXRAND(ISTART,ISTOP)

        REAL PROB,SUNIF
        INTEGER ISTART,ISTOP,IRANGE,IXRAND,IDM
        DATA IDM/-1/

        PROB=SUNIF(IDM)

        IRANGE=ISTOP-ISTART+1

        IXRAND=ISTART+IRANGE*PROB

        RETURN
        END

C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012
