pro single_test,fname

read_hr_ha_new,fname,img,header
xdim=sxpar(header,'NAXIXX')
ydim=sxpar(header,'NAXIYY')
factor=0.35
xdim1=xdim*factor
ydim1=ydim*factor
window,xs=xdim1,ys=ydim1
tvscl,congrid(img,xdim1,ydim1)

p   = { NX: 2008, NY: 2008, XC: 1000, YC: 1000, XD: 1910, YD: 1910}
	; first guest for the parameter
mask=img
lim=mean(mask)/6.
vmin=5
mask = ( SOBEL( TEMPORARY( mask ) ) GT lim )
;----- elminate problems with borders and a few bad colums
b = 4
;----- get x-coordinates and diameter
xx = TOTAL( mask, 2 ) GT vmin  ;----- horizontal sum
;----- solar image is centered in x
IF ( xx( b ) EQ 0 ) AND ( xx( p.nx - b - 2 ) EQ 0 ) THEN BEGIN
    index = WHERE( xx eq 1, n )
    x1 = index( 0 )
    x2 = index( n - 1 )
    xd = x2 - x1 + 1
    xc = xd / 2 + x1
ENDIF

;----- solar image is on the left

IF ( xx( b  ) EQ 1 ) AND ( xx( p.nx - b - 2 ) EQ 0 ) THEN BEGIN
    xx( 0 : b - 1 ) = 1    ;----- set border to 1
    index = WHERE( 1 - xx )
    x2 = index( 0 ) - 1
    xd = -1
    xc = x2 - p.xd / 2
ENDIF

;----- solar image is on the right

IF ( xx( b  ) EQ 0 ) AND ( xx( p.nx - b - 2 ) EQ 1 ) THEN BEGIN
  index = WHERE( xx )
  x1 = index( 0 )
  xd = -1
  xc = x1 + p.xd / 2
ENDIF

;----- get y-coordinates and diameter

yy = TOTAL( mask, 1 ) GT vmin    ;----- vertical sum

;----- solar image is centered in y

IF ( yy( b ) EQ 0 ) AND ( yy( p.ny - b - 2 ) EQ 0 ) THEN BEGIN
    index = WHERE( yy, n )
    y1 = index( 0 )
    y2 = index( n - 1 )
    yd = y2 - y1 + 1
    yc = yd / 2 + y1
ENDIF

;----- solar image is on the bottom

IF ( yy( b ) EQ 1 ) AND ( yy ( p.ny - b - 2 ) EQ 0 ) THEN BEGIN
  yy( 0 : b - 1 ) = 1    ;----- set border to 1
  index = WHERE( 1 - yy )
  y2 = index( 0 ) - 1
  yd = -1
  yc = y2 - p.yd / 2
  IF yc LT 0 THEN yc = y2 + yc		;---- solves some occacional problems

ENDIF

;----- solar image is on the top

IF ( yy( b ) EQ 0 ) AND ( yy( p.ny - b - 2 ) EQ 1 ) THEN BEGIN
  index = WHERE( yy )
  y1 = index( 0 )
  yd = -1
  yc = y1 + p.yd / 2
ENDIF

;;------ here starts the new part (M. Steinegger, 2000/07/20)

;----- check initial guesses for diameter
IF xd EQ -1 AND yd EQ -1 THEN xd = p.xd
IF xd EQ -1 AND yd NE -1 THEN xd = yd
IF yd EQ -1 AND xd NE -1 THEN yd = xd

;-----  eliminate center of the disk (active regions!)

r  = FIX( ROUND( ( p.xd + p.yd ) / 4. ) )
DIST_CIRCLE, d, [ p.nx, p.ny ], xc, yc
;;mask( WHERE( d LT r * 0.95 OR d GT r * 1.1 ) ) = 0
mask( WHERE( d LT r * 0.9 OR d GT r * 1.1 ) ) = 0	; for KSO images

;-----  fit a circle to obtain radius and center coordinates,
;       iterate until all points are within 2 * sigma of radius

bad = [0L]
WHILE bad( 0 ) NE -1 DO BEGIN


   wo = WHERE( mask GT 0 )
   xfit = wo MOD p.nx
   yfit = wo / FLOAT( p.nx )
   result = FIT_CIRCLE( xfit, yfit, [ xc, yc, r ] )
   radius = SQRT( ( xfit - result( 0 ) ) ^ 2. + ( yfit - result( 1 ) ) ^ 2. )
   sigma  = SQRT ( TOTAL( ( radius - result( 2 ) ) ^ 2. ) / $
                 ( N_ELEMENTS ( wo ) - 1 ) )
   bad = WHERE( radius GT ( result( 2 ) + 2. * sigma ) OR $
                radius LT ( result( 2 ) - 2. * sigma ) )
   IF bad( 0 ) NE -1 THEN BEGIN
      mask( xfit( bad ), yfit( bad) ) = 0
      xc = FIX( ROUND ( result( 0 ) ) )
      yc = FIX( ROUND ( result( 1 ) ) )
      xd = FIX( ROUND ( result( 2 ) * 2. ) )
      yd = xd
;;     r  = ( xd + yd ) / 4.
      r  = ( p.xd + p.yd ) / 4.

;;      print,sigma,bad(0),N_ELEMENTS(bad), xc, yc, xd

   ENDIF
ENDWHILE

yd = xd

help,xc,yc,xd,yd
dr_cir1,[xc,yc,xd/2.]*factor,xmax=xdim1,ymax=ydim1

END











