function nrZBrent, Func, X1, X2, Tol, arg=Arg, nar=NaR

@compile_opt.pro		; On error, return to caller

;+
; NAME:
;	nrZBrent
; PURPOSE:
;	Returns roots of function Func between X1 and X2 (with accuracy Tol)
; CATEGORY:
;	MATH: Numerical Recipes
; CALLING SEQUENCE:
;	Z0 = nrZBrent(Func,X1,X2,Tol,nar=NaR,arg=Arg)
; INPUTS:
;	Func		string		name of function for which roots are needed
;	X1,X2		arrays of same size
;					containing two coordinate values which
;					bracket the root
;	Tol		scalar		tolerance for the calculated root
; OPTIONAL INPUTS:
;	nar=NaR		scalar		Not a Root: number used to indicate a
;					root that was not bracketed by the
;					input X1,X2 values
;					default: !Values.F_NAN
;	arg=Arg		anything	is passed unmodified to all Func calls
;					For instance, Arg can be an array of
;					the same size as X1,X2, with separate
;					argument values for each of the
;					intervals [X1,X2]
; OUTPUTS:
;	nrZBrent 	array of same size as X1, X2
;					the calculated roots
; RESTRICTIONS:
;	The root must be bracketed by X1 and X2, i.e. Func(X1)*Func(X2) <= 0.
;	The function for which the root is calculated is called in the form
;	F = call_function(Func, X, Arg), where X is an array of the same
;	size as X1 and X2.
; PROCEDURE:
;	See Numerical Recipes, Press and Teukolsky
;	The main difference with the Numerical Recipes version is that
;	exact roots are intercepted right at the start of the iteration loop
;	(in Numerical Recipes they are intercepted at the convergence check).
; MODIFICATION HISTORY:
;	JAN-1998, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

ITMax	= 100		; Maximum number of iterations
EPS 	= 3.E-8		; Floating point accuracy
if n_elements(NaR) eq 0 then NaR = !Values.F_NAN

Zero = 0.
SyncArgs, X1, X2, Zero; X1, X2, Zero have same dimensions
Zero[*] = NaR		; Initialize to 'Not a Root' value

A  = X1
B  = X2
fA = call_function(Func, A, Arg)
fB = call_function(Func, B, Arg)

Tmp = fB*fA LE 0	; Are roots bracketed ?

IZ = where( Tmp )	; These roots are bracketed
IF IZ[0] EQ -1 THEN BEGIN
	message, /info, 'none of the roots are bracketed
	return, Zero
ENDIF ELSE IF (where( 1B-Tmp ))[0] NE -1 THEN		$
	message, /info, 'not all roots bracketed

A  = A [IZ]		; Pick up bracketed roots
B  = B [IZ]		; ... (unbracketed roots are ditched and remain at NaR)
fA = fA[IZ]
fB = fB[IZ]

C  = B			; Initialize third trial point to B
fC = fB

REPEAT BEGIN
	Tmp = fB EQ 0			; First check for exact roots
	IP  = where( Tmp )
	IF IP[0] NE -1 THEN BEGIN	; Exact roots found

		Zero[IZ[IP]] = B[IP]	; Set exact roots
		IP = where( 1B-Tmp )	; Remaining non-exact roots
		IF IP[0] EQ -1 THEN return, Zero ; All roots found

		A  = A [IP]		; Pick up remaining non-exact roots,
		B  = B [IP]		; ... ditching just-found exact roots
		C  = C [IP]
		fA = fA[IP]
		fB = fB[IP]
		fC = fC[IP]

		IZ = IZ[IP]		; Update Zero-indices remaining roots
	ENDIF

; On the first pass fC=fB. Since exact roots (fB=0) have already been
; intercepted, all remaining non-exact roots will pass the next test.
; So D and E are initialized here during the first pass.

	IP = where( fB*fC GT 0 )	; A and C,B on different sides of root
	IF IP[0] NE -1 THEN BEGIN
		C [IP] = A [IP]
		fC[IP] = fA[IP]
		D = B[IP]-A[IP]
		E = D
	ENDIF				; Now A,C and B on different sides of root

	IP = where( abs(fC) LT abs(fB) ); C closer to zero than B
	IF IP[0] NE -1 THEN BEGIN 	; Interchange B and C
		A [IP] = B [IP]		; A is set to new C, so A and C are
		B [IP] = C [IP]		; ... still on same side of root
		C [IP] = A [IP]
		fA[IP] = fB[IP]
		fB[IP] = fC[IP]
		fC[IP] = fA[IP]
	ENDIF

	Tol1 = 2.*EPS*abs(B)+.5*Tol	; Convergence check
	XM = .5*(C-B)

	Tmp = abs(XM) LE Tol1

	IP = where( Tmp )		; Roots located within tolerance
	IF IP[0] NE -1 THEN Zero[IZ[IP]] = B[IP]

	IP = where( 1B-Tmp )		; All roots found
	IF IP[0] EQ -1 THEN return, Zero

	A  = A [IP]			; Pick up remaining roots
	B  = B [IP]			; .. ditching the just-found roots
	C  = C [IP]
	fA = fA[IP]
	fB = fB[IP]
	fC = fC[IP]

	IZ = IZ[IP]			; Indices into Zero array

	D  = D [IP]
	E  = E [IP]
	XM = XM[IP]

	Tmp = abs(E) GE Tol1 AND abs(fA) GT abs(fB)

	IP = where( 1B-Tmp )		; Bounds decreasing too slowly, use bisection
	IF IP[0] NE -1 THEN BEGIN
		E[IP] = XM[IP]
		D[IP] = XM[IP]
	ENDIF

	IP = where( Tmp )
	IF IP[0] NE -1 THEN BEGIN
		R = fB[IP]/fC[IP]
		S = fB[IP]/fA[IP]	; Attempt inverse quadratic interpolation
		T = fA[IP]/fC[IP]

		Tmp = A[IP] EQ C[IP]
		Tmq = 1B-Tmp

		P = Tmp*( 2.*XM[IP]*S ) + Tmq*( S*(2.*XM[IP]*T*(T-R)-(B[IP]-A[IP])*(R-1.)) )
		Q = Tmp*( 1.-S ) + Tmq*( (T-1.)*(R-1.)*(S-1.) )

		Q = (2*(P LE 0)-1)*Q	; !! Don't use 2B and 1B; must be signed integer
		P = abs(P)

		Tmp = 2.*P LT ( (3.*XM[IP]*Q-abs(Tol1+Q)) < (abs(E[IP]+Q)) )
		Tmq = 1B-Tmp

		E[IP] = Tmp*D[IP] + Tmq*XM[IP]
		D[IP] = Tmp*(P/Q) + Tmq*XM[IP]

	ENDIF

	A  = B
	fA = fB
	Tmp = abs(D) GT Tol1		; Evaluate new trial root
	B  = B + Tmp*D + (1B-Tmp)*( abs(Tol1)*(2*(XM GE 0)-1) )

; The following trick was introduced to allow specification of Arg as an array
; of the form Arg[ Dim(X1), *], with separate arguments for each of the
; intervals [X1,X2] for which roots are required.

	Tmp = B[0]+0.*X1			; Same array size as input X1,X2
	Tmp[IZ] = B					; Insert B at proper positions
	fB = call_function(Func, Tmp, Arg)
	fB = fB[IZ]					; Extract fB values needed

	IF n_elements(IT) EQ 0 THEN IT = 0 ; Initialize iteration counter
	IT = IT+1					; Increment iteration counter
ENDREP UNTIL IT EQ ITMax		; At most ITMax iterations

Zero[IZ] = B
message, /info, 'not all roots converged in'+strcompress(ITMax)+' iterations'

return, Zero  &  end