;+
; NAME:
;     clean_beam
; PURPOSE:
;     determine clean beam out of a dirty beam
; CATEGORY:
;     OVSA, APC Imaging
; CALLING SEQUENCE:
;     cbm = clean_beam(dbeam)
; INPUTS:
;     dbeam: dirty beam.
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
;     none
; ROUTINES CALLED:
;
; OUTPUTS:
;     cbm[0]= a, semi-minor axis in pixel units.
;     cbm[1]= b, semi-major axis in pixel units.
;     cbm[2]= p, orientation of b w.r.t y-axis in degree
;             clockwise=-, c.c.wise=+
;
; COMMENTS:
;     The central peak region is searched downward
;     until 0. or increase is encounered.
;     For the selected peak region, logarithmic least square gaussian fit,
;     which minimizes total((B * (a[0] * x^2 + a[1] * xy + a[2] * y^2 + a[3] + lnB))^2),
;     is calculated analytically.
;     Then a is converted to cbm.
; SIDE EFFECTS:
; RESTRICTIONS:
;     Dirty beam should be square and symmetric in respect to the center
;     and the center should be a local maximum.
; MODIFICATION HISTORY:
;     Written 27-JUN-2001 by Su-Chan Bong.
;-
function clean_beam, dbeam
    sz = size(dbeam)
    m = sz[1]
    m1 = m - 1l
    m2 = m / 2l
    m1p = m + 1l

    cen = m2 * m1p
    xmax = cen  ;find the primary beam
    pri = xmax
    dbpre = dbeam[xmax]
    dn = [-m1p, -m, -m1, -1l, 1l]
    up = [m1, m, m1p]
    nbr = xmax + dn
    outflag = bytarr(m, m2 + 1l) + 1b
    outflag[pri] = 0b
    outflag[nbr] = 0b
    go = 1b
    while go do begin
        dbnow = max(dbeam[nbr], inow)
        xnow = nbr[inow]
        if dbnow le 0. or dbnow gt dbpre then begin
            go = 0b
        endif else begin    ;accept as a part of the primary beam
            dbpre = dbnow
            pri = [pri, xnow]
            nbr = shift(nbr, -inow)
            nbr = nbr[1 : *]
            nbrcan = xnow + dn
            if xnow / m lt m2 then nbrcan = [nbrcan, xnow + up]
            xnbrcan = nbrcan mod m
            edge = where(xnbrcan eq 0l or xnbrcan eq m1 or $
             nbrcan lt m, nedge)
            if nedge gt 0l then begin
                go = 0b
            endif else begin
                new = where(outflag[nbrcan], nnew)
                if nnew gt 0l then begin
                    nbrnew = nbrcan[new]
                    nbr = [nbr, nbrnew]
                    outflag[nbrnew] = 0b
                endif
            endelse
        endelse
    endwhile
    pri = pri[where(pri le cen)]
    pb = dbeam[pri]

    t2 = pb^2       ;fit the clean beam
    t2lnt = t2 * alog(pb)
    x = pri mod m - m2
    y = pri / m - m2
    xy = x * y
    x2 = x^2
    y2 = y^2
    t2x2 = t2 * x2
    t2y2 = t2 * y2
    st2x2 = total(t2x2)
    st2xy = total(t2 * xy)
    st2y2 = total(t2y2)
    st2x4 = total(t2x2 * x2)
    st2x3y = total(t2x2 * xy)
    st2x2y2 = total(t2x2 * y2)
    st2xy3 = total(t2y2 * xy)
    st2y4 = total(t2y2 * y2)

    b = [[st2x4, st2x3y, st2x2y2, st2x2], $
         [st2x3y, st2x2y2, st2xy3, st2xy], $
         [st2x2y2, st2xy3, st2y4, st2y2], $
         [st2x2, st2xy, st2y2, total(t2)]]
    c = -[total(t2lnt * x2), total(t2lnt * xy), total(t2lnt * y2), total(t2lnt)]
    a = invert(b) ## c

    co2 = a[2] - a[0]
    si2 = -a[1]
    a20 = a[2] + a[0]
    wx = sqrt(2. / (a20 - sqrt(co2^2 + si2^2)))
    wy = sqrt(2. / (a20 + sqrt(co2^2 + si2^2)))
    pa = (atan(si2, co2) / 2. * !radeg + 180.) mod 180. - 90.

    return, [wy, wx, pa]
end