FUNCTION rtbis,x1,x2,func,flg,xacc,itmax
;
; Using BISECTION, find the root of a function FUNC known to lie
; between X1 and X2. The root will be refined until its accuracy
; is +/- XACC. The calculation will finish after ITMAX iterations.
; In case of a failure (non-existing root, too many iterations),
; RTBIS returns zero and FLG=0. When root is found, FLG=1.
; FUNC: A scalar string specifying the name of a user-supplied IDL
;       function that defines the nonlinear function.
;       This function must accept an argument X.
; Parameters FLG, XACC, ITMAX are optional.  
;
; EXAMPLE: root=RTBIS(x1,x2,'func',flg,1.e-5,20)
;
; Based on the routine RTBIS from Press et al.: Numerical Recipes,
; Fortran Version (Cambridge, 1989), p. 247.
; 14.4.2004, Michal

; Defaults

if n_params() lt 6 then itmax=40
if n_params() lt 5 then xacc=1.e-4
flg=1

; Does the root exist between x1 and x2 ?

fmid=call_function(func,x2)
f=call_function(func,x1)
if (f*fmid ge 0.) then goto,bad		;root is not there

; Orient the search so that f>0 lies at x+dx

if (f lt 0) then begin
  root=x1
  dx=x2-x1
endif else begin
  root=x2
  dx=x1-x2
endelse

; Bisection loop

for i=1,itmax do begin
  dx=dx*0.5
  xmid=root+dx
  fmid=call_function(func,xmid)
  if (fmid le 0.) then root=xmid
  if (abs(dx) lt xacc) or (fmid eq 0.) then RETURN,root
endfor

print,''
print,' RTBIS: Too many iterations'

bad:
  print,''
  print,' RTBIS: Root not found'
  print,''
  flg=0
  RETURN,0

END