pro lrisautoproc, datadir=datadir, prepare=prepare, makeflat=makeflat, flatten=flatten, makefringe=makefringe, rmfringe=rmfringe, crclean=crclean, astrometry=astrometry, stack=stack, bluefirst=bluefirst, badamp=badamp, oldfileprefix=oldfileprefix, chip=chip, restart=restart, continue=continue, oldred=oldred, only=only, winok=winok, noxfix=noxfix, nofringe=nofringe, setzeal=setzeal ; Automatic processing pipeline for LRIS imaging data. All steps are automated. ; Instructions: Place all raw files in the current directory with processed flats, then run lrisautoproc (no arguments) in IDL. ; If raw files are not in working directory, specify with datadir='(directory)' ; To not do an end-to-end reduction, can do steps exclusively using the flags above. ; set a flag to 1 (or use /keyword) to do that step exclusively ; set a flag to 0 to do skip that step (i.e., do all steps *except* this one) ; set a flag to -1 to do all steps *up to* that step (stops after specified step). ; set a flag to 2 to do all steps *starting with* that step. ; any other value processes as normal ; Order of events is: ; Prepare - split multi-extension header into left and right sides, bias-subtract, crop, fix some bad pixels and columns. ; Adds _l or _r suffix. ; Flatten - Apply the flat field. Adds 'f' prefix. ; Cosmic ray cleaning - Remove (zap) cosmic rays. Adds 'z' prefix. ; Solve astrometry - Solve for astrometry. Adds 'a' prefix. ; Stack - Combine (using swarp) images of the same field. autoastrocommand = 'python ~/progs/python/autoastrometry/autoastrometry.py' if n_elements(datadir) eq 0 then datadir = '' if n_elements(bluefirst) eq 0 then bluefirst = 0 if n_elements(prepare) eq 0 then prepare = 3 if n_elements(makeflat) eq 0 then makeflat = 3 if n_elements(flatten) eq 0 then flatten = 3 if n_elements(makefringe) eq 0 then makefringe = 3 if n_elements(rmfringe) eq 0 then rmfringe = 3 if n_elements(crclean) eq 0 then crclean = 3 if n_elements(astrometry) eq 0 then astrometry = 3 if n_elements(stack) eq 0 then stack = 3 if keyword_set(nofringe) then begin makefringe = 0 rmfringe = 0 endif if crclean ne 0 then flag=1 if n_elements(chip) eq 0 then chip = 'both' chip = strmid(strlowcase(chip),0,1) ; Look for an example file and find the date to infer the camera and file format ; (MJD upgrade dates are approximate) try = (findfile(datadir+'r??????_????'+'.fits')) [0] if try[0] eq '' then try = (findfile(datadir+'b??????_????'+'.fits')) [0] if try[0] eq '' then try = (findfile(datadir+'*r??????_????'+'.fits')) [0] if try[0] eq '' then try = (findfile(datadir+'*b??????_????'+'.fits')) [0] if try[0] eq '' then try = (findfile(datadir+'lred????'+'.fits')) [0] if try[0] eq '' then try = (findfile(datadir+'lblue????'+'.fits')) [0] if try[0] eq '' then try = (findfile(datadir+'*lred????'+'.fits')) [0] if try[0] eq '' then try = (findfile(datadir+'*lblue????'+'.fits')) [0] if try[0] ne '' then begin test = mrdfits(try[0],0,h,/silent) testmjd = float(sxpar(h,'MJD-OBS')) if testmjd gt 0 and testmjd lt 54952 then lrisversion = 1 if testmjd ge 54952 and testmjd lt 55562 then lrisversion = 2 if testmjd ge 55562 then lrisversion = 3 if lrisversion eq 1 then oldred = 1 if lrisversion gt 1 then oldred = 0 print, 'Instrument version is LRIS',clip(lrisversion) endif if n_elements(lrisversion) eq 0 then begin print, 'Cannot recognize LRIS version / file format' lrisversion = -1 oldred = keyword_set(oldred) endif if chip ne 'l' and chip ne 'r' and chip ne 'b' then print, 'Chip ', chip, ' not recognized: use L, R, or both' overwrite = 0 ; current default is to not overwrite. if n_elements(restart) eq 1 then overwrite = 1 if n_elements(continue) eq 1 then overwrite = 0 if n_elements(restart) eq 1 and n_elements(continue) eq 1 then begin print, 'Cannot both restart and continue.' return endif flatfail = [''] catastrofail = [''] relastrofail = [''] fullastrofail = [''] for c = 0, 1 do begin ; camera loop (blue/red) if keyword_set(only) and c eq 1 then continue ; --------------------------- ; Configure variables for appropriate camera print print, '------------------------' if c+bluefirst eq 1 then begin print, 'Now processing BLUE SIDE' camera = 'blue' instname = 'LRISBLUE' filtkey = 'BLUFILT' ;pxscale = 0.1353 prefchar = 'b' if keyword_set(oldfileprefix) then prefchar = 'lblue' ;rcrop=[0,0,1704,-1] rcrop=[0,0,0,0] endif if c-bluefirst eq 0 then begin camera = 'red' instname = 'LRIS' filtkey = 'REDFILT' prefchar = 'r' if keyword_set(oldfileprefix) then prefchar = 'lred' if keyword_set(oldred) eq 0 then begin print, 'Now processing RED SIDE' ;pxscale = 0.135 rcrop=[3,33,-5,2400] if keyword_set(badamp) then rcrop[2] = 1015 endif else begin print, 'Now processing OLD RED SIDE' ocrop=[150,10,1660,2020] if keyword_set(winok) then ocrop[1]=0 if keyword_set(winok) then ocrop[3]=0 endelse endif wildchar = '??????_????' if keyword_set(oldfileprefix) then wildchar = '????' ;if prepare eq 1 or prepare eq 2 then goto, preparestep (this is the first step anyway) if makeflat eq 1 or makeflat eq 2 then goto, makeflatstep if flatten eq 1 or flatten eq 2 then goto, flattenstep if makefringe eq 1 or makefringe eq 2 then goto, makefringestep if rmfringe eq 1 or rmfringe eq 2 then goto, fringecorrectstep if crclean eq 1 or crclean eq 2 then goto, crcleanstep if astrometry eq 1 or astrometry eq 2 then goto, astrostep if stack eq 1 or stack eq 2 then goto, stackstep ; ----------------------------------- ; Prepare images for processing (debias, crop, fix) preparestep: files = findfile(datadir+prefchar+wildchar+'.fits') if chip eq 'b' then begin ; b means both files_l = findfile(prefchar+wildchar+'_l.fits') files_r = findfile(prefchar+wildchar+'_r.fits') endif else begin files_chip = findfile(prefchar+wildchar+'_'+chip+'.fits') endelse if oldred then files_o = findfile(prefchar+wildchar+'_o.fits') if prepare ne 0 then begin if (oldred eq 0 or camera eq 'blue') and n_elements(files) gt 1 then begin openr, 1, camera+'gain.dat' inline = '' readf, 1, inline print, 'Using gain data from ', camera+'gain.dat' close, 1 gain = float(strsplit(inline,/extract)) endif for f = 0, n_elements(files)-1 do begin if files[f] eq '' then continue filearr = (strsplit((strsplit(files[f],'.',/extract)) [0], '/', /extract)) fileroot = filearr[n_elements(filearr)-1] if chip ne 'b' then begin if oldred eq 0 or camera eq 'blue' then begin match = where(fileroot+'_'+chip+'.fits' eq files_chip, ct) ; check if file exists endif else begin match = where(fileroot+'_o.fits' eq files_o, ct) ; check if right file exists endelse endif else begin if oldred eq 0 or camera eq 'blue' then begin match = where(fileroot+'_r.fits' eq files_r, ctr) ; check if right file exists match = where(fileroot+'_l.fits' eq files_l, ctl) ; check if left file exists ct = ctr and ctl endif else begin match = where(fileroot+'_o.fits' eq files_o, ct) ; check if right file exists endelse endelse if ct eq 0 or overwrite then begin ; don't redo and overwrite the file, unless restart flag is set h = headfits(files[f], /silent) slitname = strtrim(sxpar(h, 'SLITNAME'),2) pane = strtrim(sxpar(h, 'PANE'),2) window = strtrim(sxpar(h, 'WINDOW'),2) if slitname ne 'direct' then begin print, files[f], ' is not imaging: skipping' continue endif if camera eq 'red' and ((pane ne '400,625,3400,2500' and pane ne '0') or $ (window ne '0,0,0,2048,2048' and window ne '0')) and $ keyword_set(winok) eq 0 then begin print, files[f], ' is improperly windowed: skipping' continue endif print, 'Preparing ', files[f] if oldred eq 0 then begin if lrisversion eq 3 then begin if camera eq 'red' and chip eq 'r' then splitlrisred, files[f], gain=gain, /rightonly if camera eq 'red' and chip eq 'l' then splitlrisred, files[f], gain=gain, /rightonly if camera eq 'red' and chip eq 'b' then splitlrisred, files[f], gain=gain, /rightonly endif else begin if camera eq 'red' and chip eq 'r' then splitlrisred, files[f], gain=gain, /rightonly, flag=flag, /namefix if camera eq 'red' and chip eq 'l' then splitlrisred, files[f], gain=gain, /leftonly, flag=flag, /namefix if camera eq 'red' and chip eq 'b' then splitlrisred, files[f], gain=gain, flag=flag, /namefix endelse if camera eq 'blue' and chip eq 'r' then splitlrisblue, files[f], gain=gain, /rightonly, /namefix if camera eq 'blue' and chip eq 'l' then splitlrisblue, files[f], gain=gain, /leftonly, /namefix if camera eq 'blue' and chip eq 'b' then splitlrisblue, files[f], gain=gain, /namefix endif else begin if camera eq 'red' then lrisred1prepare, files[f], prefix='', suffix='_o', /namefix if camera eq 'blue' and chip eq 'b' then splitlrisblue1, files[f], gain=gain, lcropx=[500,-1], rcropx=[1,-350], cropy=[1000,-890], /namefix if camera eq 'blue' and chip eq 'r' then splitlrisblue1, files[f], gain=gain, lcropx=[500,-1], rcropx=[1,-350], cropy=[1000,-890], /rightonly, /namefix if camera eq 'blue' and chip eq 'l' then splitlrisblue1, files[f], gain=gain, lcropx=[500,-1], rcropx=[1,-350], cropy=[1000,-890], /leftonly, /namefix endelse endif endfor endif if prepare eq -1 or prepare eq 1 then continue ; ----------------------------------------- ; Create the flatfields makeflatstep: for p = 0, 1 do begin if chip eq 'b' then begin if p eq 0 then fchip = 'r' if p eq 1 then fchip = 'l' endif else begin fchip = chip if p eq 1 then continue endelse if oldred and camera eq 'red' and p eq 1 then continue if oldred and camera eq 'red' then fchip = 'o' ;loop continues... files_r = findfile(prefchar+wildchar+'_'+fchip+'.fits') isdomeflat = intarr(n_elements(files_r)) isskyflat = intarr(n_elements(files_r)) isscience = intarr(n_elements(files_r)) filters = strarr(n_elements(files_r)) dichs = strarr(n_elements(files_r)) wins = strarr(n_elements(files_r)) counts = fltarr(n_elements(files_r)) targets = strarr(n_elements(files_r)) exps = fltarr(n_elements(files_r)) for f = 0, n_elements(files_r)-1 do begin if files_r[f] eq '' then continue fileroot = (strsplit(files_r[f],'.',/extract)) [0] h = headfits(files_r[f], /silent) filters[f] = clip(sxpar(h, filtkey)) counts[f] = sxpar(h,'COUNTS') dichs[f] = clip(sxpar(h,'DICHNAME')) wins[f] = clip(sxpar(h,'WINDOW')) if clip(sxpar(h,'PANE')) ne '0' then wins[f] = clip(sxpar(h,'PANE')) targets[f] = clip(sxpar(h,'TARGNAME')) trapdoor = clip(sxpar(h,'TRAPDOOR')) exptime = sxpar(h, 'ELAPTIME') exps[f] = exptime azimuth = sxpar(h,'AZ') elevation = sxpar(h,'EL') domeazimuth = sxpar(h,'DOMEPOSN') domerelaz = domeazimuth - azimuth while domerelaz lt -180.0 do domerelaz += 360. while domerelaz ge 180.0 do domerelaz -= 360. if trapdoor eq 'open' then begin ;abs(domerelaz-90) lt 1 or abs(domerelaz+90) lt 1)) if (abs(elevation-45) lt 0.5 and abs(domerelaz) gt 3) or strpos(targets[f],'cass') ge 0 then begin isdomeflat[f] = 1 ;if counts[f] gt 7000 and counts[f] lt 40000 then endif else begin if ((counts[f] gt 9000 and camera eq 'blue') or (counts[f] gt 12000 and camera eq 'red')) and counts[f] lt 43000 and exps[f] lt 60 then isskyflat[f] = 1 if abs(domerelaz) lt 5.0 and ((counts[f] lt 40000) and (counts[f] le 15000 or (filters[f] eq 'I' and counts[f] le 30000) or (filters[f] eq 'RG850') or (exps[f] ge 60))) then isscience[f] = 1 endelse endif ;print, files_r[f], ' ', clip(targets[f],12), string(exps[f],format='(I4)'), ' ', clip(filters[f],5), ' ', counts[f], domerelaz, isdomeflat[f], isskyflat[f], isscience[f] ;printzxx, files_r[f], counts[f], isdomeflat[f], isskyflat[f], isscience[f], elevation, domeazimuth, azimuth, domerelaz endfor sciindex = where(isscience, ct) if ct eq 0 then continue scifilterlist = unique(filters[sciindex]) for i = 0, n_elements(scifilterlist)-1 do begin filt = scifilterlist[i] dichlist = unique(dichs[sciindex[where(filt eq filters[sciindex])]]) for d = 0, n_elements(dichlist)-1 do begin dich = dichlist[d] if n_elements(dichlist) gt 1 then dichstr = dichlist[d] else dichstr = '' outflatname = 'lris'+strmid(camera,0,1)+'flat'+filt+dichstr+'_'+fchip+'.fits' if file_test(outflatname) and overwrite eq 0 then continue ; flat exists already flatsuccess = 0 ;print, filt, ' ', dich ;for f = 0, n_elements(filtskyflatlist)-1 do print, ' ', filtskyflatlist[f], ' (', strtrim(string(counts[(where(files_r eq filtskyflatlist[f])) [0]]),2), ' counts)' ; Good supersky flat. ctsupersky = 0 if (oldred eq 1 and camera eq 'red' and (filt eq 'I' or filt eq 'RG850' or filt eq 'GG570')) ne 1 then begin ; don't supersky fringed frames ; First, see if we have enough frames to do a supersky flat. poss = where(filters eq filt and dichs eq dich and isdomeflat eq 0 and counts gt 600 and exps gt 10., ctposs) if ctposs gt 0 then begin ;otherwise usually means some weird situation like the dichroic info is missing, etc. medcounts = median(counts[poss]) superskymaxcount = (medcounts * 5) < (30000. > medcounts*1.2) ; superskyminexp = 10 superskyflats = where(filters eq filt and dichs eq dich and isdomeflat eq 0 and counts gt 600 and counts gt medcounts/3. and counts lt superskymaxcount and exps ge superskyminexp, ctsupersky) if file_test('badflattargets.txt') eq 0 then begin badflattargets = [''] endif else begin ;badflattargets = grabcolumn('badflattargets.txt', 0, /str) nbf = countlines ('badflattargets.txt') openr, 6, 'badflattargets.txt' badflattargets = strarr(nbf) readf, 6, badflattargets close, 6 endelse print, badflattargets print, superskyflats for b = 0, n_elements(badflattargets)-1 do begin w = where(targets[superskyflats] ne badflattargets[b], ctsupersky) superskyflats = superskyflats[w] print, superskyflats endfor if ctsupersky gt 1 then superskyfields = unique(targets[superskyflats]) else superskyfields = [''] nsuperskyfields = n_elements(superskyfields) ; this criterion is a bit loose and depends on the number of very bright stars in the image and how big the dither steps were. ; eventually, should make flatcombine able to find bright (ultra-saturated) objects and add halo masks, and make ; lrisautoproc treat big dithers like new fields. ;print, ctsupersky, nsuperskyfields ;if (ctsupersky ge 9 and nsuperskyfields ge 5) or (ctsupersky ge 17 and nsuperskyfields ge 4) or (ctsupersky gt 17 and nsuperskyfields gt 3) or ctsupersky ge 21 then begin if (ctsupersky ge 21 and nsuperskyfields ge 5) or (ctsupersky ge 25 and nsuperskyfields ge 4) or (ctsupersky gt 29 and nsuperskyfields gt 3) or ctsupersky ge 33 then begin print, filt, '-band / ', 'D'+dich,': ', strtrim(string(ctsupersky),2), ' supersky flats (', clip(nsuperskyfields), ' unique fields)' flatcombine, files_r[superskyflats], outflatname, /removeobjects, type='supersky' flatsuccess = 1 endif endif endif ; Regular sky flat. if flatsuccess eq 0 then begin skyflats = where(isskyflat and filters eq filt and dichs eq dich, ctsky) if ctsky ge 5 then begin print, filt, '-band / ', 'D'+dich,': ', strtrim(string(ctsky),2), ' sky flats.' flatcombine, files_r[skyflats], outflatname, /removeobjects, type='sky' flatsuccess = 1 endif endif ; Dome flat. if flatsuccess eq 0 then begin domeflats = where(isdomeflat and filters eq filt and dichs eq dich and counts gt 7000. and counts lt 39000., ctdome) ;print, ctdome if ctdome ge 3 then begin print, filt, '-band / ', 'D'+dich,': ', strtrim(string(ctdome),2), ' dome flats' flatcombine, files_r[domeflats], outflatname, type='dome' flatsuccess = 1 endif endif ; marginal sky flat if flatsuccess eq 0 then begin if ctsky ge 3 then begin print, filt, '-band / ', 'D'+dich,': ', strtrim(string(ctsky),2), ' sky flats (WARNING - marginal!)' flatcombine, files_r[skyflats], outflatname, /removeobjects, type='sky' flatsuccess = 1 endif endif ; marginal supersky flat if flatsuccess eq 0 then begin if ctsupersky ge 3 then begin print, filt, '-band / ', 'D'+dich,': ', strtrim(string(ctsupersky),2), ' supersky flats (WARNING - EXTREMELY marginal)' print, 'Flat field for this filter/dichroic combination is likely to be extremely poor.' flatcombine, files_r[superskyflats], outflatname, /removeobjects, type='supersky' flatsuccess = 1 endif endif if flatsuccess eq 0 then begin print, 'Unable to produce a flat field for this filter-dichroic: ', filt, '-band / ', 'D'+dich+'.' noproc = where(isscience and filters eq filt and dichs eq dich, ctnoproc) print, 'Will not be able to further process ', clip(ctnoproc), ' images without a flat from another source:' for j = 0, ctnoproc-1 do print, ' ', files_r[noproc[j]], exps[noproc[j]], ' ', targets[noproc[j]] endif endfor endfor endfor if makeflat eq -1 or makeflat eq 1 then continue ; ----------------------------------------- ; Divide images by flatfield flattenstep: ; Search for existing flatfields noflats = 0 if flatten ne 0 then begin flats = findfile('*flat*.fits') flatfilts = strarr(n_elements(flats)) flatchips = strarr(n_elements(flats)) flatdichs = strarr(n_elements(flats)) flatwins = strarr(n_elements(flats)) if flats[0] ne '' then begin for f = 0, n_elements(flats)-1 do begin h = headfits(flats[f], /silent) inst = strtrim(sxpar(h, 'INSTRUME'),2) if inst ne instname then continue filter = clip(sxpar(h, filtkey)) dich = clip(sxpar(h, 'DICHNAME')) win = clip(sxpar(h,'WINDOW')) flatchip = clip(sxpar(h, 'CHIP')) flatfilts[f] = filter flatdichs[f] = dich flatwins[f] = strmid(win,2) flatchips[f] = flatchip endfor endif else begin print, 'No flats found! Cannot perform flatfielding.' noflats = 1 ; skip to next step - maybe the early intermediate files were deleted endelse files_r = findfile(prefchar+wildchar+'_'+'?'+'.fits') ffiles_r = findfile('f'+prefchar+wildchar+'_'+'?'+'.fits') if noflats eq 0 then begin for f = 0, n_elements(files_r)-1 do begin if files_r[f] eq '' then continue fileroot = (strsplit(files_r[f],'.',/extract)) [0] match = where('f'+fileroot+'.fits' eq ffiles_r, ct) ; check if output file exists if ct eq 0 or overwrite then begin ; don't redo step unles restart flag is set h = headfits(files_r[f], /silent) ;;fits = mrdfits(files_r[f], 0, h, /silent) counts = sxpar(h, 'SKYCTS') > sxpar(h,'COUNTS') exptime = sxpar(h, 'ELAPTIME') azimuth = sxpar(h,'AZ') elevation = sxpar(h,'EL') domeazimuth = sxpar(h,'DOMEPOSN') target = sxpar(h,'TARGNAME') domerelaz = domeazimuth - azimuth while domerelaz lt -180.0 do domerelaz += 360. while domerelaz ge 180.0 do domerelaz -= 360. ;print, files_r[f], elevation, domerelaz, target, counts, exptime, counts/exptime ; (abs(domerelaz-90) lt 15 or abs(domerelaz+90) lt 15) if abs(elevation-45) lt 0.5 and abs(domerelaz) gt 3 then continue ; silently skip dome flats if strpos(target,'cass') ge 0 then continue if strpos(target,'dome') ge 0 then continue ; sometimes the telescope is pointed wrong; we still don't want these... if strpos(target,'horizon') ge 0 then continue if counts gt 40000 then continue ; skip anything close to saturation if counts gt 30000. and exptime lt 10 then continue if counts/exptime gt 10000. and exptime gt 1 then begin print, files_r[f], ' is probably a sky flat: skipping' continue endif filter = clip(sxpar(h, filtkey)) dich = clip(sxpar(h, 'DICHNAME')) win = strmid(clip(sxpar(h,'WINDOW')),2) flatchip = clip(sxpar(h, 'CHIP')) flatfileno = where(flatfilts eq filter and flatchips eq flatchip and flatdichs eq dich and flatwins eq win, ct) ;print, flatfilts, '(', filter,')' ;print, flatchips,'(',flatchip ,')' ;print, flatdichs, '(', dich,')' ;print, flatwins, '(', win,')' if ct eq 0 then begin print, 'Flat field not found for ', files_r[f], ' (', strtrim(filter,2), '-band / ', strtrim(dich), ' / ', win, ')' flatfail = [flatfail, files_r[f]] continue endif flatfile = flats[flatfileno[0]] print, 'Flattening ', files_r[f], ' using ', flatfile if flatchip eq 'right' then flatproc, files_r[f], flatfile, crop=rcrop if flatchip eq 'left' then flatproc, files_r[f], flatfile, crop=lcrop if flatchip eq 'old' then flatproc, files_r[f], flatfile, crop=ocrop if flatchip eq 'old' and keyword_set(noxfix) eq 0 then lrisxfix, 'f'+files_r[f], prefix='' if flatchip ne 'right' and flatchip ne 'left' and flatchip ne 'old' then print, 'Chip not recognized - cannot flatten' endif endfor endif endif if flatten eq -1 or flatten eq 1 then continue makefringestep: if oldred and makefringe ne 0 and camera eq 'red' then begin ffiles_r = findfile('f'+prefchar+wildchar+'_'+'?'+'.fits') fileinfo = replicate({filename:'', filter:'', counts:0., exptime:0., airmass:0., object:''}, n_elements(ffiles_r)) fileinfo.filename = ffiles_r for f = 0, n_elements(fileinfo)-1 do begin if fileinfo[f].filename eq '' then continue h = headfits(fileinfo[f].filename, /silent) fileinfo[f].filter = strtrim(sxpar(h, filtkey),2) fileinfo[f].counts = sxpar(h, 'COUNTS') fileinfo[f].exptime = sxpar(h, 'ELAPTIME') fileinfo[f].object = sxpar(h, 'OBJECT') fileinfo[f].airmass = sxpar(h, 'AIRMASS') ;print, fileinfo[f].filename, ' ', fileinfo[f].filter, fileinfo[f].counts, fileinfo[f].exptime endfor fringablefilts = ['I','RG850'] fringeexcludeobj = [''];['GRB070412', '050412', 'GRB051109B', 'GRB070521', 'GRB100205A'] for ff = 0, n_elements(fringablefilts)-1 do begin ffilt = fringablefilts[ff] fringefile = ffilt + 'fringe.fits' if file_test(fringefile) then continue thisfilt = where(fileinfo.filter eq ffilt, ct) if ct lt 3 then continue ffileinfo = fileinfo[thisfilt] thisdark = where(ffileinfo.exptime gt 60 and ffileinfo.counts lt 40000, ct) if ct lt 3 then continue ffileinfo = ffileinfo[thisdark] ctratethreshold = median((ffileinfo.counts/ffileinfo.exptime)/ffileinfo.airmass) * 1.8 thisdark = where((ffileinfo.counts/ffileinfo.exptime)/ffileinfo.airmass le ctratethreshold) if ct lt 3 then continue ffileinfo = ffileinfo[thisdark] for nf = 0, n_elements(fringeexcludeobj)-1 do begin notforbiddenobj = where(ffileinfo.object ne fringeexcludeobj[nf], ct) if ct ge 1 then ffileinfo = ffileinfo[notforbiddenobj] endfor if ct lt 3 then continue print, 'Making ', ffilt+'-band fringe frame...' print, ffileinfo.filename mkfringe, ffileinfo.filename, outfile=fringefile, /removeobjects endfor endif if makefringe eq -1 or makefringe eq 1 then continue fringecorrectstep: if oldred and rmfringe ne 0 and camera eq 'red' then begin ffiles_r = findfile('f'+prefchar+wildchar+'_'+'?'+'.fits') ffiles_r = unmatched(ffiles_r,'i') fileinfo = replicate({filename:'', filter:'', counts:0., exptime:0., object:''}, n_elements(ffiles_r)) fileinfo.filename = ffiles_r for f = 0, n_elements(fileinfo)-1 do begin if fileinfo[f].filename eq '' then continue h = headfits(fileinfo[f].filename, /silent) fileinfo[f].filter = strtrim(sxpar(h, filtkey),2) fileinfo[f].counts = sxpar(h, 'COUNTS') fileinfo[f].exptime = sxpar(h, 'ELAPTIME') fileinfo[f].object = sxpar(h, 'OBJECT') endfor fringablefilts = ['I','RG850'] for ff = 0, n_elements(fringablefilts)-1 do begin ffilt = fringablefilts[ff] fringefile = ffilt + 'fringe.fits' if file_test(fringefile) eq 0 then continue thisfringe = where(fileinfo.filter eq ffilt and fileinfo.counts/fileinfo.exptime lt 1000., ctfringe) if ctfringe ge 1 then begin fringablefiles = fileinfo[thisfringe].filename if fringablefiles[0] ne '' then $ fringecorrect, fringablefiles, fringefile=fringefile, prefix='i';, /iterate, itersec=[50,880,200,-200] ; alternative method uses weights rather than iteration ;fringescale = filesi.counts ;continuumimg = where(filesi.counts/(filesi.exptime*filesi.airmass) gt 10, ctcontinuum) ;if ctcontinuum gt 0 then fringescale[continuumimg] = (filesi[continuumimg].exptime*filesi[continuumimg].airmass)/0.15 ;fringecorrect, filesi.name, fringefile='Ifringe.fits', /iterate, itersec=[300,700,400,600] endif endfor endif if rmfringe eq -1 or rmfringe eq 1 then continue ; ----------------------------------------------------- ; Clean cosmic rays from the image (using pzap). crcleanstep: if crclean ne 0 then begin ffiles_r = choosefiles(prefchar+wildchar+'_'+'?'+'.fits','if','f') if overwrite eq 0 then ffiles_r = unmatched(ffiles_r,'z') ;zffiles_r = findfile('zf'+prefchar+wildchar+'_'+'?'+'.fits') for f = 0, n_elements(ffiles_r)-1 do begin if ffiles_r[f] eq '' then continue h = headfits(ffiles_r[f]) counts = sxpar(h,'COUNTS') exptime = sxpar(h,'ELAPTIME') target = sxpar(h,'TARGNAME') if counts gt 40000. then continue if counts gt 30000. and exptime lt 10. then continue if counts gt 20000. and exptime gt 5. and exptime lt 10. then continue if target eq '' then continue ; blank target is certainly not a science object. if strpos(target,'flat') ge 0 then continue if strpos(target,'sky') ge 0 then continue if strpos(target,'twilight') ge 0 then continue print, 'Cleaning cosmic rays from ', ffiles_r[f] zeal = 0.85 if camera eq 'blue' then zeal = 0.75 if (camera eq 'red' and oldred) then usamp=1 else usamp=0 if camera eq 'red' and oldred eq 0 then zeal = 0.6 ;evidently should be the new default if n_elements(setzeal) eq 1 then zeal = setzeal print, zeal, usamp pzap, ffiles_r[f], /weight, zeal=zeal, usamp=usamp, /quiet endfor endif if crclean eq -1 or crclean eq 1 then continue ; ----------------------------------------------------- ; Do astrometry on images using autoastrometry.py astrostep: if astrometry ne 0 then begin for p = 0, 1 do begin if chip eq 'b' then begin if p eq 0 then achip = 'r' if p eq 1 then achip = 'l' endif else begin achip = chip if p eq 1 then continue endelse if oldred and camera eq 'red' and p eq 1 then continue if oldred and camera eq 'red' then achip = 'o' ;loop continues... if crclean eq 0 then begin zffiles_r = choosefiles(prefchar+wildchar+'_'+achip+'.fits','if','f') ; zffiles_r = findfile('f'+prefchar+wildchar+'_'+achip+'.fits') ; azffiles_r = findfile('af'+prefchar+wildchar+'_'+achip+'.fits') endif else begin zffiles_r = choosefiles(prefchar+wildchar+'_'+achip+'.fits','zif','zf') ; zffiles_r = findfile('zf'+prefchar+wildchar+'_'+achip+'.fits') ; azffiles_r = findfile('azf'+prefchar+wildchar+'_'+achip+'.fits') endelse filetargets = strarr(n_elements(zffiles_r)) fileexposures = strarr(n_elements(zffiles_r)) filecounts = fltarr(n_elements(zffiles_r)) for f = 0, n_elements(zffiles_r)-1 do begin if zffiles_r[f] eq '' then continue h = headfits(zffiles_r[f], /silent) ; print, 'autoastrometry '+zffiles_r[f]+ ' -o' ; Initially, get a "rough" solution for every image using the catalog ; spawn, 'autoastrometry '+zffiles_r[f]+ ' -o' ; (for deep images to be used in stacking we will improve this later) ; ; Overwrites initial file ; filetargets[f] = repstr(repstr(strtrim(sxpar(h,'TARGNAME'),2),' ', '_'),'/','_') ; spaces cause barfing in filenames fileexposures[f] = string(sxpar(h,'ELAPTIME')) filecounts[f] = sxpar(h,'COUNTS') endfor ; Make a reference catalog using a representative image out of an image block (several images of the same field) targets = unique(filetargets) for t = 0, n_elements(targets)-1 do begin targetfiles = where(zffiles_r eq targets[t]) refcatfile = targets[t] + '.'+camera+'_'+achip+'.ref.cat' if file_test(refcatfile) and overwrite eq 0 then continue thistarget = where(filetargets eq targets[t]) maxexp = max(fileexposures(thistarget)) minctrate = min(filecounts[thistarget]/fileexposures[thistarget]) if maxexp lt 5 or minctrate gt 5000 then begin print, targets[t], ' is a shallow field (standard or sky): not making a reference catalog.' continue endif imagesthistarg = zffiles_r[thistarget] ;for i = 0, n_elements(imagesthistarg)-1 do begin ; h = headfits(imagesthistarg[i]) ;endfor refimagename = imagesthistarg[n_elements(imagesthistarg)/2] ;h = headfits(refimagename, /silent) ;pa = strtrim(string((float(sxpar(h, 'ROTPOSN'))+0) MOD 360),2) print print, 'Making reference catalog for ', targets[t], ' using ', refimagename print, autoastrocommand+' '+refimagename+' -upa 2 -q' spawn, autoastrocommand+' '+refimagename+' -upa 2 -q' if file_test('a'+refimagename) eq 0 then begin catastrofail = [catastrofail, refimagename] print, 'WARNING - astrometry on the reference image was unsuccessful!' print endif else begin print, autoastrocommand+' '+'a'+refimagename+' -n '+refcatfile + ' -x 55000 -q' spawn, autoastrocommand+' '+'a'+refimagename+' -n '+refcatfile + ' -x 55000 -q' print endelse ; ' -px '+strtrim(string(pxscale),2)+' -pa '+pa+ ;;+' -s 6' ; (respecifying px and pa should no longer be necessary with the new splitlrisred and splitlris blue which add astrometry.) endfor zffiles_r = unmatched(zffiles_r,'a') ; catalogs should always be the same; only do this now ; Use the reference catalog to do a more precise relative astrometric solution for f = 0, n_elements(zffiles_r)-1 do begin if zffiles_r[f] eq '' then continue if file_test('a'+zffiles_r[f]) and overwrite eq 0 then continue fileroot = (strsplit(zffiles_r[f],'.',/extract)) [0] h = headfits(zffiles_r[f], /silent) pa = strtrim(string((float(sxpar(h, 'ROTPOSN'))+0) MOD 360),2) exptime = sxpar(h,'ELAPTIME') counts = sxpar(h,'COUNTS') if counts/exptime gt 5000. or exptime lt 5 then begin ; only do catalog astrometry for short exposures (no catalog). if counts gt 30000 and exptime lt 10 then begin print, zffiles_r[f], ' is in bright twilight; not solving astrometry.' continue endif print, zffiles_r[f], ' is a twilight/standard frame, solving astrometry directly against a catalog.' print, autoastrocommand+' '+zffiles_r[f]+' -upa 2 -q' spawn, autoastrocommand+' '+zffiles_r[f]+' -upa 2 -q' print if file_test('a'+zffiles_r[f]) eq 0 then fullastrofail = [fullastrofail, zffiles_r[f]] endif else begin ; use the short exposure first, but fall back on direct catalog. targname = repstr(strtrim(sxpar(h,'TARGNAME'),2),' ', '_') ; spaces cause barfing in filenames if targname eq '' then continue if strpos(targname,'flat') ge 0 then continue refcatfile = targname + '.'+camera+'_'+achip+'.ref.cat' if file_test(refcatfile) then begin print, targname print, autoastrocommand+' '+zffiles_r[f]+' -c '+refcatfile+' -upa 2' + ' -x 55000 -q' ;-upa 0.1 spawn, autoastrocommand+' '+zffiles_r[f]+' -c '+refcatfile+' -upa 2' + ' -x 55000 -q' ;-upa 0.1 endif else begin print, 'No reference catalog '+refcatfile+' exists for this field.' endelse if file_test('a'+zffiles_r[f]) eq 0 then begin print, 'Refined astrometry of ', zffiles_r[f], ' was not successful. Trying direct astrometry:' print, autoastrocommand+' '+zffiles_r[f]+' -upa 2 -q' spawn, autoastrocommand+' '+zffiles_r[f]+' -upa 2 -q' print if file_test('a'+zffiles_r[f]) then relastrofail = [relastrofail, zffiles_r[f]] $ else fullastrofail = [fullastrofail, zffiles_r[f]] endif endelse ; ' -px '+strtrim(string(pxscale),2)+' -pa '+pa+ ;;+' -s 6' endfor endfor endif if astrometry eq -1 or astrometry eq 1 then continue ; ----------------------------------------------------- ; Stack images together to produce final image. stackstep: if file_test('default.swarp') eq 0 then $ spawn, 'swarp -d > default.swarp' print print print for p = 0, 1 do begin if chip eq 'b' then begin if p eq 0 then schip = 'r' if p eq 1 then schip = 'l' endif else begin schip = chip if p eq 1 then continue endelse if oldred and camera eq 'red' and p eq 1 then continue if oldred and camera eq 'red' then schip = 'o' ;loop continues... ;schip = '?' this is a hack to combine both sides ;if crclean ne 0 then azffiles_r = findfile('azf'+prefchar+wildchar+'_'+schip+'.fits')$ ; else azffiles_r = findfile( 'af'+prefchar+wildchar+'_'+schip+'.fits') azffiles_r = findfile( 'a*'+prefchar+wildchar+'_'+schip+'.fits') realfiles = where(azffiles_r ne '', ct) if ct eq 0 then continue ; can't stack if there are no astrometry files if ct ge 1 then azffiles_r = azffiles_r[realfiles] ;schip = 'lr' ; this is a hack to combine both sides camver = camera if camver eq 'red' then camver = camver + clip(lrisversion) if file_test('autophotsummaryflux.'+camver+'.txt') then begin readcol,'autophotsummaryflux.'+camver+'.txt',pfile,pexp,pfilt,pdich,pair,dum,pdmag,pfluxratio,pseeing,format='a,i,a,a,f,a,f,f,f,f',/silent photodata = replicate({filename:'',dmag:0.,fluxratio:0.,seeing:0.},n_elements(pfile)) photodata.filename = pfile photodata.dmag = pdmag photodata.fluxratio = pfluxratio photodata.seeing = pseeing endif filetargets = strarr(n_elements(azffiles_r)) fileexposures = fltarr(n_elements(azffiles_r)) filefilters = strarr(n_elements(azffiles_r)) filesatval = fltarr(n_elements(azffiles_r)) fileskyval = fltarr(n_elements(azffiles_r)) fileairval = fltarr(n_elements(azffiles_r)) fileseeingpix = fltarr(n_elements(azffiles_r)) filefluxratio = fltarr(n_elements(azffiles_r)) datestr = '' for f = 0, n_elements(azffiles_r)-1 do begin h = headfits(azffiles_r[f], /silent) if f eq 0 then datestr = sxpar(h,'DATE') filetargets[f] = repstr(strtrim(sxpar(h,'TARGNAME'),2),' ', '_') ; spaces cause barfing in filenames fileexposures[f] = float(sxpar(h,'ELAPTIME')) filefilters[f] = string(sxpar(h,filtkey)) filesatval[f] = sxpar(h,'SATURATE') fileskyval[f] = sxpar(h,'COUNTS') fileairval[f] = sxpar(h,'AIRMASS') if n_elements(photodata) gt 0 then begin pmatch = where(azffiles_r[f] eq photodata.filename, ct) if ct gt 0 then begin fileseeingpix[f] = photodata[pmatch].seeing filefluxratio[f] = photodata[pmatch].fluxratio endif endif endfor targets = unique(filetargets) for t = 0, n_elements(targets)-1 do begin target = targets[t] thistarget = where(filetargets eq target, cttarg) if cttarg eq 0 then continue thistargetfilts = unique(filefilters[thistarget]) for l = 0, n_elements(thistargetfilts)-1 do begin filter = thistargetfilts[l] thistargfilt = where(filetargets eq target and filefilters eq filter, cttargfilt) if cttargfilt eq 0 then continue if max(fileexposures[thistargfilt]) lt 30 then continue expthresh = (median(fileexposures[thistargfilt]) * 0.8) > (max(fileexposures[thistargfilt]) * 0.25) stacki = where(filetargets eq target and filefilters eq filter and fileexposures gt expthresh, ctstack) ;if ctstack lt 2 then continue if ctstack eq 0 then continue if ctstack eq 1 then print, 'Only one deep exposure: trivial stack' stacklist = azffiles_r[stacki] stackexps = fileexposures[stacki] medianexp = median(stackexps) medair = median(fileairval[stacki]) minair = min(fileairval[stacki]) maxair = max(fileairval[stacki]) minseeingpix = min(fileseeingpix[stacki]) maxseeingpix = max(fileseeingpix[stacki]) medseeingpix = median(fileseeingpix[stacki]) minfluxratio = min(filefluxratio[stacki]) maxfluxratio = max(filefluxratio[stacki]) medfluxratio = median(filefluxratio[stacki]) totalexp = total(stackexps) nstack = n_elements(stacklist) if nstack gt 1 then stdfluxratio = stdev(filefluxratio[stacki]) else stdfluxratio = 0 outfile = 'coadd' + strtrim(target,2) +'_'+ strtrim(filter,2) +'_'+ schip + '.fits' outweightfile = 'coadd' + strtrim(target,2) +'_'+ strtrim(filter,2) +'_'+ schip + '.weight.fits' stackcmd = 'swarp ' for s = 0, n_elements(stacklist)-1 do begin if s eq 0 then stackcmd = stackcmd + stacklist[s] if s gt 0 then stackcmd = stackcmd + ',' + stacklist[s] endfor for s = 0, n_elements(stacklist)-1 do begin weightfilename = strmid(stacklist[s],0,strlen(stacklist[s])-5-0) + '.weight.fits' weightexists = file_test(weightfilename) if weightexists eq 0 then begin weightfilenameinit = strmid(weightfilename,1) weightexists = file_test(weightfilenameinit) if weightexists then begin print, 'mv '+weightfilenameinit+' '+weightfilename spawn, 'mv '+weightfilenameinit+' '+weightfilename endif endif endfor stackcmd = stackcmd + ' -IMAGEOUT_NAME ' + outfile + ' -WEIGHTOUT_NAME ' + outweightfile if nstack gt 1 then begin ; used to be 3... stackcmd = stackcmd + ' -WEIGHT_TYPE MAP_WEIGHT' endif else begin stackcmd = stackcmd + ' -WEIGHT_TYPE BACKGROUND' endelse ;stackcmd = stackcmd + ' -VERBOSE_TYPE QUIET' stackcmd = stackcmd + ' -FSCALE_DEFAULT ' for s = 0, n_elements(stacklist)-1 do stackcmd = stackcmd + strtrim(string(medianexp/stackexps[s]),2) + ',' ; multiplicative factor... stackcmd = stackcmd + ' -COPY_KEYWORDS OBJECT,TARGNAME,TELESCOP,FILTER,'+filtkey+',DICHNAME,'+$ 'SLITNAME,GRISNAME,GRANAME,GRANGLE,INSTRUME,UT,UTC,MJD_OBS,MJD-OBS,ST,'+$ 'DATE,DATE-OBS,AIRMASS,AZ,RA,DEC,EL,HA,TELFOCUS,ROTPOSN,DOMEPOSN,CURRINST,OBSERVER,FLATTYPE';+$ ;'DOMEPOSN,CURRINST,ADCWAVE0,ADCWAVE1,AMPMODE,CCDGAIN,CCDSPEED,NUMAMPS' ; 'DATE_BEG,DATE,END,UTC-END,DATE-END,BINNING,TRAPDOOR,LAMPS,FLIMAGIN,FLSPECTR' ; Too many keywords causes pyfits to be unable to load the mosaic. It's rather freakish. ;stackcmd = stackcmd + ' -BACK_SIZE 256' if (file_test(outfile) eq 0) or overwrite then begin if nstack eq 1 then $ ; used to be 3 for some bizarre reason? print, 'Warning - only ', clip(nstack), ' exposures; not flagging bad pixels.' print, 'Stacking ', target, ':' for s = 0, n_elements(stacklist)-1 do print, stacklist[s], ' ', strtrim(string(stackexps[s]),2) + 's' print, stackcmd spawn, stackcmd ; do the stack data = mrdfits(outfile, 0, h, /silent) sxaddpar, h, 'DATE', datestr sxaddpar, h, 'NSTACK', nstack sxaddpar, h, 'AIRMASS', medair, 'Median exposure airmass' sxaddpar, h, 'AIRMIN', minair, 'Minimum exposure airmass' sxaddpar, h, 'AIRMAX', maxair, 'Maximum exposure airmass' sxaddpar, h, 'EXPOSURE', medianexp, 'Effective rescaled exposure time' sxaddpar, h, 'TOTALEXP', totalexp, 'Total summed integration time' sxaddpar, h, 'MAXEXP', max(stackexps), 'Length of longest exposure' sxaddpar, h, 'MINEXP', min(stackexps), 'Length of shortest exposure' sxaddpar, h, 'SATURATE', min(filesatval[stacki]-fileskyval[stacki]) sxaddpar, h, 'MEDSKY', median(fileskyval[stacki], /even) if minseeingpix gt 0 then begin sxaddpar, h, 'SEEING', medseeingpix, 'Median seeing in pixels' sxaddpar, h, 'SEEMIN', minseeingpix sxaddpar, h, 'SEEMAX', maxseeingpix if camera eq 'red' and oldred then medseeingarcsec = medseeingpix*0.210 else medseeingarcsec = medseeingpix*0.135 sxaddpar, h, 'SEEARCSC', medseeingarcsec, 'Median seeing in arcsec' endif if minfluxratio gt 0 then begin sxaddpar, h, 'TRANSRAT', maxfluxratio/minfluxratio, 'Transmission ratio of max/min images' sxaddpar, h, 'TRANSSTD', stdfluxratio, 'Std. dev. of relative transmission' endif for s = 0, n_elements(stacklist)-1 do sxaddpar, h, 'IMAGE'+strtrim(string(s),2), stacklist[s] for s = 0, n_elements(stacklist)-1 do sxaddpar, h, 'IMEXP'+strtrim(string(s),2), stackexps[s] sxaddhist, ['Keck GRB host program','Beta release - use with caution'], h, /COMMENT mwrfits, data, outfile, h, /create ; add this header info print, autoastrocommand+' '+outfile+' -q' spawn, autoastrocommand+' '+outfile+' -q' ; do another round of precise astrometry on the stack (probably unnecessary now) endif endfor endfor endfor qq = 1 if qq eq 0 then begin; chip eq 'b' and (camera ne 'red' and oldred ne 1) then begin ; Combine left and right sides (or multiple exposure blocks, if necessary) for COADDS coaddfilesr = findfile('*coadd*_r.fits') coaddfilesl = findfile('*coadd*_l.fits') if coaddfilesr ne [''] then coaddfilesr = coaddfilesr[where(coaddfilesr ne '')] else continue if coaddfilesl ne [''] then coaddfilesl = coaddfilesl[where(coaddfilesl ne '')] else continue ; match every r with an l and coadd for f = 0, n_elements(coaddfilesr)-1 do begin rfilename = coaddfilesr[f] lfilename = repstr(rfilename, '_r.','_l.') match = where(coaddfilesl eq lfilename, ct) if ct eq 1 then begin stackcmd = 'swarp ' + rfilename + ',' + lfilename outfile = repstr(strmid(coaddfilesr,6), '_r.', '.') if file_test(outfile) eq 0 or overwrite then begin rweightname = repstr(rfilename, '.fits','.weight.fits') lweightname = repstr(lfilename, '.fits','.weight.fits') if file_test(rweightname) eq 0 then spawn, 'mv ' + strmid(rweightname,1) + ' ' + rweightname if file_test(lweightname) eq 0 then spawn, 'mv ' + strmid(lweightname,1) + ' ' + lweightname print print, 'Combining ', rfilename, ' and ', lfilename stackcmd = stackcmd + ' -IMAGEOUT_NAME ' + outfile stackcmd = stackcmd + ' -WEIGHT_TYPE MAP_WEIGHT' stackcmd = stackcmd + ' -COPY_KEYWORDS OBJECT,TARGNAME,TELESCOP,FILTER,'+filtkey+',DICHNAME,'+$ 'SLITNAME,GRISNAME,GRANAME,GRANGLE,INSTRUME,UTC,MJD_OBS,ST,'+$ 'DATE-OBS,AIRMASS,AZ,RA,DEC,EL,HA,TELFOCUS,ELAPTIME,'+$ 'DOMEPOSN,CURRINST,ADCWAVE0,ADCWAVE1,AMPMODE,CCDGAIN,CCDSPEED,NUMAMPS,'+$ 'EXPOSURE,TOTALEXP,MINEXP,MAXEXP' print, stackcmd spawn, stackcmd data = mrdfits(outfile, 0, h) hr = headfits(rfilename) hl = headfits(lfilename) for i = 0, 100 do begin rimagekeyn = sxpar(hr,'IMAGE'+clip(i)) if clip(rimagekeyn) ne '0' then sxaddpar, h, 'RIMAGE'+clip(i), clip(rimagekeyn) endfor for i = 0, 100 do begin limagekeyn = sxpar(hl,'IMAGE'+clip(i)) if clip(limagekeyn) ne '0' then sxaddpar, h, 'LIMAGE'+clip(i), clip(limagekeyn) endfor mwrfits, data, outfile, h, /create endif endif endfor ; For individual images (standards, found if exptime < 30 sec) if crclean eq 0 then begin azffiles_r = findfile('af'+prefchar+wildchar+'_r.fits') endif else begin azffiles_r = findfile('azf'+prefchar+wildchar+'_r.fits') endelse realfiles = where(azffiles_r ne '', ct) ; Combine l and r images. if ct gt 0 then begin for f = 0, n_elements(azffiles_r)-1 do begin infiler = azffiles_r[f] h = headfits(infiler) exptime = sxpar(h, 'ELAPTIME') filen = strmid(infiler,strpos(infiler,'_r.fits')-4,4) target = repstr(clip(sxpar(h,'TARGNAME')),' ', '_') if exptime le 30 then begin infilel = repstr(infiler,'_r.','_l.') outfile = target + '_' + filen + '_o.fits' if file_test(infilel) eq 0 then continue ; mate doesn't exist if file_test(outfile) eq 1 and overwrite eq 0 then continue ; already stacked these two print print, 'Preparing to combine ' + infiler + ' and ' + infilel print, autoastrocommand+' '+infiler+' -q' spawn, autoastrocommand+' '+infiler+' -q' print, autoastrocommand+' '+infilel+' -q' spawn, autoastrocommand+' '+infilel+' -q' stackcmd = 'swarp a' + infiler + ',a' + infilel print print, 'Combining ' + 'a'+infiler + ' and a' + infilel stackcmd = stackcmd + ' -IMAGEOUT_NAME ' + outfile ; no weighting for this one stackcmd = stackcmd + ' -VERBOSE_TYPE QUIET' ; fast, so hide this stackcmd = stackcmd + ' -COPY_KEYWORDS OBJECT,TARGNAME,TELESCOP,FILTER,'+filtkey+',DICHNAME,'+$ 'SLITNAME,GRISNAME,GRANAME,GRANGLE,INSTRUME,UTC,MJD_OBS,ST,'+$ 'DATE-OBS,AIRMASS,AZ,RA,DEC,EL,HA,TELFOCUS,ELAPTIME,'+$ 'DOMEPOSN,CURRINST,ADCWAVE0,ADCWAVE1,AMPMODE,CCDGAIN,CCDSPEED,NUMAMPS,' print, stackcmd spawn, stackcmd print, ' -> ', outfile print endif endfor endif endif endfor ; camera loop print print, flatfail nflatfail = n_elements(flatfail)-1 if nflatfail gt 1 then begin ;element 0 is blank print, 'Unable to flat-field the following images:' for f = 1, n_elements(flatfail)-1 do begin h = headfits(flatfail[f]) print, clip(flatfail[f],25), string(sxpar(h,'ELAPTIME'),format='(I5)'), ' ', clip(sxpar(h,'FILTER'),5), ' ', clip(sxpar(h,'DICHNAME'),7), ' ',clip(sxpar(h,'TARGNAME'),16), ' ', clip(sxpar(h,'WINDOW')), ' ', clip(sxpar(h,'PANE')) endfor endif print nafail = n_elements(relastrofail)-1 + n_elements(fullastrofail)-1 + n_elements(catastrofail)-1 if nafail gt 1 then begin ;element 0 is blank if n_elements(catastrofail) gt 1 then print, 'Unable to produce astrometric catalogs for the following reference images:' for f = 1, n_elements(catastrofail)-1 do begin h = headfits(catastrofail[f]) print, clip(catastrofail[f],25), string(sxpar(h,'ELAPTIME'),format='(I5)'), ' ', clip(sxpar(h,'FILTER'),4), ' ', clip(sxpar(h,'TARGNAME'),16), ' ', clip(sxpar(h,'COUNTS')) endfor if n_elements(relastrofail) gt 1 then print, 'Relative astrometry failed for the following images, but absolute was successful:' for f = 1, n_elements(relastrofail)-1 do begin h = headfits(relastrofail[f]) print, clip(relastrofail[f],25), string(sxpar(h,'ELAPTIME'),format='(I5)'), ' ', clip(sxpar(h,'FILTER'),4), ' ', clip(sxpar(h,'TARGNAME'),16), ' ', clip(sxpar(h,'COUNTS')) endfor if n_elements(fullastrofail) gt 1 then print, 'All astrometry failed for the following images (not stacked):' for f = 1, n_elements(fullastrofail)-1 do begin h = headfits(fullastrofail[f]) print, clip(fullastrofail[f],25), string(sxpar(h,'ELAPTIME'),format='(I5)'), ' ', clip(sxpar(h,'FILTER'),4), ' ', clip(sxpar(h,'TARGNAME'),16), ' ', clip(sxpar(h,'COUNTS')) endfor endif; else begin ;print, 'Finished with no astrometric problems.' ;endelse end