|
##############################################################################
# Reduction of CTIO 0.9m telescope images taken with the cassegrain direct
# imager on Mar 20--27 1999 by Paul Schmidtke. Rolf Jansen -- Dec 1,4-6 2006
##############################################################################
# Dec 1 2006 [RAJ]:
#==================
$> cd /data2/raj/ast598/class_project/
$> mkdir 990320 ; mkdir 990321 ; mkdir 990322 ; mkdir 990323
$> mkdir 990324 ; mkdir 990325 ; mkdir 990326 ; mkdir 990327
# Open an 'xgterm' window and fire up IRAF in that window:
$> xgterm
%> cd ~/iraf
%> cl
# Enter the working directory where we will process the data:
cl> cd /data2/raj/ast598/class_project/
# Retrieve all data and move it into a 'raw/' subdirectory:
cl> dir.maxch=36 <--- type 'help dir' to see what this does
cl> dir raw/
night1 night2 night3 night4 night5 night6 night7 night8
# For some reason, the data in the directories for nights 7 and 8 lack write
# permission -- fix this:
cl> !chmod 755 raw/night7
cl> !chmod 644 raw/night7/*.fits
cl> !chmod 755 raw/night8
cl> !chmod 644 raw/night8/*.fits
# Create copies into our working directories per night; restore original
# frame numbering:
cl> !\ls -1 raw/night1/raw*.fits > get_n1.lis
cl> hselect @get_n1.lis irafname yes > tmp.lis # view contents of this file
cl> !sed 's:.imh:.fits:g' tmp.lis | sed 's:dflat::g' | sed 's:zero::g' | \
>>> sed 's:sflat::g' | sed 's:obj::g' | sed 's:^:990320/a0:g' > put_n1.lis
cl> count get_n1.lis,put_n1.lis
147 147 3528 get_n1.lis
147 147 2499 put_n1.lis
cl> delete tmp.lis yes
cl> !\ls -1 raw/night2/raw*.fits > get_n2.lis
cl> hselect @get_n2.lis irafname yes > tmp.lis # view contents of this file
cl> !sed 's:.imh:.fits:g' tmp.lis | sed 's:dflat::g' | sed 's:zero::g' | \
>>> sed 's:sflat::g' | sed 's:obj::g' | sed 's:^:990321/b0:g' > put_n2.lis
cl> count get_n2.lis,put_n2.lis
165 165 3960 get_n2.lis
165 165 2805 put_n2.lis
cl> delete tmp.lis yes
cl> !\ls -1 raw/night3/raw*.fits > get_n3.lis
cl> hselect @get_n3.lis irafname yes > tmp.lis # view contents of this file
cl> !sed 's:.imh:.fits:g' tmp.lis | sed 's:dflat::g' | sed 's:zero::g' | \
>>> sed 's:sflat::g' | sed 's:obj::g' | sed 's:^:990322/c0:g' > put_n3.lis
cl> count get_n3.lis,put_n3.lis
163 163 3912 get_n3.lis
163 163 2771 put_n3.lis
cl> delete tmp.lis yes
cl> !\ls -1 raw/night4/raw*.fits > get_n4.lis
cl> hselect @get_n4.lis irafname yes > tmp.lis # view contents of this file
cl> !sed 's:.imh:.fits:g' tmp.lis | sed 's:dflat::g' | sed 's:zero::g' | \
>>> sed 's:sflat::g' | sed 's:obj::g' | sed 's:^:990323/d0:g' > put_n4.lis
cl> count get_n4.lis,put_n4.lis
167 167 4008 get_n4.lis
167 167 2839 put_n4.lis
cl> delete tmp.lis yes
cl> !\ls -1 raw/night5/raw*.fits > get_n5.lis
cl> hselect @get_n5.lis irafname yes > tmp.lis # view contents of this file
cl> !sed 's:.imh:.fits:g' tmp.lis | sed 's:dflat::g' | sed 's:zero::g' | \
>>> sed 's:sflat::g' | sed 's:obj::g' | sed 's:^:990324/e0:g' > put_n5.lis
cl> count get_n5.lis,put_n5.lis
136 136 3264 get_n5.lis
136 136 2312 put_n5.lis
cl> delete tmp.lis yes
cl> !\ls -1 raw/night6/raw*.fits > get_n6.lis
cl> hselect @get_n6.lis irafname yes > tmp.lis # view contents of this file
cl> !sed 's:.imh:.fits:g' tmp.lis | sed 's:dflat::g' | \
>>> sed 's:^:990325/f:g' > put_n6.lis
cl> count get_n6.lis,put_n6.lis
14 14 336 get_n6.lis
14 14 252 put_n6.lis
cl> delete tmp.lis yes
cl> !\ls -1 raw/night7/raw*.fits > get_n7.lis
cl> hselect @get_n7.lis irafname yes > tmp.lis # view contents of this file
cl> !sed 's:.imh:.fits:g' tmp.lis | sed 's:dflat::g' | sed 's:zero::g' | \
>>> sed 's:sflat::g' | sed 's:obj::g' | sed 's:^:990326/g:g' | \
>>> sed 's:g9:g09:g' > put_n7.lis
cl> count get_n7.lis,put_n7.lis
170 170 4080 get_n7.lis
170 170 3060 put_n7.lis
cl> delete tmp.lis yes
cl> !\ls -1 raw/night8/raw*.fits > get_n8.lis
cl> hselect @get_n8.lis irafname yes > tmp.lis # view contents of this file
cl> !sed 's:.imh:.fits:g' tmp.lis | sed 's:dflat::g' | sed 's:zero::g' | \
>>> sed 's:obj::g' | sed 's:^:990327/h:g' > put_n8.lis
cl> count get_n8.lis,put_n8.lis
137 137 3288 get_n8.lis
137 137 2466 put_n8.lis
cl> delete tmp.lis yes
* In hindsight, there is also a 'OPICNUM' keyword in the image header, which
* we could have used to create the 'put_n?.lis' files more easilly. Oh well...
cl> !paste -d" " get_n1.lis put_n1.lis | sed 's:^:!\\cp -p :g' > getall.cl
cl> !paste -d" " get_n2.lis put_n2.lis | sed 's:^:!\\cp -p :g' >> getall.cl
cl> !paste -d" " get_n3.lis put_n3.lis | sed 's:^:!\\cp -p :g' >> getall.cl
cl> !paste -d" " get_n4.lis put_n4.lis | sed 's:^:!\\cp -p :g' >> getall.cl
cl> !paste -d" " get_n5.lis put_n5.lis | sed 's:^:!\\cp -p :g' >> getall.cl
cl> !paste -d" " get_n6.lis put_n6.lis | sed 's:^:!\\cp -p :g' >> getall.cl
cl> !paste -d" " get_n7.lis put_n7.lis | sed 's:^:!\\cp -p :g' >> getall.cl
cl> !paste -d" " get_n8.lis put_n8.lis | sed 's:^:!\\cp -p :g' >> getall.cl
# Check that the contents of file 'getall.cl' are indeed as expected, then:
cl> cl < getall.cl
cl> !du -sh
4.7G .
# Create a list of all data in the subdirectories per night:
cl> !\ls -1 9903??/*.fits > all.lis
# Print an inventory of the unique "object" names as they occur in the
# image headers. After displaying a few image headers with 'imheader' (with
# the 'longheader=yes' or 'lo+' parameter) and some experimenting:
cl> hselect @all.lis object yes | translit "-" "\"" " " | fields "-" 1 | \
>>> sort | unique
PG0918 <-- Landolt photometric standard field
PG1323 ""
PG1633 ""
Ru149 ""
X0925 <-- science target; this X-ray binary is focus of our class project
CAL83 <-- other science target(s)
GX339 ""
USco ""
USco#1 ""
USco#2 ""
X0513 ""
X0921 ""
X1822 ""
dflats <-- dome flats
b <-- twilight sky flats
v ""
r ""
i ""
zeroes <-- bias frames
cl> hselect @all.lis imagetyp yes | sort | unique
"DOME FLAT"
"SKY FLAT"
BIAS
OBJECT
# Now we know the unique portions of the object names, and the image types as
# recorded in the headers, create seperate lists for biases, (no darks were
# taken), dome flats, sky flats, and science target frames (only the Landolt
# photometric standards and the one science target of interest for this
# project):
cl> hselect @all.lis $I 'imagetyp=="BIAS"' > bias.lis
cl> hselect @all.lis $I 'imagetyp=="DOME FLAT"' > flat.lis
cl> hselect @all.lis $I 'imagetyp=="SKY FLAT"' > sky.lis
cl> hselect @all.lis $I,object 'imagetyp=="OBJECT"' > tmp.lis
cl> count all.lis
1099 1099 19782 all.lis
cl> count bias.lis,flat.lis,sky.lis,tmp.lis
105 105 1890 bias.lis
404 404 7272 flat.lis
79 79 1422 sky.lis
511 2146 17234 tmp.lis
1099 2734 27818 Total
cl> !egrep -e '(Ru149|X0925|PG1323|PG1633|PG0918)' tmp.lis | \
>>> awk '{print $1}' > sci.lis
cl> count sci.lis
290 290 5220 sci.lis
cl> delete tmp.lis yes
# In most image headers, keyword 'FILTER2' lists the filter used (b,v,r or i).
# Some headers, however, appear incomplete with a range of keywords missing.
# Luckily, in these particular data, the 'OBJECT' keyword contains not just
# the object ID, but also the filter information.
# Add the 'FILTER2' keyword if it is missing (with arbitrary value "-"):
cl> hedit ("@all.lis", "filter2", "-", add-, addonly+, del-, ver-, show+, upd+)
# Create an inventory with the relevant information:
cl> hselect @all.lis $I,imagetyp,filter2,object yes > _inventory
# After some experimenting (make sure to print the output of 'egrep' to the
# screen, before piping it into the 'awk' command, to make sure you only
# capture the frames in a given filter!):
cl> !egrep -e '(filter2 = b| B |b sky)' _inventory | \
>>> awk '{print $1}' > _tmpB.lis
cl> !egrep -e '(filter2 = v| V |v sky)' _inventory | \
>>> awk '{print $1}' > _tmpV.lis
cl> !egrep -e '(filter2 = r| R |r sky)' _inventory | \
>>> awk '{print $1}' > _tmpR.lis
cl> !egrep -e '(filter2 = i| I |i sky)' _inventory | \
>>> awk '{print $1}' > _tmpI.lis
cl> count _tmp?.lis
226 226 4068 _tmpB.lis
179 179 3222 _tmpI.lis
208 208 3744 _tmpR.lis
381 381 6858 _tmpV.lis
994 994 17892 Total
# So we're still missing 1099 - 994 = 105 frames -- but this is exactly the
# number of bias frames (see earlier). Thus every frame is accounted for.
# Now, add a new header keyword 'FILTER' to all image headers:
cl> unlearn hedit ; hedit.add=yes ; hedit.addonly=no ; hedit.verify=no
cl> hedit @_tmpB.lis filter "B"
cl> hedit @_tmpV.lis filter "V"
cl> hedit @_tmpR.lis filter "R"
cl> hedit @_tmpI.lis filter "I"
cl> hedit @bias.lis filter "D"
# Check in what filters our science target "X0925" (i.e.,X-ray binary
# RX J0925.7-4758) was observed:
cl> hselect @all.lis object yes | match "X0925" | sort | unique
"X0925 B 450s"
"X0925 R 240s"
"X0925 R 300s"
"X0925 V 300s"
# Since there are no I-filter observations for our object of interest, we can
# remove any and all frames from the lists in that filter:
cl> hselect @flat.lis $I 'filter!="I"' > _flat.lis
cl> !\mv _flat.lis flat.lis
cl> hselect @sky.lis $I 'filter!="I"' > _sky.lis
cl> !\mv _sky.lis sky.lis
cl> hselect @sci.lis $I 'filter!="I"' > _sci.lis
cl> !\mv _sci.lis sci.lis
cl> !\cp -p all.lis all.lis~
cl> !cat bias.lis flat.lis sky.lis sci.lis | sort > all.lis
cl> count bias.lis,flat.lis,sky.lis,sci.lis
105 105 1890 bias.lis
314 314 5652 flat.lis
61 61 1098 sky.lis
264 264 4752 sci.lis
744 744 13392 Total
cl> count all.lis
744 744 13392 all.lis
cl> delete _tmp?.lis yes verify-
# Fire up our image display tool, 'ds9':
cl> !ds9 &
# Check the size and data format of the raw images in pixels:
cl> imheader @put_n1.lis --> 1098x1024 pixels; unsigned integers ("ushort")
cl> gdevices
# ALIASES NX NY DESCRIPTION
... ... ...
imt4 imt1600 1600 1600
imt5 imt2048 2048 2048
... ... ...
imt44 imt1100 1100 1100 Mosaic block averaged by 8
... ... ...
# Each of the devices displayed above would work, choose one:
cl> set stdimage = imt44
cl> set stdimage = imt44
# Load any required packages that aren't loaded by default already:
cl> stsdas
cl> imred
cl> bias
cl> rjtools
# Construct a bad pixel/bad column list using the median of a few twilight
# sky flat exposures:
cl> head sky.lis nlines=5 > tmp.lis
cl> unlearn imcombine
cl> imcombine ("@tmp.lis", "tmp.fits", combine="average", reject="avsigclip",
>>> scale="mean", zero="none", weight="none", statsec="[65:1098,1:1024]",
>>> nkeep=1, mclip=yes, lsigma=3., hsigma=3., grow=1.0)
Dec 2 0:45: IMCOMBINE
combine = average, scale = mean, zero = none, weight = none
reject = avsigclip, mclip = yes, nkeep = 1
lsigma = 3., hsigma = 3.
grow = 1.
blank = 0.
statsec = Dec 2 0:45
Images Mean Scale
990320/a0061.fits 27452. 1.000
990320/a0062.fits 30291. 0.906
990320/a0063.fits 30139. 0.911
990320/a0064.fits 27768. 0.989
990320/a0065.fits 24857. 1.104
Output image = tmp.fits, ncombine = 5
cl> unlearn getregion
cl> getregion ("tmp.fits", "tmp.reg", format="basic", chkbound+, append+,
>>> verbose+)
GETREGION: IRAF2.12.2 raj@hedgehag Dec 02 00:47:32 2006
image = tmp.fits [1098 x 1024]
output= tmp.reg (format="basic")
Displaying image "tmp.fits": z1=94.47 z2=28930.48
Wait for the cursor cross to appear in the image display area,
then mark all regions by hitting "b" once at each of two dia-
gonally oposite corners to mark an arbitrary rectangular region,
or "c" once to mark a single pixel. Hit "q" to quit.
Type b again to draw box.
Type b again to draw box.
Type b again to draw box.
# GETREGION Finished
# Overlay the defined regions onto the image to check their validity:
cl> unlearn tvmarkall
cl> tvmarkall ("tmp.fits", gal-, star-, cosmic-, del-, region+,
>>> regfil="tmp.reg")
# Save the regions file as our badpixel file, with a descriptive header line
# prepended (see also the detector information in the Observing Log):
cl> hselect tmp.fits detector yes
Tek2K_3
cl> printf("# CTIO 0.9m Tek2K 2048x2046 24x24mu pixel, ", > "badpix.Tek2K_3")
cl> printf("quad-amplifier CCD;\n# Upper-right Amp#3 ", >> "badpix.Tek2K_3")
cl> printf("[74+]1024x1024 only, unbinned, untrimmed\n",>>"badpix.Tek2K_3")
cl> !cat tmp.reg >> badpix.Tek2K_3
cl> delete tmp.lis,tmp.reg yes verify-
# Determine the location and extent of the virtual overscan region proper:
cl> implot tmp.fits
:l 1 1024
:x 1 100
:y 1 28000
:y 575 625 --> [1:64,1:1024] is only overscan; [1:54,1:1024] is well-behaved
q
cl> !head -20 bias.lis | sed 's/$/[1:54,1:1024]/g' > tmp.lis
cl> imstat @tmp.lis fields="image,mean,stddev,min,max"
# IMAGE MEAN STDDEV MIN MAX
990320/a0046.fits[1:54,1:1024] 606.7 1.588 601. 611.
990320/a0047.fits[1:54,1:1024] 606.4 1.59 601. 611.
990320/a0048.fits[1:54,1:1024] 606.6 1.594 601. 611.
990320/a0049.fits[1:54,1:1024] 606.5 1.588 601. 611.
990320/a0050.fits[1:54,1:1024] 606.6 1.583 601. 611.
990320/a0051.fits[1:54,1:1024] 606.4 1.588 601. 611.
990320/a0052.fits[1:54,1:1024] 606.3 1.585 600. 611.
990320/a0053.fits[1:54,1:1024] 606.4 1.595 601. 611.
990320/a0054.fits[1:54,1:1024] 606.4 1.583 601. 611.
990320/a0055.fits[1:54,1:1024] 606.5 1.587 601. 611.
990320/a0056.fits[1:54,1:1024] 606.4 1.584 601. 611.
990320/a0057.fits[1:54,1:1024] 606.4 1.587 601. 611.
990320/a0058.fits[1:54,1:1024] 606.1 1.588 600. 610.
990320/a0059.fits[1:54,1:1024] 606.3 1.593 600. 611.
990320/a0060.fits[1:54,1:1024] 606.3 1.585 600. 611.
990321/b0209.fits[1:54,1:1024] 607.9 1.588 602. 612.
990321/b0210.fits[1:54,1:1024] 607.8 1.595 602. 612.
990321/b0211.fits[1:54,1:1024] 607.4 1.594 602. 612.
990321/b0212.fits[1:54,1:1024] 607.8 1.591 602. 612.
990321/b0213.fits[1:54,1:1024] 607.5 1.602 602. 612.
# From which we find: mean~606.5 ADU and rms~1.6 ADU, and min/max consistent
# with mean +/- 4 sigma. So, in the following, adopt:
TRIMSEC = STATSEC = '[70:1093,1:1024]' and BIASSEC = '[1:54,1:1024]'
(trim 5 extra columns on either side to wind up with 1024 x 1024 image format)
cl> delete tmp.lis,tmp.fits yes verify-
# Measure/verify the CCD gain and read-noise using pairs of flats and biases:
cl> obsutil
cl> !grep '990320' flat.lis | wc -l
46
cl> !grep '990320' bias.lis | wc -l
15
cl> !grep '990320' flat.lis | head -14 | sed 's/.fits//g' > _f1
cl> !grep '990320' flat.lis | head -15 | tail -14 | sed 's/.fits//g' > _f2
cl> !grep '990320' bias.lis | head -14 | sed 's/.fits//g' > _b1
cl> !grep '990320' bias.lis | tail -14 | sed 's/.fits//g' > _b2
cl> !paste _f1 _f2 _b1 _b2
cl> !echo 'findgain ("$1", "$2", "$3", "$4",' > getccdpar.tem
cl> !echo ' section="[70:1093,1:1024]", verbose-)' >> getccdpar.tem
cl> type getccdpar.tem
findgain ("$1", "$2", "$3", "$4",
section="[70:1093,1:1024]", verbose-)
cl> !narg getccdpar.tem _f1 '$1' _f2 '$2' _b1 '$3' _b2 '$4' > getccdpar.cl
cl> type getccdpar.cl
findgain ("990320/a0002", "990320/a0003", "990320/a0046", "990320/a0047",
section="[70:1093,1:1024]", verbose-)
findgain ("990320/a0003", "990320/a0004", "990320/a0047", "990320/a0048",
section="[70:1093,1:1024]", verbose-)
findgain ("990320/a0004", "990320/a0005", "990320/a0048", "990320/a0049",
section="[70:1093,1:1024]", verbose-)
... (etc) ...
cl> cl < getccdpar.cl
2.94 4.78
2.96 4.81
2.96 4.81
2.95 4.78
2.95 4.78
2.95 4.78
2.92 4.74
2.95 4.79
2.96 4.78
2.95 4.78
2.95 4.81
2.96 4.80
2.96 4.81
2.06 3.35
# And similarly for subsequent nights. This yields the following results:
night1: G = 2.950 +/- 0.020 R = 4.78 +/- 0.035
night2: G = 2.955 +/- 0.005 R = 4.81 +/- 0.010
night3: G = 2.960 +/- 0.005 R = 4.82 +/- 0.010
night4: G = 2.960 +/- 0.010 R = 4.82 +/- 0.020
night5: G = 2.960 +/- 0.005 R = 4.81 +/- 0.005
night7: G = 2.950 +/- 0.005 R = 4.75 +/- 0.015
night8: G = 2.950 +/- 0.030 R = 4.73 +/- 0.055
# So, on average, we find a gain value of 2.955 +/- 0.005 e-/ADU and a
# read-noise of 4.81 +/- 0.045 e-. The image headers and the table in the
# Observing Log file listed 3.0 e-/ADU and 4.6 e-, respectively.
# Our measurements are within are within 5% of the nominal values.
cl> delete _f?,_b? yes verify-
cl> bye
# Visually inspect all bias frames (remove any biases with problems):
cl> imexamine @bias.lis 1 allframes- nframes=1
All bias frames show a pattern-noise (analog/digital cross-talk or pick-up
noise), which is particularly strong in the biases of nights 2,4 and 7.
Although, using Fourier techniques, such noise can be filtered out in the
time-domain, such corrections would be way beyond the scope of this project.
# Convert all biases from 16-bit integer ("ushort") to 32-bit floating point
# ("real") format and log processing steps to the FITS headers:
cl> chpixtype @bias.lis @bias.lis "real" oldpixtype="all" verbose+
cl> type hdrlog_init.tem
gdate() ; sysinfo() ; unlearn hedit ; hedit.add=yes ; hedit.verify=no
s1 = "CCD image processing -- "//sysinfo.username//"@"//sysinfo.host
s1 = "-------- "//s1//" --------"
hedit ("@%.lis", "HISTORY", s1)
hedit ("@%.lis", "CHPXTYPE", "COMPLETE")
hedit ("@%.lis", "CHPXDATE", gdate.fdate)
cl> !sed 's:%:bias:g' hdrlog_init.tem > hdrlog_init.cl
cl> cl < hdrlog_init.cl
# Interpolate over the bad pixel columns identified:
cl> fixpix ("@bias.lis", "badpix.Tek2K_3", verbose+, >> "ccdproc.log")
cl> type hdrlog_fxpx.tem
gdate()
hedit ("@%.lis", "FIXPIX", "COMPLETE")
hedit ("@%.lis", "FXPXFILE", "/data2/raj/ast598/badpix.Tek2K_3")
hedit ("@%.lis", "FXPXDATE", gdate.fdate)
cl> !sed 's:%:bias:g' hdrlog_fxpx.tem > hdrlog_fxpx.cl
cl> cl < hdrlog_fxpx.cl
# Subtract the bias level as measured in the virtual overscan strip.
# First, measure that level and record it for later use:
cl> !sed 's/$/[1:54,1:1024]/g' bias.lis > _bias.lis
cl> unlearn imstat
cl> imstat @_bias.lis fields="image,midpt" nclip=3 format- >> overscan.dat
cl> delete _bias.lis yes
# Then, fit and remove the overscan level from all biases:
cl> unlearn colbias ; printf("\n", >> "ccdproc.log")
cl> colbias ("@bias.lis", "@bias.lis", bias="[1:54,1:1024]", trim="",
>>> median+, interactive-, function="legendre", order=1, niterate=3,
>>> logfiles="ccdproc.log")
cl> type hdrlog_ovsc.tem
gdate()
hedit ("@&.lis", "BIASSEC,TRIMSEC,DATASEC,CCDSEC,ORIGSEC", add-, del+)
hedit ("@&.lis", "OVERSCAN", "COMPLETE")
hedit ("@&.lis", "ORIGSEC", "[1:1098,1:1024]")
hedit ("@&.lis", "BIASSEC", "[1:54,1:1024]")
hedit ("@&.lis", "TRIMSEC", "[70:1093,1:1024]")
hedit ("@&.lis", "IMAGSEC", "[1:1024,1:1024]")
print ("match \"$1\" \"overscan.dat\" | fields \"-\" 2 | scanf(\"%8f\", x)",
> "addovscmean.tem")
print ("hedit (\"$1\", \"OVSCMEAN\", x)", >> "addovscmean.tem")
!narg addovscmean.tem &.lis '$1' > addovscmean.cl
cl < addovscmean.cl
hedit ("@&.lis", "OVSCDATE", gdate.fdate)
delete addovscmean.tem,addovscmean.cl yes verify-
cl> !sed 's:&:bias:g' hdrlog_ovsc.tem > hdrlog_ovsc.cl
cl> cl < hdrlog_ovsc.cl
# But since we didn't actually trim the images:
cl> hedit @bias.lis trimsec "[1:1098,1:1024]" add-
cl> hedit @bias.lis imagsec "[1:1098,1:1024]" add-
# Average the overscan-subtracted biases together (no multiplicative scaling,
# no additive offsets allowed, and no weighting):
cl> unlearn imcombine ; printf("\n", >> "ccdproc.log")
cl> !grep '990320' bias.lis > _b1.lis
cl> !grep '990321' bias.lis > _b2.lis
cl> !grep '990322' bias.lis > _b3.lis
cl> !grep '990323' bias.lis > _b4.lis
cl> !grep '990324' bias.lis > _b5.lis
cl> !grep '990326' bias.lis > _b7.lis
cl !grep '990327' bias.lis > _b8.lis
cl> type biascomb.tem
imcombine ("@%.lis", "BIAS%.fits", logfile="ccdproc.log",
combine="average", reject="avsigclip", outtype="real", scale="none",
zero="none", weight="none", statsec="[70:1093,1:1024]",
lthreshold=INDEF, hthreshold=64000., nkeep=1, mclip+, lsigma=3.,
hsigma=3., sigscale=0.1, grow=1.0)
cl> !sed 's:%:_b1:g' biascomb.tem > biascomb.cl
cl> cl < biascomb.cl
cl> !sed 's:%:_b2:g' biascomb.tem > biascomb.cl
cl> cl < biascomb.cl
cl> !sed 's:%:_b3:g' biascomb.tem > biascomb.cl
cl> cl < biascomb.cl
cl> !sed 's:%:_b4:g' biascomb.tem > biascomb.cl
cl> cl < biascomb.cl
cl> !sed 's:%:_b5:g' biascomb.tem > biascomb.cl
cl> cl < biascomb.cl
cl> !sed 's:%:_b7:g' biascomb.tem > biascomb.cl
cl> cl < biascomb.cl
cl> !sed 's:%:_b8:g' biascomb.tem > biascomb.cl
cl> cl < biascomb.cl
cl> imexamine BIAS_b?.fits 1 allframes- nframes=1
cl> !\ls -1 BIAS_b?.fits > _b0.lis
cl> !sed 's:%:_b0:g' biascomb.tem > biascomb.cl
cl> cl < biascomb.cl
cl> !\mv BIAS_b0.fits BIAS.fits
cl> delete BIAS_b?.fits,_b?.lis yes verify-
cl> gdate()
cl> hselect BIAS.fits ncombine yes | scanf("%d", i)
cl> print ("Combined ",i,"individual bias frames.")
Combined 105 individual bias frames.
cl> hedit ("BIAS.fits", "NCOMBINE,IMCMB*", add-, del+)
cl> hedit ("BIAS.fits", "BIASCOMB", "COMPLETE")
cl> hedit ("BIAS.fits", "NCOMBINE", i)
cl> hedit ("BIAS.fits", "BSCBDATE", gdate.fdate)
# Verify that the level in the overscan strip of the combined frame is exactly
# 0. If not, subtract the residual level:
cl> imstat BIAS.fits[1:54,1:1024] fields="mean" format- | scan (s1)
cl> print ("Residual overscan level is "//s1//" ADU")
Residual overscan level is 0.04551154 ADU
cl> unlearn imarith ; imarith.pixtype="real" ; imarith.calctype="real"
cl> printf("\n", >> "ccdproc.log")
cl> imarith ("BIAS.fits", "-", s1, "MASTERBIAS.fits", title="MASTERBIAS",
>>> verbose+, >> "ccdproc.log")
# Inspect the resulting MASTERBIAS frame; does the noise agree with our
# theoretical expectations?
cl> display MASTERBIAS.fits 1 zscale- zrange- z1=-2. z2=4.
# There is some electronic structure visible in the combined bias frame,
# including some semi-repeating pattern of short (~24x1 pixel) horizontal
# stripes. These are also visible in frames in which the biases of a single
# night are combined. There is also noticeable vertical structure.
cl> implot MASTERBIAS.fits
Note the near-perfect Gaussian distribution of pixel values around the primary
peak at ~0.15 ADU (left), and the odd-even effect and slight gradient toward
higher columns (right).
cl> s1="MASTERBIAS[1:1098,1:1024],MASTERBIAS[70:1093,1:1024],"
cl> s1=s1//"MASTERBIAS[1:54,1:1024]"
cl> imstat (s1, fields="image,mean,midpt,stddev,min,max", lower=INDEF,
>>> upper=INDEF, nclip=3, lsigma=3., usigma=3.)
# IMAGE MEAN MIDPT STDDEV MIN MAX
MASTERBIAS[1:1098,1:1024] 0.1320 0.1253 0.2046 -0.5028 0.7713
MASTERBIAS[70:1093,1:1024] 0.1378 0.1316 0.1987 -0.4605 0.7370
MASTERBIAS[1:54,1:1024] -0.00291 -0.00346 0.2326 -0.7014 0.6957
# Assuming Gaussian read-noise, the expected noise in the absence of any
# genuine variations in the column-to-column bias level would be:
# 4.81 [e-] / 2.955 [e-/ADU] / sqrt(105) = 0.158852 [ADU]
# The measured standard deviation is considerably larger (in the overscan
# strip by ~46%; in the image region by ~25%; in the entire frame by ~29%).
# Inspectrion of the MASTERBIAS image shows that this is likely due to the
# pronounced column-to-column struction that is visible. This CCD has some
# issues with its electronics!
# But the signal in most features in the MASTERBIAS frame is smaller than the
# read-noise (R = 4.81 e- = 1.63 ADU).
# There is no significant offset in the bias level as measured in overscan
# and image regions.
cl> imdelete @bias.lis yes verify- default+
# Visually inspect all dome flat frames and remove any flats with problems
# (e.g., overexposed/saturated, anomalous gradients, etc..):
cl> imstat @flat.lis fields="image,mean,max" --> none are saturated
cl> imexamine @flat.lis 1 allframes- nframes=1
# Convert all dome flats from 16-bit to 32-bit floating point format:
cl> chpixtype @flat.lis @flat.lis "real" oldpixtype="all" verbose+
cl> !sed 's:%:flat:g' hdrlog_init.tem > hdrlog_init.cl
cl> cl < hdrlog_init.cl
# Interpolate over the bad pixel columns:
cl> printf("\n", >> "ccdproc.log")
cl> fixpix ("@flat.lis", "badpix.Tek2K_3", verbose+, >> "ccdproc.log")
cl> !sed 's:%:flat:g' hdrlog_fxpx.tem > hdrlog_fxpx.cl
cl> cl < hdrlog_fxpx.cl
# As proven in a previous homework problem, we don't have to correct the dome
# flats for structure in the bias level.
# Subtract the bias level as measured in the virtual overscan strip.
# First, measure that level and record it for later use:
cl> !sed 's/$/[1:54,1:1024]/g' flat.lis > _flat.lis
cl> unlearn imstat
cl> imstat @_flat.lis fields="image,midpt" nclip=3 format- >> overscan.dat
cl> delete _flat.lis yes
# Now, fit and remove the overscan level from all dome flat frames:
cl> printf("\n", >> "ccdproc.log")
cl> colbias ("@flat.lis", "@flat.lis", bias="[1:54,1:1024]", median+,
>>> trim="[70:1093,1:1024]", interactive-, function="legendre", order=1,
>>> niterate=3, logfiles="ccdproc.log")
cl> !sed 's:&:flat:g' hdrlog_ovsc.tem > hdrlog_ovsc.cl
cl> cl < hdrlog_ovsc.cl
# Check whether there is a sequence of short exposures that may be used to
# correct for shutter shading:
cl> hselect @flat.lis exptime,filter yes | sort | unique
60.000 R
60.000 V
90.000 B
--> only long exposures; the shortest exposures (see also Observing Log) of
any type are 5 sec, so shutter shading is unlikely to be an issue.
# Let's see if we can combine all flats for the entire run into a single flat
# per filter. If it works well, we'll save time -- if not, we'll just have to
# create flats per night and reprocess the science frames from scratch.
cl> hselect @flat.lis $I 'filter=="B"' > flatB.lis
cl> hselect @flat.lis $I 'filter=="V"' > flatV.lis
cl> hselect @flat.lis $I 'filter=="R"' > flatR.lis
# NB: scale the flats to the level of the first (reference) flat in the list:
cl> unlearn imcombine ; printf("\n", >> "ccdproc.log")
cl> imcombine ("@flatB.lis", "FLAT_B.fits", combine="average",
>>> logfile="ccdproc.log", reject="avsigclip", outtype="real",
>>> scale="mean", zero="none", weight="none", statsec="[1:1024,1:1024]",
>>> lthreshold=INDEF, hthreshold=64000., nkeep=1, mclip=yes, lsigma=3.,
>>> hsigma=3., sigscale=0.1, grow=1.0)
cl> gdate()
cl> hselect FLAT_B.fits ncombine yes | scanf("%d", i)
cl> print ("Combined ",i,"individual B-flat frames.")
Combined 104 individual B-flat frames.
cl> hedit ("FLAT_B.fits", "NCOMBINE", add-, del+)
cl> hedit ("FLAT_B.fits", "FLATCOMB", "COMPLETE")
cl> hedit ("FLAT_B.fits", "NCOMBINE", i)
cl> hedit ("FLAT_B.fits", "FTCBDATE", gdate.fdate)
cl> printf("\n", >> "ccdproc.log")
cl> imcombine ("@flatV.lis", "FLAT_V.fits", combine="average",
>>> logfile="ccdproc.log", reject="avsigclip", outtype="real",
>>> scale="mean", zero="none", weight="none", statsec="[1:1024,1:1024]",
>>> lthreshold=INDEF, hthreshold=64000., nkeep=1, mclip=yes, lsigma=3.,
>>> hsigma=3., sigscale=0.1, grow=1.0)
cl> gdate()
cl> hselect FLAT_V.fits ncombine yes | scanf("%d", i)
cl> print ("Combined ",i,"individual V-flat frames.")
Combined 105 individual V-flat frames.
cl> hedit ("FLAT_V.fits", "NCOMBINE", add-, del+)
cl> hedit ("FLAT_V.fits", "FLATCOMB", "COMPLETE")
cl> hedit ("FLAT_V.fits", "NCOMBINE", i)
cl> hedit ("FLAT_V.fits", "FTCBDATE", gdate.fdate)
cl> printf("\n", >> "ccdproc.log")
cl> imcombine ("@flatR.lis", "FLAT_R.fits", combine="average",
>>> logfile="ccdproc.log", reject="avsigclip", outtype="real",
>>> scale="mean", zero="none", weight="none", statsec="[1:1024,1:1024]",
>>> lthreshold=INDEF, hthreshold=64000., nkeep=1, mclip=yes, lsigma=3.,
>>> hsigma=3., sigscale=0.1, grow=1.0)
cl> gdate()
cl> hselect FLAT_R.fits ncombine yes | scanf("%d", i)
cl> print ("Combined ",i,"individual R-flat frames.")
Combined 105 individual R-flat frames.
cl> hedit ("FLAT_R.fits", "NCOMBINE", add-, del+)
cl> hedit ("FLAT_R.fits", "FLATCOMB", "COMPLETE")
cl> hedit ("FLAT_R.fits", "NCOMBINE", i)
cl> hedit ("FLAT_R.fits", "FTCBDATE", gdate.fdate)
# Normalize the combined flats such that they will have a mean value of 1.0:
cl> unlearn imsurfit
cl> imstat FLAT_B.fits fields="midpt" format- | scanf("%8f", x)
cl> imsurfit ("FLAT_B.fits", "RESP_B.fits", xorder=9, yorder=9,
>>> type_output="response", function="spline3", cross_terms=yes,
>>> lower=3., upper=3., ngrow=1, niter=5, div_min=0.2)
cl> type hdrlog_nmlz.tem
gdate()
print ("Level prior to normalization: ",x,"ADU")
hedit ("%.fits", "NORMALIZ", "COMPLETE")
hedit ("%.fits", "NMLZMEAN", x)
hedit ("%.fits", "NMLZDATE", gdate.fdate)
cl> !sed 's:%:RESP_B:g' hdrlog_nmlz.tem > hdrlog_nmlz.cl
cl> cl < hdrlog_nmlz.cl
cl> imstat RESP_B.fits fields="image,mean,midpt,stddev,min,max" nclip=0
cl> imstat RESP_B.fits fields="image,mean,midpt,stddev,min,max" nclip=3
cl> imstat RESP_B.fits fields="image,mean,midpt,stddev,min,max" nclip=5
# IMAGE MEAN MIDPT STDDEV MIN MAX
nclip=0: RESP_B.fits 0.9934 0.9968 0.02456 0.2302 1.159
nclip=3: RESP_B.fits 0.998 0.9994 0.01231 0.9572 1.037
nclip=5: RESP_B.fits 0.9984 0.9991 0.01167 0.9627 1.034
cl> imstat RESP_B.fits fields="mean" nclip=0 format- | scan (s1)
cl> imarith ("RESP_B.fits", "/", s1, "RESP_B.fits")
cl> imstat FLAT_V.fits fields="midpt" format- | scanf("%8f", x)
cl> imsurfit ("FLAT_V.fits", "RESP_V.fits", xorder=9, yorder=9,
>>> type_output="response", function="spline3", cross_terms=yes,
>>> lower=3., upper=3., ngrow=1, niter=5, div_min=0.2)
cl> !sed 's:%:RESP_V:g' hdrlog_nmlz.tem > hdrlog_nmlz.cl
cl> cl < hdrlog_nmlz.cl
cl> imstat RESP_V.fits fields="image,mean,midpt,stddev,min,max" nclip=0
cl> imstat RESP_V.fits fields="image,mean,midpt,stddev,min,max" nclip=3
cl> imstat RESP_V.fits fields="image,mean,midpt,stddev,min,max" nclip=5
# IMAGE MEAN MIDPT STDDEV MIN MAX
nclip=0: RESP_V.fits 0.9956 0.9981 0.01955 0.2229 1.165
nclip=3: RESP_V.fits 0.9987 0.9995 0.01079 0.9641 1.032
nclip=5: RESP_V.fits 0.9989 0.9995 0.01047 0.9672 1.031
cl> imstat RESP_V.fits fields="mean" nclip=0 format- | scan (s1)
cl> imarith ("RESP_V.fits", "/", s1, "RESP_V.fits")
cl> imstat FLAT_R.fits fields="midpt" format- | scanf("%8f", x)
cl> imsurfit ("FLAT_R.fits", "RESP_R.fits", xorder=9, yorder=9,
>>> type_output="response", function="spline3", cross_terms=yes,
>>> lower=3., upper=3., ngrow=1, niter=5, div_min=0.2)
cl> !sed 's:%:RESP_R:g' hdrlog_nmlz.tem > hdrlog_nmlz.cl
cl> cl < hdrlog_nmlz.cl
cl> imstat RESP_R.fits fields="image,mean,midpt,stddev,min,max" nclip=0
cl> imstat RESP_R.fits fields="image,mean,midpt,stddev,min,max" nclip=3
cl> imstat RESP_R.fits fields="image,mean,midpt,stddev,min,max" nclip=5
# IMAGE MEAN MIDPT STDDEV MIN MAX
nclip=0: RESP_R.fits 0.9966 0.998 0.01802 0.1963 1.162
nclip=3: RESP_R.fits 0.9991 0.9996 0.01017 0.967 1.031
nclip=5: RESP_R.fits 0.9992 0.9997 0.009972 0.9691 1.029
cl> imstat RESP_R.fits fields="mean" nclip=0 format- | scan (s1)
cl> imarith ("RESP_R.fits", "/", s1, "RESP_R.fits")
# Did we reach our goal of 1,000,000 e- per pixel in the combined dome flats?
cl> hselect RESP_?.fits ncombine,nmlzmean yes
104 15727.73 --> 1,635,684 ADU * 2.955 e-/ADU = 4,833,446 e- --> OK!
105 15339.42 --> 1,610,639 ADU * 2.955 e-/ADU = 4,759,439 e- --> OK!
105 18395.92 --> 1,931,572 ADU * 2.955 e-/ADU = 5,707,794 e- --> OK!
# Visually inspect all twilight flat frames and remove any flats with problems
# (e.g., overexposed/saturated, anomalous gradients, etc..):
cl> imstat @sky.lis fields="image,mean,max" --> none are saturated
cl> imexamine @sky.lis 1 allframes- nframes=1
--> Although there are a few frames in which the stars are trailed (telescope
was not tracking), none of the twilight flats has a problem
# Convert all twilight flats from 16-bit to 32-bit floating point format:
cl> chpixtype @sky.lis @sky.lis "real" oldpixtype="all" verbose+
cl> !sed 's:%:sky:g' hdrlog_init.tem > hdrlog_init.cl
cl> cl < hdrlog_init.cl
# Interpolate over the bad pixel columns:
cl> printf("\n", >> "ccdproc.log")
cl> fixpix ("@sky.lis", "badpix.Tek2K_3", verbose+, >> "ccdproc.log")
cl> !sed 's:%:sky:g' hdrlog_fxpx.tem > hdrlog_fxpx.cl
cl> cl < hdrlog_fxpx.cl
# As proven in a previous homework problem, we don't have to correct the sky
# flats for structure in the bias level. The largest error we make (for three
# frames with a signal level below 10,000 ADU) is:
# ~3.0 ADU/(6551.-602.5) ~ 0.05% (where 602.5 ADU is the overscan level).
# Subtract the bias level as measured in the virtual overscan strip.
# First, measure that level and record it for later use:
cl> !sed 's/$/[1:54,1:1024]/g' sky.lis > _sky.lis
cl> unlearn imstat
cl> imstat @_sky.lis fields="image,midpt" nclip=3 format- >> overscan.dat
cl> delete _sky.lis yes
# Now, fit and remove the overscan level from all twilight sky flats:
cl> printf("\n", >> "ccdproc.log")
cl> colbias ("@sky.lis", "@sky.lis", bias="[1:54,1:1024]", median+,
>>> trim="[70:1093,1:1024]", interactive-, function="legendre", order=1,
>>> niterate=3, logfiles="ccdproc.log")
cl> !sed 's:&:sky:g' hdrlog_ovsc.tem > hdrlog_ovsc.cl
cl> cl < hdrlog_ovsc.cl
# Let's see if we can combine all sky flats for the entire run into a single
# frame per filter:
cl> hselect @sky.lis $I 'filter=="B"' > skyB.lis
cl> hselect @sky.lis $I 'filter=="V"' > skyV.lis
cl> hselect @sky.lis $I 'filter=="R"' > skyR.lis
cl> unlearn imcombine ; printf("\n", >> "ccdproc.log")
cl> imcombine ("@skyB.lis", "SKY_B.fits", combine="average",
>>> logfile="ccdproc.log", reject="avsigclip", outtype="real",
>>> scale="mean", zero="none", weight="none", statsec="[1:1024,1:1024]",
>>> lthreshold=INDEF, hthreshold=64000., nkeep=1, mclip=yes, lsigma=3.,
>>> hsigma=3., sigscale=0.1, grow=1.0)
cl> gdate()
cl> hselect SKY_B.fits ncombine yes | scanf("%d", i)
cl> print ("Combined ",i,"individual B-filter twilight flats.")
Combined 23 individual B-filter twilight flats.
cl> hedit ("SKY_B.fits", "NCOMBINE,IMCMB", add-, del+)
cl> hedit ("SKY_B.fits", "SKYCOMB", "COMPLETE")
cl> hedit ("SKY_B.fits", "NCOMBINE", i)
cl> hedit ("SKY_B.fits", "SKCBDATE", gdate.fdate)
cl> printf("\n", >> "ccdproc.log")
cl> imcombine ("@skyV.lis", "SKY_V.fits", combine="average",
>>> logfile="ccdproc.log", reject="avsigclip", outtype="real",
>>> scale="mean", zero="none", weight="none", statsec="[1:1024,1:1024]",
>>> lthreshold=INDEF, hthreshold=64000., nkeep=1, mclip=yes, lsigma=3.,
>>> hsigma=3., sigscale=0.1, grow=1.0)
cl> gdate()
cl> hselect SKY_V.fits ncombine yes | scanf("%d", i)
cl> print ("Combined ",i,"individual V-filter twilight flats.")
Combined 20 individual V-filter twilight flats.
cl> hedit ("SKY_V.fits", "NCOMBINE,IMCMB*", add-, del+)
cl> hedit ("SKY_V.fits", "SKYCOMB", "COMPLETE")
cl> hedit ("SKY_V.fits", "NCOMBINE", i)
cl> hedit ("SKY_V.fits", "SKCBDATE", gdate.fdate)
cl> printf("\n", >> "ccdproc.log")
cl> imcombine ("@skyR.lis", "SKY_R.fits", combine="average",
>>> logfile="ccdproc.log", reject="avsigclip", outtype="real",
>>> scale="mean", zero="none", weight="none", statsec="[1:1024,1:1024]",
>>> lthreshold=INDEF, hthreshold=64000., nkeep=1, mclip=yes, lsigma=3.,
>>> hsigma=3., sigscale=0.1, grow=1.0)
cl> gdate()
cl> hselect SKY_R.fits ncombine yes | scanf("%d", i)
cl> print ("Combined ",i,"individual R-filter twilight flats.")
Combined 18 individual R-filter twilight flats.
cl> hedit ("SKY_R.fits", "NCOMBINE,IMCMB*", add-, del+)
cl> hedit ("SKY_R.fits", "SKYCOMB", "COMPLETE")
cl> hedit ("SKY_R.fits", "NCOMBINE", i)
cl> hedit ("SKY_R.fits", "SKCBDATE", gdate.fdate)
# Normalize the combined twilight flats such that they will have mean=1.0:
cl> imstat SKY_B.fits fields="mean" nclip=0 format- | scanf("%8f", x)
cl> imarith ("SKY_B.fits", "/", str(x), "SKY_B.fits")
cl> !sed 's:%:SKY_B:g' hdrlog_nmlz.tem > hdrlog_nmlz.cl
cl> cl < hdrlog_nmlz.cl
cl> imstat SKY_V.fits fields="mean" nclip=0 format- | scanf("%8f", x)
cl> imarith ("SKY_V.fits", "/", str(x), "SKY_V.fits")
cl> !sed 's:%:SKY_V:g' hdrlog_nmlz.tem > hdrlog_nmlz.cl
cl> cl < hdrlog_nmlz.cl
cl> imstat SKY_R.fits fields="mean" nclip=0 format- | scanf("%8f", x)
cl> imarith ("SKY_R.fits", "/", str(x), "SKY_R.fits")
cl> !sed 's:%:SKY_R:g' hdrlog_nmlz.tem > hdrlog_nmlz.cl
cl> cl < hdrlog_nmlz.cl
# Fit a normalized illumination frame to the combined twilight sky flats
# after correcting for the pixel-to-pixel variations in the CCD response:
cl> imarith ("SKY_B.fits", "/", "RESP_B.fits", "_ILLUM_B.fits")
cl> display _ILLUM_B.fits 1
cl> imsurfit ("_ILLUM_B.fits", "ILLUM_B.fits", xorder=9, yorder=9,
>>> type_output="fit", function="spline3", cross_terms=yes, lower=3.,
>>> upper=3., ngrow=1, niter=5)
cl> imstat ILLUM_B.fits fields="mean" format- | scanf("%8f", x)
cl> imarith ("ILLUM_B.fits", "/", str(x), "ILLUM_B.fits", title="ILLUM_B")
cl> imstat ILLUM_B.fits fields="image,mean,midpt,stddev,min,max" nclip=0
# IMAGE MEAN MIDPT STDDEV MIN MAX
ILLUM_B.fits 1. 1.001 0.005207 0.9819 1.008
cl> imarith ("SKY_V.fits", "/", "RESP_V.fits", "_ILLUM_V.fits")
cl> display _ILLUM_V.fits 1
cl> imsurfit ("_ILLUM_V.fits", "ILLUM_V.fits", xorder=9, yorder=9,
>>> type_output="fit", function="spline3", cross_terms=yes, lower=3.,
>>> upper=3., ngrow=1, niter=5)
cl> imstat ILLUM_V.fits fields="mean" format- | scanf("%8f", x)
cl> imarith ("ILLUM_V.fits", "/", str(x), "ILLUM_V.fits", title="ILLUM_V")
cl> imstat ILLUM_V.fits fields="image,mean,midpt,stddev,min,max" nclip=0
# IMAGE MEAN MIDPT STDDEV MIN MAX
ILLUM_V.fits 1. 1.001 0.005612 0.9808 1.008
cl> imarith ("SKY_R.fits", "/", "RESP_R.fits", "_ILLUM_R.fits")
cl> display _ILLUM_R.fits 1
cl> imsurfit ("_ILLUM_R.fits", "ILLUM_R.fits", xorder=9, yorder=9,
>>> type_output="fit", function="spline3", cross_terms=yes, lower=3.,
>>> upper=3., ngrow=1, niter=5)
cl> imstat ILLUM_R.fits fields="mean" format- | scanf("%8f", x)
cl> imarith ("ILLUM_R.fits", "/", str(x), "ILLUM_R.fits", title="ILLUM_R")
cl> imstat ILLUM_R.fits fields="image,mean,midpt,stddev,min,max" nclip=0
# IMAGE MEAN MIDPT STDDEV MIN MAX
ILLUM_R.fits 1. 1.001 0.005483 0.9795 1.007
top-left: B; top-right: V; bottom-left: R
# NOTE, that the '_ILLUM_?.fits' frames, produced by dividing the twilight sky
# flats by the corresponding pixel-to-pixel response frame (constructed from
# the dome flats), show still significant high-frequency structure (lots of
# "speckles", particularly in B and somewhat less in R). Only the V filter
# frame is more or less as expected. Theoretically, there should not be any
# such features in these frames, only large-scale features and noise.
# This behavior might result from, e.g., a red-leak in the B filter, and/or
# very different light paths (and possibly a contribution from scattered
# light) for the dome and sky flats. Also, all dome flats appear to have
# been taken with an additional optic installed: both keywords 'FILTERS' and
# 'OBJECT' list the use of a "cb" (Color Balance) filter. The twilight sky
# flats, on the other hand, appear to have been taken with through a "dia"
# (Diaphragm), which may or may not be the default (it may simply mean an
# empty filter slot -- any hole is a diaphragm).
# Dec 4 2006 [RAJ]:
#==================
# Inspection of the headers of the science frames shows that these were taken
# with the "dia" and not with the "cb", so let's proceed with the first steps
# of the reduction of the science frames and then test whether division of the
# science frames by the twilight sky flats results in better flatfielding than
# division by the response and illumination frames. The S/N in the combined
# twilight flats is disappointingly low, however: we typically have only 3 sky
# flats per night per filter --- far fewer than required to accumulate 1e6 e-!
# Visually inspect all science frames; remove any which appear to have
# uncorrectible problems:
cl> imstat @sci.lis fields="image,mean,max" --> none are saturated
cl> imexamine @sci.lis 1 allframes- nframes=1 --> no problems
# There a couple of additional portions of bad columns visible in the science
# frames: add these to the existing bad pixel file:
cl> getregion ("990321/b0251.fits", "badpix.Tek2K_3", format="basic",
>>> chkbound+, append+, verbose+)
cl> type badpix.Tek2K_3
# CTIO 0.9m Tek2K 2048x2046 24x24mu pixel, quad-amplifier CCD;
# Upper-right Amp#3 [74+]1024x1024 only, unbinned, untrimmed
88 88 1 636
100 100 1 656
371 371 1 841
116 116 120 148
78 78 604 623
104 104 629 637
115 115 630 664
128 128 637 670
134 134 661 672
1001 1001 162 197
1000 1001 166 197
999 1001 185 194
998 1001 190 194
401 401 861 880
371 371 834 1024
# Convert all science frames from 16-bit to 32-bit floating point format:
cl> chpixtype @sci.lis @sci.lis "real" oldpixtype="all" verbose+
cl> !sed 's:%:sci:g' hdrlog_init.tem > hdrlog_init.cl
cl> cl < hdrlog_init.cl
# Interpolate over the bad pixel columns:
cl> printf("\n", >> "ccdproc.log")
cl> fixpix ("@sci.lis", "badpix.Tek2K_3", verbose+, >> "ccdproc.log")
cl> !sed 's:%:sci:g' hdrlog_fxpx.tem > hdrlog_fxpx.cl
cl> cl < hdrlog_fxpx.cl
# Correct the science frames for structure in the bias level:
cl> printf("\n", >> "ccdproc.log")
cl> imarith ("@sci.lis", "-", "MASTERBIAS.fits", "@sci.lis",
>>> verbose+, >> "ccdproc.log")
cl> type hdrlog_zrcr.tem
gdate()
hedit ("@%.lis", "ZEROCOR", "COMPLETE")
hedit ("@%.lis", "ZRCRIMAG", "/data2/raj/ast598/MASTERBIAS.fits")
hedit ("@%.lis", "ZRCRDATE", gdate.fdate)
cl> !sed 's:%:sci:g' hdrlog_zrcr.tem > hdrlog_zrcr.cl
cl> cl < hdrlog_zrcr.cl
# Subtract the bias level as measured in the virtual overscan strip.
# First, measure that level and record it for later use:
cl> !sed 's/$/[1:54,1:1024]/g' sci.lis > _sci.lis
cl> unlearn imstat
cl> imstat @_sci.lis fields="image,midpt" nclip=3 format- >> overscan.dat
cl> delete _sci.lis yes
# Now, fit and remove the overscan level from all science frames:
cl> printf("\n", >> "ccdproc.log")
cl> colbias ("@sci.lis", "@sci.lis", bias="[1:54,1:1024]", median+,
>>> trim="[70:1093,1:1024]", interactive-, function="legendre", order=1,
>>> niterate=3, logfiles="ccdproc.log")
cl> !sed 's:&:sci:g' hdrlog_ovsc.tem > hdrlog_ovsc.cl
cl> cl < hdrlog_ovsc.cl
# Now, test which of the flat fields correct best for the pixel-to-pixel
# variations in the response:
# B-filter:
cl> hselect @sci.lis $I 'filter=="B"' > sciB.lis
cl> imstat @sciB.lis fields=image,mean
--> 990320/a0145, 990322/c0403, 990326/g1077, g1081, g1085, and 990327/h1226,
h1235, h1239, h1243, h1247, h1251, h1257 and h1261 have significant sky
background levels and are therefore useful for this test.
cl> type testflat.cl
print ("Using dome flat:")
imarith (s1, "/", "RESP_B", "tmp1.fits", pixtype="real", calctype="real")
imarith ("tmp1", "/", "ILLUM_B", "tmp1.fits", pixtype="real", calctype="real")
getsky tmp1.fits lsigma=4. usigma=3. border- label- show-
x = getsky.smidpt-3*getsky.srms ; y = getsky.smidpt+3*getsky.srms
imreplace ("tmp1.fits", getsky.smidpt, lower=INDEF, upper=x, radius=1.)
imreplace ("tmp1.fits", getsky.smidpt, lower=y, upper=INDEF, radius=1.)
unlearn xdisplay ; xdispl ("tmp1.fits", buff=1)
getsky tmp1.fits lsigma=4. usigma=3. border- label- show+ verbose+
print ("Using sky flat:")
imarith (s1, "/", "SKY_B", "tmp2.fits", pixtype="real", calctype="real")
getsky tmp2.fits lsigma=4. usigma=3. border- label- show-
x = getsky.smidpt-3*getsky.srms ; y = getsky.smidpt+3*getsky.srms
imreplace ("tmp2.fits", getsky.smidpt, lower=INDEF, upper=x, radius=1.)
imreplace ("tmp2.fits", getsky.smidpt, lower=y, upper=INDEF, radius=1.)
xdispl ("tmp2.fits", buff=2, z1=xdisplay.tz1, z2=xdisplay.tz2)
getsky tmp2.fits lsigma=4. usigma=3. border- label- show+ verbose+
delete tmp1.fits,tmp2.fits yes verify-
cl> !sed 's:%:B:g' testflat.tem | sed 's:\$:990320/a0145:g' > testflat.cl
cl> cl < testflat.cl
Using dome flat:
iter npix mean midpt mode stddev min max
0 1048576 523.1012 522.1895 521.2093 13.20809 482.2268 562.6574
1 1048111 523.0839 522.7086 522.4329 13.18527 482.2268 561.8112
Using sky flat:
iter npix mean midpt mode stddev min max
0 1048576 523.1208 522.5424 522.5943 13.16755 482.4451 562.5724
1 1048314 523.1110 522.6833 522.5591 13.15460 482.4451 562.0449
cl> !sed 's:%:B:g' testflat.tem | sed 's:\$:990322/c0403:g' > testflat.cl
cl> cl < testflat.cl
Using dome flat:
iter npix mean midpt mode stddev min max
0 1048576 216.4738 216.0077 215.8981 8.608541 189.6467 242.6903
1 1047961 216.4587 215.9406 215.8368 8.588402 189.6467 241.8273
2 1047841 216.4558 215.9627 215.8256 8.584621 189.6467 241.7054
Using sky flat:
iter npix mean midpt mode stddev min max
0 1048576 216.4816 216.0209 215.9133 8.592964 189.7106 242.6437
1 1047969 216.4666 215.9506 215.8524 8.573112 189.7106 241.7973
2 1047847 216.4637 215.9741 215.8409 8.569274 189.7106 241.6687
cl> !sed 's:%:B:g' testflat.tem | sed 's:\$:990326/g1077:g' > testflat.cl
cl> cl < testflat.cl
Using dome flat:
iter npix mean midpt mode stddev min max
0 1048576 406.5255 406.0715 405.8487 11.74235 370.0365 442.0753
1 1048131 406.5106 405.8831 405.7863 11.72247 370.0365 441.2979
2 1047981 406.5056 405.9380 405.7676 11.71596 370.0365 441.0504
Using sky flat:
iter npix mean midpt mode stddev min max
0 1048576 406.5325 406.0594 405.8045 11.69614 370.1297 441.8622
1 1048183 406.5194 405.8626 405.7465 11.67873 370.1297 441.1471
2 1048024 406.5142 405.9150 405.7275 11.67188 370.1297 440.8987
cl> !sed 's:%:B:g' testflat.tem | sed 's:\$:990327/h1226:g' > testflat.cl
cl> cl < testflat.cl
Using dome flat:
iter npix mean midpt mode stddev min max
0 1048576 517.4422 516.9583 516.6761 13.16113 476.5320 557.3641
1 1048121 517.4250 516.7634 516.6047 13.13829 476.5320 556.4378
2 1047980 517.4199 516.8188 516.5848 13.13143 476.5320 556.1782
Using sky flat:
iter npix mean midpt mode stddev min max
0 1048576 517.4478 516.9378 516.5963 13.08870 476.6731 557.0079
1 1048160 517.4323 516.7161 516.5302 13.06802 476.6731 556.2035
2 1048003 517.4265 516.7762 516.5073 13.06044 476.6731 555.9169
==> CONCLUSION: despite the low S/N of the combined twilight sky flat in B,
the science images flatfielded with the twilight sky flat show slightly
smaller residuals than those flatfielded with the response frame derived
from the much higher S/N combined dome flat.
# R-filter:
cl> hselect @sci.lis $I 'filter=="R"' > sciR.lis
cl> imstat @sciR.lis fields=image,mean
--> 990320/a0147, 990321/b0288, 990324/e0734, 990326/g1079 and g1083,
990327/h1237, h1241, h1245, h1249 and h1259 have significant background
levels and are useful for this test.
cl> !sed 's:%:R:g' testflat.tem | sed 's:\$:990320/a0147:g' > testflat.cl
cl> cl < testflat.cl
Using dome flat:
iter npix mean midpt mode stddev min max
0 1048576 1078.535 1077.646 1077.678 18.91041 1019.995 1135.212
1 1048276 1078.518 1077.865 1077.617 18.88914 1019.995 1134.377
Using sky flat:
iter npix mean midpt mode stddev min max
0 1048576 1078.535 1077.212 1075.835 18.92808 1019.995 1135.305
1 1048073 1078.508 1078.002 1077.617 18.89261 1019.995 1133.996
cl> !sed 's:%:R:g' testflat.tem | sed 's:\$:990321/b0288:g' > testflat.cl
cl> cl < testflat.cl
Using dome flat:
iter npix mean midpt mode stddev min max
0 1048576 203.4695 203.2024 203.1548 8.338694 176.8858 229.7385
3 1046804 203.4268 203.0463 202.9721 8.280639 176.8858 227.8867
Using sky flat:
iter npix mean midpt mode stddev min max
0 1048576 203.4664 203.1867 203.1311 8.338676 176.8614 229.7161
3 1046813 203.4239 203.0300 202.9495 8.280992 176.8614 227.8718
cl> !sed 's:%:R:g' testflat.tem | sed 's:\$:990324/e0734:g' > testflat.cl
cl> cl < testflat.cl
Using dome flat:
iter npix mean midpt mode stddev min max
0 1048576 362.1544 361.5642 361.4330 11.00192 327.8738 395.9093
2 1047679 362.1262 361.5432 361.3180 10.96422 327.8738 394.4147
Using sky flat:
iter npix mean midpt mode stddev min max
0 1048576 362.1492 361.5271 361.3544 11.00441 327.7884 395.8335
2 1047685 362.1211 361.5010 361.2442 10.96704 327.7884 394.3755
cl> !sed 's:%:R:g' testflat.tem | sed 's:\$:990327/h1237:g' > testflat.cl
cl> cl < testflat.cl
Using dome flat:
iter npix mean midpt mode stddev min max
0 1048576 380.8853 379.8403 379.4620 11.30364 344.9834 415.9710
1 1047152 380.8393 380.5177 380.3947 11.24196 344.9834 413.7509
Using sky flat:
iter npix mean midpt mode stddev min max
0 1048576 380.8877 379.8513 379.4826 11.31265 344.9769 416.0338
1 1047108 380.8401 380.5249 380.4093 11.24901 344.9769 413.7892
==> CONCLUSION: like for B, despite the low S/N of the combined sky flat in R,
the science images flatfielded with the twilight sky flat show slightly
smaller residuals than those flatfielded with the response frame derived
from the much higher S/N combined dome flat.
# V-filter:
cl> hselect @sci.lis $I 'filter=="V"' > sciV.lis
cl> imstat @sciV.lis fields=image,mean
--> e.g.: 990324/e0733 and e0764, 990326/g1065, g1066, and g1101, 990327/h1258
cl> !sed 's:%:V:g' testflat.tem | sed 's:\$:990324/e0733:g' > testflat.cl
cl> cl < testflat.cl
Using dome flat:
iter npix mean midpt mode stddev min max
0 1048576 297.1774 296.6868 296.5258 9.948063 266.1793 327.5798
2 1047700 297.1525 296.6386 296.4266 9.914887 266.1793 326.3595
Using sky flat:
iter npix mean midpt mode stddev min max
0 1048576 297.1747 296.6720 296.4876 9.947205 266.1455 327.5273
2 1047714 297.1503 296.6158 296.3903 9.914663 266.1455 326.3471
cl> !sed 's:%:V:g' testflat.tem | sed 's:\$:990324/e0764:g' > testflat.cl
cl> cl < testflat.cl
Using dome flat:
iter npix mean midpt mode stddev min max
0 1048576 277.4610 276.6611 276.4467 9.656654 246.9906 306.8544
2 1047399 277.4288 276.6273 276.3168 9.614150 246.9906 305.4659
Using sky flat:
iter npix mean midpt mode stddev min max
0 1048576 277.4637 276.9170 276.4789 9.656885 247.0222 306.9037
1 1047798 277.4423 277.3152 277.3499 9.628307 247.0222 305.8871
cl> !sed 's:%:V:g' testflat.tem | sed 's:\$:990326/g1065:g' > testflat.cl
cl> cl < testflat.cl
Using dome flat:
iter npix mean midpt mode stddev min max
0 1048576 4305.776 4304.534 4303.989 37.72071 4188.934 4420.656
2 1047976 4305.711 4304.204 4303.728 37.63447 4188.934 4417.018
Using sky flat:
iter npix mean midpt mode stddev min max
0 1048576 4305.769 4304.424 4303.728 37.85661 4188.255 4420.808
2 1047974 4305.704 4304.059 4303.464 37.77006 4188.255 4417.271
cl> !sed 's:%:V:g' testflat.tem | sed 's:\$:990327/h1258:g' > testflat.cl
cl> cl < testflat.cl
Using dome flat:
iter npix mean midpt mode stddev min max
0 1048576 523.2329 522.3889 522.1535 13.17266 481.9750 563.4398
2 1047694 523.1998 522.3489 522.0222 13.12872 481.9750 561.7194
Using sky flat:
iter npix mean midpt mode stddev min max
0 1048576 523.2357 522.4069 522.1796 13.18278 481.9701 563.4783
2 1047683 523.2022 522.3661 522.0477 13.13819 481.9701 561.7564
==> CONCLUSION: like for B, despite the low S/N of the combined sky flat in V,
the science images flatfielded with the twilight sky flat show slightly
smaller residuals than those flatfielded with the response frame derived
from the much higher S/N combined dome flat.
# In the following, we will therefore discard the "master" calibration frames
# derived from the dome flats, i.e., 'RESP_?.fits' and 'ILLUM_?.fits, and only
# use the combined and normalized twilight sky flats (we don't need an extra
# illumination frame in this case, since the illumination pattern in the
# twilight sky flats should resemble quite closely that in the night sky.
# Flatfield the science frames per filter, correcting for both pixel-to-pixel
# and illumination variations:
cl> type hdrlog_ftcr.tem
gdate()
hedit ("@%.lis", "FLATCOR", "COMPLETE")
hedit ("@%.lis", "FTCRIMAG", "/data2/raj/ast598/$.fits")
hedit ("@%.lis", "FTCRDATE", gdate.fdate)
cl> unlearn imarith ; imarith.pixtype="real" ; imarith.calctype="real"
cl> printf("\n", >> "ccdproc.log")
cl> imarith ("@sciB.lis", "/", "SKY_B.fits", "@sciB.lis",
>>> verbose+, >> "ccdproc.log")
cl> !sed 's:%:sciB:g' hdrlog_ftcr.tem | sed 's:\$:SKY_B:g' > hdrlog_ftcr.cl
cl> cl < hdrlog_ftcr.cl
cl> printf("\n", >> "ccdproc.log")
cl> imarith ("@sciV.lis", "/", "SKY_V.fits", "@sciV.lis",
>>> verbose+, >> "ccdproc.log")
cl> !sed 's:%:sciV:g' hdrlog_ftcr.tem | sed 's:\$:SKY_V:g' > hdrlog_ftcr.cl
cl> cl < hdrlog_ftcr.cl
cl> printf("\n", >> "ccdproc.log")
cl> imarith ("@sciR.lis", "/", "SKY_R.fits", "@sciR.lis",
>>> verbose+, >> "ccdproc.log")
cl> !sed 's:%:sciR:g' hdrlog_ftcr.tem | sed 's:\$:SKY_R:g' > hdrlog_ftcr.cl
cl> cl < hdrlog_ftcr.cl
# To test quality of the flatfielding, examine a couple of representative
# frames with a non-negligible sky background level:
cl> display 990327/h1257.fits 1
cl> implot 990327/h1257.fits
:l 1 1024
:y 515 575
:c 1 1024
q
# The residual systematic variations in the background on scales of ~100 pixels
# or more appear to be less than ~1.5 ADU (peak-to-peak) on an average level of
# ~523.5 ADU, i.e., < +/-0.15% (<~0.30% peak to peak).
# On small scales, the expected noise in the background in 5x5 pixel regions
# is: 1
# rms (ADU) = -------- * ( 25*(523.5*2.955 + 4.81**2) ) / 2.955 = 13.41 ADU
# sqrt(25)
# We observe, on the other hand:
# B-filter:
cl> imexamine
# SECTION NPIX MEAN MEDIAN STDDEV MIN MAX
[410:414,474:478] 25 521.1 521.7 15.53 492.6 544.2
[230:234,529:533] 25 525.3 526. 15.84 492.3 560.2
[264:268,567:571] 25 524.5 523.5 13.76 488.7 547.9
[314:318,601:605] 25 520.1 523.6 16.55 486.4 548.6
[375:379,653:657] 25 523.4 522.3 15.58 493. 551.6
[555:559,627:631] 25 521.4 522.3 12.1 489.3 547.
[694:698,640:644] 25 524.9 527.1 15.83 482.8 554.8
[761:765,650:654] 25 526.8 524.9 14.91 505.1 559.3
[772:776,724:728] 25 523.2 521.5 13.04 501.7 544.9
[731:735,831:835] 25 522.7 522.5 13.63 500.4 551.3
# i.e., a mean rms of 14.68 ADU. This corresponds to an excess noise of:
# sqrt( 14.68**2 - 13.41**2 ) = 5.97 ADU compared to the theoretical noise.
# Indeed, close inspection of blank regions in the images appear to show a
# slightly "grainy/salt-'n-pepper" appearance.
# V-filter:
cl> display 990327/h1272.fits 1
cl> implot 990327/h1272.fits
:l 1 1024
:y 515 575
:c 1 1024
q
cl> imexamine
# SECTION NPIX MEAN MEDIAN STDDEV MIN MAX
[275:279,571:575] 25 578.3 573.1 12.99 561.6 602.8
[484:488,743:747] 25 579.7 581.4 10.31 555.1 597.9
[662:666,757:761] 25 578.2 574.4 12.25 557.7 602.3
[702:706,670:674] 25 581.3 577.1 17.2 553.1 624.9
[700:704,600:604] 25 581.7 581.2 14.36 544.3 606.5
[688:692,510:514] 25 578.7 581.6 12.42 554.7 608.6
[587:591,511:515] 25 576.3 575. 11.3 549.5 602.8
[524:528,129:133] 25 575.6 575.6 13.16 542.8 596.6
[672:676,194:198] 25 575.5 574.8 13.01 544.6 596.9
[887:891,460:464] 25 583.2 585.8 12.32 564.5 607.1
# Similar, but no noticeable excess noise.
# Next, we will try to photometrically calibrate the images. For differential
# photometry with respect to stars with known magnitudes this may not be
# necessary, as long as the color terms of the instrumental system are
# negligible with respect to the standard (Landolt) photometric system. But,
# we have no a priori information regarding the significance of such color
# terms.
# Load the digital photometry packages:
cl> digiphot
cl> apphot
cl> photcal
# Create the lists of standard star images:
cl> hselect @sci.lis $I,object yes | match "PG" | fields "-" 1 > tmp.lis
cl> hselect @sci.lis $I,object yes | match "Ru" | fields "-" 1 >> tmp.lis
cl> !sort tmp.lis | sed 's:.fits::g' | sed 's: ::g' > standards.lis
# (the latter 'sed' strips the trailing space characters)
cl> delete tmp.lis yes
cl> hselect @standards.lis object yes | translit "-" "\"" " " | \
>>> fields "-" 1 > tmp.lis
cl> !sed 's: ::g' tmp.lis > stdnames.lis
cl> delete tmp.lis yes
# Create pixel coordinate files of each photometric standard star within the
# observed Landolt field. First we need to know what number of pixels
# corresponds to a 7" radius. The headers list a pixel scale of 0.384"/pix,
# so 7" ~ 18.2 pixels. The choice of the inner radius of the sky annulus is
# fairly arbitrary, but should be well away from the outer radius of the
# measurement aperture. I choose ~2 x 7" --> 36 pixels. The minimum width
# that is acceptable would have the annulus encompass ~1000 pixels, yielding:
# (R_out)_min = sqrt( 1000/pi + 36**2 ) = 40.18 pixels. I choose a width of
# 11 pixels (i.e., R_out = 46 pixels) giving me 2576 sky pixels: 2.5 times
# the minimum number of pixels.
cl> type getstdcoo.tem
print ("\nLandolt field &2 ...")
markstds ("&1.fits", radstr="18.2,36,46", refstar="&2")
fields ("&1.stdcoo", "3", > "&1.stdcoo1")
ctrcoo ("&1.stdcoo", "&1.stdcoo2", "&1.fits",
cntrbox=7, redispl-, verbose-)
!paste -d" " &1.stdcoo2 &1.stdcoo1 | \
awk '{printf("%.2f %.2f \"%s\"\n",$1,$2,$3)}' > &1.stdcoo
print ("CTRCOO: outab = &1.stdcoo") ; type &1.stdcoo
delete &1.stdcoo? yes verify-
cl> !narg getstdcoo.tem standards.lis '&1' stdnames.lis '&2' > getstdcoo.cl
cl> type getstdcoo.cl <-- to check that the substitution went OK
# Have your Landolt finding-charts at hand, then run the script we prepared.
# It may be convenient to flip the DS9 view in X and/or Y using the "Zoom" ->
# "Invert X","Invert Y","Invert X&Y" pull-down menu, or the "Zoom" -> "x","y",
# or "xy" buttons of DS9; doing so does not affect IRAF).
# Since the images for the class project represent the upper right quadrant
# of the full 2048x2048 pixel CCD and were read out through the amplifier
# associated with that quadrant, the read-out coordinates are flipped in Y
# with respect to the celestial coordinates. The images headers do not
# contain the appropriate WCS keywords (CRVALn,CRPIXn,CDELTm,CTYPEn) to tell
# us about this flip, here is the place we first find out (had you obtained
# the observations yourself, you probably would have noticed at the telescope
# and made a note in your observing log regarding the orientation).
cl> unlearn markstds
cl> cl < getstdcoo.cl
... (etc) ...
Landolt field PG0918 ...
MARKSTDS: NOAO/IRAF2.12.2 raj@hedgehag Dec 05 20:18:30 2006
image = 990324/e0762.fits [1024 x 1024]
Wait for the cursor cross to appear in the image display area, then mark
all (standard) stars in the image using the 'c' key. Hit 'q' to quit.
Note: start with the main (i.e., naming or reference) star in the field.
MARKSTDS: Finished.
CTRCOO: outab = 990324/e0762.stdcoo
704.61 810.27 "PG0918"
446.88 769.05 "1"
526.71 519.98 "2"
177.68 726.93 "3"
937.85 594.09 "4"
# For the standard star aperture photometry to be meaningful, we must make
# sure that each FITS header contains keyword values for the exposure time,
# airmass, filter (we already made sure of that right at the start of our
# image processing), and universal time. If not, we must add/fix such header
# keywords/values, or (if desperate) omit such images from our analysis.
# By examining a few random headers, it looks like there are two keywords
# related to the UT: UT and UTSHUT, where the latter appears to be delayed by
# couple of seconds with respect to the former (commanded versus actual time
# of opening the shutter, probably corresponding to the time it takes to clear
# the CCD). The exposure time keyword is EXPTIME, the airmass keyword is
# AIRMASS, and the filter keyword that we added ourselves is FILTER:
cl> hselect @standards.lis $I,exptime,airmass,utshut,filter yes
990320/a0079 60.000 1.155 23:52:17.0 R
990320/a0080 60.000 1.154 23:54:38.0 V
990320/a0081 90.000 1.153 23:56:58.0 B
990320/a0095 90.000 2.087 04:02:09.0 B
990320/a0096 60.000 2.124 04:04:49.0 V
990320/a0097 60.000 2.156 04:06:59.0 R
990320/a0105 90.000 1.120 05:06:41.0 B
990320/a0106 60.000 1.116 05:09:25.0 V
990320/a0107 60.000 1.113 05:11:42.0 R
990320/a0122 90.000 1.700 06:53:47.0 B
990320/a0123 60.000 1.679 06:57:06.0 V
990320/a0124 60.000 1.665 06:59:15.0 R
990320/a0142 60.000 09:44:28.0 R <--- problem: AIRMASS missing
990320/a0143 60.000 1.310 09:46:58.0 V
990320/a0144 90.000 1.311 09:49:07.0 B
... ... ... ... ..
# Luckily the Airmass for frame 990320/a0142 was noted in the Observing Logs
# as being 1.308. Add this keyword to the header:
cl> hedit 990320/a0142.fits airmass 1.308 add+
# Perform the aperture photometry to measure the instrument magnitudes for
# all standard stars in each field. For my arbitrary instrumental zeropoint,
# I will use 26 mag throughout. We use the gain and read-noise values we
# measured earlier for this CCD:
cl> type getstdphot.tem
qphot ("%.fits", cbox=5., annulus=36., dannulus=11., apertures="18.2",
coords="%.stdcoo", output="%.stdmag", zmag=26.,
exposure="exptime", airmass="airmass", filter="filter",
obstime="utshut", epadu=2.955, interactive-, verbose+)
cl> unlearn apphot ; unlearn qphot
cl> !narg getstdphot.tem standards.lis '%' > getstdphot.cl
cl> type getstdphot.cl <-- to check that the substitution went OK
cl> cl < getstdphot.cl
a0079.fits 549.04 682.69 85.17189 18.038 ok
a0079.fits 592.01 654.82 84.34295 18.410 ok
a0079.fits 429.39 688.55 84.29446 16.378 ok
a0079.fits 438.44 587.74 83.99836 18.433 ok
a0079.fits 511.12 642.34 84.87129 15.579 ok
a0079.fits 396.35 420.24 84.41537 17.502 ok
a0079.fits 558.54 469.45 84.28796 16.983 ok
a0079.fits 639.44 517.38 84.28383 16.616 ok
a0080.fits 549.63 683.10 40.96022 17.989 ok
a0080.fits 592.60 655.35 40.00893 18.601 ok
a0080.fits 430.24 688.80 39.92291 16.744 ok
a0080.fits 439.24 588.35 39.90656 18.529 ok
... ... ... ... ... ..
--> none of the instrument magnitudes seem outrageously small or large, and
no errors were issued (the final column lists "ok" for all).
# Prepare matched input and configuration files required by 'photcal' to solve
# for the absolute photometric zeropoints, atmospheric extinction and color
# terms. This has to be done per night (since not all nights need have been
# photometric and not all photometric nights have the same atmospheric
# transparency):
cl> hselect @standards.lis object,filter,utshut,$I yes | sort | \
>>> match 990320 > standards.n1
cl> hselect @standards.lis object,filter,utshut,$I yes | sort | \
>>> match 990321 > standards.n2
cl> hselect @standards.lis object,filter,utshut,$I yes | sort | \
>>> match 990322 > standards.n3
cl> hselect @standards.lis object,filter,utshut,$I yes | sort | \
>>> match 990323 > standards.n4
cl> hselect @standards.lis object,filter,utshut,$I yes | sort | \
>>> match 990324 > standards.n5
# Using the information in files 'standards.n?', create (with favorite ASCII
# text editor) a file that lists per night the matched sets of images of the
# same field in the prescribed format. In 'qphot', we specified the image
# extension ".fits" as part of the file names, so here we must do so too:
cl> copy standards.n1 990320/.
cl> type 990320/standards.sets <-- this shows the output AFTER
PG1323-086 : a0105.fits a0106.fits a0107.fits editing, not that after the
PG1323-086 : a0145.fits a0146.fits a0147.fits 'copy' command!
PG1633+099 : a0122.fits a0123.fits a0124.fits
PG1633+099 : a0144.fits a0143.fits a0142.fits
RU_149 : a0081.fits a0080.fits a0079.fits
RU_149 : a0095.fits a0096.fits a0097.fits
cl> copy standards.n2 990321/.
cl> type 990321/standards.sets
PG0918+029 : b0269.fits b0270.fits b0271.fits
PG1323-086 : b0281.fits b0282.fits b0283.fits
PG1633+099 : b0290.fits b0291.fits b0292.fits
RU_149 : b0249.fits b0250.fits b0251.fits
RU_149 : b0260.fits b0261.fits b0262.fits
cl> copy standards.n3 990322/.
cl> type 990322/standards.sets
PG0918+029 : c0422.fits c0423.fits c0424.fits
PG1323-086 : c0443.fits c0444.fits c0445.fits
PG1633+099 : c0448.fits c0449.fits c0450.fits
PG1633+099 : c0473.fits c0474.fits c0475.fits
RU_149 : c0403.fits c0404.fits c0405.fits
RU_149 : c0426.fits c0427.fits c0428.fits
cl> copy standards.n4 990323/.
cl> type 990323/standards.sets
PG0918+029 : d0581.fits d0582.fits d0583.fits
PG1323-086 : d0607.fits d0608.fits d0609.fits
PG1323-086 : d0639.fits d0640.fits d0641.fits
PG1633+099 : d0611.fits d0612.fits d0613.fits
PG1633+099 : d0635.fits d0636.fits d0637.fits
RU_149 : d0580.fits INDEF INDEF
cl> copy standards.n5 990324/.
cl> type 990324/standards.sets
PG0918+029 : e0760.fits e0761.fits e0762.fits
RU_149 : e0732.fits e0733.fits e0734.fits
RU_149 : e0755.fits e0756.fits e0757.fits
# Create the observations files, tabulating the measured instrument magnitudes
# and associated information (filter, UT, airmass). Formally, we would have
# to compute the airmasses at the **MIDDLE** of each exposure. The headers
# contain the airmass at the **START** of each exposure. Since the exposures
# are all short, I'll ignore this step here (at the celestial equator, objects
# move at most (90/3600)*15deg = 0.375 deg).
cl> cd 990320/
cl> unlearn mkobsfile
cl> mkobsfile ("*.stdmag", "B,V,R", "standards.obs", imsets="standards.sets",
>>> minmagerr=0.001, shifts="", apercors="", tolerance=5., verbose+)
Observations file: standards.obs
Image set: PG1323-086 4 stars written to the observations file
Image set: PG1323-086 4 stars written to the observations file
Image set: PG1633+099 5 stars written to the observations file
Image set: PG1633+099 5 stars written to the observations file
Image set: RU_149 8 stars written to the observations file
Image set: RU_149 8 stars written to the observations file
--> reported number of stars matches that expected, so there were no shifts
between exposures in a given sequence.
# Edit file 'standards.obs' after creating a backup copy:
cl> copy standards.obs standards.obs.orig
# edit
cl> type standards.obs
# FIELD FILTER OTIME AIRMASS XCENTER YCENTER MAG MERR
PG1323-086 B 5:06:41.0 1.120 876.840 494.882 17.743 0.002
* V 5:09:25.0 1.116 875.257 495.078 17.589 0.002
* R 5:11:42.0 1.113 875.462 495.605 17.644 0.002
PG1323-086A B 5:06:41.0 1.120 498.825 658.040 18.368 0.002
* V 5:09:25.0 1.116 497.613 658.289 17.614 0.002
* R 5:11:42.0 1.113 497.582 658.737 17.374 0.002
PG1323-086B B 5:06:41.0 1.120 464.906 736.915 18.652 0.003
* V 5:09:25.0 1.116 463.816 737.116 17.494 0.002
* R 5:11:42.0 1.113 463.708 737.638 17.081 0.002
PG1323-086C B 5:06:41.0 1.120 478.510 397.435 19.195 0.003
* V 5:09:25.0 1.116 477.381 397.690 18.092 0.003
* R 5:11:42.0 1.113 477.447 398.285 17.712 0.002
PG1323-086 B 9:57:45.0 1.783 763.012 538.284 17.920 0.004
* V 10:00:24.0 1.809 764.468 538.493 17.671 0.005
* R 10:02:36.0 1.832 765.431 537.794 17.710 0.007
PG1323-086A B 9:57:45.0 1.783 385.423 701.446 18.885 0.009
* V 10:00:24.0 1.809 386.871 701.530 18.040 0.007
* R 10:02:36.0 1.832 387.662 701.040 17.747 0.007
PG1323-086B B 9:57:45.0 1.783 351.518 780.211 18.802 0.008
* V 10:00:24.0 1.809 353.137 780.487 17.574 0.004
* R 10:02:36.0 1.832 353.918 779.841 17.141 0.004
PG1323-086C B 9:57:45.0 1.783 364.956 440.529 19.331 0.013
* V 10:00:24.0 1.809 366.513 440.999 18.182 0.007
* R 10:02:36.0 1.832 367.472 440.458 17.769 0.007
PG1633+099 B 6:53:47.0 1.700 833.544 506.665 18.736 0.003
* V 6:57:06.0 1.679 831.482 506.568 18.561 0.003
* R 6:59:15.0 1.665 830.897 506.441 18.641 0.004
PG1633+099A B 6:53:47.0 1.700 761.205 498.716 20.749 0.010
* V 6:57:06.0 1.679 759.183 498.866 19.412 0.006
* R 6:59:15.0 1.665 758.542 498.577 18.892 0.005
PG1633+099B B 6:53:47.0 1.700 492.621 731.087 18.706 0.003
* V 6:57:06.0 1.679 490.940 731.193 17.112 0.001
* R 6:59:15.0 1.665 490.429 730.985 16.525 0.001
PG1633+099C B 6:53:47.0 1.700 346.427 744.418 19.022 0.003
* V 6:57:06.0 1.679 344.532 744.442 17.362 0.002
* R 6:59:15.0 1.665 344.161 744.359 16.752 0.001
PG1633+099D B 6:53:47.0 1.700 242.430 680.807 18.810 0.003
* V 6:57:06.0 1.679 240.543 680.881 17.838 0.002
* R 6:59:15.0 1.665 239.998 680.604 17.503 0.002
PG1633+099 B 9:49:07.0 1.311 799.543 457.162 18.642 0.003
* V 9:46:58.0 1.310 799.675 456.438 18.522 0.004
* R 9:44:28.0 1.308 799.470 457.579 18.605 0.004
PG1633+099A B 9:49:07.0 1.311 727.336 449.355 20.653 0.018
* V 9:46:58.0 1.310 727.596 448.561 19.376 0.007
* R 9:44:28.0 1.308 727.178 449.956 18.866 0.005
PG1633+099B B 9:49:07.0 1.311 459.042 681.546 18.613 0.003
* V 9:46:58.0 1.310 459.520 680.943 17.070 0.002
* R 9:44:28.0 1.308 458.961 682.439 16.494 0.001
PG1633+099C B 9:49:07.0 1.311 312.566 694.832 18.945 0.004
* V 9:46:58.0 1.310 313.292 694.175 17.324 0.002
* R 9:44:28.0 1.308 312.665 695.581 16.725 0.001
PG1633+099D B 9:49:07.0 1.311 208.503 631.457 18.720 0.004
* V 9:46:58.0 1.310 209.214 630.598 17.802 0.002
* R 9:44:28.0 1.308 208.575 632.208 17.475 0.002
RU_149 B 23:56:58.0 1.153 548.802 682.364 18.164 0.002
* V 23:54:38.0 1.154 549.630 683.102 17.989 0.003
* R 23:52:17.0 1.155 549.038 682.692 18.038 0.003
RU_149A B 23:56:58.0 1.153 591.762 654.463 19.265 0.004
* V 23:54:38.0 1.154 592.604 655.347 18.601 0.004
* R 23:52:17.0 1.155 592.005 654.824 18.410 0.004
RU_149B B 23:56:58.0 1.153 429.051 687.941 17.797 0.002
* V 23:54:38.0 1.154 430.239 688.796 16.744 0.001
* R 23:52:17.0 1.155 429.387 688.549 16.378 0.001
RU_149C B 23:56:58.0 1.153 438.097 587.360 19.083 0.004
* V 23:54:38.0 1.154 439.240 588.348 18.529 0.004
* R 23:52:17.0 1.155 438.439 587.741 18.433 0.004
RU_149D B 23:56:58.0 1.153 510.882 641.702 15.881 0.001
* V 23:54:38.0 1.154 511.836 642.571 15.594 0.001
* R 23:52:17.0 1.155 511.121 642.339 15.579 0.001
RU_149E B 23:56:58.0 1.153 396.030 419.737 18.704 0.003
* V 23:54:38.0 1.154 397.162 420.769 17.810 0.002
* R 23:52:17.0 1.155 396.350 420.236 17.502 0.002
RU_149F B 23:56:58.0 1.153 558.473 468.716 19.142 0.004
* V 23:54:38.0 1.154 559.396 469.806 17.564 0.002
* R 23:52:17.0 1.155 558.541 469.446 16.983 0.002
RU_149G B 23:56:58.0 1.153 639.221 516.747 17.846 0.002
* V 23:54:38.0 1.154 640.042 517.730 16.928 0.001
* R 23:52:17.0 1.155 639.440 517.384 16.616 0.001
RU_149 B 4:02:09.0 2.087 537.754 753.515 18.422 0.002
* V 4:04:49.0 2.124 539.488 753.419 18.106 0.003
* R 4:06:59.0 2.156 538.545 754.224 18.120 0.003
RU_149A B 4:02:09.0 2.087 580.853 725.587 19.511 0.005
* V 4:04:49.0 2.124 582.400 725.467 18.736 0.004
* R 4:06:59.0 2.156 581.529 726.483 18.520 0.004
RU_149B B 4:02:09.0 2.087 418.283 759.321 18.020 0.002
* V 4:04:49.0 2.124 419.809 758.966 16.861 0.001
* R 4:06:59.0 2.156 419.134 759.886 16.466 0.001
RU_149C B 4:02:09.0 2.087 427.166 658.523 19.320 0.004
* V 4:04:49.0 2.124 428.699 658.458 18.647 0.004
* R 4:06:59.0 2.156 427.947 659.284 18.513 0.004
RU_149D B 4:02:09.0 2.087 499.800 713.041 16.129 0.001
* V 4:04:49.0 2.124 501.508 712.838 15.715 0.001
* R 4:06:59.0 2.156 500.641 713.618 15.664 0.001
RU_149E B 4:02:09.0 2.087 385.055 491.016 18.948 0.003
* V 4:04:49.0 2.124 386.590 491.118 17.934 0.002
* R 4:06:59.0 2.156 385.920 491.773 17.586 0.002
RU_149F B 4:02:09.0 2.087 547.662 539.869 19.360 0.004
* V 4:04:49.0 2.124 549.061 540.176 17.680 0.002
* R 4:06:59.0 2.156 548.338 540.750 17.071 0.002
RU_149G B 4:02:09.0 2.087 628.401 588.104 18.079 0.002
* V 4:04:49.0 2.124 629.612 588.039 17.044 0.001
* R 4:06:59.0 2.156 628.901 588.695 16.696 0.001
# Weed out measurements that might be affected by saturated stellar cores, or
# cosmic rays within the measurement aperture. This step may require visual
# inspection:
cl> type standards.stat.tem
imexamine ("%.fits", 1, "", logfile="standards.stat", keeplog+,
defkey="a", allframes-, nframes=1, ncstat=27, nlstat=27,
imagecur="%.stdcoo")
cl> !\ls -1 *.stdcoo | sed 's:.stdcoo::g' > standards.stat.lis
cl> !narg standards.stat.tem standards.stat.lis '%' > standards.stat.cl
cl> cl < standards.stat.cl
--> none have saturated cores; highest peak value ~ 39k ADU
# Create, edit and verify a configurations file listing the transformation
# equations we are trying to solve. Since we have B, V, and R, we do not
# need to edit the standard configuration as extensively as we did in the
# example in class (where we only had B and R), but we do need to remove the
# entries for the U and I filters (make sure you know the basic 'vi' commands,
# as this is the default editor in 'mkconfig'):
cl> unlearn mkconfig
cl> mkconfig ("standards.cfg", "nlandolt", "standards.obs", "nlandolt",
>>> verify+, edit+, check+, verbose+)
# Declare the new Landolt UBVRI standards catalog variables
catalog
V 4 # the V magnitude
BV 5 # the (B-V) color
UB 6 # the (U-B) color
VR 7 # the (V-R) color
RI 8 # the (R-I) color
VI 9 # the (V-I) color
error(V) 12 # the V magnitude error
error(BV) 13 # the (B-V) color error
error(UB) 14 # the (U-B) color error
error(VR) 15 # the (V-R) color error
error(RI) 16 # the (R-I) color error
error(VI) 17 # the (V-I) color error
# Declare the observations file variables
observations
TB 3 # time of observation in filter B
XB 4 # airmass in filter B
xB 5 # x coordinate in filter B
yB 6 # y coordinate in filter B
mB 7 # instrumental magnitude in filter B
error(mB) 8 # magnitude error in filter B
TV 10 # time of observation in filter V
XV 11 # airmass in filter V
xV 12 # x coordinate in filter V
yV 13 # y coordinate in filter V
mV 14 # instrumental magnitude in filter V
error(mV) 15 # magnitude error in filter V
TR 17 # time of observation in filter R
XR 18 # airmass in filter R
xR 19 # x coordinate in filter R
yR 20 # y coordinate in filter R
mR 21 # instrumental magnitude in filter R
error(mR) 22 # magnitude error in filter R
# Sample transformation section for the new Landolt UBVRI system
transformation
fit b1=0.0, b2=0.35, b3=0.000
const b4=0.0
BFIT : mB = (BV + V) + b1 + b2 * XB + b3 * BV + b4 * BV * XB
fit v1=0.0, v2=0.17, v3=0.000
const v4=0.0
VFIT : mV = V + v1 + v2 * XV + v3 * BV + v4 * BV * XV
fit r1=0.0, r2=0.08, r3=0.000
const r4=0.0
RFIT : mR = (V - VR) + r1 + r2 * XR + r3 * VR + r4 * VR * XR
:wq
"standards.cfg" 65L, 2190C written
** Beginning of compilation **
** End of compilation **
CATALOG VARIABLES, COLUMNS, AND ERROR COLUMNS:
1 V 4 12
2 BV 5 13
3 UB 6 14
4 VR 7 15
5 RI 8 16
6 VI 9 17
OBSERVATIONAL VARIABLES, COLUMNS, AND ERROR COLUMNS:
1 TB 3 INDEF
2 XB 4 INDEF
3 xB 5 INDEF
4 yB 6 INDEF
5 mB 7 8
7 TV 10 INDEF
8 XV 11 INDEF
9 xV 12 INDEF
10 yV 13 INDEF
11 mV 14 15
13 TR 17 INDEF
14 XR 18 INDEF
15 xR 19 INDEF
16 yR 20 INDEF
17 mR 21 22
FIT AND CONSTANT PARAMETER VALUES:
1 b1 0.
2 b2 0.35
3 b3 0.
4 b4 0. (constant)
5 v1 0.
6 v2 0.17
7 v3 0.
8 v4 0. (constant)
9 r1 0.
10 r2 0.08
11 r3 0.
12 r4 0. (constant)
TRANSFORMATION EQUATIONS:
1 BFIT: mB = (BV+V)+b1+b2*XB+b3*BV+b4*BV*XB
delta(BFIT, b1) = 0.1
delta(BFIT, b2) = 0.1
delta(BFIT, b3) = 0.1
delta(BFIT, b4) = 0.1
error = , min = , max =
weight = , min = , max =
plot x = , y =
2 VFIT: mV = V+v1+v2*XV+v3*BV+v4*BV*XV
delta(VFIT, v1) = 0.1
delta(VFIT, v2) = 0.1
delta(VFIT, v3) = 0.1
delta(VFIT, v4) = 0.1
error = , min = , max =
weight = , min = , max =
plot x = , y =
3 RFIT: mR = (V-VR)+r1+r2*XR+r3*VR+r4*VR*XR
delta(RFIT, r1) = 0.1
delta(RFIT, r2) = 0.1
delta(RFIT, r3) = 0.1
delta(RFIT, r4) = 0.1
error = , min = , max =
weight = , min = , max =
plot x = , y =
Catalog input variables = 12
First catalog column = 4
Last catalog column = 17
Observational input variables = 18
First observational column = 3
Last observational column = 22
Fitting parameters = 9
Constant parameters = 3
Auxiliary (set) equations = 0
Transformation equations = 3
Warnings = 0
Errors = 0
--> standards.cfg Note that I deleted the two transformations for U and I.
# Interactively solve (fit) the transformation equations listed in file
# 'standards.cfg', one at a time. Check that the solutions indeed converge.
# On photometric nights, the rms scatter should be no more than ~0.03 mag
# after rejection of outliers, and the reduced χ² should be close to 1.
cl> unlearn fitparams
cl> fitparams ("standards.obs", "nlandolt", "standards.cfg", "standards.sol",
>>> weighting="photometric", nreject=5, logfile="photcal.log", log_fit+,
>>> log_results+)
:tolerance 1e-2 <-- make fit less restrictive; allow poor fit to convergence
h <-- plot function (observed B mag + mean offset) against fit (B mag)
t <-- toggle overlay with indication of fit (indicated by open boxes)
l <-- plot residuals with respect to the fit versus function
= <-- hardcopy of screen to default printer (optional)
# most of the variance appears to be with the B~19 mag stars but there may
# also be some non-linearity: the residuals for the B~16 stars are
# systematically positive by ~0.02 mag, while the B~20 stars show residuals
# that are systematically negative by ~0.02 mag.
d <-- deleted four of the most pronounced outliers 18 < B < 19.5 mag
f <-- refit
:tolerance 1e-3 <-- reduce tolerance somewhat to see if fit still converges
f <-- yes it does
k <-- plot the residuals versus time of observation, UT (denoted by TB)
=
# there does not seem to be a lot of difference in the variance at different
# times.
g j <-- redefine the "j" key to plot:
XB,residuals <--- residuals vs. Airmass (symbolically denoted by XB)
# no clear trend, although the residuals are systematically negative for
# most of the AM~1.7--1.8 observations (but not for the AM~2.2 ones)
j
=
g j <-- redefine the "j" key to plot:
BV,residuals <--- residuals vs. B-V color
j
d # no clear trend visible; deleted two points with the largest deviations
f # RMS is down to 0.01333 mag
=
g j <-- redefine the "j" key to plot:
er(mB),residuals <--- residuals vs. error in instrumental B magnitude
j
# no clear trend visible; errors are all <0.02 mag.
:errors <-- print an overview of the errors, and fitted parameters
#Wed 01:32:41 06-Dec-2006
#mB = (BV+V)+b1+b2*XB+b3*BV+b4*BV*XB
#Solution converged
low_reject 3.
high_reject 3.
nreject 5
grow 0.
tol 0.001
maxiter 15
niterations 10
total_points 34
rejected 0
deleted 6
standard deviation 0.01410886
reduced chi 1.005472
average error 0.01403208
average scatter 0.01107532
RMS 0.01333162
parameter value error
b1 4.1318784 0.0105783 (fit)
b2 0.2610954 0.0064962 (fit)
b3 0.0766843 0.0056174 (fit)
b4 0.0000000 0.0000000 (constant)
# since the reduced χ² ~ 1 and the rms<0.03 mag, and since the fitted
# values appear to make sense (errors on fit <10% of fitted value),
# let's accept this fit as a valid solution and move on to the next filter:
q
yes
next
# and similar commands for the V and for the R filter, but with TV and TR
# instead of TB, XV and XR instead of XB.
# Result are files 'standards.sol' and 'photcal.log'.
# Check in file 'photcal.log' that nothing is listed under the header
# "#UNMATCHED OBJECT" --> all stellar identifications were indeed found in
# the 'nlandolt' standard table. Check under the header "#RESULTS: BFIT",
# "#RESULTS: VFIT" and "#RESULTS: RFIT" whether there are any stars that are
# consistently bad and hence rejected from the fit (values of INDEF in each
# column). It turns out that in every fit, star PG1323-086A is deleted. In
# a footnote to his photometry table, Landolt notes that this star was found
# to be a **variable** star by Paul Schmidtke. No wonder it didn't match its
# supposed magnitude well!
cl> type standards.sol
# Wed 01:36:27 06-Dec-2006
begin BFIT
status 0 (Solution converged)
variance 1.990599E-4
stdeviation 0.01410886
avsqerror 1.968992E-4
averror 0.01403208
avsqscatter 1.226627E-4
avscatter 0.01107532
chisqr 1.010974
msq 1.777320E-4
rms 0.01333162
reference mB
fitting (BV+V)+b1+b2*XB+b3*BV+b4*BV*XB
weights photometric
parameters 4
b1 (fit)
b2 (fit)
b3 (fit)
b4 (constant)
derivatives 4
0.1
0.1
0.1
0.1
values 4
4.131878
0.2610954
0.07668426
0.
errors 4
0.01057834
0.00649625
0.005617434
0.
# Wed 01:42:04 06-Dec-2006
begin VFIT
status 0 (Solution converged)
variance 5.818197E-5
stdeviation 0.00762771
avsqerror 5.696101E-5
averror 0.007547251
avsqscatter 1.945486E-5
avscatter 0.004410767
chisqr 1.021435
msq 5.255145E-5
rms 0.007249238
reference mV
fitting V+v1+v2*XV+v3*BV+v4*BV*XV
weights photometric
parameters 4
v1 (fit)
v2 (fit)
v3 (fit)
v4 (constant)
derivatives 4
0.1
0.1
0.1
0.1
values 4
3.969304
0.1219603
-0.02661422
0.
errors 4
0.005747955
0.003516401
0.002802127
0.
# Wed 01:44:50 06-Dec-2006
begin RFIT
status 0 (Solution converged)
variance 6.000784E-5
stdeviation 0.007746472
avsqerror 5.797408E-5
averror 0.007614071
avsqscatter 2.445702E-5
avscatter 0.004945404
chisqr 1.035081
msq 5.400705E-5
rms 0.007348949
reference mR
fitting (V-VR)+r1+r2*XR+r3*VR+r4*VR*XR
weights photometric
parameters 4
r1 (fit)
r2 (fit)
r3 (fit)
r4 (constant)
derivatives 4
0.1
0.1
0.1
0.1
values 4
4.016019
0.08676858
-0.02758684
0.
errors 4
0.00562534
0.003344508
0.005497807
0.
# So our conclusion is that the night was indeed photometric and that the
# zeropoint offsets with respect to our assumed zeropoint of 26 mag are
# as given in 'standards.sol' (the absolute zeropoint to place our instrument
# magnitudes on the Landolt photometric system is (26 - offset) mag).
# Hence, our photometric solution for night1 is:
------------------------------------------------------------------------------
coefficient: B filter: V filter: R filter: units:
------------------------------------------------------------------------------
zeropoint 21.8681 +- 0.0106 22.0307 +- 0.0057 21.9840 +- 0.0056 [mag]
extinction 0.2611 +- 0.0065 0.1220 +- 0.0035 0.0868 +- 0.0033 [mag/AM]
color_term 0.0767 +- 0.0056 -0.0266 +- 0.0028 -0.0276 +- 0.0055 [mag/mag]
------------------------------------------------------------------------------
To photometrically calibrate any instrument magnitude measured during this
night, simply do:
m = m_instr + (zeropoint-26) + extinction*airmass + color_term*color
where color is (B-V) for the B and V filters, and (V-R) for the R filter.
(To apply this, one would either have to iterate [we don't know the correct
color a priori, only the difference of the two instrument magnitudes], or to
apply it *correctly*, without loss of precision or introduction of additional
correlated errors, one would follow the procedure given in Appendix B of
Jansen et al. 2000, ApJS 126, 271).
# Now repeat the above steps for nights 2 through 5 (no standards star fields
# were observed on nights 6--8):
# Night 2:
cl> cd ../990321/
cl> mkobsfile ("*.stdmag", "B,V,R", "standards.obs", imsets="standards.sets",
>>> minmagerr=0.001, shifts="", apercors="", tolerance=5., verbose+)
Observations file: standards.obs
Image set: PG0918+029 4 stars written to the observations file
Image set: PG1323-086 4 stars written to the observations file
Image set: PG1633+099 5 stars written to the observations file
Image set: RU_149 8 stars written to the observations file
Image set: RU_149 8 stars written to the observations file
--> reported number of stars matches that expected, so there were no shifts
between exposures in a given sequence.
# Edit file 'standards.obs' after creating a backup copy:
cl> copy standards.obs standards.obs.orig
# edit identifiers in the 'standards.obs' table, and delete the entries for
# star PG1323-086A (variable star).
# Weed out measurements that might be affected by saturated stellar cores, or
# cosmic rays within the measurement aperture. This step may require visual
# inspection:
cl> copy ../990320/standards.stat.tem .
cl> !\ls -1 *.stdcoo | sed 's:.stdcoo::g' > standards.stat.lis
cl> !narg standards.stat.tem standards.stat.lis '%' > standards.stat.cl
cl> cl < standards.stat.cl
--> star RU_149D has a core in frame b0250 (V-filter) that is very close to
saturation (peak ~ 65558 ADU). Replaced the V-filter entries in the first
sequence of entries for this star in 'standards.obs' with "INDEF"s, i.e.:
RU_149D B 2:37:53.0 1.431 568.464 553.424 15.959 0.001
* V INDEF INDEF INDEF INDEF INDEF INDEF
* R 2:43:14.0 1.456 569.753 555.098 15.610 0.001
# Copy the configurations file that lists the photometric transformation
# equations we are trying to solve from the night1 (we can, of course, create
# and edit it from scratch with 'mkconfig', but the filters are the same for
# all our nights:
cl> copy ../990320/standards.cfg .
# Interactively solve (fit) the transformation equations listed in file
# 'standards.cfg', one at a time. Check that the solutions indeed converge.
# Since we know from night1 that the default tolerance of 3e-5 may be a bit
# to restrictive, increase it to 1e-3 initially:
cl> fitparams ("standards.obs", "nlandolt", "standards.cfg", "standards.sol",
>>> weighting="photometric", tolerance=1.0E-3, nreject=5,
>>> logfile="photcal.log", log_fit+, log_results+)
h <-- plot function (observed B mag + mean offset) against fit (B mag)
# solution doesn't converge
d <-- deleted some apparent outliers
f <-- refit
# solution converges
t <-- toggle overlay with indication of fit (indicated by open boxes)
l <-- plot residuals with respect to the fit versus function
# as for the previous night, we detect a hint of non-linearity (residuals
# for the B~16 are positive by ~0.01, those for B~21 negative by ~0.02 mag.
k <-- plot the residuals versus time of observation, UT (denoted by TB)
=
# residuals for UT~2.6, 3.9 and 7.0 hours are larger than those for UT~4.8
# and UT~6.1 hours, but they are still <~0.02 mag.
g j <-- redefine the "j" key to plot:
XB,residuals <--- residuals vs. Airmass (symbolically denoted by XB)
# no clear trend discernable
j
=
g j <-- redefine the "j" key to plot:
BV,residuals <--- residuals vs. B-V color (symbolically denoted by BV)
j
d # no clear trend visible; deleted one more outlier
f # RMS is down to 0.007754 mag
=
g j <-- redefine the "j" key to plot:
er(mB),residuals <--- residuals vs. error in instrumental B magnitude
# no clear trend visible; errors are all <0.02 mag.
l <-- plot residuals versus function for current best fit
:tolerance 1e-4 <-- reduce tolerance a bit
f <-- still OK
=
:errors <-- print an overview of the errors, and fitted parameters
#Wed 18:23:27 06-Dec-2006
#mB = (BV+V)+b1+b2*XB+b3*BV+b4*BV*XB
#Solution converged
low_reject 3.
high_reject 3.
nreject 5
grow 0.
tol 1.000000E-4
maxiter 15
niterations 5
total_points 28
rejected 0
deleted 5
standard deviation 0.00831557
reduced chi 1.016289
average error 0.00818229
average scatter 0.00211959
RMS 0.00775431
parameter value error
b1 4.1338058 0.0072645 (fit)
b2 0.2663105 0.0044419 (fit)
b3 0.0731637 0.0030974 (fit)
b4 0.0000000 0.0000000 (constant)
# since the reduced Chi^2 ~ 1 and the rms<<0.03 mag, and since the fitted
# values appear to make sense (errors on fit <10% of fitted value), accept
# this fit as a valid solution and move on to the next filter:
q
yes
next
# and similar commands for the V and for the R filter, but with TV and TR
# instead of TB, XV and XR instead of XB.
HINT: if autoscaling stretches the y-axis window limits to include a wildly
deviant point (which likely was automatically rejected, and displayed as
a cross symbol), then one can adjust those window limits more appropriately
by placing the cursor on or near a valid data point and hitting the "w"
key followed by the "p" key ("Pan x and y axis about cursor"), and then
the "w" key followed by the "m" key ("Autoscale x axis"). This centers
the y-range on the valid data. Next, while placing the cursor over a
valid data point, hit "w y" any number of times to zoom in, in y only.
If you subsequently switch to a different plot and no data points appear
within the graph-window, hit "w a" to autoscale the axes in x and y.
# Result are files 'standards.sol' and 'photcal.log'.
# Check in file 'photcal.log' that nothing is listed under the header
# "#UNMATCHED OBJECT" --> all stellar identifications were indeed found in
# the 'nlandolt' standard table. Check under the header "#RESULTS: BFIT",
# "#RESULTS: VFIT" and "#RESULTS: RFIT" whether there are any stars that are
# consistently bad and hence rejected from the fit (values of INDEF in each
# column). It turns out that here in every fit, star PG1633+099 is deleted.
# It is not quite clear why, in this case (perhaps because this star is one
# of the bluer stars). In any case, this star was not as significant an
# outlier as was PG1323-086A in the first night.
# From file 'standards.sol' we find (given our assumed zeropoint of 26 mag):
------------------------------------------------------------------------------
coefficient: B filter: V filter: R filter: units:
------------------------------------------------------------------------------
zeropoint 21.8662 +- 0.0073 22.0381 +- 0.0065 21.9826 +- 0.0073 [mag]
extinction 0.2663 +- 0.0044 0.1315 +- 0.0038 0.0915 +- 0.0043 [mag/AM]
color_term 0.0732 +- 0.0031 -0.0327 +- 0.0028 -0.0362 +- 0.0060 [mag/mag]
------------------------------------------------------------------------------
# Dec 6 2006 [RAJ]:
#==================
# Night 3:
cl> cd ../990322/
cl> mkobsfile ("*.stdmag", "B,V,R", "standards.obs", imsets="standards.sets",
>>> minmagerr=0.001, shifts="", apercors="", tolerance=5., verbose+)
Observations file: standards.obs
Image set: PG0918+029 5 stars written to the observations file
Image set: PG1323-086 4 stars written to the observations file
Image set: PG1633+099 5 stars written to the observations file
Image set: PG1633+099 5 stars written to the observations file
Image set: RU_149 8 stars written to the observations file
Image set: RU_149 8 stars written to the observations file
--> reported number of stars matches that expected, so there were no shifts
between exposures in a given sequence.
# Edit file 'standards.obs' after creating a backup copy:
cl> copy standards.obs standards.obs.orig
# edit identifiers in the 'standards.obs' table, and delete the entries for
# star PG1323-086A (variable star).
# Weed out measurements that might be affected by saturated stellar cores, or
# cosmic rays within the measurement aperture. This step may require visual
# inspection:
cl> copy ../990320/standards.stat.tem .
cl> !\ls -1 *.stdcoo | sed 's:.stdcoo::g' > standards.stat.lis
cl> !narg standards.stat.tem standards.stat.lis '%' > standards.stat.cl
cl> cl < standards.stat.cl
--> star RU_149D has a core in frame b0405 (R-filter) that is very close to
saturation (peak ~ 64297 ADU). Kept it anyway -- let's see if it breaks
the fit...
# Copy the configurations file that lists the photometric transformation
# equations we are trying to solve from the night 1:
cl> copy ../990320/standards.cfg .
# Interactively solve (fit) the transformation equations listed in file
# 'standards.cfg', one at a time. Check that the solutions indeed converge.
# Since we know from night1 that the default tolerance of 3e-5 may be a bit
# too restrictive, increase it to 1e-3 initially:
cl> fitparams ("standards.obs", "nlandolt", "standards.cfg", "standards.sol",
>>> weighting="photometric", tolerance=1.0E-3, nreject=5,
>>> logfile="photcal.log", log_fit+, log_results+)
# keystrokes as for the previous nights...
# Result are files 'standards.sol' and 'photcal.log'.
# Check in file 'photcal.log' that nothing is listed under the header
# "#UNMATCHED OBJECT" --> all stellar identifications were indeed found in
# the 'nlandolt' standard table. Check under the header "#RESULTS: BFIT",
# "#RESULTS: VFIT" and "#RESULTS: RFIT" whether there are any stars that are
# consistently bad and hence rejected from the fit (values of INDEF in each
# column). It turns out that stars PG1633+099 and PG1633+099D were deleted
# most consistently. Again, not quite clear why.
# From file 'standards.sol' we find (given our assumed zeropoint of 26 mag):
------------------------------------------------------------------------------
coefficient: B filter: V filter: R filter: units:
------------------------------------------------------------------------------
zeropoint 21.8634 +- 0.0096 22.0273 +- 0.0054 21.9718 +- 0.0067 [mag]
extinction 0.2683 +- 0.0066 0.1278 +- 0.0037 0.0870 +- 0.0047 [mag/AM]
color_term 0.0822 +- 0.0036 -0.0270 +- 0.0021 -0.0271 +- 0.0056 [mag/mag]
------------------------------------------------------------------------------
# Night 4:
cl> cd ../990323/
cl> mkobsfile ("*.stdmag", "B,V,R", "standards.obs", imsets="standards.sets",
>>> minmagerr=0.001, shifts="", apercors="", tolerance=5., verbose+)
Observations file: standards.obs
Image set: PG0918+029 5 stars written to the observations file
Image set: PG1323-086 4 stars written to the observations file
Image set: PG1323-086 4 stars written to the observations file
Image set: PG1633+099 5 stars written to the observations file
Image set: PG1633+099 5 stars written to the observations file
Image set: RU_149 8 stars written to the observations file
--> reported number of stars matches that expected, so there were no shifts
between exposures in a given sequence.
# Edit file 'standards.obs' after creating a backup copy:
cl> copy standards.obs standards.obs.orig
# edit identifiers in the 'standards.obs' table, and delete the entries for
# star PG1323-086A (variable star).
# Weed out measurements that might be affected by saturated stellar cores, or
# cosmic rays within the measurement aperture. This step may require visual
# inspection:
cl> copy ../990320/standards.stat.tem .
cl> !\ls -1 *.stdcoo | sed 's:.stdcoo::g' > standards.stat.lis
cl> !narg standards.stat.tem standards.stat.lis '%' > standards.stat.cl
cl> cl < standards.stat.cl
--> none have saturated cores; highest peak value ~ 24k ADU
# Copy the configurations file that lists the photometric transformation
# equations we are trying to solve from the night 1:
cl> copy ../990320/standards.cfg .
# Interactively solve (fit) the transformation equations listed in file
# 'standards.cfg', one at a time. Check that the solutions indeed converge.
# Since we know from night1 that the default tolerance of 3e-5 may be a bit
# to restrictive, increase it to 1e-3 initially:
cl> fitparams ("standards.obs", "nlandolt", "standards.cfg", "standards.sol",
>>> weighting="photometric", tolerance=1.0E-3, nreject=5,
>>> logfile="photcal.log", log_fit+, log_results+)
# keystrokes as for the previous nights...
# Result are files 'standards.sol' and 'photcal.log'.
# Check in file 'photcal.log' that nothing is listed under the header
# "#UNMATCHED OBJECT" --> all stellar identifications were indeed found in
# the 'nlandolt' standard table. Check under the header "#RESULTS: BFIT",
# "#RESULTS: VFIT" and "#RESULTS: RFIT" whether there are any stars that are
# consistently bad and hence rejected from the fit (values of INDEF in each
# column). Since there were no V and R observations of Rubin149, these
# entries are all INDEF, but that's OK.
# From file 'standards.sol' we find (given our assumed zeropoint of 26 mag):
------------------------------------------------------------------------------
coefficient: B filter: V filter: R filter: units:
------------------------------------------------------------------------------
zeropoint 21.8439 +- 0.0148 22.0297 +- 0.0118 21.9878 +- 0.0150 [mag]
extinction 0.2474 +- 0.0093 0.1259 +- 0.0085 0.0969 +- 0.0103 [mag/AM]
color_term 0.0946 +- 0.0052 -0.0179 +- 0.0038 -0.0139 +- 0.0091 [mag/mag]
------------------------------------------------------------------------------
(Fitted coefficients have somewhat larger uncertainties, but scatter was still
within the range for likely photometric conditions; transparency of the
atmosphere may have been somewhat different from the previous nights...)
# Night 5:
cl> cd ../990324/
cl> mkobsfile ("*.stdmag", "B,V,R", "standards.obs", imsets="standards.sets",
>>> minmagerr=0.001, shifts="", apercors="", tolerance=5., verbose+)
Observations file: standards.obs
Image set: PG0918+029 5 stars written to the observations file
Image set: RU_149 8 stars written to the observations file
Image set: RU_149 8 stars written to the observations file
--> reported number of stars matches that expected, so there were no shifts
between exposures in a given sequence.
# Edit file 'standards.obs' after creating a backup copy:
cl> copy standards.obs standards.obs.orig
# edit identifiers in the 'standards.obs' table.
# Weed out measurements that might be affected by saturated stellar cores, or
# cosmic rays within the measurement aperture. This step may require visual
# inspection:
cl> copy ../990320/standards.stat.tem .
cl> !\ls -1 *.stdcoo | sed 's:.stdcoo::g' > standards.stat.lis
cl> !narg standards.stat.tem standards.stat.lis '%' > standards.stat.cl
cl> cl < standards.stat.cl
--> star RU_149D has a core in frames e0733 (V filter) and e0734 (R filter)
that is very close to saturation (peak > 65000 ADU). Replaced entries
with "INDEF"s:
RU_149D B 23:34:34.0 1.157 526.917 704.955 15.923 0.001
* V INDEF INDEF INDEF INDEF INDEF INDEF
* R INDEF INDEF INDEF INDEF INDEF INDEF
# Copy the configurations file that lists the photometric transformation
# equations we are trying to solve from the night 1:
cl> copy ../990320/standards.cfg .
# Interactively solve (fit) the transformation equations listed in file
# 'standards.cfg', one at a time. Check that the solutions indeed converge.
# Since we know from night1 that the default tolerance of 3e-5 may be a bit
# to restrictive, increase it to 1e-3 initially:
cl> fitparams ("standards.obs", "nlandolt", "standards.cfg", "standards.sol",
>>> weighting="photometric", tolerance=1.0E-3, nreject=5,
>>> logfile="photcal.log", log_fit+, log_results+)
# keystrokes as for the previous nights...
--> All observations at airmass ~ 1.28 (i.e., of field PG0918) had to be
deleted in order to obtain a sensible fit. This poses concern regarding
the photometricity of this night around UT ~ 3 hours: although by
excluding PG0918 a solution is obtained, a true solution should have
fitted both the PG0918 data and the Rubin149 data well, since they were
obtained very close in time to each another.
CONCLUSION: this was most probably a non-photometric night.
# Result are files 'standards.sol' and 'photcal.log'.
# From file 'standards.sol' we find (given our assumed zeropoint of 26 mag):
------------------------------------------------------------------------------
coefficient: B filter: V filter: R filter: units:
------------------------------------------------------------------------------
zeropoint 21.8790 +- 0.0146 22.0205 +- 0.0139 21.9582 +- 0.0134 [mag]
extinction 0.3179 +- 0.0104 0.1492 +- 0.0093 0.0951 +- 0.0088 [mag/AM]
color_term 0.0708 +- 0.0065 -0.0274 +- 0.0055 -0.0194 +- 0.0101 [mag/mag]
------------------------------------------------------------------------------
But this solution cannot be transferred to data taken at any pointing other
than of these two standards star fields themselves, because there may have
been clouds (probably no thicker than cirrus) around!
NOTE: One of the conclusions of this exercise is that the color terms of the
instrumental system are not negligible for the B filter, if the (B-V) color
of the science target is redder than ~ 0.01/~0.08 = +0.125 mag or bluer than
-0.125 mag; they are not negligible for the V filter, if the (B-V) color of
the science target is redder than ~ 0.01/~0.027 = 0.37 mag or bluer than
-0.37 mag; and they are not negligible for the R filter, if the (V-R) color
of the science target is redder than ~0.01/~0.02 = +0.5 mag or bluer than
-0.5 mag. From the table on the "class project" web-page, it is clear that
the reference stars are much redder than these limits.
# So let's measure the aperture magnitudes of science target and reference
# stars in each of the observations:
cl> hselect @sci.lis $I,object yes | match "X0925" | fields "-" 1 > tmp.lis
cl> !sort tmp.lis | sed 's:.fits::g' | sed 's: ::g' > targ.lis
# (the latter 'sed' strips the trailing space characters)
cl> delete tmp.lis yes
# Create pixel coordinate files for each frame for the science target and the
# three reference stars. We can use task 'markstds' for this:
cl> type gettargcoo.tem
markstds ("&1.fits", radstr="7,18,36,44", refstar="X0925")
fields ("&1.stdcoo", "3", > "&1.targcoo1")
ctrcoo ("&1.targcoo", "&1.targcoo2", "&1.fits",
cntrbox=7, redispl-, verbose-)
!paste -d" " &1.targcoo2 &1.targcoo1 | \
awk '{printf("%.2f %.2f \"%s\"\n",$1,$2,$3)}' > &1.targcoo
print ("CTRCOO: outab = &1.targcoo") ; type &1.targcoo
delete &1.stdcoo,&1.targcoo? yes verify-
cl> !narg gettargcoo.tem targ.lis '&1' > gettargcoo.cl
cl> type gettargcoo.cl <-- to check that the substitution went OK
# Have the target finding-chart at hand, then run the script we prepared. It
# may be convenient to flip either the finding chart or the DS9 view in "Y" to
# match the view.
cl> unlearn markstds
cl> cl < gettargcoo.cl
... (etc) ...
MARKSTDS: NOAO/IRAF2.12.2 raj@hedgehag Dec 07 01:08:08 2006
image = 990327/h1272.fits [1024 x 1024]
Wait for the cursor cross to appear in the image display area, then mark
all (standard) stars in the image using the 'c' key. Hit 'q' to quit.
Note: start with the main (i.e., naming or reference) star in the field.
MARKSTDS: Finished.
CTRCOO: outab = 990327/h1272.targcoo
492.80 498.47 "X0925"
521.74 428.63 "1"
453.36 567.43 "2"
482.22 349.25 "3"
# Check that each of the coordinate files has indeed 4 objects:
cl> count 990320/a*.targcoo --> OK
cl> count 990321/b*.targcoo ""
cl> count 990322/c*.targcoo ""
cl> count 990323/d*.targcoo ""
cl> count 990324/e*.targcoo --> OK
cl> count 990326/g*.targcoo --> not OK: two files have 5 and 3 objects!
...
5 15 94 990326/g1093.targcoo
...
3 9 58 990326/g1113.targcoo
...
cl> count 990327/h*.targcoo --> OK
# Run 'imexamine' and mark the 4 objects using the "a" key. This will perform
# a centroiding algorithm to produce listed coordinates:
cl> imexamine 990326/g1093.fits 1 allframes- nframes=1
# COL LINE COORDINATES
# R MAG FLUX SKY PEAK E PA BETA ENCLOSED MOFFAT DIRECT
550.32 494.79 550.32 494.79
12.61 14.01 24984. 480.6 938. 0.03 55 4.79 4.14 4.24 4.20
579.30 424.98 579.30 424.98
12.50 13.17 54036. 480.8 2057. 0.03 -89 5.04 4.12 4.26 4.16
510.86 563.73 510.86 563.73
12.48 12.99 63866. 481.8 2446. 0.01 83 3.94 4.12 4.24 4.16
539.72 345.62 539.72 345.62
12.23 13.87 28244. 481.5 1107. 0.06 -83 2.55 4.08 4.11 4.08
cl> imexamine 990326/g1113.fits 1 allframes- nframes=1
# COL LINE COORDINATES
# R MAG FLUX SKY PEAK E PA BETA ENCLOSED MOFFAT DIRECT
550.01 497.93 550.01 497.93
14.14 14.00 25027. 437.7 827. 0.25 52 7.50 4.72 4.76 4.71
579.04 428.09 579.04 428.09
14.39 13.19 52876. 438.7 1710. 0.08 83 3.00 4.66 4.81 4.81
510.62 566.89 510.62 566.89
14.36 13.02 62131. 438.8 2049. 0.04 84 11.1 4.60 4.85 4.79
539.52 348.72 539.52 348.72
14.29 13.88 27965. 438.4 925. 0.18 81 14.2 4.79 4.81 4.77
# Manually fix the two coordinate files using the coordinates above.
# During the marking of all 188 science target frames, it became obvious that
# quite a few frames will have cosmic ray induced signal deposited on or near
# the target objects. It would be good to remove that signal to attain the
# highest quality photometry and avoid spurious scatter in the measured
# magnitudes.
# This can be done by manually marking any pixels affected by cosmic rays near
# the target objects (e.g., with task 'markall' -> 'cosmic+', the rest '-')
# and subsequently interpolating over them (e.g., with task 'imclean',
# available from http://www.asu.edu/~rjansen/) or possibly with task 'imedit'.
# A better and automated way would be to use 'la_cosmic' (van Dokkum 2001;
# particularly the IDL implementation thereof by Josh Bloom).
# Here, I'll just go ahead and check later whether there is a large fraction
# of spurious measurements or other suspect results. To do so, and to be able
# to select the highest S/N matched apertures later, I will specify a series
# of aperture radii for measurement and worry later about which ones I'll
# actually use (that way we only have to make one pass through all the data).
# I keep the radii defining sky annulus the same as specified for the
# standards (they seemed to do the job just fine).
# But first, we must make sure that each FITS header contains keyword values
# for the exposure time, airmass, filter (we already made sure of that right
# at the start of our image processing), and universal time. If not, we must
# add/fix such header keywords/values:
cl> hselect @targ.lis $I,exptime,airmass,utshut,filter yes
990320/a0086 300.000 00:45:26.0 V <-- missing AIRMASS keyword
990320/a0087 450.000 00:52:21.0 B <-- missing AIRMASS keyword
990320/a0090 300.000 1.075 03:18:59.0 V
990320/a0091 450.000 1.082 03:27:11.0 B
990320/a0092 300.000 1.091 03:36:33.0 V
... ... ... ... ..
--> only two frames with problems. From the Observing Log, we can substitute
a value that will be close enough:
cl> hedit 990320/a0086 airmass 1.102 add+
cl> hedit 990320/a0087 airmass 1.094 add+
cl> type gettargphot.tem
qphot ("%.fits", cbox=5., annulus=36., dannulus=11.,
apertures="5.,6.,7.,8.,9.,10.,12.,14.,16.,18.,20.,22.,25.",
coords="%.targcoo", output="%.targmag", zmag=26.,
exposure="exptime", airmass="airmass", filter="filter",
obstime="utshut", epadu=2.955, interactive-, verbose+)
cl> unlearn qphot
cl> !narg gettargphot.tem targ.lis '%' > gettargphot.cl
cl> type gettargphot.cl <-- to check that the substitution went OK
cl> cl < gettargphot.cl
a0086.fits 580.835 531.615 113.7867 21.337 21.293 21.270 21.254 21.246 \
21.240 21.229 21.218 21.198 21.181 21.179 21.187 ok
... (etc) ...
h1272.fits 492.81 498.50 577.9795 21.379 21.313 21.277 21.259 21.252 \
21.248 21.246 21.242 21.230 21.213 21.210 21.210 ok
h1272.fits 521.76 428.56 578.4391 20.498 20.435 20.399 20.377 20.362 \
20.350 20.338 20.327 20.326 20.322 20.315 20.312 ok
h1272.fits 453.50 567.43 578.1895 20.309 20.241 20.204 20.184 20.173 \
20.164 20.153 20.151 20.154 20.155 20.152 20.149 ok
h1272.fits 482.39 349.23 579.1058 21.199 21.133 21.089 21.068 21.055 \
21.047 21.034 21.022 21.020 21.030 21.039 21.034 ok
.... todo: extract magnitudes per filter, compute relative magnitudes
wrt reference stars, calibrate on absolute scale (incl. color terms),
plot results.
|