;+
; NAME:
;     clean_bijm
; PURPOSE:
;     determine clean beam out of a dirty beam
; CATEGORY:
;     OVSA, APC Imaging
; CALLING SEQUENCE:
;     cbm = clean_bijm(uv_ij,dbeam)
; INPUTS:
;     uv_ij: uv points
;     dbeam: dirty beam.
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
;     none
; ROUTINES CALLED:
;
; OUTPUTS:
;     cbm[0]= a, semi-major/minor axis closer to x-axis in pixel units
;     cbm[1]= b, semi-minor/major axis closer to y-axis
;     cbm[2]= p, orientation of b w.r.t y-axis in degree
;    		  clockwise=-, c.c.wise=+
;
; COMMENTS:
;     Instead of a widely-used general gaussian fitting algorithms
;     here use is made of u,v data in such a way that the uv points with
;     the largest distance determines the orientation of the ellipse.
; SIDE EFFECTS:
; RESTRICTIONS:
;     assume 0.1 <= a, b <= 10.1
;     accuracy within 0.01 pixel units
;     use all but 1-5, 1-6 bsls for determining the semi-minor axis
; MODIFICATION HISTORY:
;     Written 05-JAN-2001 by JL
;-

function clean_bijm, uv_ij, dbeam
 m  = n_elements(dbeam(*,0))
 nt = n_elements(uv_ij(*,0))
 nb = n_elements(uv_ij(0,*))

;--------determine the position angle of semi-major axis

bsl =[0,1,4,5,6] ; meant to exclude 1-5 & 1-6 baselines in judging semi-minor axis
uv  =[uv_ij[*,bsl],-uv_ij[*,bsl]]
up  = where(imaginary(uv) ge 0.)

uv_ = uv(up)

pp_ = -atan(float(uv_),imaginary(uv_)) / !dtor
; atan itself give +ve value at the 1st quadrant and -ve at the 2nd quad.
; to follow the sign convention in the above we flip the sign
; the next line is then put all numbers in the 2nd & 3rd quad.
nok=where(pp_ lt 0.,nbad) & if (nbad ne 0) then pp_(nok)=pp_(nok)+180.
ww_=abs(uv_)	; weights according to the uv-distance from (0,0)

pp = total(pp_*ww_)/total(ww_)

;--------Actually we want p angle of the semi-major axis, thus

pp = pp - 90.

;--------Determine semi-major and minor axes

sb = rot(dbeam(m/2-20:m/2+20, m/2-20:m/2+20),pp,/cubic)

al = 0.1+(10.1-0.1)*findgen(1001)/1000.

; along x-direction - usually semi-minor axis
yy = reform(sb(20:*,20),21) > 0.
   ds  = shift(yy,1)-yy & ds[0]=1.0
  bnd  = min(where (ds le 0.))-1
;  yy[bnd:*]=0. ; this is worse than the following
yy=yy[0:bnd-1]
ss=findgen(bnd)
  chi2 = fltarr(1001)
  for i=0,1000 do chi2(i)=total( $
	(yy-1.0/exp((ss/al(i))^2<55.))^2 )
  aa=al(where(chi2 eq min(chi2)))

; along y-direction - usually semi-major axis
yy = reform(sb(20,20:*),21) > 0.
   ds  = shift(yy,1)-yy & ds[0]=1.0
  bnd  = min(where (ds le 0.))-1
;  yy[bnd:*]=0.
yy=yy[0:bnd-1]
ss=findgen(bnd)
  chi2 = fltarr(1001)
  for i=0,1000 do chi2(i)=total( $
	(yy-1.0/exp((ss/al(i))^2<55.))^2 )
  bb=al(where(chi2 eq min(chi2)))

return,[aa[0],bb[0],pp]

end