;+
; NAME: mkdomeflat
; 	based on mktwiflat, but simpler
; PURPOSE:
; 	makes a dome flat by loading a bunch of images, summing them, 
; 	and then flattening. Optionally deals with polarization issues.
;
;
; INPUTS:
; KEYWORDS:
;   infile	name of a text file containing a list of files.
;   		if present, this overrides istart and iend
; 	fdir	file directory
; 	noflag	don't try to flag negative pixels with BADVAL.
; 	/polarimetry	handle polarimetry subregions
; 	/auto	run automatically without asking for user intervention or
; 			confirmation
; 	
; OUTPUTS:
;
; HISTORY:
; 
; 2002-04-11 18:35:45 started based on mktwiflat.pro
; 2003-02-24 		  polarimetry code added. Converted to use fitsloader.
;
;----------------------------------------
PRO mkdomeflat,nstart,nend,infile=infile,$
	fdir=fdir,outdir=outdir,noflag=noflag,$
	polarimetry=polarimetry,outname=outnam,images=images,header=header,$
	auto=auto

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

;----------------------------------------

BADVAL = -1e6
;goto, SPT


if not(keyword_set(images)) then begin
	; this loads the specified images and handles bias subtraction.
	fitsloader,transpose([[nstart],[nend]]),instr='ircal',fdir=fdir,caldir=caldir,$
            dome,hdrs,/dobias,imlist=imlist
	images = dome
	head0=hdrs[*,0]
	header=head0
endif else begin
	dome = temporary(images)
	head0=temporary(header)
endelse


if not(keyword_set(outnam)) then outnam = 'domeflat-'+strc(sxpar(head0,"FILTER1"))

;------------------------------------------------------------
sz=size(dome)
imgsize=sz[1]
lincoeff = fltarr(imgsize, imgsize) 
; history string
sxaddhist, 'MKDOMEFLAT.PRO: '+systime()+ userid(), head0

; 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, dome, dome_lin, lincoeff, satmap, /nan
  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, dome, dome_lin, fltarr(imgsize, imgsize), satmap, /nocorrect,/nan $
 else dome_lin = temporary(dome_nb)
  sxaddhist, '  LINCORR: no linearity correction', head0
endelse
dome =  0

SPT:

;deal with multiple polarization regions. We have to do this -after-
; the step above so we can sort the images first.
if (keyword_set(polarimetry)) then begin
	medarr,dome_lin,dome_med
	pol_getregions,dome_med,regionmask,badregion,leftregion,rightregion, $
		badpix=badmask,middleregion=middleregion,/parallel
endif


; initialize output arrays
;dome_flat = fltarr(sz(1), sz(2))
;dome_const = fltarr(sz(1), sz(2))
;dome_masks = fltarr(sz(1), sz(2), sz(3))+1B

;==================================================

; mean? Or median?
;dome_flat = total(dome_lin,3)
medarr,dome_lin, dome_flat
    
;==================================================


; normalize and flag bad pixels (<0.001)
dome_flatraw = dome_flat
if keyword_set(polarimetry) then begin
	dome_flat[leftregion] = normal(dome_flat[leftregion])
	dome_flat[rightregion] = normal(dome_flat[rightregion])
	dome_flat[middleregion] = normal(dome_flat[middleregion])
	edgeregion = cmset_op(badregion, "and", /NOT2, middleregion)
	;dome_flat[edgeregion] = normal(dome_flat[edgeregion])
	dome_flat[edgeregion] = !values.f_nan ; these are just bad pixels

endif else begin ; regular mode
	dome_flat = normal(dome_flat)
endelse
w = where(dome_flat le 0.01, nw)
if (nw gt 0) then begin
    dome_flat(w) = !values.f_nan
    sxaddhist, '  setting <=0.01 pixels to NaN, # = '+strc(nw), head0
    print, '  setting <=0.01 pixels to NaN, # = '+strc(nw)
endif


badmask = bytarr(sz[1],sz[2])+1B
nanpixels = where(finite(dome_flat,/NaN),nanct)
if nanct gt 0 then badmask[nanpixels]=0B

; display results
loadct, 15
win,0
display2,regionmask,title="Polarimetry Regions"
win, 1
display2, dome_flat, tit = 'final domeflat'
win,2
display2, badmask,title= "bad pixel mask (1=good,0=bad)"
;win, 1
;display2, dome_const, tit = 'constant emission'
;iterstat, dome_const, istat, /silent
;print
;print, '----------------------------------------'
;print, '   output statistics'
;print, '----------------------------------------'
;print, 'constant frame median = ', strc(median(dome_const))
;print, 'constant frame (iter)sigma = ', strc(istat(3))
;print, 'constant frame sigma = ',  $
  ;strc(stdev(dome_const(where(dome_const ne BADVAL))))
;print
;stack, dome_masks, masktot
;win, 2
;display2, /tv, masktot eq 1, tit = 'total of fitting masks'



; 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

	i0 = nstart[0]
    outfile = outdir+outnam+"."+strc(i0)+'.fits'
    objname = sxpar(head0, 'OBJECT')

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

    outfile = outdir+outnam+"."+strc(i0)+'_badpix.fits'
    sxaddpar, head0, 'OBJECT',"BAD PIXEL MASK"
    if filecheck(outfile) then  $
      writefits, outfile, badmask, head0

	if keyword_set(polarimetry) then begin
	    outfile = outdir+outnam+"."+strc(i0)+'_regions.fits'
	    sxaddpar, head0, 'OBJECT',"POLARIMETRY REGION MASK"
	    if filecheck(outfile) then  $
	      writefits, outfile, regionmask, head0
	endif
;
;    outfile = outdir+outnam+strc(i0)+'_const.fits'
;    sxaddpar, head0, 'OBJECT', objname+' - CONSTANT'
;    if filecheck(outfile) then  $
;      writefits, outfile, dome_const, head0
;

endif

print
print, '** MKDOMEFLAT.PRO done **'
print

end