PRO Amoeba_x, p, y, func_name, iter, ftol=ftol0, alpha=alpha0, beta=beta0, $
              gammx=gammx0, itmax=itmax0, quiet=quiet

   iter = 0

;set keyword defaults
   IF(KEYWORD_SET(quiet)) THEN qq = 1 ELSE qq = 0
   IF(N_ELEMENTS(alpha0) GT 0) THEN alpha = alpha0 ELSE alpha = 1.0
   IF(N_ELEMENTS(beta0) GT 0) THEN beta = beta0 ELSE beta = 0.5
   IF(N_ELEMENTS(gammx0) GT 0) THEN gammx = gammx0 ELSE gammx = 2.0
   IF(N_ELEMENTS(ftol0) GT 0) THEN ftol = ftol0 ELSE ftol = 1.0e-5
   IF(N_ELEMENTS(itmax0) GT 0) THEN itmax = itmax0 ELSE itmax = 500

   typ = size(p) & typ = typ(typ(0)+1)
   
   n = N_ELEMENTS(p(0, *))
   IF(n LT 2) THEN BEGIN
      message, /info, 'THIS IS SUPPOSED TO BE MULTIDIMENSIONAL, '
      print, 'IT WON''T WORK WITHOUT AT LEAST N=2'
      RETURN
   ENDIF
   
   m = N_ELEMENTS(p(*, 0))
   IF(m NE n+1) THEN BEGIN
      message, /info, 'WRONG NUMBER OF ELEMENTS IN P'
      RETURN
   ENDIF

   ny = N_ELEMENTS(y)
   IF(ny NE m) THEN BEGIN
      message, /info, 'WRONG NUMBER OF ELEMENTS IN Y'
      RETURN
   ENDIF

;loop
   iter = 0
   loop1: ys = sort(y)
;get max and min of y, and also the second highest point, the easiest way to
;do this is sort it.
      ihi = ys(ny-1)
      inhi = ys(ny-2)
      ilo = ys(0)

;compute the fractional range from the highest to lowest, if it's ok, return
      rtol =  2.0*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo)))
      IF(rtol LT ftol) THEN RETURN
      IF(iter EQ itmax) THEN BEGIN
         IF(qq EQ 0) THEN message, 'EXCEEDING MAX. ITERATIONS', /cont
         RETURN
      ENDIF
      iter = iter+1

;New iteration, compute the vector average of all points except the highest.
;i.e., the center of the "face" of the simplex across from the high point.
;We will subsequently explore explore along the ray from the high point
;to the center.
      pbar = make_array(n, type = typ) + 0 ;to assure same type as p
      FOR j = 0, n-1 DO pbar(j) = total(p(ys(0:ny-2), j))/n
;reflect the high point through the face, extrapolate through the face
;by a factor of alpha.
      pr = (1.0+alpha)*pbar-alpha*p(ihi, *)
;Get the function at the reflected point
      ypr = call_function(func_name, pr) ; hey, ypr is a scalar

      IF(ypr LE y(ilo)) THEN BEGIN ;Result is better than low pt.
         prr = gammx*pr+(1.0-gammx)*pbar ;so try an additional gammx's worth
         yprr = call_function(func_name, prr) ;of extrapolation,
         IF(yprr LT y(ilo)) THEN BEGIN ;which works, so replace y(ihi)
            p(ihi, *) = prr
            y(ihi) = yprr
         ENDIF ELSE BEGIN       ;additional extrap failed, but point prr is ok
            p(ihi, *) = pr
            y(ihi) = ypr
         ENDELSE
      ENDIF ELSE IF(ypr GE y(inhi)) THEN BEGIN ;highest pt. is worse than 2nd
         IF(ypr LT y(ihi)) THEN BEGIN ;but better than the highest, replace it
            p(ihi, *) = pr
            y(ihi) = ypr
         ENDIF                  ;but look for an intermediate pt, i.e., do a
         prr = beta*p(ihi, *)+(1.0-beta)*pbar ;contraction of the simplex in 1d
         yprr = call_function(func_name, prr)
         IF(yprr LT y(ihi)) THEN BEGIN ;contraction is an improvement
            p(ihi, *) = prr
            y(ihi) = yprr
         ENDIF ELSE BEGIN       ;Cant' seem to get rid of the hi pt,
            FOR i = 1, ny-1 DO BEGIN ;So we contract around the low pt
               pr = 0.5*(p(ys(i), *)+p(ys(0), *))
               p(ys(i), *) = pr
               y(ys(i)) = call_function(func_name, pr)
            ENDFOR
         ENDELSE
      ENDIF ELSE BEGIN          ;We're here if we got a middling point,
         p(ihi, *) = pr         ;Replace the high point and continue.
         y(ihi) = ypr
      ENDELSE
      
   GOTO, loop1

   RETURN
END