; NAME: ; PROC4k ; ; PURPOSE: ; Does the basic reductions (xtalk,overscan,flat,fringe,wcs) for a list ; of MDM 4K images. ; ; PROCEDURE: ; 1) Corrects the pixels in each quadrant that are contaminated by ; crosstalk signals from the other three quadrants. This is done by ; subtracting the rotated source quadrant times a crosstalk coefficient ; for each quadrant in the image. A file filename.x.fits is created ; after this step. ; NOTE: Since the xtalk is proportional to the number of electrons ; in each well, it is unreliable to fix pixels affected by ; ADC-saturated (65535) sources. By default, these pixels are set to ; the saturation level, but this routine allows the user to ignore ; them or attempt a correction. The method should be chosen ; carefully by the science goals of the data. ; ; 2) Subracts the median value of the overscan for each row, flips ; image to be N up, E left, adds header keywords EPOCH=2000, ; SECPIX=0.315*CCDXBIN for every file in filelist and writes new ; file, filename.xo.fits. Skips and issues a warning if ; filename.xo.fits already exists. ; ; 3) Creates a normalized masterflat for each CCDXBIN and MISFLTID by ; scaling each by the sky value (calculated by MMM) and median ; combining all IMAGETYP=FLAT of a given bin/filter. Outputs ; masterflat[bin][filter].fits (eg for a 2x2 binned I band master ; flat, the file will be masterflat2I.fits. Skips and issues a ; warning if the masterflat already exists. Assumes square binning. ; ; 4) Applies the appropriate flat to all IMAGETYP=OBJECT ; frames. Writes to filename.xof.fits. Skips and issues a warning if ; if filename.xof.fits already exists or if a masterflat of the same ; bin/filter does not exist. ; ; 5) Applies a fringe correction to all IMAGETYP=OBJECT, MISFLTID=I ; frames using fringe[bin].fits in the current working directory. If ; this is not created beforehand, and FORCEFRINGE is set, it will ; retrieve a fringe image from ; http://www.astronomy.ohio-state.edu/~jdeast/4k/fringe[bin].fits. ; These images should be stable for months (?? - feedback ; appreciated on that). Outputs to filename.xoff.fits. Skips and ; issues a warning if there is no fringe[bin].fits and FORCEFRINGE ; is not set, or if filename.xoff.fits already exists. If anyone ; creates a more recent, higher SN fringe image, please send it to me. ; ; 6) Solves for the coordinate solution using SExtractor and ; WCSTools. Outputs to filename.xofw.fits and creates a DS9 region ; file to check it, which should circle the 50 brightest stars in the ; frame. Skips and issues a warning if filename.xofw.fits already exists. ; NOTE: You need the environment variable UB1_PATH to be defined and ; point to the location of the USNOb catalog for this to work. ; ; CALLING SEQUENCE: ; proc4k, path [,/noflip, /nowcs ,/forcefringe, xtcoefflo=, xtcoeffhi=,$ ; /xtalksat, /ignorehigh, /fixhigh, /mask] ; ; EXAMPLE: ; prock4k, path ; ; INPUTS: ; PATH: A string containing the path of files to process, a ; string array of files to process, or the name of file with a list ; of image to process. It will distinguish between the 3 scenerios ; like this: ; ; If PATH has more than one element, it assumes it's a string array ; If PATH ends in '.fits' or '.fits.fz', it assumes it's a path of ; one or more images Otherwise, PATH will be interpreted as the name ; of a file containing a list name image filenames, one per line. ; ; These images must have the correct IMAGETYP ,MISFLTID, and CCDXBIN ; keywords, as well as square binning for proper processing. There ; is no quality control, so exclude bad images (especially flats) ; from the list before hand. It does not currently handle images of ; different sizes; they should be put in separate filelists and ; processed individually. ; ; NOTE: The list should include all images, including raw flats, but ; not the products of this code (e.g., filename.xo.fits) or the ; master calibration frames (e.g., fringe1.fits, masterflat1I.fits). ; ; NOFLIP: Normally the program will flip every image in the input ; list to be oriented N up E left. Set this keyword to avoid ; this. Note this step is required for the WCS solution. WCS will ; be skipped if this keyword is set. ; ; NOWCS: Skip the World Coordinate Solution. Set this keyword if you ; don't have sextractor and WCSTools installed, or you don't want to ; wait for the solution (a few seconds to 15 minutes per image, ; depending on how good the pointing of the telescope was -- with ; the new control system, the pointing is generally accurate and the ; coordinate solution should be quick). ; ; FORCEFRINGE: If no fringe correction image is available and this ; keyword is set, a stock fringe correction from ; http://www.astronomy.ohio-state.edu/~jdeast/4k/fringe[bin].fits ; will be used. This fringe image was created in Sept 2007. Caution ; should be used if using this image long afterward. ; ; XTCOEFFLO: A 4x4 array specifying the xtalk coefficient, ; xtcoeff[i,j], of jth source quadrant on ith victim quadrant in the ; source regime ADU <= 65534. By default, the routine will use a set ; of crosstalk coefficients measured by Calen Henderson (OSU) in ; 2007 and 2009. ; NOTE: While these coefficients were stable between 2007 and 2009, ; if you want to calculate your own set based on your data, ; xtcoeffs.pro has been provided here: ; http://www.astronomy.ohio-state.edu/~jdeast/4k/xtcoeffs.pro ; ; XTALKSAT: Pixels affected by xtalk from saturated sources can be ; extremely unreliable. By default, they are set to XTALKSAT ; (default = 65535) so they won't be trusted for later processing. ; If /IGNOREHIGH or /FIXHIGH is specified, XTALKSAT will be ignored. ; NOTE: Subsquent processing will mean that these pixels will not ; remain 65535 in the final calibrated data, but should still fall ; above your saturation threshhold. Use MASK to be able to identify ; these pixels individually. ; ; MASK: If set, PROC4K will output filename.mask.fits, a byte image ; where 0s are pixels affected by xtalk from saturated sources and ; 1s are good. ; ; IGNOREHIGH: Rather than setting pixels affected by saturated ; sources to XTALKSAT, the user can set /IGNOREHIGH and apply only ; the low regime xtalk coefficient. Caution should be used when ; using pixels affected by saturated sources, as these pixels are ; guaranteed to be >= their true values, and may differ ; substantially from their true value. ; ; FIXHIGH: If set, pixels contaminated by saturated pixels will be ; corrected using XTCOEFFHIGH. Caution should be used if attempting ; to use these pixels, as these coefficients are an ensemble average ; of all saturated pixels in a given quadrant, which is field ; dependent, and not necessarily applicable to any given pixel. ; ; XTCOEFFHI: A 4x4 array specifying the xtalk coefficient, ; xtcoeffhi[i,j], of jth source quadrant on ith victim quadrant in the ; source regime ADU = 65535. This routine has a set of built-in ; coefficients derived by Calen Henderson from the Cygnus OB field ; that will be used if none are given, however, it is **HIGHLY ; RECOMMENDED** that you calculate your own coefficients based on ; your field, using xtcoeff.pro, if attempting to correct these pixels. ; ; DEPENDENCIES: ; IDL Astronomy library ; SExtractor (if WCS solution desired) ; WCSTools (if WCS solution desired) ; UB1_PATH environment variable set (if WCS solution desired) ; -- This should point to your local copy of the USNOb catalog ; -- If you don't have a copy of the catalog, use Harvard's (slow): ; setenv UB1_PATH http://tdc-www.harvard.edu/cgi-bin/scat ; ; LIMITATIONS: ; 1) Requires square bins ; 2) All images taken with different ROI settings should be processed ; separately. ; ; Modification History: ; Written by Jason Eastman OSU ; 2008/08/01 Fixed bug that failed when flats contained 0s ; (0s are set to 1) ; 2008/08/07 Fixed bug that didn't allow paths in the filelist ; 2009/04/07 Fixed bug that didn't properly take into account BZERO keyword ; 2009/04/09 Simplified DS9 region file to 50 stars labeled by magnitude ; Crop y overscan (no correction applied) ; 2009/04/17 Forced mode to be greater than 0 ; 2009/11/17 Fixed bug in unrotated NAXIS keywords ; use faster file_test instead of file.exist ; changed 'filelist' to 'path' to simplify calling ; procedure (backward compatible) ; changed output file structure to filebase.xofw.fits ; used dimension keyword in median to get rid of deprecated medarr ; Added xtalk correction (Calen Henderson, OSU) ; Replaced mode with MMM for sky background measurement ; Got rid of naxis1, naxis2, xover, yover arrays, now only ; carry them for each image ; Optionally output mask image, which identifies pixels ; affected by xtalk from saturated sources ; Use MMM to normalize flat instead of mean ; 2010/09/23 Changed ambigous keywords xtcoeff and xtcoeffhigh to ; xtcoefflo and xtcoeffhi ; Added /even keyword to all median calls (less biased) ; Changed all for loop indices to type long ; Works directly with fpacked images now ; 2011/04/13 Now works for OSMOS imaging data (uses "INSTRUME" header key) ; Will xtalk and overscan correct all images, only flat ; field images where SLITID = 'Open' and DISPID = 'Open' ; Fixed bug using "*.fits" as path ; 2011/05/02 Fixed bug in NAXIS1/NAXIS2 assignments for OSMOS (caught ; when writing the FITS file, this fix prevents the warning). pro proc4k, path, noflip=noflip, nowcs=nowcs, forcefringe=forcefringe, $ xtcoefflo=xtcoefflo, xtcoeffhi=xtcoeffhi, xtalksat=xtalksat, $ ignorehigh=ignorehigh, fixhigh=fixhigh, mask=mask if n_elements(path) gt 1 then begin ;; path is a string array of filenames files = path nfiles = n_elements(files) endif else if (strpos(path, '.fits') ne -1 and $ strpos(path,'.fits') eq strlen(path) - 5) or $ (strpos(path, '.fz') ne -1 and $ strpos(path,'.fz') eq strlen(path) - 3) then begin ;; path is a path of images files = file_search(path) if files[0] eq '' then message, 'ERROR: No files found' nfiles = n_elements(files) endif else begin ;; assumes a listname nfiles = file_lines(path) files = strarr(nfiles) openr, list, path, /get_lun for i=0L, nfiles-1 do begin filename = '' readf, list, filename files[i] = filename endfor free_lun, list endelse nfiles = n_elements(files) xbin = intarr(nfiles) filter = strarr(nfiles) imgtype = strarr(nfiles) filebase = strarr(nfiles) instrument = strarr(nfiles) ;; for OSMOS slitid = strarr(nfiles) + 'Open' dispid = strarr(nfiles) + 'Open' ;; Crosstalk coefficients, where first index is the victim frame and ;; second number is the source frame. ;; Low is the count regime x < 65534.0 if n_elements(xtcoefflo) ne 16 then begin xtcoefflo = [[0.000000000, 0.000512182, 0.00182688, 0.00123051],$ [0.000559781, 0.000000000, 0.00037473, 0.00112503],$ [0.001783810, 0.000417698, 0.00000000, 0.00174711],$ [0.001234770, 0.001176940, 0.00183285, 0.00000000]] endif ;; High is x = 65535.d0 -- this is really just a false sense of ;; security and should probably be removed if n_elements(xtcoeffhi) ne 16 then begin xtcoeffhi = [[0.000000000, 0.000665058, 0.00232165, 0.00172752],$ [0.000582549, 0.000000000, 0.00038356, 0.00119812],$ [0.001565830, 0.000380083, 0.00000000, 0.00157540],$ [0.000823900, 0.000669940, 0.00144859, 0.00000000]] endif if n_elements(xtalksat) eq 0 then xtalksat=65535 thresh = 65534.0 rot = [[0,5,7,2], $ [5,0,2,7], $ [7,2,0,5], $ [2,7,5,0]] quad = fltarr(2032,2032,4) quadxtfix = fltarr(2032,2032,4) for n=0L, nfiles-1 do begin filebase[n] = (strsplit(files[n],'.fits',/extract,/regex))(0) ;; get the path of the images if n eq 0 then begin filepath = strpos(filebase[0], '/', /REVERSE_SEARCH) if filepath eq -1 then filepath = '' $ else filepath = strmid(filebase[0],0,filepath+1) endif if file_test(filebase[n] + '.x.fits') then begin hdr = headfits(files[n]) print, 'Warning: ' + filebase[n] +'.x.fits' + $ ' exists. To reprocess, delete current file' endif else image = float(readfits(files[n], hdr, /silent)) ;; get the useful header entries xbin[n] = sxpar(hdr,'CCDXBIN') instrument[n] = strtrim(sxpar(hdr, 'INSTRUME'),2) if instrument[n] eq 'MDM4K' then begin filter[n] = strtrim(sxpar(hdr, 'MISFLTID'),2) endif else if instrument[n] eq 'OSMOS' then begin slitid[n] = strtrim(sxpar(hdr,'SLITID'),2) dispid[n] = strtrim(sxpar(hdr,'DISPID'),2) ;; use the first non-open filter as the filter ID if strtrim(sxpar(hdr,'FILTID1'),2) eq 'Open' then begin filter[n] = strtrim(sxpar(hdr,'FILTID2'),2) endif else begin filter[n] = strtrim(sxpar(hdr,'FILTID1'),2) endelse endif else message, 'ERROR: undefined instrument (' + instrument[n] + ')' imgtype[n] = strtrim(sxpar(hdr, 'IMAGETYP'),2) ;; correct pixels with cross talk if not file_test(filebase[n] + '.x.fits') then begin naxis1 = sxpar(hdr,'NAXIS1') naxis2 = sxpar(hdr,'NAXIS2') rawimage = image x1 = [0,naxis1/2,0,naxis1/2] x2 = [naxis1/2-1,naxis1-1,naxis1/2-1,naxis1-1] y1 = [naxis2/2,naxis2/2,0,0] y2 = [naxis2-1,naxis2-1,naxis2/2-1,naxis2/2-1] nbad = 0 if keyword_set(mask) then maskimg = bytarr(naxis1,naxis2) + 1 for i=0L,3 do begin victim = rawimage[x1[i]:x2[i],y1[i]:y2[i]] for j=0L,3 do begin if i ne j then begin source = rotate(rawimage[x1[j]:x2[j],y1[j]:y2[j]],rot[i,j]) victim -= source*xtcoefflo[i,j] if not keyword_set(ignorehigh) then begin high = where(source gt thresh, nhigh) nbad += nhigh if high[0] ne -1 then begin if keyword_set(mask) then begin y = long(2*high/naxis1) x = long(high - y*naxis1/2) maskimg[x+x1[i],y+y1[i]] = 0 endif if keyword_set(fixhigh) then begin victim[high] -= source[high]*xtcoeffhi[i,j] endif else begin victim[high] = xtalksat endelse endif endif endif endfor image[x1[i]:x2[i],y1[i]:y2[i]] = victim endfor ;; update header keywords sxaddpar, hdr, 'BITPIX', -32 sxaddpar, hdr, 'BZERO', 0 writefits, filebase[n] + '.x.fits', image, hdr if keyword_set(mask) then $ writefits, filebase[n] + '.mask.fits', maskimg if not keyword_set(ignorehigh) then begin print, 'Crosstalk Corrected: ' + filebase[n] + '.fits, ' + $ strtrim(nbad,2) + ' pixels were affected by saturated sources' endif else begin print, 'Crosstalk Corrected: ' + filebase[n] + '.fits' endelse endif endfor ;; do the overscan subtraction for n=0L, nfiles-1 do begin ;; if the overscan subtracted image already exists, don't re-process if not file_test(filebase[n] + '.xo.fits') then begin image = readfits(filebase[n] + '.x.fits', hdr, /silent) xover = sxpar(hdr,'OVERSCNX')/xbin[n] yover = sxpar(hdr,'OVERSCNY')/xbin[n] naxis1 = sxpar(hdr,'NAXIS1') naxis2 = sxpar(hdr,'NAXIS2') xtrim = naxis1-2*xover ytrim = naxis2-2*yover trimmed = fltarr(xtrim, ytrim) trimmed2 = fltarr(xtrim, ytrim) ;; trim the x axis and correct overscan if xover ne 0 then begin ;; for loop ~2x faster than matrix math! for i=0L,ytrim-1 do begin ;; subtract the median overscan value from the row trimmed[0:xtrim/2-1,i] = $ image[xover:xtrim/2+xover-1,i+yover] - $ median(image[0:xover-1,i+yover],/even) trimmed[xtrim/2:xtrim-1,i] = $ image[xtrim/2+xover:naxis1-xover-1,i+yover] - $ median(image[naxis1-xover:naxis1-1,i+yover],/even) endfor endif else trimmed = image ;;orient north up east left if keyword_set(noflip) then begin naxis1 = xtrim naxis2 = ytrim endif else begin if instrument[n] eq 'MDM4K' then begin naxis1 = ytrim naxis2 = xtrim trimmed = rotate(trimmed, 4) endif else if instrument[n] eq 'OSMOS' then begin naxis1 = xtrim naxis2 = ytrim trimmed = rotate(trimmed, 7) endif else message,'ERROR: undefined instrument ('+instrument[n]+')' endelse ;;add header values for coordinate solution sxaddpar, hdr, 'NAXIS1', naxis1 sxaddpar, hdr, 'NAXIS2', naxis2 sxaddpar, hdr, 'EPOCH', 2000 if instrument[n] eq 'MDM4K' then begin sxaddpar, hdr, 'SECPIX', 0.315*xbin[n] endif else if instrument[n] eq 'OSMOS' then begin sxaddpar, hdr, 'SECPIX', 0.273*xbin[n] endif else message,'ERROR: undefined instrument ('+instrument[n]+')' writefits, filebase[n] + '.xo.fits', trimmed, hdr print, 'Overscan Subtracted: ' + filebase[n] + '.x.fits' endif endfor ;; median combine sky-scaled flats by filter and bin ;;determine the unique filters sorted = sort(filter) filterndx = uniq(filter[sorted]) filters = filter[sorted[filterndx]] ;; determine the unique bins (assumes square binning) sorted = sort(xbin) binndx = uniq(xbin[sorted]) bins = xbin[sorted[binndx]] ;; find flats in each filter and bin for i=0L, n_elements(filters) - 1 do begin for j=0L, n_elements(bins) - 1 do begin flatlist = where(imgtype eq 'FLAT' and filter eq filters[i] and $ xbin eq bins[j] and slitid eq 'Open' and dispid eq 'Open') if n_elements(flatlist) ge 3 then begin if file_test(filepath + 'masterflat' + strtrim(bins[j],2) + $ strtrim(filters[i],2) + '.fits') then begin print, 'Warning: ' + 'masterflat' + strtrim(bins[j],2) + $ strtrim(filters[i],2) + $ '.fits exists. To reprocess, delete current file' endif else begin for k=0L, n_elements(flatlist)-1 do begin flat = float(readfits(filebase[flatlist[k]] + $ '.xo.fits',hdr, /silent)) ;; find sky value if instrument[0] eq 'OSMOS' then begin ;; OSMOS is not fully illuminated ;; determine scaling from central 1/4 of image naxis1 = sxpar(hdr, 'NAXIS1') naxis2 = sxpar(hdr, 'NAXIS2') mmm,flat[naxis1/2:naxis1*.75,naxis2/2:naxis2*.75],mode endif else if instrument[0] eq 'MDM4K' then begin mmm, flat, mode endif else message, 'ERROR: undefined instrument (' + $ instrument[0] + ')' mode = float(mode) if k eq 0 then begin ;; initialize image cube masterflats = flat/mode endif else begin ;; append to data cube masterflats = [[[masterflats]],[[flat/mode]]] endelse endfor ;;median combine all flats masterflat = median(masterflats, dimension=3,/even) ;; normalize to 1 mmm, masterflat, norm masterflat = masterflat / norm ;; make 0s into 1s (which cause infinities) bad = where(masterflat eq 0) if bad[0] ne -1 then masterflat[bad] = 1 writefits, filepath + 'masterflat' + strtrim(bins[j],2) + $ strtrim(filters[i],2) + '.fits', float(masterflat), hdr print, 'Master flat created: ' + filepath + 'masterflat' + $ strtrim(bins[j],2) + strtrim(filters[i],2) + '.fits' endelse endif endfor endfor ;; apply flats to all images for i=0L, n_elements(filters) - 1 do begin for j=0L, n_elements(bins) - 1 do begin images = where(imgtype eq 'OBJECT' and filter eq filters[i] and $ xbin eq bins[j] and slitid eq 'Open' and dispid eq 'Open') if images[0] ne - 1 then begin flatname = filepath + 'masterflat' + strtrim(bins[j],2) + $ strtrim(filters[i],2) + '.fits' if not file_test(flatname) then begin print, 'WARNING: No flat for '+ transpose(filebase[images]) + $ '.fits' ; stop endif else begin flat = readfits(filepath+'masterflat'+strtrim(bins[j],2) + $ strtrim(filters[i],2) + '.fits', /silent) for k=0L, n_elements(images)-1 do begin if file_test(filebase[images[k]] + '.xof.fits') then begin print, 'Warning: ' + filebase[images[k]] + $ 'of.fits exists. To reprocess, delete current file' endif else begin image = readfits(filebase[images[k]] + '.xo.fits', $ hdr, /silent) image = image / flat writefits, filebase[images[k]]+'.xof.fits', image, hdr print, 'Flat fielded: ' + filebase[images[k]] + $ '.xo.fits' endelse endfor endelse endif endfor endfor ;; apply fringe correction to I band images images = where(imgtype eq 'OBJECT' and filter eq 'I' and slitid eq 'Open' and dispid eq 'Open') if images[0] ne -1 then begin for i=0L, n_elements(images)-1 do begin if file_test(filebase[images[i]] + '.xoff.fits') then begin print, 'Warning: ' + filebase[images[i]] + $ 'off.fits exists. To reprocess, delete current file' endif else begin ;;fringe correction... fringename = 'fringe' + strtrim(xbin[images[i]],2) + '.fits' if not file_test(fringename) then begin if keyword_set(forcefringe) then begin print, 'No fringe image present; retrieving from ' + $ 'http://www.astronomy.ohio-state.edu/~jdeast/4k' spawn, 'wget http://www.astronomy.ohio-state.edu/' + $ '~jdeast/4k/fringe' + strtrim(xbin[images[i]],2)+'.fits' endif else begin print, 'Warning: No fringe image available. ' + $ 'Use /forcefringe to download a pre-made fringe image' endelse endif if file_test(fringename) then begin image = readfits(filebase[images[i]]+'.xof.fits',hdr, /silent) fringe = readfits('fringe' + strtrim(xbin[images[i]],2) + $ '.fits', /silent) fringe = fringe - mean(fringe) ;; find background ;; I think a better background fit would make for a ;; better fringe subtraction... mmm, image, background background = float(background) ;; find scale that solves ;; <((image-background) - scale*(fringe-))* ;; (fringe-)> = 0 scale = $ mean((image-background)*fringe,/nan)/mean(fringe^2,/nan) ;; subtract the fringe image = image - scale*fringe ;; add the scale of the fringe image to the header sxaddpar, hdr, 'fringecor', scale ;; write the fringe corrected image writefits, filebase[images[i]] + '.xoff.fits', image, hdr print, 'Fringe Corrected: ' + filebase[images[i]] + $ '.xof.fits with scale = ' + strtrim(scale,2) endif endelse endfor endif ;; coordinate solution reqfiles = ['daofind.sex','daofind.param','default.nnw','default.conv'] images = where(imgtype eq 'OBJECT' and slitid eq 'Open' and dispid eq 'Open') if images[0] ne -1 and not keyword_set(nowcs) and $ not keyword_set(noflip) then begin ;; check for required files, retrieve them if non-existant for j=0L,n_elements(reqfiles)-1 do begin if not file_test(reqfiles[j]) then begin print, 'Warning: ' + reqfiles[j] + ' required; retrieving from web' spawn, 'wget http://www.astronomy.ohio-state.edu/~jdeast/4k/' + $ reqfiles[j] if not file_test(reqfiles[j]) then $ message, 'Cannot retrieve ' + reqfiles[j] endif endfor ;; for each image for i=0L, n_elements(images)-1 do begin ;; make sure to analyze the fringe corrected image if exists if filter[images[i]] eq 'I' then begin ext = '.xoff' if not file_test(filebase[images[i]] + ext + '.fits') then $ ext = '.xof' endif else ext = '.xof' if file_test(filebase[images[i]] + ext + 'w.fits') then begin print, 'Warning: ' + filebase[images[i]] + ext + $ 'w.fits exists. To reprocess, delete current file' endif else begin ;; find stars for solution spawn, 'sex -c daofind.sex ' + filebase[images[i]] + ext + '.fits' ;; get 500 brightest stars for solution bright = read_ascii('bright.cat',data_start=4) openw, brightest, 'brightest.cat', /get_lun last = 500 < (file_lines('bright.cat')-4) printf, brightest, $ bright.field1[*,(sort(bright.field1[2,*]))[0:last-1]] free_lun, brightest ;; match brightest in image with brightest in USNOb catalog spawn, 'imwcs -vdw brightest.cat -c ub1 -h ' + strtrim(last,2) + $ ' ' + filebase[images[i]] + ext + '.fits -o ' + $ filebase[images[i]] + ext + 'w.fits' ;; creates DS9 region file for error checking ;; Open image->Region->Load Regions... ;; The 50 brightest stars should be circled spawn, 'imcat -c ub1 -n 50 -s mI ' + filebase[images[i]] + ext + $ 'w.fits | awk ' + "'{print " + $ '"fk5; circle("$2","$3",0.005) # color=blue text={"$8"}"}' + $ "' > " + filebase[images[i]] + ext + 'w.ub1.reg' print, 'Solved Coordinates: ' + filebase[images[i]] + ext + '.fits' endelse endfor endif end