CCD data reduction

Basic data reduction using IRAF

Below is a bash script that creates lists of all images: biases, flats, darks and targets. Special case of dark frames of 600sec that Milankovic generates is used Last lines create ascii files with median values in flat frames that will be used for their normalization.

The script creates reduce.cl IRAF script that should be run from iraf:
cl < reduce.cl

 #!/bin/bash

rm *list
for frame in $(ls *fit); do
	filterType=$(gethead filter $frame)
	imageType=$(gethead imagetyp $frame)
	outFileName=$(echo $imageType | tr '[:upper:]' '[:lower:]')
	# test if DARK
	if [ ${imageType} = "DARK" ]; then
 		exptime=$(gethead exptime $frame | cut -d. -f1) 
 		ls $frame >> ${outFileName}${exptime}"s".list 
		echo "b"$frame >> "b"${outFileName}${exptime}"s".list
	else ls $frame >> ${outFileName}${filterType}.list
	fi
	# test if LIGHT and create additional lists
	if [ ${imageType} = "LIGHT" ]; then
	 	exptime=$(gethead exptime $frame | cut -d. -f1)
		echo $exptime >> exptime${filterType}.list
		echo "b"$frame >> "b"${outFileName}${filterType}.list
		echo "bd"$frame >> "bd"${outFileName}${filterType}.list
		echo "bdf"$frame >> "bdf"${outFileName}${filterType}.list
	 else if [ ${imageType} = "FLAT" ]; then
                echo "b"$frame >> "b"${outFileName}${filterType}.list
                echo "bd"$frame >> "bd"${outFileName}${filterType}.list
		echo "nbd"$frame >> "nbd"${outFileName}${filterType}.list
	      fi
	 fi
done

# Create dark current corresponding to each light frame using exposure time
for exptime in $(ls exptime*.list); do 
	filter=$(echo $exptime | sed -e s/exptime//g -e s/\.list//g);
	awk -v filter=${filter} '{print "dc"NR"_"filter".fits"}' exptime${filter}.list  > dc${filter}.list
done

if [ -e reduce.cl ]; then rm reduce.cl; fi
echo "imdelete master_bias,dc_600s,dc_5s,dc_1s" >> reduce.cl
echo "imcombine @bias.list master_bias comb=med" >> reduce.cl
echo "imarith @dark600s.list - master_bias @bdark600s.list" >> reduce.cl
echo "imarith @dark5s.list - master_bias @bdark5s.list" >> reduce.cl
echo "imcombine @bdark5s.list dc_5s comb=med" >> reduce.cl
echo "imcombine @bdark600s.list dc_600s comb=med" >> reduce.cl 
echo "imarith dc_600s / 600 dc_1s" >> reduce.cl

for flats in $(ls flat*list); do
        filterType=$(echo $flats | sed s/flat//g | cut -d. -f1)
        #echo $filterType
	echo "imarith @flat${filterType}.list - master_bias @bflat${filterType}.list" >> reduce.cl
	echo "imarith @bflat${filterType}.list - dc_5s @bdflat${filterType}.list" >> reduce.cl
	echo "imstat @bdflat${filterType}.list fields=\"image,midpt\" format- > norm${filterType}.txt" >> reduce.cl
done

Flat field normalization

From terminal (bash):

  
for norm in $(ls norm*txt); do while read line; do echo $line | awk '! /^#/{print "imarith "$1" / "$2,"n"$1}'; done < $norm ; done > doflatnorm.cl

From terminal (IRAF):

cl < doflatnorm.cl

Calibration of science frames

#!/bin/bash/

if [ -e reduce_images.cl ]; then rm reduce_images.cl; fi

for flat_list in $(ls nbdflat*list); do
	filter=$(echo ${flat_list} | sed -e s/nbdflat//g -e s/\.list//g);
	echo "imcomb @nbdflat${filter}.list master_flat_${filter} comb=med" >> reduce_images.cl
	echo "imarith @light${filter}.list - master_bias @blight${filter}.list" >> reduce_images.cl
	echo "imarith dc_1s * @exptime${filter}.list @dc${filter}.list" >> reduce_images.cl	
	echo "imarith @blight${filter}.list - @dc${filter}.list @bdlight${filter}.list" >> reduce_images.cl
	echo "imarith @bdlight${filter}.list / master_flat_${filter} @bdflight${filter}.list" >> reduce_images.cl

done

cat blight*list bdlight*list bflat*list bdflat*list nbdflat*list bdark*list dc*.list > tmp_files.list

while read line; do rm $line; done < tmp_files.list

This script creates IRAF's script reduce_images.cl that should be started from IRAF's terminal. At the same time do some cleaning (removes intermediate files that are no longer necessary).

BPM & object mask creation

# cat bpmask.cl
imstat dc_600s field="image,mean,stddev" for- > hot.txt
imstat master_flat_FILTER fields="image,mean" for- > cold.txt
cl < bpmask.cl
# awk -v filter="L" '{print "imexpr \"((a > b) ? 1 : 0)\" hot"filter".pl a="$1" b="($2+3.*$3)*0.75}' hot.txt > hot.cl
# awk -v filter="L" '{print "imexpr \"((a < b) ? 1 : 0)\" cold"filter".pl a="$1" b="$2*0.9}' cold.txt > cold.cl
cl < hot.cl
cl < cold.cl
# cat bpm.cl
imexpr "max(a,b)" bpmFILTER.pl a=hotFILTER.pl b=coldFILTER.pl
hedit @bdfFILTER.list BPM bpmFILTER.pl add+ ver- 
# sed 's/FILTER/'${filter}'/g'  bpm.cl > dobpm.cl
cl < dobpm.cl
  
# Object masks creation
nproto
objmask @bdfFILTER.list @objbdfFILTER.list masks="!BPM" hsig=2 lsig=2

Astrometry

# cat astrometry.sh
gethead OBJCTRA @bdfFILTER.list > name_ra.txt 
gethead OBJCTDEC @bdfFILTER.list | awk '{print $2,$3,$4}' > dec.txt
paste -d" " name_ra.txt dec.txt | awk '{print $1,($2+$3/60+$4/3600)*15,$5+$6/60+$7/3600}' > astrometryFILTER.txt
# sed s/FILTER/${filter}/g astrometry.sh  > astrometry${filter}.sh

sh astrometry${filter}.sh

 # This could be used:
 # cat komanda.sh
 solve-field img --ra imagera --dec imagedec --radius 0.2 --scale-units arcsecperpix --scale-low 0.38 --scale-high 0.40 --crpix-center --wcs img.wcs 
 ### The next command creates new files instead of just wcs binaries that can later be added to the header. It's a stand-alone command. And should be ignored.
 ### solve-field img --ra imagera --dec imagedec --radius 0.2 --scale-units arcsecperpix --scale-low 0.38 --scale-high 0.40 --crpix-center -p -N wcsimg.fit  \ 
 ###  -O -I noneI.fits -M none.fits -R none.fits -B none.fits -P none -k none -U none -y -S none --axy noneaxy 
 #while read line; do img=$(echo $line | awk '{print $1}'); imagera=$(echo $line | awk '{print $2}'); imagedec=$(echo $line | awk '{print $3}'); sed -e 's/img/'$img'/g' \ 
-e 's/imagera/'${imagera}'/g' -e 's/imagedec/'${imagedec}'/g' komanda.sh; done < astrometry${filter}.txt > doAstrometry.sh 
# Or this:
 awk '{print "solve-field "$1" --ra "$2" --dec "$3" --radius 0.2 --scale-units arcsecperpix --scale-low 0.38 --scale-high 0.40 --crpix-center -p -N wcs"$1" \ 
 -O -I noneI.fits -M none.fits -R none.fits -B none.fits -P none -k none -U none -y -S none --axy noneaxy --wcs bla.wcs"}' astrometry${filter}.txt > doAstrometry.sh

# /bin/bash doAstrometry.sh
 #ls *wcs > wcs${filter}.list
 # echo "wcscopy @bdf${filter}.list @wcs${filter}.list ver-" > dowcs.cl
cl < dowcs.cl

Supersky flat

# Creation of super sky flat, normalization to each individual image median background and final stacking
# cat stat.cl
real mini
real sky
imstat("img",fields="min",lower=INDEF,upper=INDEF,for-) | scanf("%g", mini)
imcalc(input="img,msk",output="skyimage",equals="if(im2==0.) then im1 else "//mini-0.001,ver-)
imstat("skyimage",fields="midpt",lower=mini-0.001,upper=INDEF,for-) | scanf("%g", sky)
print("img"," ",sky)
imarith("img","/",sky,"norm")
imdelete("skyimage",yes)
# while read line; do img=$(echo $line | awk '{print $1}'); msk="obj"${img%.*}.pl; rsky="norm"${img}; cat stat.cl | sed -e 's/img/'${img}'/g' \
-e 's/msk/'${msk}'/g'  -e 's/norm/'${rsky}'/g'; done < bdf${filter}.list >  bla.cl
cl < bla.cl > med_values.txt
# ls normbdf*fit > norm.list
imcomb @norm.list norm.fits comb=med masktype="!OBJMASK" maskval=0
# awk '{print "imarith norm.fits * "$2," n"$1}' med_values.txt > donorm.cl
cl < donorm.cl
# awk '{print "imarith",$1,"-","n"$1,"s"$1}' med_values.txt > do.cl
cl < do.cl
   
imcomb s@bdfFILTER.list final.fits comb=med offset="wcs"

Surface brightness limit

1. Calibrate the image taking non-saturated stars and not too faint stars: intercept and slope (TOPCAT with SDSS DR12 from Vizier)
2. Calculate statistics inside at least 20 (empty, i.e. without objects) boxes of the size of 10" (for example if 1pix=0.396" than 10"~25 pix so boxes = 25 pix × 25 pix)
3. Take the mean value of the standard deviation inside boxes and it's standard deviation
4. Calculate SBmin and SBmax as: 
   SBmax =-2.5*log10(<σ>-stdev(<σ>))*slope+intercept 
   SBmin =-2.5*log10(<σ>+stdev(<σ>))*slope+intercept 
where <σ> is the average standard deviation inside boxes and stdev(<σ>) its standard deviation.
5. Finally, 1 sigma surface brightness limit (SBlim) is given as:
   SBlim = (SBmin+SBmax)/2 +/- (SBmax-SBmin)/2