kenkehoe · October 1, 2014 20:43
diff --git a/zbrent.pro b/zbrent.pro
 function ZBRENT, x1, x2, FUNC_NAME=func_name,    $
                         MAX_ITERATIONS=maxit, TOLERANCE=TOL
 ;+
 ; NAME:
 ;     ZBRENT
 ; PURPOSE:
 ;     Find the zero of a 1-D function up to specified tolerance with IDL.
 ; EXPLANTION:
 ;     This routine assumes that the function is known to have a zero.
 ;     Adapted from procedure of the same name in "Numerical Recipes" by
 ;     Press et al. (1992), Section 9.3
 ;
 ; CALLING:
 ;       x_zero = ZBRENT( x1, x2, FUNC_NAME="name", MaX_Iter=, Tolerance= )
 ;
 ; INPUTS:
 ;       x1, x2 = scalars, 2 points which bracket location of function zero,
 ;                                               that is, F(x1) < 0 < F(x2).
 ;       Note: computations are performed with
 ;       same precision (single/double) as the inputs and user supplied function.
 ;
 ; REQUIRED INPUT KEYWORD:
 ;       FUNC_NAME = function name (string)
 ;               Calling mechanism should be:  F = func_name( px )
 ;               where:  px = scalar independent variable, input.
 ;                       F = scalar value of function at px,
 ;                           should be same precision (single/double) as input.
 ;
 ; OPTIONAL INPUT KEYWORDS:
 ;       MAX_ITER = maximum allowed number iterations, default=100.
 ;       TOLERANCE = desired accuracy of minimum location, default = 1.e-3.
 ;
 ; OUTPUTS:
 ;       Returns the location of zero, with accuracy of specified tolerance.
 ;
 ; PROCEDURE:
 ;       Brent's method to find zero of a function by using bracketing,
 ;       bisection, and inverse quadratic interpolation,
 ;
 ; EXAMPLE:
 ;       Find the root of the COSINE function between 1. and 2.  radians
 ;
 ;        IDL> print, zbrent( 1, 2, FUNC = 'COS')
 ;
 ;       and the result will be !PI/2 within the specified tolerance
 ; MODIFICATION HISTORY:
 ;       Written, Frank Varosi NASA/GSFC 1992.
 ;       FV.1994, mod to check for single/double prec. and set zeps accordingly.
 ;       Converted to IDL V5.0   W. Landsman   September 1997
 ;       Use MACHAR() to define machine precision   W. Landsman September 2002
 ; kenkehoe
 ;-
        if N_params() LT 2 then begin
             print,'Syntax - result = ZBRENT( x1, x2, FUNC_NAME = ,'
             print,'                  [ MAX_ITER = , TOLERANCE = ])'
             return, -1
        endif

        if N_elements( TOL ) NE 1 then TOL = 1.e-3
        if N_elements( maxit ) NE 1 then maxit = 100

        if size(x1,/TNAME) EQ 'DOUBLE' OR size(x2,/TNAME) EQ 'DOUBLE' then begin
                xa = double( x1 )
                xb = double( x2 )
                zeps = (machar(/DOUBLE)).eps   ;machine epsilon in double.
          endif else begin
                xa = x1
                xb = x2
                zeps = (machar(/DOUBLE)).eps   ;machine epsilon, in single 
           endelse

        fa = call_function( func_name, xa )
        fb = call_function( func_name, xb )
        fc = fb

        if (fb*fa GT 0) then begin
                message,"root must be bracketed by the 2 inputs",/INFO
                return,xa
           endif

        for iter = 1,maxit do begin

                if (fb*fc GT 0) then begin
                        xc = xa
                        fc = fa
                        Din = xb - xa
                        Dold = Din
                   endif

                if (abs( fc ) LT abs( fb )) then begin
                        xa = xb   &   xb = xc   &   xc = xa
                        fa = fb   &   fb = fc   &   fc = fa
                   endif

                TOL1 = 0.5*TOL + 2*abs( xb ) * zeps     ;Convergence check
                xm = (xc - xb)/2.

                if (abs( xm ) LE TOL1) OR (fb EQ 0) then return,xb

                if (abs( Dold ) GE TOL1) AND (abs( fa ) GT abs( fb )) then begin

                        S = fb/fa       ;attempt inverse quadratic interpolation

                        if (xa EQ xc) then begin
                                p = 2 * xm * S
                                q = 1-S
                          endif else begin
                                T = fa/fc
                                R = fb/fc
                                p = S * (2*xm*T*(T-R) - (xb-xa)*(R-1) )
                                q = (T-1)*(R-1)*(S-1)
                           endelse

                        if (p GT 0) then q = -q
                        p = abs( p )
                        test = ( 3*xm*q - abs( q*TOL1 ) ) < abs( Dold*q )

                        if (2*p LT test)  then begin
                                Dold = Din              ;accept interpolation
                                Din = p/q
                          endif else begin
                                Din = xm                ;use bisection instead
                                Dold = xm
                           endelse

                  endif else begin

                        Din = xm    ;Bounds decreasing to slowly, use bisection
                        Dold = xm
                   endelse

                xa = xb
                fa = fb         ;evaluate new trial root.

                if (abs( Din ) GT TOL1) then xb = xb + Din $
                                        else xb = xb + TOL1 * (1-2*(xm LT 0))

                fb = call_function( func_name, xb )
          endfor

        message,"exceeded maximum number of iterations: "+strtrim(iter,2),/INFO

 return, xb
 end
	function ZBRENT, x1, x2, FUNC_NAME=func_name, $
	MAX_ITERATIONS=maxit, TOLERANCE=TOL
	;+
	; NAME:
	; ZBRENT
	; PURPOSE:
	; Find the zero of a 1-D function up to specified tolerance with IDL.
	; EXPLANTION:
	; This routine assumes that the function is known to have a zero.
	; Adapted from procedure of the same name in "Numerical Recipes" by
	; Press et al. (1992), Section 9.3
	;
	; CALLING:
	; x_zero = ZBRENT( x1, x2, FUNC_NAME="name", MaX_Iter=, Tolerance= )
	;
	; INPUTS:
	; x1, x2 = scalars, 2 points which bracket location of function zero,
	; that is, F(x1) < 0 < F(x2).
	; Note: computations are performed with
	; same precision (single/double) as the inputs and user supplied function.
	;
	; REQUIRED INPUT KEYWORD:
	; FUNC_NAME = function name (string)
	; Calling mechanism should be: F = func_name( px )
	; where: px = scalar independent variable, input.
	; F = scalar value of function at px,
	; should be same precision (single/double) as input.
	;
	; OPTIONAL INPUT KEYWORDS:
	; MAX_ITER = maximum allowed number iterations, default=100.
	; TOLERANCE = desired accuracy of minimum location, default = 1.e-3.
	;
	; OUTPUTS:
	; Returns the location of zero, with accuracy of specified tolerance.
	;
	; PROCEDURE:
	; Brent's method to find zero of a function by using bracketing,
	; bisection, and inverse quadratic interpolation,
	;
	; EXAMPLE:
	; Find the root of the COSINE function between 1. and 2. radians
	;
	; IDL> print, zbrent( 1, 2, FUNC = 'COS')
	;
	; and the result will be !PI/2 within the specified tolerance
	; MODIFICATION HISTORY:
	; Written, Frank Varosi NASA/GSFC 1992.
	; FV.1994, mod to check for single/double prec. and set zeps accordingly.
	; Converted to IDL V5.0 W. Landsman September 1997
	; Use MACHAR() to define machine precision W. Landsman September 2002
	; kenkehoe
	;-
	if N_params() LT 2 then begin
	print,'Syntax - result = ZBRENT( x1, x2, FUNC_NAME = ,'
	print,' [ MAX_ITER = , TOLERANCE = ])'
	return, -1
	endif

	if N_elements( TOL ) NE 1 then TOL = 1.e-3
	if N_elements( maxit ) NE 1 then maxit = 100

	if size(x1,/TNAME) EQ 'DOUBLE' OR size(x2,/TNAME) EQ 'DOUBLE' then begin
	xa = double( x1 )
	xb = double( x2 )
	zeps = (machar(/DOUBLE)).eps ;machine epsilon in double.
	endif else begin
	xa = x1
	xb = x2
	zeps = (machar(/DOUBLE)).eps ;machine epsilon, in single
	endelse

	fa = call_function( func_name, xa )
	fb = call_function( func_name, xb )
	fc = fb

	if (fb*fa GT 0) then begin
	message,"root must be bracketed by the 2 inputs",/INFO
	return,xa
	endif

	for iter = 1,maxit do begin

	if (fb*fc GT 0) then begin
	xc = xa
	fc = fa
	Din = xb - xa
	Dold = Din
	endif

	if (abs( fc ) LT abs( fb )) then begin
	xa = xb & xb = xc & xc = xa
	fa = fb & fb = fc & fc = fa
	endif

	TOL1 = 0.5TOL + 2abs( xb ) * zeps ;Convergence check
	xm = (xc - xb)/2.

	if (abs( xm ) LE TOL1) OR (fb EQ 0) then return,xb

	if (abs( Dold ) GE TOL1) AND (abs( fa ) GT abs( fb )) then begin

	S = fb/fa ;attempt inverse quadratic interpolation

	if (xa EQ xc) then begin
	p = 2 * xm * S
	q = 1-S
	endif else begin
	T = fa/fc
	R = fb/fc
	p = S * (2xmT(T-R) - (xb-xa)(R-1) )
	q = (T-1)(R-1)(S-1)
	endelse

	if (p GT 0) then q = -q
	p = abs( p )
	test = ( 3xmq - abs( qTOL1 ) ) < abs( Doldq )

	if (2*p LT test) then begin
	Dold = Din ;accept interpolation
	Din = p/q
	endif else begin
	Din = xm ;use bisection instead
	Dold = xm
	endelse

	endif else begin

	Din = xm ;Bounds decreasing to slowly, use bisection
	Dold = xm
	endelse

	xa = xb
	fa = fb ;evaluate new trial root.

	if (abs( Din ) GT TOL1) then xb = xb + Din $
	else xb = xb + TOL1 * (1-2*(xm LT 0))

	fb = call_function( func_name, xb )
	endfor

	message,"exceeded maximum number of iterations: "+strtrim(iter,2),/INFO

	return, xb
	end
No results found