;+ ; NAME: ; comb_spe_zero ; ; PURPOSE: (one line) ; combine spe_zero files ; ; DESCRIPTION: ; combine spe_zero files ; ; CATEGORY: ; SPE ; ; CALLING SEQUENCE: ; comb_spe_zero, zfilenames, zero, zerr, zbad, znoise, $ ; silent=silent, $ ; aexptime=aexptime, ainfo=ainfo ; INPUTS: ; zfilenames : filenames of SPE zero-flux (dark+bias) average images ; ; OPTIONAL INPUT PARAMETERS: ; -------- behavior of the routine ------------ ; silent : 1 if the routine is to perform silently. Default = 0, ; print out statements ; ; -------- description of the file the zero file will be applied to ------ ; aexptime : exposure time in seconds ; ainfo: string (filename), array of strings (FITS header) or ; structure (SPE header) ; ; OUTPUTS: ; zero - the average zero-flux image ; zerr - the error in the average zero-flux image (NOT error per ; image per pixel) ; zbad - bad pixel image for avg (0 = good) ; znoise - read noise (or read noise plus dark noise) scalar, in DN ; ; SIDE EFFECTS: ; ; EXAMPLE: ; ; MODIFICATION HISTORY: ; Written 14 Jun 2007 Leslie Ann Young SwRI ; Modified 13 Sep 2007 LAY ;- pro comb_spe_zero,zfilenames, zero, zerr, zbad, znoise, $ silent=silent, aexptime=aexptime, ainfo=ainfo, $ exptime=exptime, $ method=method, bias=bias, dark=dark, $ caldir=caldir, $ outfilename = outfilename, $ overridefilename = overridefilename, $ nowrite = nowrite n = n_elements(zfilenames) if n eq 1 then zfilenames = reform(zfilenames, 1) ; read the first one and make the larger arrays i = 0 fn = zfilenames[i] a = get_spe_zero(fn, e, b, fitshead=h, $ silent=silent, aexptime=aexptime, ainfo=ainfo) nxnyn, a, nx0, ny0 aarr = fltarr(nx0,ny0,n) earr = fltarr(nx0,ny0,n) barr = fltarr(nx0,ny0,n) harr = strarr( n_elements(h), n) narr = fltarr(n) ; read noise per pixel exparr = fltarr(n) if n eq 1 then begin aarr = reform(aarr, nx0,ny0,n) earr = reform(earr, nx0,ny0,n) barr = reform(barr, nx0,ny0,n) harr = reform(harr, n_elements(h), n) endif for i = 0, n-1 do begin fn = zfilenames[i] if not isfile(fn) then stop a = get_spe_zero(fn, e, b, fitsh=h, exp=exp, silent=silent) nxnyn, a, nx, ny if nx eq nx0 and ny eq ny0 then begin bd = where(e le 0., nbd) if nbd gt 0 then begin e[bd] = median(e) b[bd] = 1 endif harr[*,i] = h narr[i] = sxpar(harr[*,i],'READN_DN') aarr[*,*,i] = a barr[*,*,i] = b earr[*,*,i] = e exparr[i] = exp endif else begin print, 'comb_spe_zero : ', fn, ' has different dimensions than ', zfilenames[0] endelse end if keyword_set(aexptime) then begin exptime = aexptime endif else begin valid = 0 if keyword_set(ainfo) then begin head = spe_head(ainfo, valid, speh, exptime=exptime) endif if valid eq 0 then begin exptime = mean(exparr) endif endelse ; if not keyword_set(aexptime) and not keyword_set(ainfo) then begin ; exptime = mean(exparr) ; endif else begin ; head = spe_head(ainfo, valid, speh, exptime=aexptime) ; exptime = aexptime ; endelse if not keyword_set(method) then method = 'average' case method of 'average': begin zero = total(aarr,3)/n zerr = sqrt(total( reform(earr^2,nx,ny,n), 3)/n) zbad = ( total(barr,3) gt 0 ) znoise = mean(narr) end 'weight_by_nimage': begin nimages = lonarr(n) for i = 0, n-1 do nimages[i] = sxpar(harr[*,i],'NIMAGES') wght = fltarr(nx,ny,n) for i = 0, n-1 do wght[*,*,i] = nimages[i] / total(nimages) zero = total( reform(aarr * wght,nx,ny,n) ,3) zerr = sqrt(total( reform(wght^2 * earr^2,nx,ny,n) , 3) ) zbad = ( total(barr,3) gt 0 ) znoise = total( nimages * narr) / total(nimages) end 'weight_by_error': begin wght = 1./earr^2 wghttot = total(wght, 3) for i = 0, n-1 do wght[*,*,i] = wght[*,*,i] / wghttot zero = total( reform(aarr * wght,nx,ny,n) ,3) zerr = sqrt(total( reform(wght^2 * earr^2,nx,ny,n), 3) ) zbad = ( total(barr,3) gt 0 ) wgt_per_image = fltarr(n) for i=0,n-1 do wgt_per_image[i] = median(wght[*,*,i]) znoise = total( reform(wgt_per_image * narr,nx,ny,n) ) / total(wgt_per_image) end 'bias+dark': begin i = 0 w = 1/earr[*,*,i]^2. s = w sx = exparr[i] * w sy = aarr[*,*,i] * w sxx = exparr[i]^2. * w sxy = exparr[i] * aarr[*,*,i] * w for i = 1, n-1 do begin w = 1/e^2. s = s + w sx = sx + exparr[i] * w sy = sy + aarr[*,*,i] * w sxx = sxx + exparr[i]^2. * w sxy = sxy + exparr[i] * aarr[*,*,i] * w endfor del = s * sxx - sx^2 if min(del) gt 0. then begin ; ----------------------- ; normal equations for weighted fit to a line ; ----------------------- bias = (sxx * sy - sx * sxy)/del dark = (s * sxy - sx * sy)/del biaserr = sqrt( sxx/del ) darkerr = sqrt( s/del ) cov = -sx/del endif else begin ; ----------------------- ; del = 0 if one file or all times are equal. ; assume all of the signal is due to the bias ; ----------------------- bias = sy/s dark = fltarr(nx,ny) biaserr = sqrt(1/s) darkerr = fltarr(nx,ny) cov = fltarr(nx,ny) endelse zero = bias + dark * exptime zerr = sqrt( biaserr^2 + exptime^2 * dark^2 + exptime * cov) zbad = ( total(barr,3) gt 0 ) znoise = poly(exptime,poly_fit(exparr, narr,1)) end endcase if not keyword_set(nowrite) then begin ;-------------------------------- ; Get or create the output filename ;-------------------------------- if not keyword_set(caldir) then caldir = 'cal/' if keyword_set(overridefilename) then begin zerofn = outfilename endif else begin zerofn = caldir + strtrim(string(round(exptime*1000)),2)+'_zero.fits' outfilename = zerofn endelse if not silent then print, zerofn ;--------------------------------------------- ; write it out ;--------------------------------------------- mkhdr, hdr, zero, /EXTEND iend = where( strmid(hdr,0,8) eq 'END ') header = spe_head(h, /justspe) ; the lines between FIRSTSPE and LASTSPE header = [ hdr[0:iend-1], header, hdr[iend] ] sxaddpar, header, 'OBJECT', $ 'median image of zero-flux files', before='FIRSTSPE' sxaddpar, header, $ 'NFILES', n, 'Number of files combined', $ before='FIRSTSPE' for i = 0, n-1 do begin keyword = 'INFILE' + strtrim(string(i),2) sxaddpar, header, keyword, zfilenames[i], before='FIRSTSPE' endfor sxaddpar, header, $ 'EXPTIME', exptime, 'Exposure time in seconds', $ before='FIRSTSPE' sxaddpar, header, $ 'READN_DN', znoise, $ 'Derived noise per pixel, read noise + dark noise (DN)', $ before='FIRSTSPE' sxaddpar, header, $ 'COMBMETH', method, 'Method for combining the zero-flux files', $ before='FIRSTSPE' WRITEFITS, zerofn, zero, header mkhdr, header, zerr, /EXTEND sxaddpar, header, 'OBJECT', 'error in median image' WRITEFITS, zerofn, zerr, /append mkhdr, header, zbad, /EXTEND sxaddpar, header, 'OBJECT', 'bad pixel list (0=good)' WRITEFITS, zerofn, zbad, /append endif end