make_pairs.py

File contents

#!/usr/bin/env python

'''This python script looks at the irp files in the current directory and
makes a first guess at the dflat_pairs, tflat_pairs, unc_pairs and std_pairs
files.'''

import sys,os,string,re
from glob import glob
from pyraf import iraf
from myiraf import get_header
from Numeric import *
import string

time_factor = 3.0

if len(sys.argv) > 1:
   if sys.argv[1] == "idf":
      filepat = "idf*.fits"
      filetemp = "idf_%05d.fits"
      filtn = 'FILTER2'
   elif sys.argv[1] == 'irx':
      filepat = "irx*c1_001.fits"
      filetemp = "irx_%05d_c1_001.fits"
      filtn = 'FILTER'
   else:
      print "unknown option"
      sys.exit(1)
else:
   filepat = "irp*_001.fits"
   filetemp = "irp_%05d_001.fits"
   filtn = 'FILTER2'

pat = re.compile(r'([0-9]+)T([0-9]+):([0-9]+):([0-9]+)')
def date2day(string):
   res = pat.search(string)
   if not res:
      return -1
   day,hour,minute,second = map(int, res.groups())
   return(day + hour/24.0 + minute/1440.0 + second/86400.0)

# Get a list of filters we need calibrations for:
result = iraf.hselect(filepat, filtn, 
      'OBSTYPE="astro"||OBSTYPE="standard"', Stdout=1)
u = {}
for x in result:  u[x] = 1
filters = u.keys()
print "Found astro and standards in following filters:", filters

# Do the domeflats:
result = iraf.hselect(filepat, "FILE", 'OBSTYPE="dflat"', Stdout=1)
if len(result) != 2:
   print "Sorry, didn't get only two dlfat pairs... you're going to have to do"
   print "it by hand or put a mask in the defaults file."
else:
   f = open('dflat_pairs','w')
   print >>f, result[0],result[1]
   f.close()

# Do the twilight flats:
f = open('tflat_pairs', 'w')
for filter in filters:
   result = iraf.hselect(filepat, "FILE", 
         'OBSTYPE="tflat"&&%s="%s"' % (filtn,filter), Stdout=1)
   if len(result) == 0:
      print "You don't have tflats for filter %s, add it to defaults file" \
            % (filter)
      continue

   result = map(int, result)
   result.sort()
   print >> f, result[0], result[-1]
f.close()

# Now do the darks:
result = iraf.hselect(filepat, "EXPTIME", 'OBSTYPE!="dark"', Stdout=1)
exptimes = {}.fromkeys(result).keys()
result = iraf.hselect(filepat, "EXPTIME", 'OBSTYPE=="dark"', Stdout=1)
for exptime in exptimes:
   if exptime not in result:
      print "No darks for exposure time %s s, add it to defaults file" % exptime

# Do the uncrowded pairs:
astros = iraf.hselect(filepat, "FILE,DITHER", 'OBSTYPE="astro"', Stdout=1)
astros = map(lambda x: string.split(x, '\t'), astros)
files = array(map(int,[astros[i][0] for i in range(len(astros))]))
dithers = array(map(int, [astros[i][1] for i in range(len(astros))]))
begins = nonzero(equal(dithers, 1))
astro1 = take(files, begins)
astro2 = take(files, begins[1:] - 1)
astro2 = concatenate([astro2, files[-1:]])
astro1 = astro1.tolist()
astro2 = astro2.tolist()
#astro1 = iraf.hselect(filepat, "FILE", 'OBSTYPE="astro"&&DITHER=1', 
#      Stdout=1)
#astro2 = iraf.hselect(filepat, "FILE", 
#      'OBSTYPE="astro"&&DITHER=NDITHERS', Stdout=1)
#astro1 = map(int, astro1)
#astro2 = map(int, astro2)
# make a dictionary of airmass/times:
astro_obst = {}
astro_airm = {}
astro_filt = {}
for i in astro1:
   astro_obst[i]=date2day(get_header(filetemp % i, "DATE-OBS", str))
   astro_airm[i] = get_header(filetemp % i, "AIRMASS", float)
   astro_filt[i] = get_header(filetemp % i, filtn, str)

if len(astro1) != len(astro2):
   print "Sorry, couldn't figure out unc_pairs, do it by hand"
else:
   f = open('unc_pairs','w')
   for i in range(len(astro1)):
      print >> f, astro1[i], astro2[i]
   f.close()

# Do standard pairs:
f = open('std_pairs','w')
standard1 = iraf.hselect(filepat, "FILE", 'OBSTYPE="standard"&&DITHER=1', 
      Stdout=1)
standard2 = iraf.hselect(filepat, "FILE", 
      'OBSTYPE="standard"&&DITHER=NDITHERS', Stdout=1)
standard1 = map(int, standard1)
standard2 = map(int, standard2)
if len(standard1) != len(standard2):
   print "Sorry, couldn't figure out std_pairs, do it by hand"
else:
   for i in range(len(standard1)):
      # Find the airmass and obstime
      #print "doing standard sequence %d-%d" % (standard1[i],standard2[i])
      thisfile = filetemp % (standard1[i])
      dday = date2day(get_header(thisfile, "DATE-OBS", str))
      airm = get_header(thisfile, "AIRMASS", float)
      filt = get_header(thisfile, filtn, str)
      #print "filter=%s, airmass=%.2f, obstime=%.2f" % (filtn,airm,dday)
      mind = -1.0
      minastro = -1
      for j in astro1:
         if astro_filt[j] == filt:
            dist = (dday - astro_obst[j])**2 + \
                  (airm - astro_airm[j])**2/time_factor**2
            #print "   seq:  %d:  airmass=%.2f  obstime=%.2f  d=%.4f" % \
            #      (j,astro_airm[j],astro_obst[j],dist)
            if mind < 0 or dist < mind:
               mind = dist
               minastro = j
      #print "   chose %d" % minastro
      if minastro == -1:
         print "Sorry:  coulnd't find astro pass for standard p = %d,%d" % \
            (standard1[i],standard2[i])
         print >> f, standard1[i],standard2[i],-1,-1
      else:
         index = astro1.index(minastro)
         print >> f, standard1[i],standard2[i],astro1[index],astro2[index]
   f.close()