C FAMOEB -- amoeba
C This is taken from the AMOEBA subroutine from Numerical Recipes by
C Press, Flannery, Teukolsky and Vetterling.  The calling sequence of
C this routine and of FUNK have been modified, and the PAUSE statement
C was replaced by RETURN.
C This routine minimizes the function FUNK.

	SUBROUTINE FAMOEB (FNL, P, Y, NPAR, FTOL, FUNK, ITER)

	PARAMETER (NMAX = 20, ITMAX = 500)
	PARAMETER (ALPHA = 1.0, BETA = 0.5, GAMMA = 2.0)
	INTEGER FNL
	DIMENSION P(NPAR+1,NPAR), Y(NPAR+1)
	DIMENSION PR(NMAX), PRR(NMAX), PBAR(NMAX)
        EXTERNAL FUNK

	MPTS = NPAR + 1
	ITER = 0
1	ILO = 1
	IF (Y(1) .GT. Y(2)) THEN
	    IHI = 1
	    INHI = 2
	ELSE
	    IHI = 2
	    INHI = 1
	ENDIF
	DO 10 I = 1, MPTS
	    IF (Y(I) .LT. Y(ILO)) ILO = I
	    IF (Y(I) .GT. Y(IHI)) THEN
		INHI = IHI
		IHI = I
	    ELSE IF (Y(I) .GT. Y(INHI)) THEN
		IF (I .NE. IHI) INHI = I
	    ENDIF
 10     CONTINUE
	RTOL = 2. * ABS(Y(IHI) - Y(ILO)) / (ABS(Y(IHI)) + ABS(Y(ILO)))
C
C	Here is where we return, either because the parameters are close
C	enough to the same set of values at each vertex or because we
C	have reached the maximum number of iterations.
C
	IF (RTOL .LT. FTOL) RETURN
	IF (ITER .EQ. ITMAX) RETURN
	ITER = ITER + 1
	DO 20 J = 1, NPAR
	    PBAR(J) = 0.
 20     CONTINUE
	DO 40 I = 1, MPTS
	    IF (I .NE. IHI) THEN
		DO 30 J = 1, NPAR
		    PBAR(J) = PBAR(J) + P(I,J)
 30             CONTINUE
	    ENDIF
 40     CONTINUE
	DO 50 J = 1, NPAR
	    PBAR(J) = PBAR(J) / NPAR
	    PR(J) = (1. + ALPHA) * PBAR(J) - ALPHA * P(IHI,J)
 50     CONTINUE
	YPR = FUNK (FNL, PR)
	IF (YPR .LE. Y(ILO)) THEN
	    DO 60 J = 1, NPAR
		PRR(J) = GAMMA * PR(J) + (1. - GAMMA) * PBAR(J)
 60         CONTINUE
	    YPRR = FUNK (FNL, PRR)
	    IF (YPRR .LT. Y(ILO)) THEN
		DO 70 J = 1, NPAR
		    P(IHI,J) = PRR(J)
 70             CONTINUE
		Y(IHI) = YPRR
	    ELSE
		DO 80 J = 1, NPAR
		    P(IHI,J) = PR(J)
 80             CONTINUE
		Y(IHI) = YPR
	    ENDIF
	ELSE IF (YPR .GE. Y(INHI)) THEN
	    IF (YPR .LT. Y(IHI)) THEN
		DO 90 J = 1, NPAR
		    P(IHI,J) = PR(J)
 90             CONTINUE
		Y(IHI) = YPR
	    ENDIF
	    DO 100 J = 1, NPAR
		PRR(J) = BETA * P(IHI,J) + (1. - BETA) * PBAR(J)
 100        CONTINUE
	    YPRR = FUNK (FNL, PRR)
	    IF (YPRR .LT. Y(IHI)) THEN
		DO 110 J = 1, NPAR
		    P(IHI,J) = PRR(J)
 110            CONTINUE
		Y(IHI) = YPRR
	    ELSE
		DO 130 I = 1, MPTS
		    IF (I .NE. ILO) THEN
			DO 120 J = 1, NPAR
			    PR(J) = 0.5 * (P(I,J) + P(ILO,J))
			    P(I,J) = PR(J)
 120                    CONTINUE
			Y(I) = FUNK (FNL, PR)
		    ENDIF
 130            CONTINUE
	    ENDIF
	ELSE
	    DO 140 J = 1, NPAR
		P(IHI,J) = PR(J)
 140        CONTINUE
	    Y(IHI) = YPR
	ENDIF
	GO TO 1
	END