pro minF_bracket, xa,xb,xc, fa,fb,fc, FUNC_NAME=func_name, \$ POINT_NDIM=pn, DIRECTION=dirn ;+ ; NAME: ; MINF_BRACKET ; PURPOSE: ; Bracket a local minimum of a 1-D function with 3 points, ; EXPLANATION: ; Brackets a local minimum of a 1-d function with 3 points, ; thus ensuring that a minimum exists somewhere in the interval. ; This routine assumes that the function has a minimum somewhere.... ; Routine can also be applied to a scalar function of many variables, ; for such case the local minimum in a specified direction is bracketed, ; This routine is called by minF_conj_grad, to bracket minimum in the ; direction of the conjugate gradient of function of many variables ; CALLING EXAMPLE: ; xa=0 & xb=1 ; minF_bracket, xa,xb,xc, fa,fb,fc, FUNC_NAME="name" ;for 1-D func. ; or: ; minF_bracket, xa,xb,xc, fa,fb,fc, FUNC="name", \$ ; POINT=[0,1,1], \$ ; DIRECTION=[2,1,1] ;for 3-D func. ; INPUTS: ; xa = scalar, guess for point bracketing location of minimum. ; xb = scalar, second guess for point bracketing location of minimum. ; KEYWORDS: ; FUNC_NAME = function name (string) ; Calling mechanism should be: F = func_name( px ) ; where: ; px = scalar or vector of independent variables, input. ; F = scalar value of function at px. ; POINT_NDIM = when working with function of N variables, ; use this keyword to specify the starting point in N-dim space. ; Default = 0, which assumes function is 1-D. ; DIRECTION = when working with function of N variables, ; use this keyword to specify the direction in N-dim space ; along which to bracket the local minimum, (default=1 for 1-D). ; (xa,xb,xc) are then relative distances from POINT_NDIM. ; OUTPUTS: ; xa,xb,xc = scalars, 3 points which bracket location of minimum, ; that is, f(xb) < f(xa) and f(xb) < f(xc), so minimum exists. ; When working with function of N variables ; (xa,xb,xc) are then relative distances from POINT_NDIM, ; in the direction specified by keyword DIRECTION, ; with scale factor given by magnitude of DIRECTION. ; OPTIONAL OUTPUT: ; fa,fb,fc = value of function at 3 points which bracket the minimum, ; again note that fb < fa and fb < fc if minimum exists. ; PROCEDURE: ; algorithm from Numerical Recipes (by Press, et al.), sec.10.1 (p.281). ; MODIFICATION HISTORY: ; Written, Frank Varosi NASA/GSFC 1992. ; Converted to IDL V5.0 W. Landsman September 1997 ;- goldm = (sqrt(5)+1)/2 ;golden mean factor to march with. glimit = 100 ;maximum factor to try. tiny = 1.e-19 ;a tiny number to avoid divide by zero. if N_elements( pn ) LE 0 then begin pn = 0 dirn = 1 endif if (xa EQ xb) then xb = xa + 1 fa = call_function( func_name, pn + xa * dirn ) fb = call_function( func_name, pn + xb * dirn ) if (fb GT fa) then begin x = xa & xa = xb & xb = x f = fa & fa = fb & fb = f endif xc = xb + goldm * (xb-xa) fc = call_function( func_name, pn + xc * dirn ) while (fb GE fc) do begin zba = xb-xa zbc = xb-xc r = zba * (fb-fc) q = zbc * (fb-fa) delta = q-r sign = 1 - 2 * (delta LT 0) xu = xb - (zbc * q - zba * r)/(2* sign * (abs( delta ) > tiny) ) ulim = xb + glimit * (xc-xb) if ( (xb-xu)*(xu-xc) GT 0 ) then begin fu = call_function( func_name, pn + xu * dirn ) if (fu LT fc) then begin xa = xb & xb = xu fa = fb & fb = fu return endif else if (fu GT fb) then begin xc = xu fc = fu return endif xu = xc - goldm * zbc fu = call_function( func_name, pn + xu * dirn ) endif else if ( (xc-xu)*(xu-ulim) GT 0 ) then begin fu = call_function( func_name, pn + xu * dirn ) if (fu LT fc) then begin xb = xc & fb = fc xc = xu & fc = fu xu = xc + goldm * (xc-xb) fu = call_function( func_name, pn + xu * dirn ) endif endif else if ( (ulim-xc)*(xu-ulim) GE 0 ) then begin xu = ulim fu = call_function( func_name, pn + xu * dirn ) endif else begin xu = xc + goldm * (xc-xb) fu = call_function( func_name, pn + xu * dirn ) endelse xa = xb & xb = xc & xc = xu fa = fb & fb = fc & fc = fu endwhile return end