;+ ; NAME: ; make_spe_flat1 ; ; PURPOSE: (one line) ; make flat averaged files from spe files ; ; DESCRIPTION: ; Given a filenames of SPE zero-flux files, make averaged flat-field file ; ; CATEGORY: ; IMAGES ; ; CALLING SEQUENCE: ; make_spe_flat1,filename, avg, err, bad, d = d, $ ; silent=silent, median = median, $ ; imbad = imbad, gain=gain, $ ; zfilename = zfilename, $ ; caldir=caldir, outfilename = outfilename, $ ; overrideoutfilename = overrideoutfilename ; INPUTS: ; filename : filename of SPE image cube with zero-flux (dark+bias) images ; ; OPTIONAL INPUT PARAMETERS: ; -------- behavior of the routine ------------ ; silent : 1 if the routine is to perform silently. Default = 0, ; print out statements ; median : 1 if using median averaging, 0 if using avgclip algorithm ; ; -------- description of the input file ------ ; imbad : if set, an array of bad image indices (eg ; imbad=[0] to drop 1st image ; gain : electron/dn, used for calculating expected noise. Default = ; estimate from scatter. ; inbadmat : array of bad pixels for the flat (good for masking stars) ; zero : array or scalar for zero flux (bias + dark) ; zerr : array or scalar for error in zero ; zbad : array of bad pixels (0 = good) ; znoise : noise per pixek in a zero-flux image (eg readnoise for ; negligible dark), in DN ; exptime : if set, the exposure time in seconds, otherwise ; exposure time is taken from the SPE header (needed for ; triggered zero-flux files with PhotonMax systems) ; ; -------- output file ------------ ; caldir : data directory for cal files (with terminal '/'). Default = 'cal/' ; overrideoutfilename : if set (1), use outfilename as output file name. ; Default = 0, construct the output filename from the input file name ; outputfilename : if overrideoutfilename is set, this is the output ; filename used ; nowrite : if set, don't write the output file ; exptime : if set, use passed value for EXPTIME keyword ; ; OUTPUTS: ; avg - the average zero-flux image ; err - the error in the average zero-flux image (NOT error per ; image per pixel) ; bad - bad pixel image for avg (0 = good) ; ; OPTIONAL OUTPUTS: ; d - the input file (after deselecting with imbad) ; ; SIDE EFFECTS: ; For each file, writes ; cal/_.fits ; ; For each combination of gain, adc, and binning, write ; /_zero.fits ; with ; primary HDU - ; - INFILE - input file used to make average ; - OBJECT - 'median image of zero-flux files' ; - NIMAGES - number of images used to create average ; - EXPTIME - Exposure time in seconds ; - READN_DN - Derived noise per pixel, read noise + dark noise (DN) ; - AVG_METH - Method for making the average ; - header converted from original SPE header ; extention 1 - zero-flux error DN ; extention 2 - bad pixel map (0=good) ; ; EXAMPLE: ; ; MODIFICATION HISTORY: ; Written 14 Jun 2007 Leslie Ann Young SwRI ; Modified 13 Sep 2007 LAY ; Modified 18 Sep 2007 LAY added exptime keyword ; Modified 20 Sep 2007 LAY print input file ; Modified 21 Sep 2007 LAY median only the middle 1000 frames ; Modified 22 Sep 2007 CBO fixed two minor bugs from median'ing the middle 1000 frames ;- pro make_spe_flat1,filename, avg, err, bad, d = d, $ silent=silent, median = median, $ imbad = imbad, gain=gain, $ zfilename = zfilename, $ inbadmat = inbadmat, $ caldir=caldir, outfilename = outfilename, $ overrideoutfilename = overrideoutfilename, $ nowrite = nowrite, $ exptime = exptime if not keyword_set(silent) then silent = 0b if not keyword_set(caldir) then caldir = 'cal/' ;-------------------------------- ; Read the file using filename ;-------------------------------- fullfn = filename[0] if not silent then begin print, 'make_spe_flat1: averaging ', fullfn endif fn = file_basename(fullfn) len = strlen(fn) if strmatch(fn,'*.gz', /fold) then begin basefn = strmid(fn,0,len-7) endif else begin basefn = strmid(fn,0,len-4) endelse read_princeton_gz,fullfn,d,head=h ; n_elements rather than keyword_set in case ; user passes exptime = 0 if n_elements(exptime) eq 0 then begin exptime = h.exp_sec endif if not silent then begin print, 'EXPTIME = ', float(exptime) end nxnyn, d, nx, ny, nimages ;-------------------------------- ; Filter out the bad images ;-------------------------------- if n_elements(imbad) gt 0 then begin if not silent then begin print, 'DESELECTING IMAGES ', imbad end isgood = bytarr(nimages) + 1b isgood[imbad] = 0 igood = where(isgood eq 1b, ngood) if ngood eq 0 then begin print, 'make_spe_flat1(', filename,'...): all images deselected' return endif else begin d = d[*,*, igood] endelse nxnyn, d, nx, ny, nimages end ;-------------------------------- ; Subtract the zero-flux ;-------------------------------- if keyword_set(zfilename) then begin if zfilename eq 'auto' then begin zfn = caldir + strtrim(string(round(exptime*1000)),2)+'_zero.fits' endif else begin zfn = zfilename endelse zero = get_spe_zero(zfn, zerr, zbad, znoise, $ silent=silent, ainfo=h, aexp=exptime) end else begin zfn = '' zero = fltarr(nx, ny) zerr = fltarr(nx, ny) zbad = bytarr(nx,ny) znoise = 0. endelse d = float(d) for i = 0, nimages-1 do begin d[*,*,i] = d[*,*,i] - zero endfor ;-------------------------------- ; Set gain ; If the gain is not specified, estimate it from the ; first image ;-------------------------------- if not keyword_set(znoise) then znoise = 0. if keyword_set(gain) then begin g = gain endif else begin if nimages eq 1 then begin d0 = float(d[*,*,0]) robomean,d0,5.,0.1,avg,avgdev,stddev,var,skew,kurt,nfinal err = stddev endif else begin d0 = float(d[*,*,1]) - float(d[*,*,0]) robomean,d0,5.,0.1,avg,avgdev,stddev,var,skew,kurt,nfinal ; get noise err = stddev / sqrt(2.) d0 = float(d[*,*,1]) + float(d[*,*,0]) robomean,d0,5.,0.1,avg,avgdev,stddev,var,skew,kurt,nfinal ; get level endelse g = avg/(err^2 - znoise^2) endelse if not silent then begin print, 'GAIN (1st pass) = ', g, 'e/DN' end ;-------------------------------- ; Get the flat field frame ; ;-------------------------------- if keyword_set(median) then method = 'median' else method = 'avgclipsim' means = fltarr( nimages ) for i=0,nimages-1 do begin skysclim,d[*,*,i],lowval,hival,meanval,sigma,npts=20000 means[i] = meanval d[*,*,i]=d[*,*,i]/means[i] if not silent then begin if i lt 20 then print,'Frame ',i,' scaled by ',means[i] endif endfor normfac = median(means) case nimages of ;-------------------------------- ; 1 image: ; avg = that image ; err = passed gain, or estimated gain from robomean ; bad = 1 where points are more than 5 sigma from a 10x10 median-smoothed ; version of the average ;-------------------------------- 1: begin thresh = 5.0 avg = d[*,*,0] err = replicate( sqrt(means[0]/g), nx, ny) ; error in the flat smavg = median(avg, 10) bad = bytarr(nx,ny) if keyword_set(inbadmat) then begin bad = inbadmat[*,*,0] endif indxbad = where( abs(avg - smavg)/err gt thresh, nbad) if nbad gt 0 then bad[indxbad] = 1b end ;-------------------------------- ; 2 images: ; avg = mean of the two images ; err = passed gain, or estimated gain from robomean ; bad = 1 where the difference is > 5 sigma (5 * rn * sqrt(2)) ;-------------------------------- 2: begin thresh = 5. avg = (d[*,*,0] + d[*,*,1])/2. err = replicate(sqrt((means[0]+means[1])/g), nx, ny)/2. ; error in the zero bad = bytarr(nx,ny) if keyword_set(inbadmat) then begin d0 = d[*,*,0] d1 = d[*,*,1] b0 = inbadmat[*,*,0] b1 = inbadmat[*,*,1] ; the points that are bad in frame 0 and good in frame 1 bad0 = where(b0 eq 1 and b1 eq 0, nbad0) if nbad0 gt 0 then begin avg[bad0] = d1[bad0] err[bad0] = sqrt(2.) * err[bad0] endif ; the points that are good in frame 0 and bad in frame 1 bad1 = where(b0 eq 0 and b1 eq 1, nbad1) if nbad1 gt 0 then begin avg[bad1] = d0[bad1] err[bad1] = sqrt(2.) * err[bad1] endif ; the points that are bad in both bad01 = where(b0 eq 1 and b1 eq 1, nbad1, compl=good01, ncompl=ngood01) if ngood01 eq 0 then begin print,' Serious Warning! all pixels in imge cube were flagged bad.' bad[*,*] = 1 endif else begin if nbad1 gt 0 then begin avg[bad1] = median(avg[good01]) bad[bad01] = 1 endif endelse endif indxbad = where (abs(d[*,*,0] - d[*,*,1])/(sqrt(2.)*rn) gt thresh, nbad) if nbad gt 0 then bad[indxbad] = 1b end else: begin case method of ;-------------------------------- ; >=3 images, median: ; avg = median ; err = std deviation of mean sqrt( total(resid^2)/( (n-1)*n) ) ; bad = 1 where error is 5 sigma larger than other errors ;-------------------------------- 'median': begin ; call the (slow) medarr_mwb rather ; than risk malloc error if nimages gt 1000 then begin medarr_mwb, d, avg endif else begin avg = median(d, dim=3, /even) endelse ; create the error array res = d for i = 0, nimages-1 do res[*,*,i] = d[*,*,i]-avg ; err1 is the error in one pixel (i.e., the readnoise) err1 = fltarr(nx, ny) for i=0, nx-1 do begin for j = 0, ny-1 do begin err1[i,j] = total(res[i,j,*]^2) endfor endfor err1 = sqrt( err1/(nimages-1) ) ; error in a single pixel err = err1/sqrt(nimages) ; error in the average robomean,err1,5.,0.1,avgerr,avgerrdev,stddeverr gnew = (avgerr^2 - (znoise/normfac)^2)/normfac if gnew ge 0. then g = gnew if not keyword_set(silent) then begin print, 'GAIN = ', g, 'e/DN' endif ; create the bad-pixel array bad = bytarr(nx, ny) ibad = where(err gt (avgerr + 5.0 * stddeverr), nbad) if nbad gt 0 then bad[ibad] = 1 end ; median ;-------------------------------- ; >=3 images, median: ; avg = average, clipping 3-sigma outliers ; err = std deviation of mean excluding marked outliers ; bad = 1 where all pixels are outliers ;-------------------------------- 'avgclipsim': begin thresh = 3.0 ; 2007 Sept 21: malloc errors when ; median'ing 0.5-Gb file if nimages gt 1000 then begin first = (nimages / 2 - 500) > 0 last = (first + 999) < (nimages - 1) avg = median(d[*,*, first:last], dim=3, /even) endif else begin avg = median(d, dim=3, /even) endelse ; approximate noise sigma = sqrt( ( (avg>0)*normfac/g) + znoise^2 + zerr^2 ) $ / normfac if keyword_set(inbadmat) then begin badmat = inbadmat endif else begin badmat = bytarr(nx,ny,nimages) endelse asum = fltarr(nx,ny) acnt = fltarr(nx,ny) npts = long(nx) * long(ny) npixclipped = 0 for i=0,nimages-1 do begin new = d[*,*,i] newavg = avg * median(new)/median(avg) resid = (new - newavg)/sigma z = where(abs(resid) ge thresh,count) if count ne 0 then begin badi = badmat[*,*,i] badi[z] = 1B badmat[*,*,i] = badi endif z = where(badmat[*,*,i] eq 0B,count) if count ne 0 then begin asum[z] = asum[z]+new[z] acnt[z] = acnt[z] + 1.0 endif else begin print,' Serious Warning! all pixels in image ',strn(i), $ ' were considered bad.' ;stop endelse npixclipped = npixclipped + (npts - count) endfor ; refine the average where it's possible to z = where(acnt ne 0,count) if count ne 0 then avg[z]= asum[z]/acnt[z] ; create the bad pixel array bad = (acnt le 1) ; create the array with the error in avg err1 = sigma for i = 0, nx-1 do begin for j=0, ny-1 do begin n = acnt[i,j] if n ge 2 then begin err1[i,j] = $ sqrt( total( $ (d[i,j,*] - avg[i,j])^2 * (1.-badmat[i,j,*])) / $ (n-1.) ) endif endfor endfor err = err1 / sqrt(nimages) robomean,err1,5.,0.1,avgerr,avgerrdev,stddeverr ; gnew = sqrt(avgerr^2 - (znoise/normfac)^2)/normfac ; gnew = (1/normfac)/(avgerr^2 - (znoise/normfac)^2 ; - (median(zerr)/normfac)^2 ) gnew = normfac/( (avgerr*normfac)^2 - znoise^2 - median(zerr)^2 ) if gnew ge 0. then g = gnew if not keyword_set(silent) then begin print, 'GAIN = ', g, 'e/DN' endif end ; avgclipsim endcase end end if not keyword_set(nowrite) then begin ;-------------------------------- ; Get or create the output filename ;-------------------------------- if keyword_set(overrideoutfilename) then begin flatfn = outfilename endif else begin flatfn = caldir + basefn+'_flat.fits' outfilename = flatfn endelse if not silent then print, flatfn ;--------------------------------------------- ; write it out ;--------------------------------------------- mkhdr, header, avg, /EXTEND header = spe_head2fits(h, infits = header) sxaddpar, header, 'INFILE', fullfn, before='FIRSTSPE' sxaddpar, header, 'ZEROFILE', zfn, before='FIRSTSPE' sxaddpar, header, 'OBJECT', $ 'median image of flat-field files', before='FIRSTSPE' sxaddpar, header, $ 'NIMAGES', nimages, 'Number of images in original file', $ before='FIRSTSPE' sxaddpar, header, $ 'EXPTIME', exptime, 'Exposure time in seconds', $ before='FIRSTSPE' sxaddpar, header, $ 'NORMFAC', normfac, 'Normalizing factor' sxaddpar, header, $ 'READN_DN', znoise, $ 'Derived noise per pixel, read noise + dark noise (DN)', $ before='FIRSTSPE' sxaddpar, header, $ 'GAIN', g, $ 'Derived gain, electrons per DN', $ before='FIRSTSPE' sxaddpar, header, $ 'AVG_METH', method, 'Method for making the average', $ before='FIRSTSPE' WRITEFITS, flatfn, avg, header mkhdr, header, err, /EXTEND sxaddpar, header, 'OBJECT', 'error in median image' WRITEFITS, flatfn, err, /append mkhdr, header, bad, /EXTEND sxaddpar, header, 'OBJECT', 'bad pixel list (0=good)' WRITEFITS, flatfn, bad, /append endif end