pro find_limb, img1, x0, y0, r0, r_err, oblateness, ob_angle, bias, 	$
			brightness,sig_bright,sxt=sxt,qtest=qtest
;+
; NAME:
;	FIND_LIMB
; PURPOSE:
;	Find the solar coordinates from an aspect camera image. 
;	Uses the IDL function SOBEL to differentiate, then fits to circle.
;	Return sun center coordinates, solar radius, error estimate,
;	and oblateness estimate determined from the Fourier spectrum
;	of the limb location. Outputs are in (real) pixel units.
; CALLING SEQUENCE:
;	find_limb, img, 	$
;	   [x0, y0, r0, r_err, oblateness, ob_angle, bias, $
;	   brightness, sig_bright,/sxt,qtest=qtest]
; INPUTS:
;	img 	= input image 
; OUTPUTS:
;	x0	= pixel location of sun center, x (1st harmonic of radius)
;	y0	= pixel location of sun center, y
;	r0	= radius in pixel units
;	r_err	= uncertainty in r0 determination from scatter
;	oblateness 	= second harmonic of radius
;	ob_angle	= phase of 2nd harmonic
;	bias		= distortion of limb due higher harmonics
;	brightness 	= most probable signal at x0,y0
;	sig_bright	= sigma of brightness
; OPTIONAL INPUT PARAMETERS:
;	qtest	= 1 for messages+pause, = 2 for messages+plots+pause,
;		= 3 for messages+plots (no pause)
;	sxt	= 0 for results in units of the pixels of the input image
;		= 1 for results in units of SXT 1x1 pixels (default)
; RESTRICTIONS:
;	1. The input image is assumed to have an oblateness that is < 5%.
;	2. Missing data in the middle of the image will cause large values
;	   in the derivative which are not compensated for.
;	3. The data is assumed to be ge 0.  (Values less than zero will not
;	   be handled correctly by the histogram function).
;	4. In determining summation mode for SXT=1, uses n_elements(img(*,0)).
;	5. To compute brightness, uses a box that is 1/3 * r0 centered at x0,y0.
;	   The image must contain this box (i.e., small crescents may be 
;	    be absent, but a large portion missing at the middle of the
;	    image will cause problems).
;		
; SIDE EFFECTS
;	None
; COMMON BLOCKS
;	None
; MODIFICATION HISTORY
;	HSH written in IDL version 1, Feb. 1991
;	HSH updated with V.2 on real orbital data, Sep. 1991
;	19-Nov-91 MDM - Added a correction factor for changing the
;			resolution back to 1x1 (because of un-summed columns)
;	19-aug-92 JRL+HSH  V3.0 Numerous changes to fix/improve algorithm.
;				Much more robust for SAA data.
;				Brightness is now most probable signal at x0,y0.
;	 1-sep-92 JRL	V3.1	Brightness calculation taken from mean of 
;				histogram
;-
if n_params() eq 0 then begin			; Print information
   print,'find_limb V3.1 (sep-92)'
   print,'find_limb,img,					$'
   print,'	[x0, y0, r0, r_err, oblateness, ob_angle, bias, $'
   print,'	   brightness, sig_bright,/sxt,qtest=qtest]'
   return
endif

; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;  The following parameters control data selection
rad_limit = 0.95		; Radial limit on interior to exclude
bright_limit = .5		; Fraction of peak limb signal to cut
hist_min = .05			; (HSH to explain this)
limb_sharp = 0.01		; fraction width of annulus containing limb at
				;-  final iteration step.
bright_box = 1/3.		; Size of box at x0,y0 for brightness calculation
; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

if n_elements(sxt) eq 0 then sxt=1	; This is the default.  How did MDM let
					;- this definition turn out like this????
; Check the dimensions


dim = size(img1) 
if dim(0) lt 2 then begin
  print,'*** Error in find_limb: Input array must be 2-d'
  x0 = 0. & y0 = 0. & r0 = 0. & r_err = 1000.
  return
endif

numx = dim(1)		; Number of pixels in the x direction

;
; Do Sobel filter. Find limb pixels and centroid.
;
img = img1(*,*,0)
img = sobel(img)

; Try to eliminate the very large Sobel values which might be the result
; of some kind of noise.

tot_pix = n_elements(img)
hh = histogram(img)			; Histogram of Sobel image
j = 1 & not_finished = 1
while( not_finished ) do begin
  ii = where(hh gt j)
  max_sobel = max(  (indgen(n_elements(hh)))(ii)    )
  limb_pixels = where((img gt bright_limit*max_sobel) and (img lt max_sobel))
  num_pix = n_elements(limb_pixels)
  if (j gt 1) then begin	; Assume each loop will result in more pixels
     ccc_percent = (float(num_pix) - prev_num_pix) / prev_num_pix 
     if (ccc_percent lt hist_min) and 			$
		(float(num_pix)/tot_pix gt .005) then not_finished = 0
  endif
  if keyword_set(qtest) then begin
       if j eq 1 then begin
	  ccc_percent = -9999
	  prev_num_pix = num_pix
       endif
       print,'j,prev,num,ccc=',j,prev_num_pix,num_pix,ccc_percent
       if qtest ge 2 then begin
	  xx = indgen(n_elements(hh))
          plot,xx(ii),hh(ii),yrange=[0,200],title='j='+strtrim(j,2)+'; change='+strtrim(ccc_percent,2)
	  bell = string(7b)
	  ans = '' & if qtest lt 3 then read,'Pause '+bell+bell,ans
       endif
  endif
  j=j+1
  prev_num_pix = num_pix
  if j gt 10 then begin		; Prevent an inifinite search situation
     print,'*** Warning:  find_limb did not find a good cutoff in Sobel space'
     print,'***           The resulting fit may be bad'
     not_finished = 0
  endif
endwhile
;
; Calculate positions of limb pixels:

x = limb_pixels mod numx
y = limb_pixels  / numx

; Compute weighted average centroids for first guesses of x0,y0.
totimg = total(img) 	
y0 = total(img#findgen(n_elements(img(0,*))))/totimg	
x0 = total(findgen(n_elements(img(*,0)))#img)/totimg	
r0 = total(sqrt((x-x0)^2 + (y-y0)^2)) / n_elements(x)	; First guess to radius


for i=0,3 do begin		; Fit for solution 3 times

  if i eq 0 then begin	
;    Initial guess
     if keyword_set(qtest) then print,' ** Initial Fit ** '
  endif else if i eq 1 then begin
; Restrict the calculation to pixels near the limb
     true_limb = where(rad gt rad_limit*(total(rad))/n_elements(rad))
     x = x(true_limb) & y=y(true_limb)
     if keyword_set(qtest) then print,' ** 2nd Fit: Eliminate inside of',rad_limit
  endif else begin
; Tigher restriction on the pixels near the limb
     rad2 = (rad - rr(2))/rr(2)		; Residual
     if i eq 2 then xxx = 1-rad_limit else xxx = limb_sharp
     not_spots = where(abs(rad2) lt xxx)
     x = x(not_spots) & y=y(not_spots)
     if keyword_set(qtest) then print,' ** 3rd/4th Fit: Accept within ',xxx
  endelse
  if keyword_set(qtest) then begin
    npix = n_elements(x)
    if i eq 0 then n_start = npix else	 			$
       tot_discard = n_start - npix
    print,'Number of pixels in fit =',strtrim(npix,2),		$
		';  Number discarded=',strtrim(n_start-npix,2)
    if qtest ge 2 then begin
       erase & plots,x,y,psym=1,/dev
    endif
  endif
  rr=fit_circle(x,y,tol=.00001)
  x0 = rr(0) & y0 = rr(1) & r0 = rr(2)		; solutions
  dx=x-x0
  dy=y-y0
  angle = atan(dy,dx)		; = tan(dy/dx)
  rad = sqrt(dx^2+dy^2)
  if keyword_set(qtest) then begin
	print,'Iter:',i,';  x0,y0,r0=',x0,y0,r0,format='(a,i2,a,3f10.2)'
	bell = string(7b)
	ans = '' & if qtest lt 3 then read,'Pause '+bell+bell,ans
  endif
endfor

;
; Analyze limb shape by fitting harmonics (1,2,3,4 * theta).
;
;r0 = total(rad)/n_elements(rad)
r0 = rr(2)
rad = rad - r0
sin_fun = sin(angle)
cos_fun = cos(angle)
sin_1 = poly_fit(sin_fun, rad, 1, sin_fit, yband, ssigma, a)
cos_1 = poly_fit(cos_fun, rad, 1, cos_fit, yband, csigma, a)


r_err = sqrt(ssigma^2+csigma^2)/sqrt(n_elements(rad))
;print, sigma, ' Residuals after fundamental'
; *** What do to ***
rad = rad - (sin_1(0) + cos_1(0))/2 - sin_1(1)*sin(angle)- cos_1(1)*cos(angle)
sin_fun = sin(2*angle)
cos_fun = cos(2*angle)
sin_2 = poly_fit(sin_fun, rad, 1, sin_fit, yband, sigma, a)
cos_2 = poly_fit(cos_fun, rad, 1, cos_fit, yband, sigma, a)
;print, sigma, ' Residuals after 1st harmonic'
rad = rad-(sin_2(0)+cos_2(0))/2 - sin_2(1)*sin(2*angle) - cos_2(1)*cos(2*angle)
sin_fun = sin(3*angle)
cos_fun = cos(3*angle)
sin_3 = poly_fit(sin_fun, rad, 1, sin_fit, yband, sigma, a)
cos_3 = poly_fit(cos_fun, rad, 1, cos_fit, yband, sigma, a)
;print, sigma, ' Residuals after 2nd harmonic'
rad = rad-(sin_3(0)+cos_3(0))/2 - sin_3(1)*sin(3*angle) - cos_3(1)*cos(3*angle)
sin_fun = sin(4*angle)
cos_fun = cos(4*angle)
sin_4 = poly_fit(sin_fun, rad, 1, sin_fit, yband, sigma, a)
cos_4 = poly_fit(cos_fun, rad, 1, cos_fit, yband, sigma, a)
;print, sigma, ' Residuals after 3rd harmonic'
rad = rad-(sin_4(0)+cos_4(0))/2 - sin_4(1)*sin(4*angle) - cos_4(1)*cos(4*angle)
;
; Generate outputs.
;
x0 = x0 + cos_1(1) 
y0 = y0 + sin_1(1) 
ob_angle = atan(cos_2(1)/sin_2(1))
oblateness = sqrt(cos_2(1)^2 + sin_2(1)^2)
bias = sqrt(sin_3(1)^2 + cos_3(1)^2 + sin_4(1)^2 + cos_4(1)^2)

; --------- Compute the brightness at the center --------------------
;	    Take the mean and standard deviation of the distribution

if n_params() ge 9 then begin
  nnx = fix(bright_box * r0 + .5)	; Number of pixels
  img2 = img1(x0-nnx/2:x0+nnx/2,y0-nnx/2:y0+nnx/2)
  hh = histogram(img2)
  xx = min(img2) + indgen(n_elements(hh))
  ii = where(hh gt 0)
  cc = fltarr(4)
  fg = fltarr(4)
  sdev = stdev(img2,mean)
  fg(0) = [max(hh),mean,sdev,0.]	; 1st guesses: peak, centroid, width, background
; yfit = gaussfit2(xx(ii),hh(ii),cc,fg)
  cc = fg				;*** No gaussfit2 call
  gauss_funct2,xx(ii),cc,yfit		;*** Bec. no gaussfit2 call
  brightness = cc(1) 			; Centroid of the histogram
  sig_bright = cc(2)			; 1/e Gaussian width of histogram 
  if keyword_set(qtest) then begin
     print,'Gaussfit results=',cc
     if qtest ge 2 then begin
	plot,xx,hh,psym=10
     	oplot,xx(ii),yfit		; Show the fitted result
     endif
     bell = string(7b)
     ans = '' & if qtest lt 3 then read,'Pause '+bell+bell,ans
   endif				; keyword_set(qtest)
endif					; n_params() ge 9

; --------- Convert output to SXT 1x1 pixels ------------------------

if keyword_set(sxt) then begin
  off_corr = [-1, 0, 1, -1, 3]		;MDM added 19-Nov-91
  ;changing from 1x1 to 1x1, no offset
  ;changing from 2x2 to 1x1, offset = 1 (1x1 pixels 1&2 are 2x2 pixel 0)
  ;changing from 4x4 to 1x1, offset = 3 (1x1 pixels 3,4,5&6 are 4x4 pixel 0)

  sum = 1024/numx		; Summation for SXT (use number of x pixels)
  x0 = x0 * sum + off_corr(sum)
  y0 = y0 * sum
  r0 = r0 * sum
  r_err = r_err * sum
  oblateness = oblateness * sum
  bias = bias * sum
endif

return
end