pro circle_fit, xdata, ydata, x0i, y0i, r0i, x0o, y0o, ro, check, verbose=verbose,quiet=quiet,norfit=norfit,noiterate=noiterate


   iter_limit = 2000  ; Limit to the number of iterations
   x0o = double(x0i)
   y0o = double(y0i)
   ro = double(r0i)
   x = double(xdata)  ; protect the input
   y = double(ydata)  ; protect the input
   nx = n_elements(x)
   count = 0

   REPEAT begin    ; Iterate to the best fit

      count = count+1
      x0o_old = x0o
      y0o_old = y0o

      theta = atan((y-y0o),(x-x0o))
      fit = poly_fit(cos(theta),x,1)
      x0o = fit(0)
      fit = poly_fit(sin(theta),y,1)
      y0o = fit(0)

      ; Compute a new value for the radius, consistent with the new values
      ; of the x and y offsets.

      if NOT keyword_set(norfit) then ro=total(sqrt((x-x0o)^2 + (y-y0o)^2))/nx

      if keyword_set(verbose) then print,x0o,y0o

   endrep UNTIL (abs(x0o-x0o_old) LT 1.0e-5 AND abs(y0o-y0o_old) LT 1.e-5) $
                OR count GT iter_limit OR keyword_set(noiterate)

   if count GT iter_limit then check = 0 else check = 1

   if NOT keyword_set(norfit) then ro = total(sqrt((x-x0o)^2 + (y-y0o)^2))/nx
   step = sqrt(total(ABS((x-x0o)^2 + (y-y0o)^2 - ro^2))/nx)

   if (NOT check) and (NOT keyword_set(quiet)) then $
       print,"WARNING: circle_fit: No Convergence."

   if keyword_set(verbose) then begin
      print,'count=',count
      print,'Fit parameters:'
      print,x0o,y0o,ro,step/nx, $
      format="('x0 = ',f10.4,' y0 = ',f10.4,' r = ',f10.4,' residual = ',f10.4)"
      if check then begin
         erase
         psave = !P
         xsave = !X
         ysave = !Y
         plot,x,y,/nodata
         plots,x,y,psym=1
         theta = findgen(361)*!dtor
         oplot,(ro*cos(theta)+x0o),(ro*sin(theta)+y0o)
         !Y = ysave
         !X = xsave
         !P = psave
      endif
   endif

   x0o = float(x0o)
   y0o = float(y0o)
   ro = float(ro)

end