CCD data reduction

Basic data reduction

# 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

 #!/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

BPM & object mask creation

 # cat bpmask.cl
 imstat dc_300s 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