; 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
;          (principal limitation is in registering the images)


objn = 'psf.hd181655.brg+nd2'
istart = 662
iend =   671
nstack = 5   ; # of images per pointing
CBOX =   20    ; centering box size in pixels
satflag= 0  ; if image is saturated, use diff centroiding parms

fdir = '~/data/hokupaa.01jun/jun24/'  ; data directory
fname = '01jun24_' ; file prefix
flatf = ''         ; flatfield
badpixf ='~/idl/masks/hokupaa-badpix.kp.fits' ; bad pixel mask

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

goto, SPT


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


loadct, 15

nimgs = iend-istart+1

; history preface
hstr = ['<<'+strupcase(SNAME)+': '+ systime()+ userid()+ '>>',  $
        '  fdir = '+fdir, $
        '  fname = '+fname, $
        '  istart = '+strc(istart), $
        '  iend = '+strc(iend), $
        '  nstack = '+strc(nstack)]
; hstr = strarr(5)
; hstr(0) = '<<'+strupcase(SNAME)+': '+ systime()+ userid()+ '>>'
; hstr(1) = '  fdir = '+fdir
; hstr(2) = '  fname = '+fname
; hstr(3) = '  istart = '+strc(istart)
; hstr(4) = '  iend = '+strc(iend)
; hstr(5) = '  nstack = '+strc(nstack)


; get images
listims, istart, iend, out = objn+'-in.lst',  $
  fdir = fdir, fname = fname, /silent
loadimf, objn+'-in.lst', imgs, head0, histstart = hstr, /coadd
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


; load 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
    ;if not(getyn('using dummy flat field (all pix='$
    ;             +strc(avg(flat))+'). continue?^G')) then stop
  sxaddhist, 'no flat field loaded', head0
endif else $
  sxaddhist, 'loaded flat field = "'+flatf+'"', head0


;;load flat field
; if (flatf ne '')  then begin
;     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
;         if not(getyn('using dummy flat field (all pix='$
;                      +strc(avg(flat))+'). continue?^G')) then stop
;         sxaddhist, '  not flat fielding', head0
;     endif else $
;       sxaddhist, '  loaded flat field = "'+flatf+'"', head0

;     ; flatten
;     print, '========================='
;     print, '=== flattening frames ==='
;     print, '========================='
;     flatdiv, imgs, flat, imgs, badpix = badpix
;     if stdev(flat) eq 0 then $
;       sxaddhist, '  ** no flat field applied **', head0 $
;     else $
;       sxaddhist, '  FLATDIV: flat-fielded images', head0
; endif else  $
;   sxaddhist, '  not flat fielding', head0
  


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



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


SPT:

; reg + mosaic
ireg, imgsub, dx, dy, mos, mosexp, mossky, /setsky, $
  log = objn+'_shifts.dat', badpix = badpix,  $
  cbox = CBOX, nomaxpix = satflag, /tv
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)

    ; 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, imgs0(*, *, i), dx(j), dy(j), x, y,  $
          cbox = CBOX, /nodisp, /silent
        dx2(i) = x
        dy2(i) = y
        ; compute FWHM
        radplotf, imgs0(*, *, i), x, y, radout, ff, $
          sky = 1e-20, outrad = 3*CBOX > 15 < 30, /silent
        fwhm(i) = ff
        ; remove best guess sky level
        imgs0(*, *, i) = imsub(imgs0(*, *, i), -mossky(j))
    endfor
    print

    ; subtract the sky
    skysub, imgs0, sky, imgsub

    ; 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, 'filtered by FWHM, thrshold = '+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) then begin
        openw, unit0, outfile, /get_lun
        printf, unit0, '# RED.IDL: ', 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


; 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'


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

    outfile = objn+'_mos.fits'
    if filecheck(outfile) 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) then $
        spawn, 'cp '+SDIR+SNAME+' '+outfile
      print, 'wrote "', outfile, '" to current directory'

  endif
endif



end