function rtbis,x1,x2,y,func,xacc,_extra=extra
  JMAX=40 ; Maximum allowed number of bisections.
          ; Using bisection, find the root of a function func known to lie between x1 and x2. The root,
          ; returned as rtbis, will be refined until its accuracy is +/-xacc.
  j=0
  dx=0d
  f=call_function(func,x1,_extra=extra)-y
  fmid=call_function(func,x2,_extra=extra)-y
  if (f*fmid ge 0.0) then BEGIN
;    print, get_routine_name() + "Root must be bracketed for bisection in rtbis: func = " + string(func)
    return, 0
  ENDIF
  if f lt 0 then begin ; Orient the search so that f>0 lies at x+dx.
    dx=x2-x1
    rtb = x1
  end else begin
    dx=x1-x2
    rtb=x2
  end
  for j=1,jmax do begin
    dx *= 0.5
    xmid=rtb+dx
    fmid=call_function(func,xmid,_extra=extra)-y ; Bisection loop.
    if (fmid le 0.0) then rtb=xmid;
    if (abs(dx) lt xacc or fmid eq 0.0) then return, rtb
  end
  print, get_routine_name() + "Too many bisections in rtbis: func = " + string(func)
  return, 0
end