;------------------------------------------------------------------
;+
; NAME:
;	MCLEANUP
; PURPOSE:
;       This program "cleans up" a BBSO VMG polar mosaic image.
;       Also, pixels beyond the limb in the original mosaic are set to 0.
;
; CATEGORY:
;	POLAR
; CALLING SEQUENCE:
;	MCLEANUP,inimage,imh,h,k,imgscl
; INPUTS:
;       h       = x-coordinate of center of Sun
;       k       = y-coordinage of center of Sun
;       inimage = input VMG mosaic image.
;	imh     = BBSO VMG image FITS header (used to get date and
;	          time).
;       imgscl  = image scale (pixels / arcsec).
;
;       h and k are measured with (0,0) at the lower left corner of
;       the image.  Correct values are around (1500,0) for north pole
;       mosaics and (1500,1500) for south pole mosaics.
;       This function will work with mosaics of any size, however.
; KEYWORD PARAMETERS:
;       CUTOFF - if present, use this as factor above std. dev. for
;                cutoff instead of 3.0
; OUTPUTS:
;       Modifies input image; points beyond limb are set to -32000,
;       points within image are set to 0 if abs. value is within
;       3 * std dev. from 0.  Finally the remaining points are
;       scaled by 1/cos(r).
; COMMON BLOCKS:
; NOTES:
;       This function uses array operations for efficiency!
; MODIFICATION HISTORY:
;       J. Varsik,      29 Mar. 1999
;       J. Varsik,      19 Apr. 1999 Changed comments to refer to
;                                    3000x1500 mosaic.  No change to
;                                    code.
;       J. Varsik,      22 May  1999 Add CUTOFF keyword.
;-
;-------------------------------------------------------------------

PRO MCLEANUP,inimage,imh,h,k,imgscl,help=help,cutoff=cutoff

    ; Display IDL header if help is required.
    IF (KEYWORD_SET(help)) THEN BEGIN
      GET_IDLHDR,'mcleanup.pro'
      GOTO,finishup
    ENDIF

    ; Get input image size
    isz = SIZE(inimage)
    ixz = isz(1)
    iyz = isz(2)
    ; Create output image
    outimage = FLTARR(ixz,iyz)
    ; Read input image info.
    ; Find center, radius, p angle, b angle, etc.
    datec = FXPAR(imh,'DATE-OBS')
    timec = FXPAR(imh,'TIME-OBS')
    IF (STRMID(datec, 2, 1) EQ '/') THEN BEGIN
        READS,datec,day,mon,yrs,FORMAT='(i2,1x,i2,1x,i2)'
        yrs = yrs + 1900
    ENDIF ELSE BEGIN
        READS,datec,yrs,mon,day,FORMAT='(i4,1x,i2,1x,i2)'
    ENDELSE
    READS,timec,hour,min,sec,FORMAT='(i2,1x,i2,1x,i2)'
    hr = hour + ((min + (sec / 60.0)) / 60.0)
;   Use solar ephemeris from Johns Hopkins library
    sun,yrs,mon,day,hr,pa=p,lat0=b,sd=rsun

    rad = imgscl * rsun

;   find the part of the mosaic that is filled
    filled = WHERE(inimage NE -32000, fc)
    IF fc EQ 0 THEN BEGIN
        PRINT,'empty mosaic'
        GOTO,finishup
    ENDIF
    PRINT,fc
    ; reconstruct array subscripts and put *them* in arrays
    i = filled MOD ixz
    j = filled / ixz
    xfc = (i+1)-h
    yfc = (j+1)-k
    r1 = sqrt((xfc)^2 + (yfc)^2)
    outimage[filled] = r1  ;makes image with radial distance for each point
    big = WHERE(outimage GE rad)
    ; outimage[big] = 0.0    ;remove points beyond limb
    inimage[big] = -32000  ;remove noise beyond limb of mosaic

    ; measure noise level
    used = WHERE(inimage GT -32000)
    ondisk = inimage(used)
    tmask = ondisk
    low = WHERE(ondisk LT -30)
    ondisk(low) = 300
    low = WHERE(ondisk GT 30)
    tmask(low) = 0
    magstat = MOMENT(tmask)
    stddevmag = SQRT(magstat[1])
    PRINT,'mean background ',magstat[0],'std dev ',stddevmag

    ; set points on disk with abs value less than 3 * std dev to 0
    ; (or less than value of cutoff if set).
    IF (KEYWORD_SET(cutoff)) THEN BEGIN
        noise = WHERE(abs(inimage) LT (cutoff * stddevmag))
    ENDIF ELSE BEGIN
        noise = WHERE(abs(inimage) LT (3 * stddevmag))
    ENDELSE
    inimage[noise] = 0

    ; Now scale the magnetic field by 1/cos(r)
    usedr = FLOAT(outimage(used))
    usedr = 1.0 / SQRT(1.0 - (usedr^2 / rad^2))
    inimage(used) = FIX(FLOAT(inimage(used)) * usedr)

    finishup:

    RETURN
END