;+
; NAME:
;	PH_REGISTER
;
; PURPOSE:
;	Reduction of dithered PHARO images of faint targets.
;
; EXPLANATION:
;	Determines telescope pointing from image FITS headers, then ajusts
;	relative positions of images by centroiding a point source selected
;	by user. Outputs a stack or registered images.
;
; CALLING SEQUENCE:
;	PH_REGISTER, stack, KWD=KWD, ...
;
; INPUTS:
;       objfr=  Vector of frame numbers for object.  If not given, will
;               prompt for first and last frame numbers.
;
;       skyfr=  Vector of frame numbers for sky frames.  If not given, will
;               use "off" pointings in object frames.
;
;	flatname= Filename of flat to use.  If not specified, uses
;		'flat_X_YYmas.fits' where X is filter and YY is platescale.
;
;	bpixname= Filename of bad pixel map to use.  If not specified, uses
;		'bpix_X_YYmas.fits' where X is filter and YY is platescale.
;
;       datadir= Full path to data directory.  If not specified, determined
;               from date.
;
;       procdir= Full path to processing directory.  If not specified,
;               determined from date.
;
;       caldir= Full path to calibration file directory.  If not specified,
;               assumed to be same as processing directory
;
;	mask=	1024x1024 pixel mask image (or weight) to give each image
;		in the mosaic.  Use this to mask out bad regions of the
;		detector.
;
;	/autocray Auto-filter cosmic rays (watch out - can chop top off PSF!)
;
;       /noshift Do no cross-correlate and shift images in final stack.
;
;	skyrej=	Fraction of highest pixel values to reject from median stack.
;		Default = 0.0
;
;       /silent Supress informational messages
;
; OUTPUTS:
;	stack = stacked image (XxYxN) where the layers are:
;		0	Average of sky-subtracted frames.
;		1	Number of frames which contributed to each pixel
;
;       img = Mean image over entire stack.
;
;	Writes to directory 'procdir' the following files:
;		ph????r1.fits	Reduced individual frames.
;		ph????r2d.fits	Mean image of overlap region only.
;
; MODIFICATION HISTORY:
;       Original, based on PH_DITHER, written 09/08
;					A. Bouchez, Caltech Optical Obs.
;-
PRO     PH_REGISTER, im6, im5, objfr=objfr, skyfr=skyfr, $
          flatname=flatname, bpixname=bpixname, $
          datadir=datadir, procdir=procdir, caldir=caldir, $
          mask=mask, autocray=autocray, noshift=noshift, skyrej=skyrej, $
          silent=silent, mag=mag, best=best

;;; Define constants

  fmt = '("ph",I4.4,".fits.gz")'        ; PHARO filename format
  fmt_r1 = '("ph",I4.4,"r1.fits")'	; Filename format for saved frames
  fmt_r2s = '("ph",I4.4,"r2m.fits")'	; Filename format for saved image stack
  fmt_r2d = '("ph",I4.4,"r2d.fits")'	; Filename format for saved coadd image
  cr_zpt = 334.28                       ; Cass ring zeropoint (N up) (9/08)
  drd = 1.0				; Dither position radius (")
  csig = 10.0				; Cosmic ray filter threshold

;;; Parse keywords

  if not KEYWORD_SET(noshift) then noshift = 0
  if not KEYWORD_SET(silent) then silent = 0
  if not KEYWORD_SET(autocray) then autocray = 0
  if not KEYWORD_SET(skyrej) then skyrej = 0
  if not KEYWORD_SET(mag) then mag = 4.0	; x-correlation magnification
  if not KEYWORD_SET(best) then best = 1.0	; fraction to retain

  if not KEYWORD_SET(datadir) then begin
    str = ''
    READ, 'Please enter data directory: ', str
    datadir = str
  endif
  if not KEYWORD_SET(procdir) then procdir = './'
  if not KEYWORD_SET(caldir) then caldir = procdir

  if not KEYWORD_SET(objfr) then begin
    fstr= ''
    READ, 'Please enter OBJECT first frame number: ', fstr
    lstr = ''
    READ, 'Please enter OBJECT last frame number: ', lstr
    flno = FIX([fstr, lstr])
    objfr = flno[0] + LINDGEN(flno[1]-flno[0]+1)
  endif

  if not KEYWORD_SET(skyfr) then begin
    if not silent then $
      MESSAGE, /info, 'Using dithered object frames for sky.'
    skyfr=objfr
    dsky = 1B
  endif else dsky = 0B

;;; Determine image size, filter, platscale, and PA from first image header

  hd0 = HEADFITS(datadir + STRING(objfr[0],f=fmt), /silent)
  sz0 = [2*SXPAR(hd0,'NAXIS1'), 2*SXPAR(hd0,'NAXIS2'), N_ELEMENTS(objfr)]

  filt = STRTRIM(SXPAR(hd0,'FILTER'),2)		; USED ONLY TO CHOOSE FLAT
  case filt of
    'K_short': filt='Ks'
    'Br-gamma': filt='Ks'
    else: foo=0					; Add more lines as necessary!
  endcase

  camname = STRTRIM(SXPAR(hd0,'CAROUSEL'),2)
  case camname of
    '25 mas': pscl = 0.025
    '40 mas': pscl = 0.040
  endcase

  pa = (SXPAR(hd0,'CR_ANGLE') - cr_zpt + 360) mod 360		; PA of N

;;; Read flat and bad pixel maps (does not yet deal correctly with subarrays!)

  if not KEYWORD_SET(flatname) then $
    flatname = 'flat_' + filt + '_' + $
      STRING(ROUND(pscl*1e3),f='(I0.0)') + 'mas.fits'
  tmp = FILE_SEARCH(caldir + flatname)
  while (STRLEN(tmp[0]) eq 0) do begin
    MESSAGE, /info, 'Flat file ' + flatname + ' not found.'
    READ, 'Please enter a new file name: ', flatname
    tmp = FILE_SEARCH(caldir + flatname)
  endwhile
  if not silent then $
    MESSAGE, /info, 'Using flat: ' + flatname
  flat = READFITS(caldir + flatname, /silent)

  if not KEYWORD_SET(bpixname) then $
    bpixname = 'bpix_' + filt + '_' + $
      STRING(ROUND(pscl*1e3),f='(I0.0)') + 'mas.fits'
  tmp = FILE_SEARCH(caldir + bpixname)
  while (STRLEN(tmp[0]) eq 0) do begin
    MESSAGE, /info, 'Bad pixel file ' + bpixname + ' not found.'
    READ, 'Please enter a new file name: ', bpixname
    tmp = FILE_SEARCH(caldir + bpixname)
  endwhile
  if not silent then $
    MESSAGE, /info, 'Using bpix: ' + bpixname
  bpix = READFITS(caldir + bpixname, /silent)

;;; Sort all object frames into distinct pointings based on RA/Dec

  nfr = N_ELEMENTS(objfr)
  hd = STRARR(N_ELEMENTS(hd0), nfr)
  rade = DBLARR(2, nfr)
  for n=0,nfr-1 do begin
    hd[*,n] = HEADFITS(datadir + STRING(objfr[n],f=fmt), /silent)
    rade[*,n] = [SXPAR(hd[*,n],'CRVAL1'), SXPAR(hd[*,n],'CRVAL2')]	; deg
  endfor

  mrade = [TOTAL(MINMAX(rade[0,*])), $			; Midpoint RA/Dec (deg)
           TOTAL(MINMAX(rade[1,*]))] / 2.
  rdoff = (rade - REBIN(mrade,2,nfr)) * $		; Offset RA/Dec (")
    [COS(rade[1,*]*!dtor), REPLICATE(1.,1,nfr)] * 3600
  dgroup = ROUND(rdoff[0,*]/drd) + 100L * ROUND(rdoff[1,*]/drd)
  duniq = dgroup(UNIQ(dgroup))
  dpos = N_ELEMENTS(duniq)
  if not silent then $
    MESSAGE, /info, 'Found ' + STRING(dpos,f='(I0.0)') + $
    ' distinct object pointings'
  didx = [REFORM(objfr,1,nfr), LONARR(1,nfr)]		; [frno, pointing]
  for i=0,dpos-1 do $
    didx[1, WHERE(dgroup eq duniq[i])] = i

;;; Compute median sky frame for each pointing

  if dsky then begin				; if using dithered sky frames
    sky = LONARR(sz0[0], sz0[1], dpos)
    skydat = LONARR(3, dpos)
    for n=0,dpos-1 do begin
      offfr = objfr[WHERE(didx[1,*] ne n)]
      skydat[*,n] = [N_ELEMENTS(offfr), MINMAX(offfr)]
      if not silent then $
        MESSAGE, /info, 'Computing sky ' + $
        STRING([n+1,dpos],f='(I0.0," of ",I0.0)') + ' using ' + $
        STRING(skydat[0,n],f='(I0.0)') + ' frames.'
      sky[*,*,n] = MEDIM(datadir + STRING(offfr,f=fmt), /pharo, reject=skyrej)
    endfor

  endif else begin				; if using indep. sky frames
    nskyfr = N_ELEMENTS(skyfr)
    if not silent then $
      MESSAGE, /info, 'Computing median of ' + $
      STRING(nskyfr,f='(I0.0)') + ' independent sky frames.'
    onesky = MEDIM(datadir + STRING(skyfr,f=fmt), /pharo, reject=skyrej)
    sky = REBIN(onesky, sz0[0], sz0[1], dpos)
    skydat = REBIN([nskyfr, MINMAX(skyfr)], 3, dpos)
  endelse

;;; Determine pixel offsets of obj frames, set up image arrays, mosaic weights

  xyoff = FLTARR(2,sz0[2])
  rot = [[   COS(pa*!dtor), SIN(pa*!dtor)], $
         [-1*SIN(pa*!dtor), COS(pa*!dtor)]]
  for n=0,sz0[2]-1 do $
    xyoff[*,n] = (([-1,1]*rdoff[*,WHERE(didx[0,*] eq objfr[n])]) # rot)/pscl

  if not silent then $
    MESSAGE, /info, 'Creating image arrays (up to ' + $
    STRING(sz0,f='(I0.0,"x",I0.0,"x3x",I0.0)') + ')'
  im4 = FLTARR(sz0[0],sz0[1],sz0[2])
  im4sky = FLTARR(sz0[2])
  im5 = FLTARR(sz0[0],sz0[1],3,sz0[2])		; 0=img, 1=sky-sub, 2=activ pix
  im6 = FLTARR(sz0[0],sz0[1],3)			; same
  xx = REBIN(INDGEN(sz0[0],1),sz0[0],sz0[1])
  yy = REBIN(INDGEN(1,sz0[1]),sz0[0],sz0[1])

  if not KEYWORD_SET(mask) then begin
    wt = REPLICATE(1.0, sz0[0], sz0[1])
    ;DIST_CIRCLE,r,sz0[0],sz0[0]/2-0.5,sz0[1]/2-0.5	; dist from img center
    ;wt = 1-(SQRT(2)*r/sz0[0])^2.0
    ;wt[0,*] = 0.				; FOR 25" FOV ONLY!
    ;wt[0:30,*] = 0.				; FOR 40" FOV ONLY!
  endif else wt = mask

  c0 = [0, 0, sz0[0:1]-1]
  c1 = FLTARR(4, sz0[2])
  c2 = FLTARR(4, sz0[2])

;;; Reduce individual object frames.

  dsz = [512, 512]
  WINDOW, 10, xs=dsz[0], ys=dsz[1]
  for n=0,sz0[2]-1 do begin
    im0 = RFITS(datadir + STRING(objfr[n],f=fmt), hd, /pharo, /silent)
    idx = didx[1,(WHERE(didx[0,*] eq objfr[n]))[0]]
    im1 = (im0 - sky[*,*,idx]) / flat
    im2 = BPIXFIX(im1, bpix, neighbors=3)
    if autocray then $
      im3 = CRAYFIX(im2, t=csig) else $
      im3 = im2

    SXADDPAR, hd, 'BITPIX', -32
    SXADDPAR, hd, 'NAXIS', 2
    SXADDPAR, hd, 'NAXIS1', sz0[0]
    SXADDPAR, hd, 'NAXIS2', sz0[1]
    SXDELPAR, hd, 'NAXIS3'
    SXADDPAR, hd, 'IMGLEVEL', 1, ' Image processing level (pre-reduced)'
    SXADDPAR, hd, 'IMGTYPE', 'frame', ' Image type (frame/mosaic/dither)'
    SXADDPAR, hd, 'DITHPOS', idx, ' Dither position'
    SXADDPAR, hd, 'SKYNFR', skydat[0,idx], ' Number of sky frames'
    SXADDPAR, hd, 'SKYFRNO0', skydat[1,idx], ' First sky frame'
    SXADDPAR, hd, 'SKYFRNO1', skydat[2,idx], ' Last sky frame'
    SXADDPAR, hd, 'FLATNAME', flatname
    SXADDPAR, hd, 'BPIXNAME', bpixname
    if autocray then $
      SXADDPAR, hd, 'CRTHRESH', csig, ' Cosmic ray threshold (sigma)' else $
      SXADDPAR, hd, 'CRTHRESH', 0, ' Cosmic ray filtering disabled'
    sname = STRING(objfr[n],f=fmt_r1)
    TVSCL, REBIN(im3>0<1e3, dsz[0], dsz[1])^0.5
    XYOUTS, .5, .95, sname, chars=2, /norm, align=0.5
    if not silent then $
      MESSAGE, /info, STRING([n+1,sz0[2]],f='(I2.2,"/",I2.2)') + $
      ' Saving ' + sname
    WRITEFITS, procdir + sname, im3, hd
    if n eq 0 then hd0 = hd

;;; Insert at estimated location in stack

    im4[*,*,n] = im3			; NO DISTORTION CORRECTION YET!

    c1[*,n] = ROUND(c0 + [xyoff[[0,1,0,1],n]])
    c2[*,n] = ROUND(c0 - [xyoff[[0,1,0,1],n]])
    SKY, im4[*,*,n], skyval, /silent
    im4sky[n] = skyval
    im5[c1[0,n]>0:c1[2,n]<(sz0[0]-1), c1[1,n]>0:c1[3,n]<(sz0[1]-1), 0, n] = $
      im4[c2[0,n]>0:c2[2,n]<(sz0[0]-1), c2[1,n]>0:c2[3,n]<(sz0[1]-1), n]
    im5[c1[0,n]>0:c1[2,n]<(sz0[0]-1), c1[1,n]>0:c1[3,n]<(sz0[1]-1), 1, n] = $
      im4[c2[0,n]>0:c2[2,n]<(sz0[0]-1), c2[1,n]>0:c2[3,n]<(sz0[1]-1), n] - $
      im4sky[n]
    im5[c1[0,n]>0:c1[2,n]<(sz0[0]-1), c1[1,n]>0:c1[3,n]<(sz0[1]-1), 2, n] = $
      wt[c2[0,n]>0:c2[2,n]<(sz0[0]-1), c2[1,n]>0:c2[3,n]<(sz0[1]-1)]

;;; Compute mean image

    im6[*,*,2] = REBIN(im5[*,*,2,0:n], sz0[0], sz0[1])	; normalized weight
    wt6 = im6[*,*,2]
    w = WHERE(wt6 eq 0)
    if (w[0] ne -1) then wt6(w) = 1.0
    im6[*,*,0] = REBIN(im5[*,*,0,0:n]*im5[*,*,2,0:n], sz0[0], sz0[1]) / wt6
    im6[*,*,1] = REBIN(im5[*,*,1,0:n]*im5[*,*,2,0:n], sz0[0], sz0[1]) / wt6
  endfor

;;; Display and request user input for registration point

  sz1 = [256, 256]
  r = FLOAT(dsz) / sz0[0:1]
  TVSCL, REBIN(im6[*,*,1]>0<1e3, dsz[0], dsz[1])^0.5
  XYOUTS, .5, .95, 'Click on bright source', chars=2, /norm, align=0.5
  CURSOR, x, y, /down, /device
  TVBOX, sz1*r, x, y, /device
  xy0 = [x,y]/r

;;; Cross-correlate subimages with first image in stack

  c3 = [xy0-sz1/2, xy0+(sz1/2)-1]
  MESSAGE, /info, 'Cross-correlating with first image over ' + $
    STRING(sz1,f='(I0.0,"x",I0.0)') + ' region.'
  for n=0,sz0[2]-1 do begin
    im5a = im5[c3[0]:c3[2], c3[1]:c3[3], 1, n]
    im5b = im5[c3[0]:c3[2], c3[1]:c3[3], 1, 0]
    u = [FINDGEN(sz1[0]/2),FINDGEN(sz1[0]/2)-(sz1[0]/2)]
    v = [[FINDGEN(1,sz1[1]/2)], $
         [FINDGEN(1,sz1[1]/2)-(sz1[1]/2)]]
    uu = REBIN(u,sz1[0],sz1[1])
    vv = REBIN(v,sz1[0],sz1[1])
    cc = FFT( PFFT(im5a,1) * CONJ(PFFT(im5b,1)), 1)
    maxcc = MAX(cc, mp)
    mshift = MAX(ABS([uu[mp],vv[mp]]))

    dxy = [uu[mp], vv[mp]]
    MESSAGE, /info, STRING([n+1,sz0[2]],f='(I2.2,"/",I2.2)') + $
      ' Shift = ' + STRING(dxy,f='("(",F4.0,",",F4.0,")")')

    c1[*,n] = c1[*,n] - dxy[[0,1,0,1]]
    c2[*,n] = c2[*,n] + dxy[[0,1,0,1]]
    im5[c1[0,n]>0:c1[2,n]<(sz0[0]-1), c1[1,n]>0:c1[3,n]<(sz0[1]-1), 0, n] = $
      im4[c2[0,n]>0:c2[2,n]<(sz0[0]-1), c2[1,n]>0:c2[3,n]<(sz0[1]-1), n]
    im5[c1[0,n]>0:c1[2,n]<(sz0[0]-1), c1[1,n]>0:c1[3,n]<(sz0[1]-1), 1, n] = $
      im4[c2[0,n]>0:c2[2,n]<(sz0[0]-1), c2[1,n]>0:c2[3,n]<(sz0[1]-1), n] - $
      im4sky[n]
    im5[c1[0,n]>0:c1[2,n]<(sz0[0]-1), c1[1,n]>0:c1[3,n]<(sz0[1]-1), 2, n] = $
      wt[c2[0,n]>0:c2[2,n]<(sz0[0]-1), c2[1,n]>0:c2[3,n]<(sz0[1]-1)]

;;; Compute mean image over stack, display

    im6[*,*,2] = REBIN(im5[*,*,2,0:n], sz0[0], sz0[1])	; normalized weight
    wt6 = im6[*,*,2]
    w = WHERE(wt6 eq 0)
    if (w[0] ne -1) then wt6(w) = 1.0
    im6[*,*,0] = REBIN(im5[*,*,0,0:n]*im5[*,*,2,0:n], sz0[0], sz0[1]) / wt6
    im6[*,*,1] = REBIN(im5[*,*,1,0:n]*im5[*,*,2,0:n], sz0[0], sz0[1]) / wt6
    TVSCL, REBIN(im6[*,*,1]>0<1e3, dsz[0], dsz[1])^0.5
  endfor
  im6_v1 = im6

;;; Display and request user input for registration point

  sz1 = [64, 64]
  c4 = xy0[[0,1,0,1]] + [dsz/(-2), (dsz/2)-1]
  TVSCL, (im6[c4[0]:c4[2],c4[1]:c4[3],1]>0<1e3)^0.5
  XYOUTS, .5, .95, 'Click on point source', chars=2, /norm, align=0.5
  CURSOR, x, y, /down, /device
  TVBOX, sz1, x, y, /device
  xy0 = c4[0:1] + [x,y]

;;; Repeat, cross-correlating subimages with integer-shift mean

  c3 = [xy0-sz1/2, xy0+(sz1/2)-1]
  peak = FLTARR(sz0[2])
  MESSAGE, /info, 'Cross-correlating with mean image over ' + $
    STRING(sz1,f='(I0.0,"x",I0.0)') + ' region.'
  for n=0,sz0[2]-1 do begin
    im5a = im6_v1[c3[0]:c3[2], c3[1]:c3[3], 1]
    im5b = im5[c3[0]:c3[2], c3[1]:c3[3], 1, n]
    u = [FINDGEN(sz1[0]*mag/2),FINDGEN(sz1[0]*mag/2)-(sz1[0]*mag/2)]
    v = [[FINDGEN(1,sz1[1]*mag/2)], $
         [FINDGEN(1,sz1[1]*mag/2)-(sz1[1]*mag/2)]]
    uu = REBIN(u,sz1[0]*mag,sz1[1]*mag)/mag
    vv = REBIN(v,sz1[0]*mag,sz1[1]*mag)/mag
    cc = FFT( PFFT(im5a,mag) * CONJ(PFFT(im5b,mag)), 1)
    maxcc = MAX(cc, mp)
    mshift = MAX(ABS([uu[mp],vv[mp]]))

    dxy = [uu[mp], vv[mp]]
    MESSAGE, /info, STRING([n+1,sz0[2]],f='(I2.2,"/",I2.2)') + $
      ' Fourier shift = ' + STRING(dxy,f='("(",F5.2,",",F5.2,")")')
    im5[*,*,0,n] = FSHIFT(REFORM(im5[*,*,0,n]), dxy[0], dxy[1])
    im5[*,*,1,n] = FSHIFT(REFORM(im5[*,*,1,n]), dxy[0], dxy[1])
    ;im5[*,*,2,n] = FSHIFT(REFORM(im5[*,*,2,n]), dxy[0], dxy[1])
    peak[n] = MAX(im5[c3[0]:c3[2], c3[1]:c3[3], 1, n])

;;; Compute mean image over stack, display

    im6[*,*,2] = REBIN(im5[*,*,2,0:n], sz0[0], sz0[1])	; normalized weight
    wt6 = im6[*,*,2]
    w = WHERE(wt6 eq 0)
    if (w[0] ne -1) then wt6(w) = 1.0
    im6[*,*,0] = REBIN(im5[*,*,0,0:n]*im5[*,*,2,0:n], sz0[0], sz0[1]) / wt6
    im6[*,*,1] = REBIN(im5[*,*,1,0:n]*im5[*,*,2,0:n], sz0[0], sz0[1]) / wt6
    TVSCL, (im6[c4[0]:c4[2],c4[1]:c4[3],1]>0<1e3)^0.5
  endfor

;;; Retain only best images in mean

  MESSAGE, /info, 'Retaining best ' + STRING(best*100,f='(I0.0,"%")')
  idx = REVERSE(SORT(peak))
  keep = idx[0:ROUND(sz0[2]*best)-1]
  junk = idx[ROUND(sz0[2]*best):sz0[2]-1]
  fmt = '(I0.0,' + STRTRIM(N_ELEMENTS(junk)-1,2) + '(", ",I0.0))'
  MESSAGE, /info, 'Throwing out frames: '
  MESSAGE, /info, '  ' + STRING(objfr(junk),f=fmt)

  im6[*,*,2] = REBIN(im5[*,*,2,keep], sz0[0], sz0[1])
  wt6 = im6[*,*,2]
  w = WHERE(wt6 eq 0)
  if (w[0] ne -1) then wt6[w] = 1.
  im6[*,*,0] = REBIN(im5[*,*,0,keep]*im5[*,*,2,keep], sz0[0], sz0[1]) / wt6
  im6[*,*,1] = REBIN(im5[*,*,1,keep]*im5[*,*,2,keep], sz0[0], sz0[1]) / wt6
  TVSCL, REBIN(im6[*,*,1]>0<1e3, dsz[0], dsz[1])^0.5
  TVSCL, CONGRID(im6[c3[0]:c3[2],c3[1]:c3[3],1]>0,sz1[0]*2,sz1[1]*2)^0.5
  TVBOX, sz1[0:1]*2, sz1[0], sz1[1]
  t_eff = SXPAR(hd0,'T_EFF')*N_ELEMENTS(keep)
  MESSAGE, /info, 'Total integration time = ' + $
    STRING(t_eff*1e-3, f='(F6.1)') + ' s'

  WINDOW,11,xs=512,ys=384
  PLOT, peak, /xsty, xtit='Frame', ytit='Peak DN', chars=1.5, /noerase
  OPLOT, peak, ps=4
  OPLOT, [0,sz0[2]], peak[idx[ROUND(sz0[2]*best)]] * [1,1], line=2

;;; Save mean image

  im6 = im6[*,*,1:2]
  SXADDPAR, hd0, 'NAXIS', 3
  SXADDPAR, hd0, 'NAXIS1', sz0[0]
  SXADDPAR, hd0, 'NAXIS2', sz0[1]
  SXADDPAR, hd0, 'NAXIS3', 2, after='NAXIS2'
  SXADDPAR, hd0, 'IMGLEVEL', 2, ' Image processing level (mean image)'
  SXADDPAR, hd0, 'IMGTYPE', 'mosaic', ' Image type (frame/mosaic/dither)'
  SxADDPAR, hd0, 'T_EFF', t_eff
  SXADDPAR, hd0, 'OBSNFR', N_ELEMENTS(keep), ' Frames in dither pattern'
  SXADDPAR, hd0, 'FRNO0', MIN(objfr), ' First object frame'
  SXADDPAR, hd0, 'FRNO1', MAX(objfr), ' Last object frame'
  sname = STRING(objfr[0], f=fmt_r2d)
  if not silent then $
    MESSAGE, /info, 'Saving ' + STRING([sz0[0:1],2],f='(I0.0,2("x",I0.0))') + $
      ' mean as ' + sname
  WRITEFITS, procdir + sname, im6, hd0

END