; edited version of PTSOURCE.IDL to re-reduce images 
; using existing object masks and measured offsets to save effort,
; typically used to emply diff. biases and/or flats in the reduction
;
; 02/25/98 MCL
;
; option to avoid re-saving image offsets file, 
; avoiding erasing of the original comments
; 03/31/98 MCL


;--------------------------------------------------
; please set these to the appropriate values:
;
fdir = 'fits-11nov97/'    ; image directory
fnam = 'ctio'       ; image name prefix
objn = 's301d'  ; object name string
setn = '_n4-1'            ; set name string (user must enter "_" if desired)
ISTART = 4281          ; 1st image
NF =      10      ; number of images
NC =      1        ; step size for image numbers
biasf = 'calib/bias4231_8s.fits'       ; bias frame
flatf = 'calib/kstflat4530-nosub.fits'          ; flat field
badpix = badpix    ; bad pixel mask (0=bad, !0=good)
sflag = 0          ; 0=don't examine images, !0=examine images
;lincoeff = fltarr(256, 256)
gemini = 0         ; if Gemini files
lick = 0           ; if Lick files 
ctio = 1
CBOX = 5        ; box size for centroiding

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


BADVAL = -1e6
NSKY = 3        ; value for AUTOGAL sky subtraction
MASKRAD = 20    ; radius for masks
MASKSIG = 0.7   ; sigmascale for image masking
dlim = 50       ; for sensible displaying
filen = objn+setn    ; string for files names

shiftfile = 'standards/'+filen+'_shifts.dat'
maskfile = 'standards/'+filen+'_masks.fits'

;goto, G_WRITE
;G_START:


;; check for calibration frames
;if (n_elements(biasf) eq 0) or (n_elements(badpix) eq 0) or $
;  (n_elements(flatf) eq 0) then begin 
;    print, '** calibration frames not loaded !! **' 
;    retall
;endif


; load calibration frames
print, '=== loading bias frame: "', biasf, '" ==='
bias =  readfits(biasf, /silent)
if (n_elements(bias) eq 1) then begin
    print, '** bias file does not exist **'
    stop
endif
print, '=== loading flat field: "', flatf, '" ==='
flat =  readfits(flatf, /silent)
;flat = fltarr(256, 256)+1.0
if stdev(flat) eq 0 then $
  if not(getyn('using dummy flat field (all pix='+strc(avg(flat))+'). continue?')) then stop
if (n_elements(flat) le 0) then begin
    print, '** flatfield file does not exist **'
    stop
endif
if (n_elements(badpix) eq 0) then begin
    print, '** bad pixel array does not exist !! **' 
    retall
endif
;if (n_elements(satmap) eq 0) then begin
;    print, '**  saturation map does not exist !! **' 
;    retall
;endif


; check file numbers
print, '-- ISTART = ', strc(ISTART), ',  NCOADD = ', strc(NC), ', ' + $
  'NFILES = ', strc(NF), ' --'



; make image list
listims, istart, istart+nc*(nf-1), step = nc, $
  out = filen+'-in.lst', fname = fdir+fnam, $
  gemini = gemini, lick = lick, ctio = ctio, $
  /silent


; get images
loadimf, filen+'-in.lst', obj, head0, gemini = gemini, lick = lick, ctio = ctio


; add history info
spawn, 'date', date
spawn, 'pwd', pwd
spawn, 'echo $USER@$HOST.berkeley.edu', id
sxaddhist, 'REPROC.IDL: re-running pt source processing', head0
sxaddhist, 'PTSOURCE.IDL: '+date+', '+id, head0
sxaddhist, '  starting directory = '+pwd, head0
sxaddhist, '  LOADIMF: loaded images '+strc(istart)+$
  ' to '+strc(istart+nc*(nf-1))+', step='+strc(nc), head0
sxaddhist, '    divided by number of coadds', head0
sxaddhist, '    image directory = "'+fdir+'"', head0
sxaddhist, '    image prefix = "'+fnam+'"', head0


; subtract bias
print, '=============================='
print, '=== subtracting bias frame ==='
print, '=============================='
skysub,obj,bias,obj_nb
sxaddhist, '  SKYSUB: subtracted bias frame "'+biasf+'"', head0
obj = 0


; linearize
print, '==================================================='
print, '=== linearity correction and bad pixel flagging ==='
print, '==================================================='
if getyn('apply linearity correction?', 0) then begin
  lincorr, obj_nb, obj_lin, lincoeff, satmap
  if total(lincoeff) ne 0 then $
    sxaddhist, '  LINCORR: applied linearity correction', head0 $
  else $
    sxaddhist, '  LINCORR: no lin correction, flagged pix<0 as BADVAL', head0
endif else begin
  lincorr, obj_nb, obj_lin, lincoeff, /nocorrect 
  sxaddhist, '  LINCORR: no lin correction, flagged pix<0 as BADVAL', head0
endelse
obj_nb =  0



; examine bias-subtracted images
if (sflag ne 0) then begin
    loadct, 15
    win,0,title='IDL (0): scan mode'
    print, '=============================================================='
    print, '=== look at bias-sub., linearized images (self-subtracted) ==='
    print, '=============================================================='
    scan,obj_lin>0,/sub, const = badpix
endif


;; make sky flat
;print, '========================'
;print, '=== making sky flat  ==='
;print, '========================'
;avgmed, obj_skyflat, obj_lin
;w = where(obj_skyflat eq BADVAL, nbad)
;obj_skyflat = normal(obj_skyflat)
;if (nbad gt 0) then obj_skyflat(w) = BADVAL
;display2, obj_skyflat



SPT:

; flatten
print, '========================='
print, '=== flattening frames ==='
print, '========================='
w =  where(obj_lin le BADVAL)
flatdiv, obj_lin, flat, obj_linfl, badpix = (flat gt 0)
if stdev(flat) eq 0 then $
  sxaddhist, '  no flat field applied', head0 $
else $
  sxaddhist, '  FLATDIV: flat-fielded with "'+flatf+'"', head0
obj_linfl(w) =  BADVAL
obj_lin =  0


; examine the results
if (sflag ne 0) then begin
    print, '=================================='
    print, '=== look at flatfielded images ==='
    print, '=================================='
    scan, obj_linfl > 0, const = badpix
endif



read2, shiftfile, dat
obj_dx0 = dat(*, 2)
obj_dy0 = dat(*, 3)
if (n_elements(obj_dx0) ne NF) then begin
    print, 'shifts file does not agree with number of images!'
    stop
endif


obj_masks = readfits(maskfile)
if ((size(obj_masks))(3) ne NF) then begin
    print, 'obj masks do not agree with number of images!'
    stop
endif



; do 2nd pass sky subtraction, using masks
; note that pixels in the subtracted images which had
; no good pixels in the sky frame are set to BADVAL
;
; NOTE: uses ALL frames, even ones flagged during IREG as bad!
print, '================================================'
print, '=== doing 2nd pass sky subtraction (w/masks) ==='
print, '================================================'
if getyn('for sky subtraction, NSKY='+strc(nsky)+', change this?', 0) then $
  read, 'enter new NSKY = ', nsky
autogal, obj_linfl, obj_sub, obj_sky,  $
  nsky = NSKY, mask = obj_masks * (obj_linfl gt BADVAL), $
  badpix = badpix
sxaddhist, '  AUTOGAL: sky subtraction with nsky='+strc(NSKY), head0
sxaddhist, '    subtraction scaled by image medians', head0


;; cross correlate & mosaic
;print, '==============================================='
;print, '=== cross correlating images to find shifts ===
;print, '==============================================='
;imgcor, obj_sub(*, *, objlist), obj_dx, obj_dy, mask = $
;  (obj_masks* (obj_sub gt BADVAL))(*, *, objlist),  badpix = badpix



G_RECEN:

; register frames interactively and mosaic
; IREG allows user to flag frames to be discarded - 
;   these have offsets set to -999
print, '=========================='
print, '=== registering frames ===
print, '=========================='
wset, 0
if getyn('interactively re-register and exclude frames?',0) then begin
    ireg, obj_sub, obj_dx, obj_dy, $
      badpix = badpix, masks = (obj_sub gt BADVAL),  $
      /nomos, lscl = -dlim, uscl = dlim
    sxaddhist, '  IREG: centroiding to find offsets on 2nd pass mosaic', head0
endif else begin
    recenter, obj_sub, obj_dx0, obj_dy0, obj_dx, obj_dy, $
      badpix = badpix, mask = (obj_sub ne BADVAL), cbox = CBOX
    sxaddhist, $
      '  RECENTER: auto-register using 1st pass offsets as starting guesses', $
      head0
    if getyn('do you want to exclude any frames?', 0) then begin
        repeat begin
            wgood = where(obj_dx ne -999)
            wbad = where(obj_dx eq -999)
            print, 'good images = ', commalist(wgood)
            print, 'bad images = ', commalist(wbad)
            read, 'enter number of image to exclude (<0 to quit): ', aa
            w = where(wgood eq aa)
            if (aa ge 0) and (w(0) ne -1) then begin 
                obj_dx(aa) = -999
                obj_dy(aa) = -999
            endif
        endrep until (aa lt 0)
    endif
endelse


G_MOS:

; make mosaic
print, '======================'
print, '=== making mosaic ==='
print, '======================'
wgood = where(obj_dx ne -999, ngood)
nimgs = n_elements(obj_dx)
print, '=== using '+strc(ngood)+' out of '+strc(nimgs)+' images ==='
sxaddhist, '    '+strc(ngood)+' out of '+strc(nimgs)+' are good', head0
if (ngood eq nimgs) then begin
    mosf, obj_sub, obj_dx, obj_dy, $
      obj_mos, obj_mosexp, obj_mossky,  $
      badpix = badpix, $
      mask = (obj_sub gt BADVAL),/setsky
endif else begin
    mosf, obj_sub(*, *, wgood), obj_dx(wgood), obj_dy(wgood), $
      obj_mos, obj_mosexp, obj_mossky,  $
      badpix = badpix, $
      mask = (obj_sub(*, *, wgood) gt BADVAL),/setsky
endelse
sxaddhist, '  MOSF: using /setsky', head0



; show resulting mosaic
window,4,xs=450,ys=450
loadct,15
wset,4
display2,obj_mos


G_WRITE:



;--------------------------------------------------
; write results from reduction scripts 
;--------------------------------------------------
if not(getyn('do you want to save the files?')) then stop
if getyn('obj="'+objn+'", set="'+setn+'": do you want to change ' + $
         'these?', 0) then begin
    read, 'enter object name: ', objn
    read, 'enter set name: ', setn
;   filen = objn+'_'+setn
   filen = objn+setn
endif
;if not(getyn('did you remember to change ISTART = '+strc(istart)+'?')) $
;  then read, 'enter value for ISTART: ', istart  
writeims, obj_sub, head = head0, numlist = indgen(NF)*NC + istart, $
  objn = strupcase(filen), fname = objn+'_', ext = '.fits', /lick
listims, istart, istart+nc*(nf-1), step = nc, $
  out = filen+'-out.lst', fname = objn+'_', ext = '.fits', /lick, /silent
if getyn('save the individual sky frames?', 0) then $
  writeims, obj_sky, head = head0, numlist = indgen(NF)*NC + istart, $
  objn = strupcase(filen+' sky frames'), $
  fname = objn+'_', ext = 'sky.fits', /lick


; write mosaic files
print, '=== writing mosaic files for <', filen, '> ==='
sxaddpar, head0, 'OBJECT', strupcase(filen+' mosaic   ')
writefits, filen+'_mos.fits', obj_mos, head0
sxaddpar, head0, 'OBJECT', strupcase(filen+' sky masks')
writefits, filen+'_masks.fits', obj_masks, head0
sxaddpar, head0, 'OBJECT', strupcase(filen+' mosaic exposure map')
writefits, filen+'_mosexp.fits', obj_mosexp, head0
sxaddpar, head0, 'OBJECT', strupcase(filen+' mosaic object masks')


; write image offsets file
spawn,'date',date
if getyn('re-save image offsets file? (will erase original comments)', 0) then begin
    openw, unit0, filen+'_shifts.dat', /get_lun
    printf,unit0, '# ', filen, '_shifts.dat'
    printf,unit0, '# image offsets (dx,dy) for ', filen, ' from centroiding'
    printf, unit0, '#   and initial guess of image offsets (dx0,dy0)'
    printf,unit0, '# '+strc(ngood)+' out of '+strc(nimgs)+' images used'
    printf,unit0, '#'
    printf,unit0, '# ', date
    printf,unit0, '# '
    printf,unit0, '#     dx           dy           dx0           dy0'
    ;printf,unit0, '#     dx           dy'
    printf,unit0, ' '
    free_lun, unit0
    printarray, [[obj_dx], [obj_dy], [obj_dx0], [obj_dy0]], $
      out = filen+'_shifts.dat', /silent
    ;printarray, [[obj_dx], [obj_dy]], $
    ;  out = filen+'_shifts.dat', /silent
endif


; write sky offsets file
openw, unit0, filen+'_mossky.dat', /get_lun
printf,unit0, '# ', filen, '_mossky.dat'
printf,unit0, '# sky offsets for ', filen, ' frames from mosaicing routine'
printf,unit0, '#'
printf,unit0, '# ', date
printf,unit0, '# '
printf,unit0, '# frame number    sky offset'
printf,unit0, ' '
free_lun, unit0
printarray, [[wgood], [obj_mossky]],  $
  out =  filen+'_mossky.dat', /silent


; make a PS file of the final mosaic
if getyn('make a PS file of the final mosaic?', 0) then $
  ps_image, obj_mos, filen


end