pro scc_mk_daily_med,tel,sc,Date,polar=polar, DOREBIN=dorebin, SAVEDIR=savedir,FILES=files, VERBOSE=verbose
;
;+
; $Id: scc_mk_daily_med.pro,v 1.19 2008/12/11 19:18:51 nathan Exp $
;
; Project   : STEREO SECCHI
;                   
; Name      : SCC_MK_DAILY_MED
;               
; Purpose   : This procedure generates an image by taking all the files of a given
;	     type (up to 40) for one day and finding the median value for each pixel.
;               
; Explanation: 
;               
; Use       : IDL> scc_mk_daily_med,'CAM','SC','YYYYMMDD',polar='POL'
;    
; Inputs    :CAM = cor1, cor2, hi1, hi2
;            SC = 'A', 'B'
;            YYYYMMDD = date to be processed
;            polar= '0', '120', '240', 'tbr'  (tbr = total brightness)
;               
; Outputs   : FITS file in $SECCHI_BKG/[ab]/daily_med/YYYYMM/
;
; Keywords  : SAVEDIR: specify your own directory to save file to
;             FILES:   specify your own file list (this may not work properly yet!!)
;   	    /DOREBIN	output 512x512 instead of  1024x1024
;
;
; Calls from LASCO : 
;
; Common    : 
;               
; Restrictions: Need $SECCHI_BKG and appropriate permissions if not using savedir option
;               
; Side effects: 
;               
; Category    : DAILY
;               
; Prev. Hist. : None.
;
; Written     : Karl Battams, NRL/I2, Jan 2007
;               
; $Log: scc_mk_daily_med.pro,v $
; Revision 1.19  2008/12/11 19:18:51  nathan
; do not include zero values in median
;
; Revision 1.18  2008/08/06 19:16:58  nathan
; correct cor1 inp directory; fix bug if savedir=.
;
; Revision 1.17  2008/02/12 14:55:18  secchia
; nr - added /new to scc_putin_array call
;
; Revision 1.16  2008/02/11 15:31:06  secchia
; nr - fixed case where ng=0
;
; Revision 1.15  2007/11/19 19:22:05  nathan
; use normal SECCHI_BKG directories
;
; Revision 1.14  2007/11/01 21:49:14  nathan
; added image roll check and dstart/stop correction
;
; Revision 1.13  2007/10/26 22:24:12  nathan
; Added /DOREBIN; added cadence, readtime, cleartim, crval, nmissing, cosmics
; to output header; tried to get dsatval right; use start of first image for
; DATE-OBS and midpoint is DATE-AVG; scale COR2 images to type INT before
; saving; save HI as type FLOAT; saves results in $SECCHI_BKG/../newbkg
;
; Revision 1.12  2007/09/25 22:55:52  nathan
; mostly done--still need to update scc_monthly_min.pro
;
; Revision 1.11  2007/09/25 16:00:11  nathan
; This has been re-written to utilize secchi_prep and the full SECCHI FITS header
; and to resolve Bug 224
;
; Revision 1.10  2007/08/16 15:08:12  nathan
; made sure default output dir is group writable
;
; Revision 1.9  2007/07/25 14:15:05  reduce
; needed uppercase S/C identifier. Karl B
;
; Revision 1.8  2007/07/05 19:03:36  reduce
; Fixed 0-deg polarizer filename error. Karl B.
;
; Revision 1.7  2007/07/02 19:28:26  reduce
; Added a 'cd,orig' cmd. Karl B.
;
; Revision 1.6  2007/06/27 14:01:15  reduce
; Exits gracefully if too few files found. Karl B.
;
; Revision 1.5  2007/06/22 20:46:25  nathan
; cd,orig before return in Error case
;
; Revision 1.4  2007/06/21 19:20:00  reduce
; No longer crashes on missing images. Karl B.
;
; Revision 1.3  2007/06/18 20:48:57  reduce
; Couple of bug fixes. Should be ready for release now. Karl B.
;
; Revision 1.2  2007/02/01 20:43:52  reduce
; More mods and bug fixes. KB
;
; Revision 1.1  2007/01/30 20:29:21  reduce
; Initial Release by KB
;
;-
	
;
; PROCEDURE:
;	For each image that satifies the selection conditions, (default naxis1=1024,
;       filter and polarizer as requested), the median image is computed of 
;	the median value of all the images for a single day after being
;	normalized to the median exposure time.  
;
;	If the number of images is less than 7, there is a second pass. 
;  	In the second pass, images within +/- 2 days of the given day are used, 
;	up to 15 per day.  
;

;
;   Convert the telescope number into lower case
;   And select the standard size parameters
;

IF (STRLOWCASE(tel) EQ 'cor1') THEN BEGIN
    PRINT,'
    PRINT,'**********************************************************'
    PRINT,' Creation of median images for COR1 is not yet implemented
    PRINT,' as per the wishes of the COR1 instrument team. 
    PRINT,'**********************************************************'
    PRINT,'
    return
ENDIF

IF (n_params() LT 3) then BEGIN
    PRINT,'Incorrect number of input parameters.'
    PRINT,'
    PRINT,'Usage:
    PRINT,'   scc_mk_daily_med,CAM,S/C,DATE,POLAR
    PRINT,'EXAMPLE:
    PRINT,'   IDL> scc_mk_daily_med,"hi1","a","20070115",polar="tbr"'
    PRINT,'    --- OR ---
    PRINT,'   IDL> scc_mk_daily_med,"cor2","b","20070222",polar="120"'
    PRINT,'        
ENDIF

version='$Id: scc_mk_daily_med.pro,v 1.19 2008/12/11 19:18:51 nathan Exp $'
len=strlen(version)
version=strmid(version,1,len-2)

n_median_min = 5		; number of images needed to generate median
cam = STRLOWCASE(tel)
xaxis=1
yaxis=1
cor_flag=0 ; set to 1 if CORs are chosen

if not keyword_set(polar) then BEGIN
    PRINT,'Keyword "polar" needs to be defined.
    PRINT,'Options are as follows:'
    PRINT,'     "tbr" -- total brightness (COR or HI)
    PRINT,'     "0", "120" or "240" -- polarized (COR only)
ENDIF

pol=strlowcase(polar)

IF keyword_set(VERBOSE) THEN strn = ' Using:' ELSE strn = ''

if (strlen(date) ne 8) THEN BEGIN
    PRINT,'Date must be of the form YYYYMMDD'
    RETURN
ENDIF

CASE cam OF
	'cor1':  BEGIN
		if pol EQ 'tbr' THEN BEGIN
                    s=GETENV('SECCHI_P0')+'/'+strcompress(strlowcase(sc),/remove_all)+'/cor1/' 
                    pol_str='pTBr' 
                    cor_flag=1
                ENDIF ELSE BEGIN
                    s=GETENV('secchi')+'/lz/L0/'+strcompress(strlowcase(sc),/remove_all)+'/seq/cor1/'
                    pp=strcompress(pol,/remove_all)
                    IF (float(pp) EQ 0) THEN pol_str='p000' ELSE pol_str='p'+pp
                ENDELSE
                pref='c1'  ; filename prefix

                IF not keyword_set(anysize) THEN BEGIN 
                    xaxis=2048
                    yaxis=2048
                ENDIF
	END
	'cor2':  BEGIN
		stdexptime=4. ; am I using this?
                if pol EQ 'tbr' THEN BEGIN
                    s=GETENV('SECCHI_P0')+'/'+strcompress(strlowcase(sc),/remove_all)+'/cor2/' 
                    pol_str='pTBr' 
                    cor_flag=1
                ENDIF ELSE BEGIN
                    s=GETENV('secchi')+'/lz/L0/'+strcompress(strlowcase(sc),/remove_all)+'/seq/cor2/'
                    pp=strcompress(pol,/remove_all)
                    IF (float(pp) EQ 0) THEN pol_str='p000' ELSE pol_str='p'+pp
                ENDELSE
                pref='c2'  ; filename prefix

                IF not keyword_set(anysize) THEN BEGIN 
                    xaxis=2048
                    yaxis=2048
                ENDIF

	END
	'hi1':  BEGIN
		IF KEYWORD_SET(polar)   THEN pol=strupcase(polar) ELSE pol=''
		s=GETENV('secchi')+'/lz/L0/'+strcompress(strlowcase(sc),/remove_all)+'/img/hi_1/'
                pref='h1'
                pol='tbr'
                pol_str='pTBr'
                IF not keyword_set(anysize) THEN BEGIN 
                    xaxis=1024
                    yaxis=1024
                ENDIF
	END
	'hi2':  BEGIN
                pol=''
                s=GETENV('secchi')+'/lz/L0/'+strcompress(strlowcase(sc),/remove_all)+'/img/hi_2/'
                pref='h2'
                pol='tbr'
                pol_str='pTBr'
                IF not keyword_set(anysize) THEN BEGIN 
                    xaxis=1024
                    yaxis=1024
                ENDIF
	END
	ELSE:  BEGIN
		PRINT,'Unrecognized telescope code: '+cam
		RETURN
	END
ENDCASE
img_type=''

IF (datatype(date) NE 'STR') THEN BEGIN
    PRINT,'Input date must be a string E.G. 20061223'
    RETURN
ENDIF

dt2=strmid(date,0,4)+'-'+strmid(date,4,2)+'-'+strmid(date,6,2)

;  stop
IF keyword_set(SAVEDIR) THEN BEGIN
    ; in case savedir= '.'
    cd,savedir
    spawn,'pwd',outp,/sh
    savedir=outp[0]
ENDIF
; Check if data exists yet...
; s=GETENV('secchi')+'/lz/L0/'+strcompress(strlowcase(sc),/remove_all)+'/'+imgdir ; got here **************
CD,s,curr=orig
f=file_search(date+'/*fts')
sz=size(f)
IF (sz(0) EQ 0)  THEN BEGIN
	PRINT,'No directory for '+cam+' telescope on '+date
	CD,orig
	RETURN
ENDIF

spawn,'pwd',/SH


; Find appropriate files...
IF keyword_set(FILES) THEN BEGIN
	fnames = files 
	ng=n_elements(files)
	good = indgen(ng) 
	print,'Using ',trim(ng),' files for median.'
ENDIF ELSE BEGIN

        if (not cor_flag) then summary=scc_read_summary(DATE=dt2,SPACECRAFT=sc,TELESCOPE=strupcase(tel)) $
            else summary=scc_read_summary(DATE=dt2,/totalb,SPACECRAFT=sc)
;stop
        PRINT,'Found a total of ',strcompress(n_elements(summary)),' files for the day in the summary file.'

        IF (n_elements(summary) LE 5) THEN BEGIN
            PRINT,'Based on the summary file, there do not appear to be enough files for making a median for today on this telescope.'
	    CD,orig
            return
        ENDIF

	IF datatype(summary) EQ 'INT' THEN BEGIN
		CD,orig
		return
	ENDIF

        ; Either an individual COR polarizer angle is picked, or a HI image (no polarizer) or total brightness is picked. 

        if (pol NE 'tbr') then good = where(summary.VALUE EQ pol and summary.XSIZE GE xaxis and summary.YSIZE EQ yaxis and summary.PROG NE 'Dark' and summary.PROG NE 'Doub', ng) $
            else if not cor_flag then good = where(summary.XSIZE GE xaxis and summary.YSIZE EQ yaxis and summary.DEST EQ 'SSR1'and summary.PROG NE 'Dark', ng) $
                else good = where(summary.VALUE EQ 1001 and summary.XSIZE GE xaxis and summary.YSIZE EQ yaxis and summary.PROG NE 'Dark', ng)

    	print,'Found ',trim(ng),' appropriate images for ',dt2,' from summary file.'
	IF ng EQ 0 THEN BEGIN
    	    print,'ERROR: There are not enough acceptable images to make a good daily median image. Returning.'
    	    wait,5
    	    CD,orig
    	    return
    	ENDIF ELSE files=summary[good].FILENAME  

ENDELSE

path=s+date+'/'
CD,path,curr=old

print,'Screening file headers for acceptability,,,'
;  - Figure out what roll should be
ut=anytim2utc(dt2)
uta=replicate(ut,5)
uta.mjd=uta.mjd + indgen(5)-2
ra=-1*get_stereo_roll(uta,sc,/verbose)
minr=min(ra,max=maxr)
help, minr,maxr
;
i=1
f=0
;stop
WHILE i LE 40 and f LT ng DO BEGIN
    ; There are more files than we need so lets take the first and last 20 of the day.
    ;tmp1=files[0:19]
    ;tmp2=files[n-20:n-1]
    ;files=[tmp1,tmp2]
    ;print,'reading ',files[f]
    jnk=sccreadfits(files[f],h,/nodata)
    IF h.sc_roll GT maxr or h.sc_roll LT minr THEN $
    	print,'Not using ',files[f],' because SC_ROLL=',h.sc_roll $
    ELSE BEGIN
    	IF i EQ 1 THEN gfiles=files[f] ELSE gfiles=[gfiles,files[f]]
	i=i+1
    ENDELSE
    f=f+1
ENDWHILE

files=gfiles
n=n_elements(files)
if n LT n_median_min then begin
    print,'ERROR: There are not enough acceptable images to make a good daily median image. Returning.'
    wait,5
    CD,orig
    return
endif ELSE PRINT,'Using ',trim(n),' files for median.'

datacube=   fltarr(1024,1024,n)
crotas=     fltarr(n) 
biases=     fltarr(n)
exptimes=   dblarr(n)
offsetcrs=  fltarr(n)
cadences=   fltarr(n-1)
readtimes=  dblarr(n)
cleartimes= dblarr(n)
pitches=    fltarr(n)
yaws=	    fltarr(n)

nmissings=  lonarr(n)
nsats=	    lonarr(n)
cosmicss=   lonarr(n)

midpoint=n/2

n_used=n ; used later

; This is a temporary bug fix for the COR total-B summary files
last_char=strmid(files[0],strlen(files[0])-1,1)
if (last_char) NE 's' then files=files+'s'
; end bug fix

IF keyword_set(VERBOSE) THEN quiet=0 ELSE quiet=1

FOR m=0,n-1 DO BEGIN

    if file_exist(files[m])   THEN fn=files[m]   ELSE $
    if file_exist(files[m-1]) THEN BEGIN
    	fn=files[m-1] 
    	PRINT,'Could not find file: ',files[m]
    	PRINT,'Will use previous file instead...'
    ENDIF ELSE BEGIN
    	PRINT,'OK, I can not find the previous file either. I give up.'
    	CD,orig
    	return
    ENDELSE
    print,fn
    IF (cor_flag) THEN BEGIN
    	im = sccreadfits(fn,hdr) 
	; im is type float
	IF hdr.naxis1 EQ 2048 THEN BEGIN
	    hdr.dstart1=1
	    hdr.dstop1 =2048
	    ; some early TBr images are incorrect
	ENDIF
    	im = SCC_PUTIN_ARRAY(im,hdr,1024,/new)
    	ss=''
    	k=0
    	WHILE ss NE 'EXPTIMEs' DO BEGIN
	    k=k+1
	    ss=strmid(hdr.comment[k],0,8)
	ENDWHILE
	print,hdr.comment[k]
	parts=strsplit(hdr.comment[k],' ,',/extract)
	recexp=avg(float(parts[1:3]))
    	effexp=1.0
	; TBr from SECCHI_P0 are effectively 1 second
    	dsatval=30000.	; approximate for 6sec exposure, from 1 image example (20071005_025230_0B4c2B.fts)
    ENDIF ELSE BEGIN
    	secchi_prep,fn, hdr,im, /NOCALFAC, /calimg_off, /update_hdr_off, /EXPTIME_OFF, OUTSIZE=1024, SILENT=quiet
    	; image is corrected for IP including IPSUM
	; Do not correct for EXPTIME in secchi_prep because of weighting done in hi_correction.pro
    	recexp=hdr.exptime
	;help,recexp
    ;wait,2
	im=im/recexp
    	effexp=1.0
	dsatval=hdr.dsatval
	IF dsatval LE 0 THEN BEGIN
	    IF strmid(pref,0,1) EQ 'h' THEN ccdsat=14000. ELSE ccdsat=15000.
	    dsatval = long(ccdsat * hdr.n_images * (2^(hdr.IPSUM - 1))^2)
	    IF strmid(pref,0,1) NE 'h' THEN dsatval = dsatval<60000
	ENDIF
    ENDELSE
    
    thistai=anytim2tai(hdr.date_obs)
    IF (m GT 0) THEN cadences[m-1]=   thistai-lasttai
    lasttai=thistai
    readtimes[m]=  hdr.readtime
    cleartimes[m]= hdr.cleartim
    pitches[m]=    hdr.crval1
    yaws[m]=	   hdr.crval2

    nmissings[m]=  hdr.nmissing
    cosmicss[m]=   hdr.cosmics
    
    crotas[m]=hdr.crota
    biases[m]=hdr.biasmean
    exptimes[m]=recexp
    offsetcrs[m]=hdr.offsetcr	; bias that has been subtracted
    nsats[m]=hdr.datasat
    
    datacube[*,*,m]=im
    ;--Datacube is normalized to 1 second
    ;  
   ; stop
    if (m EQ midpoint) THEN outhdr0 = hdr
    if (m EQ 0) THEN dateobstai=thistai
    if (m EQ 0) THEN date_obs=hdr.date_obs

ENDFOR

; This is experimental but what I'm trying to do here is remove from the datacube any images
; that would not make a good median image. So I start by removing images that are rolled more
; than 5-sigma of the median roll value for that day. Later I will remove pixels that are 
; freakishly bright compared to the others. Note that 5-sigma might not be tight enough... we'll see...

;crota_stdv=stddev(hdrs.crota)
;crota_med=median(hdrs.crota)
;good=where( abs(hdrs.crota - crota_med) LT (5*crota_stdv) )

<<<<<<< scc_mk_daily_med.pro
		IF datatype(summary) EQ 'INT' THEN BEGIN
			CD,orig
			return
		ENDIF
                
                ; Either an individual COR polarizer angle is picked, or a HI image (no polarizer) or total brightness is picked. 
                
                if (pol NE 'tbr') then good = where(summary.VALUE EQ pol and summary.XSIZE GE xaxis and summary.YSIZE EQ yaxis and summary.PROG NE 'Dark' and summary.PROG NE 'Doub') $
                    else if not cor_flag then good = where(summary.XSIZE GE xaxis and summary.YSIZE EQ yaxis and summary.DEST EQ 'SSR1'and summary.PROG NE 'Dark') $
                        else good = where(summary.VALUE EQ 1001 and summary.XSIZE GE xaxis and summary.YSIZE EQ yaxis and summary.PROG NE 'Dark')

                 ; Before we continue, just make sure we haven't screwed up and dropped too many images...
                if n_elements(good) LT n_median_min then begin
                    print,'ERROR: There are not enough images to make a good daily median image.'
		    CD,orig
                    return
                endif else files=summary[good].FILENAME  
                 
        ENDELSE
	
        n=n_elements(files)
        PRINT,'Found a total of ',strcompress(n),' appropriate files for the day.'          
        path=s+date+'/'
        CD,path,curr=old
        
       
       if (n gt 40) THEN BEGIN
            ; There are more files than we need so lets take the first and last 20 of the day.
            tmp1=files[0:19]
            tmp2=files[n-20:n-1]
            files=[tmp1,tmp2]
            n=n_elements(files)
            PRINT,'Trimming down number of files to ',strcompress(n)
       ENDIF
       
       datacube=lonarr(xaxis,yaxis,n)
       crotas=fltarr(n) 
       
       midpoint=n/2
       
       n_used=n ; used later
       
       ; This is a temporary bug fix for the COR total-B summary files
       last_char=strmid(files[0],strlen(files[0])-1,1)
       if (last_char) NE 's' then files=files+'s'
       ; end bug fix
        
       FOR m=0,n-1 DO BEGIN
            next:
            if file_exist(files[m]) THEN im=readfits(files[m],hdr) $
                ELSE BEGIN
                    ;maillist='karl.battams@nrl.navy.mil sungrazer@nrl.navy.mil'
	            ;spawn,'mailx -s SCC_MK_DAILY_MED: Missing image found '+maillist+' < '+files[m],/sh
                    PRINT,'Could not find file: ',files[m]
                    PRINT,'Will use previous file instead...'
                    if file_exist(files[m-1]) THEN im=readfits(files[m-1],hdr) $
                        ELSE BEGIN
                            PRINT,'OK, I can not find the previous file either. I give up.'
    	    	    	    CD,orig
                            return
                        ENDELSE
                ENDELSE
            crotas[m]=float(sxpar(hdr,'CROTA'))
            im = scc_img_trim(im, hdr)
            datacube[*,*,m]=im
           ; stop
            if (m EQ midpoint) THEN BEGIN
            ;stop
                if datatype(hdr) NE 'STC' THEN h=fitshdr2struct(hdr) else h = hdr
                mid_date=hdr.DATE_OBS
                mid_time=strmid(mid_date,11,12)
                exptime=abs(hdr.EXPTIME)
                filter=hdr.FILTER
                polar=hdr.P1COL
                p1col=hdr.P1COL
                p1row=hdr.P1ROW
                p2col=hdr.P2COL
                p2row=hdr.P2ROW
            ENDIF    
        ENDFOR
        
        ; This is experimental but what I'm trying to do here is remove from the datacube any images
        ; that would not make a good median image. So I start by removing images that are rolled more
        ; than 5-sigma of the median roll value for that day. Later I will remove pixels that are 
        ; freakishly bright compared to the others. Note that 5-sigma might not be tight enough... we'll see...
                
        ;crota_stdv=stddev(hdrs.crota)
        ;crota_med=median(hdrs.crota)
        ;good=where( abs(hdrs.crota - crota_med) LT (5*crota_stdv) )
        
        crota_stdv=stddev(crotas)
        crota_med=median(crotas)
        good=where( abs(crotas - crota_med) LT (5*crota_stdv) )
        
       ; stop        
        if (n_elements(good) NE n) THEN BEGIN
            if (n_elements(good) GT 0) THEN datacube=datacube[*,*,good] ELSE BEGIN
                PRINT,'ERROR! Something went wrong when trying to eliminate rolled images from the datacube...'
                stop
            ENDELSE
        ENDIF
       ; stop
       
        
        
        median_array=lonarr(xaxis,yaxis)
        
        for j=0,xaxis-1 DO BEGIN
            for k=0,yaxis-1 DO BEGIN
                
                pixel_row = datacube[j,k,*]
                ; Get rid of zero-value pixels
               ; stop
                w = where(pixel_row[0,0,*] GT 0)
                if (n_elements(w) GT 1 ) THEN good_pix = pixel_row[0,0,*] ELSE good_pix = pixel_row
                ; Now try to eliminate extreme-value pixels.
                ; Note that 10-sigma might be too high...
                ; pix_med = median(pixel_row)
                ; pix_sd = stddev(pixel_row)
                ; ok=where( abs(pixel_row - pix_med) LT (10*pix_sd) )
                
                median_array[j,k]=median(good_pix)
            
            endfor  
        endfor
=======
crota_stdv=stddev(crotas)
crota_med=median(crotas)
good=where( abs(crotas - crota_med) LT (5*crota_stdv) )
help,good
>>>>>>> 1.17
         
if (n_elements(good) NE n) THEN BEGIN
    if (n_elements(good) GT 0) THEN datacube=datacube[*,*,good] ELSE BEGIN
        message,'ERROR! Something went wrong with range of CROTA in selected images...'
    ENDELSE
ENDIF
; stop
exptime_stdv 	=stdev(exptimes[good],exptime_mean)
bias_stdv    	=stdev(biases[good],bias_mean)
offsetcr_stdv	=stdev(offsetcrs[good],offsetcr_mean)
crota_stdv   	=stdev(crotas[good],crota_mean)
nsat	    	=min(nsats)
help,exptime_mean
dsatval     	=float(dsatval/exptime_mean)	; normalize this also to 1 second

datendtai   	=anytim2tai(hdr.date_end)
date_avg    	=utc2str(tai2utc(dateobstai + (datendtai-dateobstai)/2))
cadence_stdv	=stdev(cadences[good],cadence_mean)
readtime_stdv	=stdev(readtimes[good],readtime_mean)
cleartime_stdv	=stdev(cleartimes[good], cleartime_mean)
crval1_stdv 	=stdev(pitches[good],	crval1_mean)
crval2_stdv 	=stdev(yaws[good], crval2_mean)
cosmics_stdv	=stdev(cosmicss[good], cosmics_mean)

outhdr0.nmissing=total(nmissings[good])

median_array=fltarr(1024,1024)

for j=0,1024-1 DO BEGIN
    for k=0,1024-1 DO BEGIN

        pixel_row = datacube[j,k,*]
        ; Get rid of zero-value pixels
       ; stop
        w = where(pixel_row[0,0,*] GT 0)
        if (n_elements(w) GT 1 ) THEN median_array[j,k]=median(pixel_row[0,0,w])
        ; Now try to eliminate extreme-value pixels.
        ; Note that 10-sigma might be too high...
        ; pix_med = median(pixel_row)
        ; pix_sd = stddev(pixel_row)
        ; ok=where( abs(pixel_row - pix_med) LT (10*pix_sd) )

        

    endfor  
endfor

; stop
window,xsize=512,ysize=512,retain=2

;if (xaxis GT 1024) THEN median_array=rebin(median_array,1024,1024)
maxmin,median_array
maxmedian=max(median_array)
IF dsatval LT 32000/exptime_mean THEN BEGIN
; COR
    bscale= 1./exptime_mean  
    median_array_out=    fix(round( ((median_array>0<dsatval)/bscale) ))
ENDIF ELSE BEGIN
; HI
    bscale=1.0
    median_array_out=median_array>0<dsatval
ENDELSE
maxmin,median_array_out
IF keyword_set(DOREBIN) THEN median_array_out=rebin(temporary(median_array_out),512,512)
tv,hist_equal(rebin(median_array,512,512))
;stop
xyouts,20,20,sc+' '+date_avg,charsiz=2, /dev

outhdr0=scc_update_hdr(median_array>0,outhdr0,verbose=verbose,satmax=dsatval)
outhdr1 = rem_tag2(outhdr0,'TIME_OBS')
outhdr = struct2fitshead(outhdr1,/allow_crota,/dateunderscore)


;IF nsat LE 0 THEN it=where(median_array_out GE dsatval, nsat)
;fxaddpar,outhdr,'DSATVAL',dsatval
;fxaddpar,outhdr,'DATASAT',nsatload
help,dsatval,nsat,bscale

get_utc,dte,/ecs
fxaddpar,outhdr,'DATE-OBS',date_obs,' Start of first exposure in N_IMAGES'
fxaddpar,outhdr,'DATE-END',hdr.date_end,' End of last exposure in N_IMAGES'
fxaddpar,outhdr,'DATE-AVG',date_avg,' Midpoint between OBS and END'
fxaddpar,outhdr,'CADENCE', median(cadences),' sec (median); stdev='+trim(cadence_stdv)
fxaddpar,outhdr,'READTIME',readtime_mean,' sec; stdev='+trim(readtime_stdv)
fxaddpar,outhdr,'CLEARTIM',cleartime_mean,' sec; stdev='+trim(cleartime_stdv)
fxaddpar,outhdr,'CRVAL1',crval1_mean,' arcsec; stdev='+trim(crval1_stdv)
fxaddpar,outhdr,'CRVAL2',crval2_mean,' arcsec; stdev='+trim(crval2_stdv)
fxaddpar,outhdr,'COSMICS',cosmics_mean,' stdev='+trim(cosmics_stdv)
fxaddpar,outhdr,'DATE',dte
fxaddpar,outhdr,'BUNIT','DN/sec'
fxaddpar,outhdr,'BSCALE',bscale
fxaddpar,outhdr,'HISTORY','Data values normalized to 1 sec'
fxaddpar,outhdr,'EXPTIME', 1.0,'effective exptime; actual avg='+trim(exptime_mean)+'; stdev='+trim(exptime_stdv)
fxaddpar,outhdr,'BIASMEAN',bias_mean,'=avg over N_IMAGES'
fxaddpar,outhdr,'BIASSDEV',bias_stdv,' Over N_IMAGES'
fxaddpar,outhdr,'OFFSETCR',offsetcr_mean,'=avg over N_IMAGES; stdev='+trim(offsetcr_stdv)
fxaddpar,outhdr,'CROTA',crota_mean,'deg; =avg over N_IMAGES; stdev='+trim(crota_stdv)
fxaddpar,outhdr,'N_IMAGES',n_used,' scc_mk_daily_med.pro'

;stop
sdir1 = GETENV('SECCHI_BKG')+'/'+strlowcase(sc)+'/daily_med/'
sdir2 = sdir1+strmid(date,0,6)
if not file_exist(sdir2) THEN BEGIN
    PRINT,'Directory ',sdir2, ' does not exist.'
    cmd='mkdir -p -m 775 '+sdir2
    PRINT,cmd
    spawn,cmd,/sh
ENDIF

IF keyword_set(SAVEDIR) THEN sdir = savedir ELSE sdir = sdir2

CD,sdir
fname0 = 'd'+pref+strupcase(sc)+'_'+pol_str+'_'+strmid(date,2,6)+'.fts'

PRINT,'Writing daily fits file: '+sdir+'/'+fname0
print,' using '+string(n_used)+' images'
fxaddpar,outhdr,'FILENAME',fname0
fxaddpar,outhdr,'HISTORY',version
fxaddpar,outhdr,'COMMENT','%SCC_MK_DAILY_MED: header values from middle image unless noted'

WRITEFITS,fname0,median_array_out,outhdr

;stop

CD,orig
;IF floor(r) NE last_r THEN GOTO, beginning
RETURN
END