;+
; NAME:
;	dbrent
; PURPOSE:
;	minimum finding with derivatives using Brent's method
; DESCRIPTION:
;  Given a function f and its derivative function df, and given a
;  bracketing triplet of abscissas ax, bx, cx [such that bx is between ax
;  and cx, and f(bx) is less than both f(ax) and f(cx)], this routine
;  isolates the minimum to a fractional precision of about tol using a modi
;  cation of Brent's method that uses derivatives. The abscissa of the
;  minimum is returned as xmin, and the minimum function value is returned
;  as dbrent, the returned function value.
; CALLING SEQUENCE:
;	min = dbrent(func, ax,bx,cx, xmin, tol=tol, extras=extras)
;
; INPUT PARAMETERS
;   func: a string containing the name of the function to be minimized. 
;     func, x, f, df,   OR  func,x,f,df, p
;   ax, bx, cx: a bracketing triplet of abscissas ax, bx, cx 
;    [such that bx is between ax and cx, and f(bx) is less than both f(ax) and f(cx)]
;
; OPTIONAL INPUT PARAMETERS
;   tol - tolerance
;   funcp - optional other parameters to be passed to f
;     
;
; OUTPUT PARAMETERS
;	xmin abscissa of the minimum
;   dbrent - minimum function value, returned
;
; PROCEDURE:
;   Based on dbrent in Numerical Recipes
;
; REVISION HISTORY
;-


function dbrent, func, ax,bx,cx, xmin, tol=tol, p=p

ITMAX = 100
ZEPS = 1.0e-10
pq = keyword_set(funcp)

;#define MOV3(a,b,c, d,e,f) (a)=(d);(b)=(e);(c)=(f);
;#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))

a = (ax < cx)                  ; a=(ax < cx ? ax : cx);
b = (ax > cx)                  ; b=(ax > cx ? ax : cx);
x = bx & w = bx & v = bx       ; x=w=v=bx
if pq then call_function(func,x,f,df,p) else call_function(func,x,f,df)
fw=f & fw=f & fx=f             ; fw=fv=fx=(*f)(x);
dw=df & dv=df & dx=df          ; dw=dv=dx=(*df)(x)
for iter = 1, itmax do begin   ; for (iter=1;iter<=ITMAX;iter++) {
  xm=0.5*(a+b)                 ; xm=0.5*(a+b);
  tol1=tol*abs(x)+ZEPS         ; tol1=tol*fabs(x)+ZEPS;
  tol2=2.0*tol1                ; tol2=2.0*tol1;
if (abs(x-xm) le (tol2-0.5*(b-a))) do begin ; if (fabs(x-xm) <= (tol2-0.5*(b-a))) {
   xmin=x                      ; *xmin=x
   return, fx                  ; return fx
endif                          ; }
if (abs(e) gt tol1) do begin         ; if (fabs(e) > tol1) {
  d1=2.0*(b-a)                 ;d1=2.0*(b-a); Initialize these d's to an out-of-bracket
  d2=d1                        ;d2=d1       ; value.
  if (dw != dx) then d1=(w-x)*dx/(dx-dw);;if (dw != dx) d1=(w-x)*dx/(dx-dw); Secant method with one point.
  if (dv != dx) d2=(v-x)*dx/(dx-dv)      ;if (dv != dx) d2=(v-x)*dx/(dx-dv); And the other.
                                         ; Which of these two estimates of d shall we take? We will insist that they be within
                                         ; the bracket, and on the side pointed to by the derivative at x:
  u1=x+d1               ; u1=x+d1
  u2=x+d2               ; u2=x+d2
  ok1 = (a-u1)*(u1-b) gt 0.0 and dx*d1 le 0.0; ok1 = (a-u1)*(u1-b) > 0.0 && dx*d1 <= 0.0
  ok2 = (a-u2)*(u2-b) gt 0.0 abd dx*d2 le 0.0; ok2 = (a-u2)*(u2-b) > 0.0 && dx*d2 <= 0.0
  olde=e;                 olde=e; Movement on the step before last
  e=d;    e=d;
  if (ok1 || ok2) do begin     ; if (ok1 || ok2) { Take only an acceptable d, and if
                               ; both are acceptable, then take
                               ; the smallest one.
     if (ok1 && ok2) do begin  ; if (ok1 && ok2)
        if (abs(d1) lt abs(d2) then d = d1 else d = d2  ;d=(fabs(d1) < fabs(d2) ? d1 : d2);
     endif else begin
        if (ok1) do begin   
           d=d1                       
        endif else begin
           d=d2;
        endelse
     endelse
     if (abs(d) le abs(0.5*olde)) do begin
         u=x+d;
        if (u-a < tol2 || b-u < tol2) do begin
           d=SIGN(tol1,xm-x);
        endif
     endif else begin
         d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
     endelse
  endif else begin
      d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
  endelse
endif else begin
    d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
endelse
if (fabs(d) >= tol1) {
u=x+d;
fu=(*f)(u);
} else {
u=x+SIGN(tol1,d);
fu=(*f)(u);
if (fu > fx) { If the minimum step in the downhill
direction takes us uphill, then
we are done.
*xmin=x;
return fx;
}
}
du=(*df)(u); Now all the housekeeping,
if (fu <= fx) {
if (u >= x) a=x; else b=x;
MOV3(v,fv,dv, w,fw,dw)
MOV3(w,fw,dw, x,fx,dx)
MOV3(x,fx,dx, u,fu,du)
} else {
if (u < x) a=u; else b=u;
if (fu <= fw || w == x) {
MOV3(v,fv,dv, w,fw,dw)
MOV3(w,fw,dw, u,fu,du)
} else if (fu < fv || v == x || v == w) {
MOV3(v,fv,dv, u,fu,du)
}
}
}
nrerror("Too many iterations in routine dbrent");
return 0.0; Never get here.
}