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