; quick redux of dithered Lick AO data
; 06/11/01 M. Liu
; 06/13/01 mods by JPL to add flatfield
; 06/23/01 hacked for Gemini AO 
;          (diff file names, no coadds, flexible image size, more parms)
; 07/03/01 parameter CBOX now passed to RECENTER.PRO
;          subtracts DC sky found from 1st pass mosaic from indiv images
;          bug fix: wasn't subtracting sky frame on 2nd pass mosaic!
; 07/04/01 added 'satflag' so can crudely reduce saturated images
;          (main limit is registering the images -> position mouse carefully)
;          lots of little improvements
; 09/16/01 slight tweak: for individual frames, does sky subtraction
;            before doing the recentering (matters for faint sources)
;          passes bad pixel mask to RECENTER
;          bug fix: divides by flat field if passed
; 09/19/01 for saturated images, use cross correlation to register
;            (IMGCOR.PRO) -> much better
; 09/29/01 write results from AORADNORM to file, also use /nofloor flag
;          **now divides by number of coadds**


goto, SPT
;goto, G_SAVE

objn = 'psf-hd130322.ns5'
istart = 272
iend =   280
nstack = 1   ; # of images per pointing
CBOX =   15    ; centering box size in pixels
satflag= 1  ; if image is saturated, use diff centroiding parms

fdir = 'raw/15may01/scam/'  ; data directory
fname = '15mai0' ; file prefix
extn = '.fits'
flatf = ''         ; flatfield
badpixf ='scam-badpix+xslit2.aug01.fits' ; bad pixel mask

pclip = 1.0  ; percentage of images to keep
setsky = 1   ; calc sky shifts for mosaic (slow for big images)
sznorm = 200 ; region size for AORADNORM (slow for big images), 0=full image

scam = 1

;--------------------------------------------------
SDIR = '~/idl/reduct/scripts/ao/'
SNAME = 'redircal.idl'


loadct, 15

nimgs = iend-istart+1
if (nimgs/nstack) ne nimgs/float(nstack) then  $
  message, 'nstack does not divide evenly into images!'


; history preface
hstr = ['<<'+strupcase(SNAME)+': '+ systime()+ userid()+ '>>',  $
        '  fdir = '+fdir, $
        '  fname = '+fname, $
        '  extn = '+extn, $
        '  istart = '+strc(istart), $
        '  iend = '+strc(iend), $
        '  nstack = '+strc(nstack)]


; get images
listims, istart, iend, out = objn+'-in.lst',  $
  fdir = fdir, fname = fname, ext = extn, /silent
loadimf, objn+'-in.lst', imgs, head0, histstart = hstr, scam = scam
imsz = (size(imgs))(1)


; load bad pixel mask
print, '=== loading bad pixel mask: "', badpixf, '" ==='
if (badpixf ne '') then badpix =  readfits(badpixf, /silent) $
else badpix = bytarr(imsz, imsz)+1B
if (n_elements(badpix) eq 1) then begin
    print, '** bad pixel image does not exist !! **^G'
    stop
endif else $
  sxaddhist, '  loaded bad pixel mask = "'+badpixf+'"', head0


; flat field
print, '=== loading flat field: "', flatf, '" ==='
if (flatf ne '') then flat =  readfits(flatf, /silent) $
else flat = fltarr(imsz, imsz)+1.0
if stdev(flat) eq 0 then begin
    print, 'using dummy flat field (all pix=1.0)'
    sxaddhist, 'no flat field loaded', head0
endif else begin
    sxaddhist, 'loaded flat field = "'+flatf+'"', head0
    flatdiv, imgs, flat, imgs, badpix = badpix
endelse  


; stack images at same position
print, '=== stacking images at each position ==='
if (nstack gt 1) then begin
    npos = (iend-istart+1)/float(nstack)
    if (npos ne fix(npos)) then  $
      message, 'non-integer positions!', /info
    imgsum = fltarr(imsz, imsz, npos)
    print, form = '($,"position: ")'
    for i = 0, npos-1 do begin
        print, form = '($,A," ")', strc(i)
        imgsum(*, *, i) =  $
          total(imgs(*, *, nstack*i:nstack*(i+1)-1), 3)/nstack
    endfor
endif else $
  imgsum = imgs


stop

; sky sub
avgmed, sky, imgsum, /noweight
skysub, imgsum, sky, imgsumsub
sxaddhist, 'making median sky frame & subtracting', head0

SPT:

; reg + mosaic
if (satflag ne 1) then begin
    ireg, imgsumsub, dx, dy, mos, mosexp, mossky, /setsky, $
      log = objn+'_shifts.dat', badpix = badpix,  $
      cbox = CBOX, /tv
endif else begin
    print, '=== registering saturated images via cross-correlation ==='
    imgcor, imgsumsub, dx, dy, badpix = badpix
    mosf, imgsumsub, dx, dy, mos, mosexp, mossky,  $
      setsky = setsky, badpix = badpix
    sxaddhist, 'images are saturated (satflag=1)', head0
    sxaddhist, 'using IMGCOR.PRO to find shifts via x-correlation', head0
    ; still need to identify image center for later work (aoradnorm)
    ireg, mos, x, y, cbox = CBOX, /nomaxpix, /tv
    dx = dx+(x-max(dx))
    dy = dy+(y-max(dy))
    ; write shifts
    outfile = objn+'_shifts.dat'
    if filecheck(outfile, /over) then begin
        forprint2, [dx], [dy], outfile = outfile,  $
          comments = [strupcase(SNAME)+': '+systime()+userid(),  $
                      'used IMGCOR.PRO to find shifts via cross correlation', $
                      'then used IREG.PRO to identify image center', ' ',  $
                      '      dx        dy']
    endif
endelse
xstar = max(dx)
ystar = max(dy)
win, 0
display2, mos, tit = '1st pass mosaic'


; if there are multiple images per position
; we can do image selection based on FWHM
if (nstack ne 0) and (satflag ne 1) then begin
    fwhm = fltarr(nimgs)
    dx2 = fltarr(nimgs)
    dy2 = fltarr(nimgs)

    ; subtract the sky
    skysub, imgs, sky, imgsub

    ; loop through images to adjust centers & get FWHMs
    print, form = '($,A," ")', 'FWHM computation: '
    for i = 0, nimgs-1 do begin
        if (i mod 10) eq 0 then  $
          print, form = '($,A," ")', strc(i)
        j = i/nstack
        ; touch up centering
        recenter, imgsub(*, *, i), dx(j), dy(j), x, y,  $
          cbox = CBOX, badpix = badpix, /nodisp, /silent
        dx2(i) = x
        dy2(i) = y
        ; compute FWHM
        radplotf, imgsub(*, *, i), x, y, radout, ff, $
          sky = 1e-20, outrad = 3*CBOX > 15 < 30,  $
          badpix = badpix, /silent
        fwhm(i) = ff
        ; remove best guess sky level
        imgsub(*, *, i) = imsub(imgsub(*, *, i), -mossky(j))
    endfor
    print

    ; plot results
    win, 0
    w = where(fwhm ne 999, nw)
    plot, fwhm, ps = 6, yr = [0, max(fwhm(w))], $
      xtit = 'image number', ytit = 'FWHM (pix)'
    wbad = where(fwhm eq 999, nbad)
    if (nbad ge 1) then begin
        oplot, wbad, intarr(nbad), ps = 8, sym = 2; FWHM=999 objects
        reg2dev, 0.1, 0.1, xd, yd
        xyouts, xd, yd, /dev, $
          '# of images with bad FWHM (999) = '+strc(nbad)
    endif

    ; draw line for best subset as a guide
    ss = sort(fwhm)
    cc = (fwhm(ss))(round(nimgs*pclip)-1)
    oplot, !x.crange, [cc, cc]
    print
    print, 'keeping '+strc(pclip*100)+'% of the images'
    print, 'FWHM threshold = ', strc(cc)
    print, 'median FWHM = ', strc(median(fwhm))
    print, 'average FWHM = ', strc(avg(fwhm))
    win, 1
    plothist, fwhm < (max(fwhm(w)) < 15), bin = 0.5
    oplot, [cc, cc], !y.crange, line = 5

    ; select cut
    fwhmcut = getnum('enter FWHM threshold:', cc)
    plots, !x.crange, fwhmcut+[0, 0], line = 5

    ; make a new mosaic
    ; if don't use /setsky, at least remove the median value
    goodlist = fwhm le fwhmcut
    wgood = where(fwhm le fwhmcut)
    ngood = total(goodlist)
    mosf, imgsub, dx2, dy2, mosnew, mosexpnew,  $
      good = goodlist, badpix = badpix,  $
      setsky = setsky, submed = (setsky eq 0)
    sxaddhist, 'images selected by FWHM <= '+strc(fwhmcut)+' pix', head0
    sxaddhist, 'number of good images = '+strc(ngood), head0

    ; compare with old mosaic
    radplotf, mos, max(dx), max(dy), radout, f0, sky = 1e-20, /silent
    radplotf, mosnew, max(dx2(wgood)), max(dy2(wgood)), radout, f1,  $
      sky = 1e-20, /silent
    whereismax, mos, x, y, maxv0, /silent
    whereismax, mosnew, x, y, maxv1, /silent
    print, 'old mosaic: FWHM = ', strc(f0), ', peak = ', strc(maxv0)
    print, 'new mosaic: FWHM = ', strc(f1), ', peak = ', strc(maxv1)
    print, 'FWHM improves by ', strc(f0/f1)
    print, 'peak improves by ', strc(maxv1/maxv0)
    if not(getyn('satisfied?')) then stop

    ; keep improved mosaic
    mos = mosnew
    mosexp = mosexpnew
    xstar = max(dx2(wgood))
    ystar = max(dy2(wgood))

    ; write list of new shifts, FWHM, and good images
    outfile = objn+'_shifts.good.dat'
    if filecheck(outfile, /over) then begin
        openw, unit0, outfile, /get_lun
        printf, unit0, '# ', strupcase(SNAME), ': ', systime(), userid()
        printf, unit0, '# cut used for FWHM = ', strc(fwhmcut)
        printf, unit0, '#'
        printf, unit0, '#     xpos       ypos      FWHM(pix)      1=good'
        printf, unit0, ''
        free_lun, unit0
        printarray, [[dx2], [dy2], [fwhm], [goodlist]],  $
          out = outfile, /silent
    endif

endif


; AO radial normalize
print, '--- normalizing image by radial noise profile ("aoradnorm") ---'
if (sznorm eq 0) then begin
  tmp = mos 
  x = xstar
  y = ystar
endif else begin
  tmp = imcut(mos, sznorm, xstar, ystar)
  x = sznorm/2
  y = sznorm/2
endelse
aoradnorm, tmp, aonorm, xcen = x, ycen = y,  $
  out = objn+'_aonorm.dat', /nofloor


; display results
win, 0 & display2, mos, tit = 'mosaic'
win, 1 & display2, aonorm > (-1) < 5, /tv, tit = 'masked & filtered'
;win, 2 & display2, imsub(mos, median(mos, 15)) < 20, tit = 'unsharp masked'


G_SAVE:

; write files
if getyn('save files?') then begin

    outfile = objn+'_mos.fits'
    if filecheck(outfile, /over) then begin
      writefits, outfile, mos, head0
      writefits, objn+'_mosexp.fits', mosexp, head0
      writefits, objn+'_aonorm.fits', aonorm, head0

      ; write parm file
      outfile = objn+'.parm'
      openw, unit0, outfile, /get_lun
      printf, unit0, '# ', strupcase(SNAME), ': ', systime(), userid()
      printf, unit0, ''
      printf, unit0, 'objn    ', objn
      printf, unit0, 'istart  ', strc(istart)
      printf, unit0, 'iend    ', strc(iend)
      printf, unit0, 'nstack  ', strc(nstack)
      printf, unit0, 'pclip   ', strc(pclip)
      printf, unit0, 'fdir    ', fdir
      printf, unit0, 'fname   ', fname
      printf, unit0, 'badpixf ', badpixf
      printf, unit0, 'flatf   ', flatf
      printf, unit0, 'cbox    ', cbox
      printf, unit0, 'satflag ', satflag
      printf, unit0, 'setsky  ', setsky
      printf, unit0, 'sznorm  ', sznorm
      free_lun, unit0

      ; save copy of reduction script
      outfile = SNAME+'.'+objn
      print, '  current script dir: "', SDIR, '"'
      print, '  current IDL script: "', SNAME, '"'
      print, '  output file name: "', outfile, '"'
      if not(getyn('  is this correct?')) then stop
      if filecheck(outfile, /over) then $
        spawn, '\cp '+SDIR+SNAME+' '+outfile
      print, 'wrote "', outfile, '" to current directory'

  endif
endif



end