
!************************************************************
!***Nonlinear Function Optimization (Minimization) Routine***
!************************************************************
	SUBROUTINE COMPLEXOPTIMIZE(NPOINT,NPARM,PARMINIT,PARMDELTA,PARMMIN,PARMMAX,PARMFINAL,FUNCFINAL,ERROR,NITER,NHIST,HISTITER,HISTFUNC,HISTPARM,PRINTFLAG)
	USE IFPORT
	INTEGER NPARM,IPARM						!Number of parameters and parameter index
	INTEGER NPOINT,IPOINT					!Number of points in complex and point index
	REAL*8 PARMINIT(6)						!Initial parameter guesses
	REAL*8 PARMFINAL(6)						!Final best parameter estimates
	REAL*8 FUNCFINAL						!Final optimization function value
	REAL*8 PARMDELTA(6)						!Deviations from PARMINIT in which initial complex is generated
	REAL*8 PARMMIN(6),PARMMAX(6)				!Limits for parameter values
	INTEGER HIPOINT(6),LOPOINT(6)			!Index of complex point having highest and lowest parameter value
	INTEGER BEST,WORST,MOVED				!Index of complex point having best and worst function value
	REAL*8 PARMSEARCH(6,200)				!Parameter values for complex
	REAL*8 FUNC(200)						!Optimization function values for complex
	REAL*8 CENTROID(6)						!Centroid of complex
	REAL*8 ALPHA/1.3/						!Multiplier for moving worst point 
	INTEGER MAXIT/2000/,ITER				!Maximum iterations and iteration number
	LOGICAL CONVERGE						!Convergence flag
	INTEGER ERROR(8)						!Error codes
	INTEGER NITER							!# of iterations used
	INTEGER NHIST							!Info saved every 10th iteration - NHIST is number of such saves
	INTEGER HISTITER(500)					!Iteration number (ITER) at each save 
	REAL*8 HISTFUNC(500,2)					!Best/worst function values at each save
	REAL*8 HISTPARM(500,6,2)				!Low/high parameter values at each save
	REAL*8 FUNCCALC						!Name for optimization function
	LOGICAL PRINTFLAG						!Flag for printing iteration information
	INTEGER ISEED,HR,MN,SEC,FSEC				!Seed for random number generator, and time variables used to pick seed
	REAL*8 XRAN							!Random number from uniform distribution for generating initial comples
    EXTERNAL FUNCCALC
!***
!***Initialize error and history variables 
!***
	DO IPARM=1,8
		ERROR(IPARM)=0
	ENDDO
	DO NHIST=1,500
		HISTITER(ITER)=0
		HISTFUNC(ITER,1)=0
		HISTFUNC(ITER,2)=0
		DO IPARM=1,6
			HISTPARM(ITER,IPARM,1)=0
			HISTPARM(ITER,IPARM,2)=0
		ENDDO
	ENDDO
	NHIST=0
!***
!***Set parameter values for first complex point to initial guesses 
!***Set parameter values for other initial complex points to random values within specified deviations
!***Compute function values for initial complex 
!***Record parameter and function values for initial complex
!***Determine index values for complex points with best and worst function values
!***Determine index values for complex points with lowest and highest parameter values
!***
	CALL GETTIM(HR,MN,SEC,FSEC)
	ISEED=SEC*100+FSEC
	XRAN=RAND(ISEED)
	DO IPARM=1,NPARM
		PARMSEARCH(IPARM,1)=PARMINIT(IPARM)
		DO IPOINT=2,NPOINT
			XRAN=RAND(0)
			PARMSEARCH(IPARM,IPOINT)=PARMINIT(IPARM)-2*(XRAN-0.5)*PARMDELTA(IPARM)*(PARMMAX(IPARM)-PARMMIN(IPARM))
		ENDDO
		IF(PRINTFLAG) TYPE "(I5,5F10.5)",IPARM,PARMINIT(IPARM),PARMDELTA(IPARM),PARMMIN(IPARM),PARMMAX(IPARM)
	ENDDO
	DO IPOINT=1,NPOINT
		FUNC(IPOINT)=FUNCCALC(PARMSEARCH(1,IPOINT))
		IF(PRINTFLAG) TYPE "(I8,F12.6,6F10.6)",IPOINT,FUNC(IPOINT),(PARMSEARCH(IPARM,IPOINT),IPARM=1,NPARM)
	ENDDO
	WORST=1
	BEST=1
	DO IPOINT=2,NPOINT
		IF(FUNC(IPOINT).GT.FUNC(WORST)) WORST=IPOINT
		IF(FUNC(IPOINT).LT.FUNC(BEST)) BEST=IPOINT
	ENDDO
	DO IPARM=1,NPARM
		LOPOINT(IPARM)=1
		HIPOINT(IPARM)=1
		DO IPOINT=2,NPOINT
			IF(PARMSEARCH(IPARM,IPOINT).LT.PARMSEARCH(IPARM,LOPOINT(IPARM))) LOPOINT(IPARM)=IPOINT
			IF(PARMSEARCH(IPARM,IPOINT).GT.PARMSEARCH(IPARM,HIPOINT(IPARM))) HIPOINT(IPARM)=IPOINT
		ENDDO
	ENDDO
!***
!***Start of iteration loop
!***At start of each iteration, set index of MOVED point to the current WORST point
!***
	DO ITER=1,MAXIT
		IF(PRINTFLAG) TYPE "(2I5,F15.8,7F12.8)",ITER,BEST,FUNC(BEST),(PARMSEARCH(IPARM,BEST),IPARM=1,NPARM)
!		IF(PRINTFLAG) ACCEPT *,IGO
		MOVED=WORST
!***
!***Determine parameter value centroid for complex, excluding worst point
!***
        	DO IPARM=1,NPARM
        		CENTROID(IPARM)=0.0
        		DO IPOINT=1,NPOINT
        			IF(IPOINT.EQ.WORST) CYCLE
        			CENTROID(IPARM)=CENTROID(IPARM)+PARMSEARCH(IPARM,IPOINT)
        		ENDDO
        		CENTROID(IPARM)=CENTROID(IPARM)/(NPOINT-1)
        	ENDDO
!***
!***Move worst point through centroid to opposite side of complex, 
!***increasing distance from centroid by factor ALPHA, but satisfying min/max limits
!***
        	DO IPARM=1,NPARM
        		PARMSEARCH(IPARM,MOVED)=CENTROID(IPARM)+ALPHA*(CENTROID(IPARM)-PARMSEARCH(IPARM,WORST))
        		PARMSEARCH(IPARM,MOVED)=MIN(PARMSEARCH(IPARM,MOVED),PARMMAX(IPARM)-1D-6)
        		PARMSEARCH(IPARM,MOVED)=MAX(PARMSEARCH(IPARM,MOVED),PARMMIN(IPARM)+1D-6)
        	ENDDO
!***
!***Evaluate function value of moved point and determine which point is now worst
!***
        	FUNC(MOVED)=FUNCCALC(PARMSEARCH(1,MOVED))
        	DO IPOINT=1,NPOINT
        		IF(FUNC(IPOINT).GT.FUNC(WORST)) WORST=IPOINT
        	ENDDO
!***
!***If moved point is still the worst point, reduce distance from centroid by half,
!***reevaluate function, and again determine which point is now worst
!***
        	IF(WORST.EQ.MOVED) THEN
        		DO IPARM=1,NPARM
        			PARMSEARCH(IPARM,MOVED)=0.5*CENTROID(IPARM)+0.5*PARMSEARCH(IPARM,MOVED)
        		ENDDO
        		FUNC(MOVED)=FUNCCALC(PARMSEARCH(1,MOVED))
        		DO IPOINT=1,NPOINT
        			IF(FUNC(IPOINT).GT.FUNC(WORST)) WORST=IPOINT
        		ENDDO
        	ENDIF
!***
!***If moved point is still the worst point, contract all points to best point
!***
        	IF(WORST.EQ.MOVED) THEN
        		DO IPOINT=1,NPOINT
        			IF(IPOINT.EQ.BEST) CYCLE
        			DO IPARM=1,NPARM
        				PARMSEARCH(IPARM,IPOINT)=0.5*PARMSEARCH(IPARM,IPOINT)+0.5*PARMSEARCH(IPARM,BEST)
        			ENDDO
				FUNC(IPOINT)=FUNCCALC(PARMSEARCH(1,IPOINT))
        		ENDDO
        	ENDIF
!***
!***Reevaluate points 
!***
        	WORST=1
        	BEST=1
        	DO IPOINT=2,NPOINT
        		IF(FUNC(IPOINT).GT.FUNC(WORST)) WORST=IPOINT
        		IF(FUNC(IPOINT).LT.FUNC(BEST)) BEST=IPOINT
        	ENDDO
        	DO IPARM=1,NPARM
        		LOPOINT(IPARM)=1
        		HIPOINT(IPARM)=1
        		DO IPOINT=2,NPOINT
        			IF(PARMSEARCH(IPARM,IPOINT).LT.PARMSEARCH(IPARM,LOPOINT(IPARM))) LOPOINT(IPARM)=IPOINT
        			IF(PARMSEARCH(IPARM,IPOINT).GT.PARMSEARCH(IPARM,HIPOINT(IPARM))) HIPOINT(IPARM)=IPOINT
        		ENDDO
        	ENDDO
!***
!***Test whether convergence has occurred
!***
        	CONVERGE=.TRUE.
        	IF(ABS((FUNC(BEST)-FUNC(WORST))/FUNC(BEST)).GT.1D-10) THEN
        		CONVERGE=.FALSE.
        	ELSE
        		DO IPARM=1,NPARM
        			IF(ABS(PARMSEARCH(IPARM,BEST)).LT.1D-3) THEN
        				IF(ABS(PARMSEARCH(IPARM,HIPOINT(IPARM))-PARMSEARCH(IPARM,LOPOINT(IPARM))).GT.1D-9) CONVERGE=.FALSE.
        			ELSE
        				IF(ABS((PARMSEARCH(IPARM,HIPOINT(IPARM))-PARMSEARCH(IPARM,LOPOINT(IPARM)))/PARMSEARCH(IPARM,BEST)).GT.1D-6) CONVERGE=.FALSE.
        			ENDIF
        		ENDDO
        	ENDIF
        	IF(CONVERGE) EXIT
!***
!***End of iteration loop - record history every 10 iterations
!***
        	IF(MOD(ITER,10).EQ.1) THEN
			NHIST=NHIST+1
        		HISTITER(NHIST)=ITER
        		HISTFUNC(NHIST,1)=FUNC(BEST)
        		HISTFUNC(NHIST,2)=FUNC(WORST)
        		DO IPARM=1,NPARM
        			HISTPARM(NHIST,IPARM,1)=PARMSEARCH(IPARM,LOPOINT(IPARM))
        			HISTPARM(NHIST,IPARM,2)=PARMSEARCH(IPARM,HIPOINT(IPARM))		
        		ENDDO
        	ENDIF
	ENDDO
!***
!***Set final parameter estimates to parameter values of best point and record final history info
!***
	DO IPARM=1,NPARM
		PARMFINAL(IPARM)=PARMSEARCH(IPARM,BEST)
	ENDDO
	FUNCFINAL=FUNC(BEST)
	NHIST=NHIST+1
	HISTITER(NHIST)=ITER
	HISTFUNC(NHIST,1)=FUNC(BEST)
	HISTFUNC(NHIST,2)=FUNC(WORST)
	DO IPARM=1,NPARM
		HISTPARM(NHIST,IPARM,1)=PARMSEARCH(IPARM,LOPOINT(IPARM))
		HISTPARM(NHIST,IPARM,2)=PARMSEARCH(IPARM,HIPOINT(IPARM))		
	ENDDO
	NITER=ITER
!***
!***Set error flags and return
!***
	IF(CONVERGE) THEN
		ERROR(1)=0
	ELSE
		ERROR(1)=1
	ENDIF
	DO IPARM=1,NPARM
		ERROR(IPARM+1)=0
		IF(ABS((PARMFINAL(IPARM)-PARMMAX(IPARM))/PARMMAX(IPARM)).LT.1D-5) ERROR(IPARM+1)=1
		IF(ABS((PARMFINAL(IPARM)-PARMMIN(IPARM))/PARMMIN(IPARM)).LT.1D-5) ERROR(IPARM+1)=-1
	ENDDO
	RETURN

	END SUBROUTINE COMPLEXOPTIMIZE