;+
; NAME: mktwiflat
;PURPOSE:
;  Median-combine a number of sky flat images into one master sky flat.
;
; INPUTS:
; KEYWORDS:
; 	datadir	file directory
; 	
; OUTPUTS:
;
; HISTORY:
; 
; 2002-04-29 14:40:55 split off from mktwiflat.pro by MDP
; 2004-05-12  Merged in some updates from mktwiflat.pro
; 2005-12-14  Don't set bad pixels in the flat to NaN; just
; 				mark the appropriate pixel in the bad pixel mask.
; 2005-12-22	Split off as 'newskyflat' for development branch.
;----------------------------------------
PRO newskyflat,istart,iend,$
	datadir=datadir,outdir=outdir,noflag=noflag,$
	images=images,header=header,auto=auto,$
	caldir=caldir,instrument=instrument,hdrs=hdrs,flatout=flatout

if not(keyword_set(outdir)) then outdir = './'
if not(keyword_set(datadir)) then datadir = './data'
checkdir,outdir
checkdir,datadir

if ~(keyword_set(intstrument)) then instrument="ircal"

;------------------------------------------------------------
i0 = ISTART
i1 = IEND

	; get images, either from disk or from supplied variables
	print, '=== loading images ==='
	if not keyword_set(images) then begin
		fitsloader,transpose([[i0],[i1]]),instrument=instrument,datadir=datadir,caldir=caldir,$
			sky,hdrs,/dobias,imlist=imlist
		head0 = hdrs[*,0]
	endif else begin
		sky = images
		head0= hdrs[*,0]
		sxaddhist,"MKSKYFLAT: data passed via images/hdrs arguments.",head0
	endelse
;------------------------------------------------------------
sz=size(sky)
lincoeff = fltarr(sz[1],sz[2]) 
sxaddhist, 'MKSKYFLAT.IDL: '+systime()+ userid(), head0

if ~(keyword_set(auto)) then begin
	atv,sky,/bl,names=imlist
	stop
endif

; select a subset of images to use
;
; TODO: automatically detect saturated images here.
; 		Also throw out images which have too few counts to be useful
; 		(or will those be effectively thrown out by the error masking?)
nz = sz[3]
print, '=== image statistics (bias-subtracted) ==='
stats, sky
print, 'there are ', strc(nz), ' twi frames'
sxaddhist, 'there are '+ strc(nz)+ ' twi frames',head0
while getyn('use only a subset of the images?', 0,noquery=auto) do begin
    sz = size(sky)
    nn = sz(3)
    a0 = -1  &  a1 = -1
    repeat begin
        read, 'enter starting image: (0-',+strc(nn-1)+') ', a0
        read, 'enter ending image: (0-',+strc(nn-1)+') ', a1
    endrep until ((a0 ge 0) and (a0 lt nn) and $
                  (a1 ge 0) and (a1 lt nn) and (a1 gt a0)) 
    sky = sky[*, *, a0:a1]
    print, '=== examine image stats (w/o bias) ==='
    stats, sky
    i0 = i0 + a0
    i1 = i0 + (a1-a0)
    sxaddhist, '  selected images '+strc(i0)+' to '+strc(i1), head0
    nz = a1-a0+1
    print, 'there are ', strc(nz), ' twi frames'
endwhile

; 2005-12-15
;  Don't use images with more than ~18,000 counts per coadd median.
;  These may have saturated pixels toward the right hand side of the array
;  which would mess up the flat computation.
median_thresh=18000
message,/info,"Rejecting frames with median above "+strc(median_thresh)
sxaddhist,"Rejecting frames with median above "+strc(median_thresh),head0
med_per_coadd = meds(sky)  / sxpararr(hdrs,"NCOADDS")
wgood = where(med_per_coadd le median_thresh,ngood,ncomp=nrej)
if nrej gt 0 then begin
	sky = sky[*,*,wgood]
	message,/info,strc(nrej)+ " frames rejected."
	sxaddhist,strc(nrej)+ " frames rejected.",head0
	nz =ngood
	sz = size(sky)
endif
stop

; apply linearity correction
; (use /nocorrect flag if no lin coeffs - will still check for saturation)
if getyn('apply linearity correction?', (avg(lincoeff) ne 0),noquery=auto) then begin
  lincorr, sky, sky_lin, lincoeff, satmap 
  if total(lincoeff) ne 0 then $
    sxaddhist, '  LINCORR: applied linearity correction', head0 $
  else $
    sxaddhist, '  LINCORR: no linearity correction', head0
endif else begin
 if not(keyword_set(noflag)) then $
	  lincorr, sky, sky_lin, fltarr(sz[1],sz[2]), satmap, /nocorrect $
 else sky_lin = temporary(sky)
  sxaddhist, '  LINCORR: no linearity correction', head0
endelse


; convention is 1 = good, 0 = bad
badpixmask = bytarr(sz[1],sz[2])+1B

; bad pixel mask based on what LINCORR flags as bad.
print,"Marking bad pixels base on LINCORR results"
sky_bad = bytarr(sz[1],sz[2],sz[3])+1
wbad = where(sky_lin eq BADVAL,badcount)
if badcount gt 0 then sky_bad[wbad]=0
badmask = total(sky_bad,3)<1
badpixmask *= badmask

;;  print,"Now normalizing individual sky flats."
;;  sky_norm = fltarr(sz[1],sz[2],sz[3])
;;  for x=0,sz[3]-1 do begin
;;  	sky_norm[*,*,x] = normal(sky_lin[*,*,x])
;;  endfor 
;;  
;;  print,"Now median combining all flats."
;;  medarr,sky_norm,sky_flat,bytarr(sz[1],sz[2],sz[3])+1B
;;


; 2005-12-14: Don't mark bad pixels in the flat, just in the badpixmask

; lincorr flags pixels which are negative as NaNs.
; but now we don't want NaNs in the flat, just the badmask.
wnan = where(sky_flat ne sky_flat,nanct)
if nanct gt 0 then begin
    sxaddhist, '  setting negative pixels to bad, # = '+strc(wnan), head0
    print, '  setting negative pixels to bad, # = '+strc(wnan)
	sky_flat[wnan]=1
	badpixmask[wnan]=0
endif

wbad = where(sky_bad eq 0,badcount)
;if badcount gt 0 then sky_flat[wbad]=!values.f_nan
if badcount gt 0 then badpixmask[wbad] = 0
	
; normalize and flag bad pixels (<0.001)
w = where(sky_flat le 0.001, nw)
if (nw gt 0) then begin
    ;sky_flat[w] = !values.f_nan
    sxaddhist, '  setting <=0.001 pixels to bad, # = '+strc(nw), head0
    print,'  setting <=0.001 pixels to bad, # = '+strc(nw)
	badpixmask[w] = 1
endif

;----------------------------------------
; guess a bad pixel mask
; IRCAL traditionally has about 1.5% of the pixels low due to the two blocked
; corners and other misc defects
;
; Let's throw away everything more than 3 standard deviations down from the mean
; 
; grab cold pixels based on slopes
;
print, "Now marking bad pixels "
;binsize = 0.001
;h = histogram(sky_flat > 0,bin=binsize,omin=omin,omax=omax)
;hi = total(h,/cumulative)
;cutoff = binsize*(value_locate(hi,0.02*n_elements(sky_flat))+1)
wg = where(sky_flat ne BADVAL and not finite(sky_flat,/NaN))
cutoff = mean(sky_flat[wg]) - 3*stddev(sky_flat[wg])

if instrument eq "NIRC2" then cutoff = 0.6

sz = size(sky_flat)
cold = where(sky_flat le cutoff,wcold)
if wcold gt 0 then badpixmask[cold] = 0B

wnan = where(finite(sky_flat,/NAN) eq 1,nanct)
if nanct gt 0 then badpixmask[wnan] = 0B

;wbad = where(badpixmask eq 0,badct)
;if badct gt 0 then sky_flat[wbad] = !values.f_nan

print, "Cold pixel cutoff: slope lt "+strc(cutoff)+",	 "+strc(wcold)+" pixels rejected."

; display results
loadct, 15
win, 0
display2, sky_flat, tit = 'final twiflat'
stat,sky_flat
win,1
display2,badpixmask,title="Bad Pixel Mask (dark = bad)"


if not(keyword_set(auto)) then stop

	outnam = 'skyflat-'+strc(sxpar(head0,"FILTER1"))

; save images if desired
if getyn('save the flat?', 1,noquery=auto) then begin
	
	if keyword_set(aobir) then begin
		; this is needed to change the FITS header since the
		; input stuff from CFHT has BZERO=32768 but our output here doesn't.
		print,"Updating BZERO keyword in output file"
		sxaddpar,head0,"BZERO",0,"offset"
	endif

    if getyn('outdir="'+outdir+'", outnam="'+outnam+'": '+ $
             'do you want to change these?', 0,noquery=auto) then begin
        read, 'enter outdir: ', outdir
		checkdir,outdir
        read, 'enter outnam: ', outnam
    endif

	; just use first image to name flat.
 	if n_elements(i0) gt 1 then i0 = i0[0]; 
	
    outfile = outdir+outnam+"."+strc(i0)+'.fits'
    objname = sxpar(head0, 'OBJECT')

    if filecheck(outfile) then  $
      writefits, outfile, sky_flat, head0

    outfile = outdir+outnam+"."+strc(i0)+"_badmask.fits"
    sxaddpar, head0, 'OBJECT', objname+' - BAD PIXELS'
    if filecheck(outfile) then  $
      writefits, outfile, badpixmask, head0


endif

flatout = sky_flat
print
print, '** MKSKYFLAT.IDL done **'
print

end