; script for making twilight flat: note that flats are made with
; *arithmatic* averaging, since routine assumes you are able to mask
; out all objects (and it ignores cosmic rays: this is ok provided the
; bad pixel mask is made by clipping high pixels in the flat, though
; we do lose a few pixels this way)
; 2/26/96 MCL
;
; added second pass to check for stars in frames
; and also subtracts a frame to account for telescope emission
; (has everything in 'twi_2pass.idl' and more)
; output is three flats (all have obj-masking, med-wt., arith.-averaging)
;   1. avg. of all the frames ('twi_flat_all')
;   2. avg. of "telescope-subtracted" subset of the frames ('twi_flat')
;   3. avg. of the subset in (2), w/o tel.-subtraction ('twi_flat_nosub')
; 'twi_flat' should be the 'best' one to use.
; extensive checking -> this thing works fine.
; 8/18/96 MCL
;
; corrected bug in making twiflat w/o telescope-subtraction
; from the subset of frames - wasn't using the right set
; of images with the masks
; 11/08/96 MCL
;
; bug fix: forgot when doing telescope subtraction to account for
;   subtracted frame's mask - so we have in the past been
;   adding negative objects (if there were any in the frame)
; now pass badpix to AVGMASK to ignore the bad pixels when
;  calculating median weights (now default for AVGMASK)
; 11/09/96 MCL
;
; bug from last time had *not* been fixed!!!
; now it finally does account for obj masks in the frame used
; to subtract telescope emission
; 04/09/98 MCL
;
; creation of MKTWIFLAT.IDL routine makes this obsolete,
; since the former uses a new & better algorithm
; 03/10/99 MCL

;----------------------------------------
; set these to the appropriate values:
ISTART =  145
NF =  6
NC =  1
bias =  readfits('calib-15apr/bias170b_2s.fits')
fdir = 'fits-15apr98/'
fnam = '15apb'
outdir = 'calib-15apr/'
outnam = 'kptflat'
badpix = badpix
CBOX = 3
gemini = 1
;----------------------------------------


BADVAL =  -1e6
MASKRAD = 7       ; for object masking
dlim =  100        ; for sensible displaying
if (!d.window lt 0) then begin
    win, 0
endif

;goto, SPT

; get images
print, '=== loading images ==='
loadims, ISTART, ISTART+NC*(NF-1), step = NC, twi, $
  fdir = fdir, fnam = fnam, gemini = gemini
;loadims, ISTART, ISTART+NC*(NF-1), step = NC, twi, fdir = fdir, /lick

; subtract bias
print, '=== subtracting bias ==='
skysub, twi, bias, twi_nb, /silent
twi =  0

; check counts
print, '=== examine image stats (w/o bias) ==='
stats, twi_nb
print

; apply linearity correction
; (use /nocorrect flag if no lin coeffs - will still check for saturation)
if getyn('apply linearity correction?') then $
  lincorr, twi_nb, twi_lin, lincoeff, satmap $
else $
  lincorr, twi_nb, twi_lin, fltarr(256, 256), satmap, /nocorrect
twi_nb =  0



; median average all frames with median-weighting
print, '=== generating average sky frame (be patient!) ==='
avgmask, twi_avg, twi_lin, twi_lin gt 0, twi_map, /silent, /median, badpix=badpix


; check the histogram & map of the contributing pixels
print, '=== check median averaging map & histogram ==='
scan, twi_map, /tvscl
plothist, twi_map, bin = 1, $
  title = 'contributions to average frame', $
  xtitle = 'frame number',  ytitle = 'number of pixels'


; do preliminary sky subtraction
print, '======================================'
print, '=== doing 1st pass sky subtraction ==='
print, '======================================'
if getyn('use AUTOGAL (NSKY=2) for 1st pass sky subtraction?') then begin
    autogal, twi_lin, twi_sub,  $
      masks = temporary(twi_lin ne BADVAL), $
      nsky = 2, /median, badpix = badpix
endif else begin
    ;subtract med-scaled average frame from every frame to check 
    ;to check for objects 
    skysub, twi_lin, twi_avg, twi_sub, /scale, /silent, $
      badpix =  twi_avg gt BADVAL
    w =  where(twi_lin eq badval, nw)
    if (nw gt 0) then twi_sub(w) =  BADVAL
endelse
print, '=== examine first pass sky-subtracted frames for objects ==='
scan, twi_sub > (-dlim) < dlim
twi_avg =  0 


; mask stars from the images
; and combine mask with the saturated pixel mask
print, '========================================='
print, '=== please mask stars from sky frames ==='
print, '========================================='
ans = 1
if n_elements(twi_masks) ne 0 then begin
    ans =  getyn('masks already defined - reinitialize them?', 0)
    if not(ans) then $
      ans2 = getyn('masks already defined - add objects to them?')
endif 
if (ans) then begin
    imask, twi_sub, twi_masks, $
      badpix = badpix, radius = MASKRAD, /new, $
      cbox = CBOX
;      lscl = -dlim, uscl = dlim
endif else if (ans2) then begin
;endif else if (ans2 eq 'y') then begin
    imask, twi_sub, twi_masks, $
      badpix = badpix, radius = MASKRAD, $
      cbox = CBOX
;      lscl = -dlim, uscl = dlim
endif



; review masked images and do more masking if desired
print, '================================='
print, '=== look at masked sky images ==='
print, '================================='
scan, twi_sub * twi_masks > (-dlim) < dlim
while getyn('do you want to do more masking?', 0) do begin
    imask, twi_sub, twi_masks, $ 
      radius = MASKRAD, badpix = badpix
;      lscl = -dlim/2., uscl = dlim/2.
    print, '================================='
    print, '=== look at masked sky images ==='
    print, '================================='
    scan, twi_sub * twi_masks > (-dlim) < dlim
endwhile

SPT:

; arithmatically average frames w/median-weighting, 
; incorporating masks to exclude the objects
print, '=== making twiflat using object masks ==='
avgmask, twi_avg, twi_lin, (twi_lin gt 0) * twi_masks, badpix=badpix

; normalize to make the flat
w =  where(twi_avg eq BADVAL)
twi_flat_all =  normal(twi_avg)
twi_flat_all(w) =  BADVAL
twi_avg = 0




print, '=================================================='
print, '      Telescope-subtraction process '
print, '=================================================='

; determine which frame to subtract
print, '=== examine bias-subtracted image stats ==='
stats, twi_lin
sz =  size(twi_lin)
ii =  sz(3)-1
repeat $
  read, 'which frame to subtract? (0=first, 1=last) ', ns $
until (ns eq 0) or (ns eq 1)

; do telescope subtraction
; bad pixels in the subtracted frame get passed to all other frames
print, 'subtracting from frames ', strc(indgen(ii)+(ns eq 0))
skysub, twi_lin(*, *, indgen(ii)+(ns eq 0)), twi_lin(*, *, (ns eq 1)*ii), $
  twi_tsub, /silent 
ww =  (twi_lin(*, *, (ns eq 1)*ii) eq BADVAL) or $
  (twi_masks(*, *, (ns eq 1)*ii) eq 0)
for i = 0, ii-1 do $
    twi_tsub(*, *, i) =  twi_tsub(*, *, i) - (ww * twi_tsub(*, *, i))
    

; select which frames to use for twiflat
print, '=== examine telescope-subtracted image stats ==='
stats, twi_tsub
repeat begin
    read, 'enter starting frame number: (0-',+strc(ii-1)+') ', n1
    read, 'enter ending frame number: (0-',+strc(ii-1)+') ', n2
endrep until (n1 ge 0) and (n2 ge 0) and (n1 le ii-1) and (n2 le ii-1) and $
    (n1 lt n2) 


; make a flat from the same subset except w/o
; telescope-subtraction so can compare the difference between 2 flats
print, '=== making NON-telescope-subtracted flat, using object masks ==='
avgmask, twi_nsavg, twi_lin(*, *, indgen(n2-n1+1)+n1+(ns eq 0)), $
  (twi_lin(*, *, n1:n2) gt 0) * $
  (twi_masks(*, *, indgen(n2-n1+1)+n1+(ns eq 0))), $
;  badpix = badpix
  badpix = badpix * twi_masks(*,*,ii*(ns eq 0))


; do initial med-wt. median averaging to see if
; telescope-subtraction helps by studying the median histogram & map
print, '=== median averaging telescope-subtracted images, using object ' + $
  'masks === ' 
avgmask, twi_tsavg, twi_tsub(*, *, n1:n2), $
  (twi_tsub(*, *, n1:n2) gt 0) * $
  (twi_masks(*, *, indgen(n2-n1+1)+n1+(ns eq 0))), twi_tsmap, /median, $
;  badpix= badpix
  badpix= badpix * twi_masks(*,*,ii*(ns eq 0))
print, '=== check median averaging map & histogram ==='
scan, twi_tsmap, /tvscl
plothist, twi_tsmap, bin = 1, $
  title = 'contributions to average frame', $
  xtitle = 'frame number',  ytitle = 'number of pixels'


; arithmatically average frames w/median-weighting, 
; incorporating masks to exclude the objects
print, '=== making final telescope-subtracted flat, using object masks ==='
avgmask, twi_tsavg, twi_tsub(*, *, n1:n2), $
  (twi_tsub(*, *, n1:n2) gt 0) * $
  (twi_masks(*, *, indgen(n2-n1+1)+n1+(ns eq 0))), $
;  badpix= badpix 
  badpix= badpix * twi_masks(*,*,(n2-n1+1)*(ns eq 0))

    
; normalize to make the flat
w =  where(twi_tsavg eq BADVAL)
twi_flat =  normal(twi_tsavg)
twi_flat(w) =  BADVAL


; normalize to make the non tel.-subtracted flat
wbad =  where(twi_nsavg eq BADVAL, nbad)
twi_flat_nosub =  normal(twi_nsavg)
if (nbad gt 0) then twi_flat_nosub(w) =  BADVAL



; save images if desired
if getyn('save the flat and masks?') then begin

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

    filen = outdir+outnam+strc(ISTART)
    if (gemini) then $
        filen = filen + strlowcase(strc(sxpar(head, 'CHANNEL'))) 

;    sxaddhist, history, head
    mkhdr, head, twi_flat
    spawn, 'date', date
    sxaddhist, 'TWI_TSUB.IDL: '+date, head
    writefits, filen+'-all.fits', twi_flat_all, head
    writefits, filen+'-nosub.fits', twi_flat_nosub, head
    writefits, filen+'-tsub.fits', twi_flat, head
;    sxaddpar, head, 'OBJECT', strupcase(filen+' sky masks')
    writefits, filen+'-masks.fits', twi_masks, head

endif

end