moved hstwcs from general development repository to stsci_python/development

git-svn-id: http://svn.stsci.edu/svn/ssb/stsci_python/development/trunk/hstwcs@6931 fe389314-cf27-0410-b35b-8c050e845b92

3 files changed, 965 insertions, 0 deletions

diff --git a/distortion/__init__.py b/distortion/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/distortion/__init__.py
diff --git a/distortion/models.py b/distortion/models.py
new file mode 100644
index 0000000..cff9a46
--- /dev/null
+++ b/distortion/models.py
@@ -0,0 +1,360 @@
+import types
+
+# Import PyDrizzle utility modules
+import mutil
+import numpy as N
+import mutil
+from mutil import combin
+
+yes = True
+no = False
+
+#################
+#
+#
+#               Geometry/Distortion Classes
+#
+#
+#################
+
+class GeometryModel:
+    """
+    Base class for Distortion model.
+    There will be a separate class for each type of
+    model/filetype used with drizzle, i.e., IDCModel and
+    DrizzleModel.
+
+    Each class will know how to apply the distortion to a
+    single point and how to convert coefficients to an input table
+    suitable for the drizzle task.
+
+    Coefficients will be stored in CX,CY arrays.
+    """
+    #
+    #
+    #
+    #
+    #
+    #
+    #
+    NORDER = 3
+
+    def __init__(self):
+        "  This will open the given file and determine its type and norder."
+
+        #       Method to read in coefficients from given table and
+        #       populate the n arrays 'cx' and 'cy'.
+        #       This will be different for each type of input file,
+        #       IDCTAB vs. drizzle table.
+
+        # Set these up here for all sub-classes to use...
+        # But, calculate norder and cx,cy arrays in detector specific classes.
+        self.cx = None
+        self.cy = None
+        self.refpix = None
+        self.norder = self.NORDER
+        # Keep track of computed zero-point for distortion coeffs
+        self.x0 = None
+        self.y0 = None
+
+        # default values for these attributes
+        self.direction = 'forward'
+
+        self.pscale = 1.0
+
+    def shift(self,cx,cy,xs,ys):
+        """
+        Shift reference position of coefficients to new center
+        where (xs,ys) = old-reference-position - subarray/image center.
+        This will support creating coeffs files for drizzle which will
+        be applied relative to the center of the image, rather than relative
+        to the reference position of the chip.
+        """
+
+        _cxs = N.zeros(shape=cx.shape,dtype=cx.dtype)
+        _cys = N.zeros(shape=cy.shape,dtype=cy.dtype)
+        _k = self.norder + 1
+        # loop over each input coefficient
+        for m in xrange(_k):
+            for n in xrange(_k):
+                if m >= n:
+                    # For this coefficient, shift by xs/ys.
+                    _ilist = range(m, _k)
+                    # sum from m to k
+                    for i in _ilist:
+                        _jlist = range(n, i - (m-n)+1)
+                        # sum from n to i-(m-n)
+                        for j in _jlist:
+                            _cxs[m,n] += cx[i,j]*combin(j,n)*combin((i-j),(m-n))*pow(xs,(j-n))*pow(ys,((i-j)-(m-n)))
+                            _cys[m,n] += cy[i,j]*combin(j,n)*combin((i-j),(m-n))*pow(xs,(j-n))*pow(ys,((i-j)-(m-n)))
+        return _cxs,_cys
+
+    def convert(self, tmpname, xref=None,yref=None,delta=yes):
+        """
+         Open up an ASCII file, output coefficients in drizzle
+          format after converting them as necessary.
+        First, normalize these coefficients to what drizzle expects
+        Normalize the coefficients by the MODEL/output plate scale.
+
+        16-May-2002:
+        Revised to work with higher order polynomials by John Blakeslee.
+        27-June-2002:
+            Added ability to shift coefficients to new center for support
+                of subarrays.
+        """
+        cx = self.cx/self.pscale
+        cy = self.cy/self.pscale
+        x0 = self.refpix['XDELTA'] + cx[0,0]
+        y0 = self.refpix['YDELTA'] + cy[0,0]
+        #xr = self.refpix['XREF']
+        #yr = self.refpix['YREF']
+        xr = self.refpix['CHIP_XREF']
+        yr = self.refpix['CHIP_YREF']
+
+
+
+        '''
+        if xref != None:
+            # Shift coefficients for use with drizzle
+            _xs = xref - self.refpix['XREF'] + 1.0
+            _ys = yref - self.refpix['YREF'] + 1.0
+
+
+            if _xs != 0 or _ys != 0:
+                cxs,cys= self.shift(cx, cy, _xs, _ys)
+                cx = cxs
+                cy = cys
+
+                # We only want to apply this shift to coeffs
+                # for subarray images.
+                if delta == no:
+                    cxs[0,0] = cxs[0,0] - _xs
+                    cys[0,0] = cys[0,0] - _ys
+
+                # Now, apply only the difference introduced by the distortion..
+                # i.e., (undistorted - original) shift.
+                x0 += cxs[0,0]
+                y0 += cys[0,0]
+        '''
+        self.x0 = x0 #+ 1.0
+        self.y0 = y0 #+ 1.0
+
+        # Now, write out the coefficients into an ASCII
+        # file in 'drizzle' format.
+        lines = []
+
+
+        lines.append('# Polynomial distortion coefficients\n')
+        lines.append('# Extracted from "%s" \n'%self.name)
+        lines.append('refpix %f %f \n'%(xr,yr))
+        if self.norder==3:
+            lines.append('cubic\n')
+        elif self.norder==4:
+            lines.append('quartic\n')
+        elif self.norder==5:
+            lines.append('quintic\n')
+        else:
+            raise ValueError, "Drizzle cannot handle poly distortions of order %d"%self.norder
+
+        str = '%16.8f %16.8g %16.8g %16.8g %16.8g \n'% (x0,cx[1,1],cx[1,0],cx[2,2],cx[2,1])
+        lines.append(str)
+        str = '%16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cx[2,0],cx[3,3],cx[3,2],cx[3,1],cx[3,0])
+        lines.append(str)
+        if self.norder>3:
+            str = '%16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cx[4,4],cx[4,3],cx[4,2],cx[4,1],cx[4,0])
+            lines.append(str)
+        if self.norder>4:
+            str = '%16.8g %16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cx[5,5],cx[5,4],cx[5,3],cx[5,2],cx[5,1],cx[5,0])
+            lines.append(str)
+        lines.append("\n")
+
+        str = '%16.8f %16.8g %16.8g %16.8g %16.8g \n'% (y0,cy[1,1],cy[1,0],cy[2,2],cy[2,1])
+        lines.append(str)
+        str = '%16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cy[2,0],cy[3,3],cy[3,2],cy[3,1],cy[3,0])
+        lines.append(str)
+        if self.norder>3:
+            str = '%16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cy[4,4],cy[4,3],cy[4,2],cy[4,1],cy[4,0])
+            lines.append(str)
+        if self.norder>4:
+            str = '%16.8g %16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cy[5,5],cy[5,4],cy[5,3],cy[5,2],cy[5,1],cy[5,0])
+            lines.append(str)
+
+        output = open(tmpname,'w')
+        output.writelines(lines)
+        output.close()
+
+
+    def apply(self, pixpos,scale=1.0,order=None):
+        """
+         Apply coefficients to a pixel position or a list of positions.
+          This should be the same for all coefficients tables.
+        Return the geometrically-adjusted position
+        in arcseconds from the reference position as a tuple (x,y).
+
+        Compute delta from reference position
+        """
+        
+        """
+        scale actually is a ratio of pscale/self.model.pscale
+        what is pscale?
+        """
+        if self.cx == None:
+            return pixpos[:,0],pixpos[:,1]
+
+        if order is None:
+            order = self.norder
+
+        # Apply in the same way that 'drizzle' would...
+        _cx = self.cx / (self.pscale * scale)
+        _cy = self.cy / (self.pscale * scale)
+        _convert = no
+        _p = pixpos
+
+        # Do NOT include any zero-point terms in CX,CY here
+        # as they should not be scaled by plate-scale like rest
+        # of coeffs...  This makes the computations consistent
+        # with 'drizzle'.  WJH 17-Feb-2004
+        _cx[0,0] = 0.
+        _cy[0,0] = 0.
+
+        if isinstance(_p,types.ListType) or isinstance(_p,types.TupleType):
+            _p = N.array(_p,dtype=N.float64)
+            _convert = yes
+
+        dxy = _p - (self.refpix['XREF'],self.refpix['YREF'])
+        # Apply coefficients from distortion model here...
+        c = _p * 0.
+        for i in range(order+1):
+            for j in range(i+1):
+                c[:,0] = c[:,0] + _cx[i][j] * pow(dxy[:,0],j) * pow(dxy[:,1],(i-j))
+                c[:,1] = c[:,1] + _cy[i][j] * pow(dxy[:,0],j) * pow(dxy[:,1],(i-j))
+        xc = c[:,0]
+        yc = c[:,1]
+
+        # Convert results back to same form as original input
+        if _convert:
+            xc = xc.tolist()
+            yc = yc.tolist()
+            # If a single tuple was input, return just a single tuple
+            if len(xc) == 1:
+                xc = xc[0]
+                yc = yc[0]
+
+        return xc,yc
+
+    def setPScaleCoeffs(self,pscale):
+        self.cx[1,1] = pscale
+        self.cy[1,0] = pscale
+
+        self.refpix['PSCALE'] = pscale
+        self.pscale = pscale
+
+
+class IDCModel(GeometryModel):
+    """
+    This class will open the IDCTAB, select proper row based on
+    chip/direction and populate cx,cy arrays.
+    We also need to read in SCALE, XCOM,YCOM, XREF,YREF as well.
+    """
+    def __init__(self, idcfile, date=None, chip=1, direction='forward',
+                filter1='CLEAR1',filter2='CLEAR2',offtab=None,  binned=1):
+        GeometryModel.__init__(self)
+        #
+        # Norder must be derived from the coeffs file itself,
+        # then the arrays can be setup. Thus, it needs to be
+        # done in the sub-class, not in the base class.
+        # Read in table.
+        # Populate cx,cy,scale, and other variables here.
+        #
+        self.name = idcfile
+        self.cx,self.cy,self.refpix,self.norder = mutil.readIDCtab(idcfile,
+                        chip=chip,direction=direction,filter1=filter1,filter2=filter2,
+                date=date, offtab=offtab)
+
+        if self.refpix.has_key('empty_model') and self.refpix['empty_model']:
+            pass
+        else:
+            self.refpix['PSCALE'] = self.refpix['PSCALE'] * binned
+            self.cx = self.cx * binned
+            self.cy = self.cy * binned
+            self.refpix['XREF'] = self.refpix['XREF'] / binned
+            self.refpix['YREF'] = self.refpix['YREF'] / binned
+            self.refpix['XSIZE'] = self.refpix['XSIZE'] / binned
+            self.refpix['YSIZE'] = self.refpix['YSIZE'] / binned
+
+        self.pscale = self.refpix['PSCALE']
+        
+
+class WCSModel(GeometryModel):
+    """
+    This class sets up a distortion model based on coefficients
+    found in the image header.
+    """
+    def __init__(self,header,rootname):
+        GeometryModel.__init__(self)
+
+
+        if header.has_key('rootname'):
+            self.name = header['rootname']
+        else:
+            self.name = rootname
+        # Initialize all necessary distortion arrays with
+        # default model...
+        #self.cx,self.cy,self.refpix,self.order = mutil.defaultModel()
+
+        # Read in values from header, and update distortion arrays.
+        self.cx,self.cy,self.refpix,self.norder = mutil.readWCSCoeffs(header)
+
+        self.pscale = self.refpix['PSCALE']
+
+
+
+class DrizzleModel(GeometryModel):
+    """
+    This class will read in an ASCII Cubic
+    drizzle coeffs file and populate the cx,cy arrays.
+    """
+
+    def __init__(self, idcfile, scale = None):
+        GeometryModel.__init__(self)
+        #
+        # We now need to read in the file, populate cx,cy, and
+        # other variables as necessary.
+        #
+        self.name = idcfile
+        self.cx,self.cy,self.refpix,self.norder = mutil.readCubicTable(idcfile)
+
+        # scale is the ratio wcs.pscale/model.pscale. 
+        # model.pscale for WFPC2 is passed from REFDATA.
+        # This is needed for WFPC2 binned data.
+        
+        if scale != None:
+            self.pscale = scale
+        else:
+            self.pscale = self.refpix['PSCALE']
+
+        """
+        The above definition looks wrong.
+        In one case it's a ratio in the other it's pscale.
+        
+        """
+
+class TraugerModel(GeometryModel):
+    """
+    This class will read in the ASCII Trauger coeffs
+    file, convert them to SIAF coefficients, then populate
+    the cx,cy arrays.
+    """
+    NORDER = 3
+
+    def __init__(self, idcfile,lam):
+        GeometryModel.__init__(self)
+        self.name = idcfile
+        self.cx,self.cy,self.refpix,self.norder = mutil.readTraugerTable(idcfile,lam)
+        self.pscale = self.refpix['PSCALE']
+        #
+        # Read in file here.
+        # Populate cx,cy, and other variables.
+        #
+
+
diff --git a/distortion/mutil.py b/distortion/mutil.py
new file mode 100644
index 0000000..c3f83f2
--- /dev/null
+++ b/distortion/mutil.py
@@ -0,0 +1,605 @@
+from pytools import fileutil
+import numpy as N
+import string
+import calendar
+
+# Set up IRAF-compatible Boolean values    
+yes = True
+no = False
+
+# This function read the IDC table and generates the two matrices with
+# the geometric correction coefficients.
+#
+#       INPUT: FITS object of open IDC table
+#       OUTPUT: coefficient matrices for Fx and Fy
+#
+#### If 'tabname' == None: This should return a default, undistorted
+####                        solution.
+#
+
+def readIDCtab (tabname, chip=1, date=None, direction='forward',
+                filter1=None,filter2=None, offtab=None):
+
+    """
+        Read IDCTAB, and optional OFFTAB if sepcified, and generate
+        the two matrices with the geometric correction coefficients.
+
+        If tabname == None, then return a default, undistorted solution.
+        If offtab is specified, dateobs also needs to be given.
+
+    """
+  
+ # Return a default geometry model if no IDCTAB filename
+    # is given.  This model will not distort the data in any way.
+    if tabname == None:
+        print 'Warning: No IDCTAB specified! No distortion correction will be applied.'
+        return defaultModel()
+
+    # Implement default values for filters here to avoid the default
+    # being overwritten by values of None passed by user.
+    if filter1 == None or filter1.find('CLEAR') == 0:
+        filter1 = 'CLEAR'
+    if filter2 == None or filter2.find('CLEAR') == 0:
+        filter2 = 'CLEAR'
+
+    # Insure that tabname is full filename with fully expanded
+    # IRAF variables; i.e. 'jref$mc41442gj_idc.fits' should get
+    # expanded to '/data/cdbs7/jref/mc41442gj_idc.fits' before
+    # being used here.
+    # Open up IDC table now...
+    try:
+        ftab = fileutil.openImage(tabname)
+    except:
+        err_str =  "------------------------------------------------------------------------ \n"
+        err_str += "WARNING: the IDCTAB geometric distortion file specified in the image     \n"
+        err_str += "header was not found on disk. Please verify that your environment        \n"
+        err_str += "variable ('jref'/'uref'/'oref'/'nref') has been correctly defined. If    \n"
+        err_str += "you do not have the IDCTAB file, you may obtain the latest version       \n"
+        err_str += "of it from the relevant instrument page on the STScI HST website:        \n"
+        err_str += "http://www.stsci.edu/hst/ For WFPC2, STIS and NICMOS data, the           \n"
+        err_str += "present run will continue using the old coefficients provided in         \n"
+        err_str += "the Dither Package (ca. 1995-1998).                                      \n"
+        err_str += "------------------------------------------------------------------------ \n"
+        raise IOError,err_str
+
+    #First thing we need, is to read in the coefficients from the IDC
+    # table and populate the Fx and Fy matrices.
+
+    if ftab['PRIMARY'].header.has_key('DETECTOR'):
+        detector = ftab['PRIMARY'].header['DETECTOR']
+    else:
+        if ftab['PRIMARY'].header.has_key('CAMERA'):
+            detector = str(ftab['PRIMARY'].header['CAMERA'])
+        else:
+            detector = 1
+
+    # Set default filters for SBC
+    if detector == 'SBC':
+        if filter1 == 'CLEAR':
+            filter1 = 'F115LP'
+            filter2 = 'N/A'
+        if filter2 == 'CLEAR':
+            filter2 = 'N/A'
+
+    # Read FITS header to determine order of fit, i.e. k
+    norder = ftab['PRIMARY'].header['NORDER']
+    if norder < 3:
+        order = 3
+    else:
+        order = norder
+
+    fx = N.zeros(shape=(order+1,order+1),dtype=N.float64)
+    fy = N.zeros(shape=(order+1,order+1),dtype=N.float64)
+
+    #Determine row from which to get the coefficients.
+    # How many rows do we have in the table...
+    fshape = ftab[1].data.shape
+    colnames = ftab[1].data._names
+    row = -1
+
+    # Loop over all the rows looking for the one which corresponds
+    # to the value of CCDCHIP we are working on...
+    for i in xrange(fshape[0]):
+
+        try:
+            # Match FILTER combo to appropriate row,
+            #if there is a filter column in the IDCTAB...
+            if 'FILTER1' in colnames and 'FILTER2' in colnames:
+
+                filt1 = ftab[1].data.field('FILTER1')[i]
+                if filt1.find('CLEAR') > -1: filt1 = filt1[:5]
+
+                filt2 = ftab[1].data.field('FILTER2')[i]
+                if filt2.find('CLEAR') > -1: filt2 = filt2[:5]
+            else:
+                if 'OPT_ELEM' in colnames:
+                    filt1 = ftab[1].data.field('OPT_ELEM')
+                    if filt1.find('CLEAR') > -1: filt1 = filt1[:5]
+                else:
+                    filt1 = filter1
+
+                if 'FILTER' in colnames:
+                    _filt = ftab[1].data.field('FILTER')[i]
+                    if _filt.find('CLEAR') > -1: _filt = _filt[:5]
+                    if 'OPT_ELEM' in colnames:
+                        filt2 = _filt
+                    else:
+                        filt1 = _filt
+                        filt2 = 'CLEAR'
+                else:
+                    filt2 = filter2
+        except:
+            # Otherwise assume all rows apply and compare to input filters...
+            filt1 = filter1
+            filt2 = filter2
+
+        if 'DETCHIP' in colnames:
+            detchip = ftab[1].data.field('DETCHIP')[i]
+            if not str(detchip).isdigit():
+                detchip = 1
+        else:
+            detchip = 1
+
+        if 'DIRECTION' in colnames:
+            direct = string.strip(string.lower(ftab[1].data.field('DIRECTION')[i]))
+        else:
+            direct = 'forward'
+
+        if filt1 == filter1.strip() and filt2 == filter2.strip():
+            if direct == direction.strip():
+                if int(detchip) == int(chip) or int(detchip) == -999:
+                    row = i
+                    break
+    if row < 0:
+        err_str = '\nProblem finding row in IDCTAB! Could not find row matching:\n'
+        err_str += '        CHIP: '+str(detchip)+'\n'
+        err_str += '     FILTERS: '+filter1+','+filter2+'\n'
+        ftab.close()
+        del ftab
+        raise LookupError,err_str
+    else:
+        print '- IDCTAB: Distortion model from row',str(row+1),'for chip',detchip,':',filter1.strip(),'and',filter2.strip()
+
+    # Read in V2REF and V3REF: this can either come from current table,
+    # or from an OFFTAB if time-dependent (i.e., for WFPC2)
+    theta = None
+    if 'V2REF' in colnames:
+        v2ref = ftab[1].data.field('V2REF')[row]
+        v3ref = ftab[1].data.field('V3REF')[row]
+    else:
+        # Read V2REF/V3REF from offset table (OFFTAB)
+        if offtab:
+            v2ref,v3ref,theta = readOfftab(offtab, date, chip=detchip)
+        else:
+            v2ref = 0.0
+            v3ref = 0.0
+
+    if theta == None:
+        if 'THETA' in colnames:
+            theta = ftab[1].data.field('THETA')[row]
+        else:
+            theta = 0.0
+
+    refpix = {}
+    refpix['XREF'] = ftab[1].data.field('XREF')[row]
+    refpix['YREF'] = ftab[1].data.field('YREF')[row]
+    refpix['XSIZE'] = ftab[1].data.field('XSIZE')[row]
+    refpix['YSIZE'] = ftab[1].data.field('YSIZE')[row]
+    refpix['PSCALE'] = round(ftab[1].data.field('SCALE')[row],8)
+    refpix['V2REF'] = v2ref
+    refpix['V3REF'] = v3ref
+    refpix['THETA'] = theta
+    refpix['XDELTA'] = 0.0
+    refpix['YDELTA'] = 0.0
+    refpix['DEFAULT_SCALE'] = yes
+    refpix['centered'] = no
+
+    # Now that we know which row to look at, read coefficients into the
+    #   numeric arrays we have set up...
+    # Setup which column name convention the IDCTAB follows
+    # either: A,B or CX,CY
+    if 'CX10' in ftab[1].data._names:
+        cxstr = 'CX'
+        cystr = 'CY'
+    else:
+        cxstr = 'A'
+        cystr = 'B'
+
+    for i in xrange(norder+1):
+        if i > 0:
+            for j in xrange(i+1):
+                xcname = cxstr+str(i)+str(j)
+                ycname = cystr+str(i)+str(j)
+                fx[i,j] = ftab[1].data.field(xcname)[row]
+                fy[i,j] = ftab[1].data.field(ycname)[row]
+
+    ftab.close()
+    del ftab
+
+    # If CX11 is 1.0 and not equal to the PSCALE, then the
+    # coeffs need to be scaled
+    
+    if fx[1,1] == 1.0 and abs(fx[1,1]) != refpix['PSCALE']:
+        fx *= refpix['PSCALE']
+        fy *= refpix['PSCALE']
+
+    # Return arrays and polynomial order read in from table.
+    # NOTE: XREF and YREF are stored in Fx,Fy arrays respectively.
+    return fx,fy,refpix,order
+
+def readOfftab(offtab, date, chip=None):
+
+ 
+#Read V2REF,V3REF from a specified offset table (OFFTAB). 
+# Return a default geometry model if no IDCTAB filenam  e
+# is given.  This model will not distort the data in any way.
+
+    if offtab == None:
+        return 0.,0.
+
+    # Provide a default value for chip
+    if chip:
+        detchip = chip
+    else:
+        detchip = 1
+
+    # Open up IDC table now...
+    try:
+        ftab = fileutil.openImage(offtab)
+    except:
+        raise IOError,"Offset table '%s' not valid as specified!" % offtab
+
+    #Determine row from which to get the coefficients.
+    # How many rows do we have in the table...
+    fshape = ftab[1].data.shape
+    colnames = ftab[1].data._names
+    row = -1
+
+    row_start = None
+    row_end = None
+
+    v2end = None
+    v3end = None
+    date_end = None
+    theta_end = None
+
+    num_date = convertDate(date)
+    # Loop over all the rows looking for the one which corresponds
+    # to the value of CCDCHIP we are working on...
+    for ri in xrange(fshape[0]):
+        i = fshape[0] - ri - 1
+        if 'DETCHIP' in colnames:
+            detchip = ftab[1].data.field('DETCHIP')[i]
+        else:
+            detchip = 1
+
+        obsdate = convertDate(ftab[1].data.field('OBSDATE')[i])
+
+        # If the row is appropriate for the chip...
+            # Interpolate between dates
+        if int(detchip) == int(chip) or int(detchip) == -999:
+            if num_date <= obsdate:
+                date_end = obsdate
+                v2end = ftab[1].data.field('V2REF')[i]
+                v3end = ftab[1].data.field('V3REF')[i]
+                theta_end = ftab[1].data.field('THETA')[i]
+                row_end = i
+                continue
+
+            if row_end == None and (num_date > obsdate):
+                date_end = obsdate
+                v2end = ftab[1].data.field('V2REF')[i]
+                v3end = ftab[1].data.field('V3REF')[i]
+                theta_end = ftab[1].data.field('THETA')[i]
+                row_end = i
+                continue
+
+            if num_date > obsdate:
+                date_start = obsdate
+                v2start = ftab[1].data.field('V2REF')[i]
+                v3start = ftab[1].data.field('V3REF')[i]
+                theta_start = ftab[1].data.field('THETA')[i]
+                row_start = i
+                break
+
+    ftab.close()
+    del ftab
+
+    if row_start == None and row_end == None:
+        print 'Row corresponding to DETCHIP of ',detchip,' was not found!'
+        raise LookupError
+    elif row_start == None:
+        print '- OFFTAB: Offset defined by row',str(row_end+1)
+    else:
+        print '- OFFTAB: Offset interpolated from rows',str(row_start+1),'and',str(row_end+1)
+
+    # Now, do the interpolation for v2ref, v3ref, and theta
+    if row_start == None or row_end == row_start:
+        # We are processing an observation taken after the last calibration
+        date_start = date_end
+        v2start = v2end
+        v3start = v3end
+        _fraction = 0.
+        theta_start = theta_end
+    else:
+        _fraction = float((num_date - date_start)) / float((date_end - date_start))
+
+    v2ref = _fraction * (v2end - v2start) + v2start
+    v3ref = _fraction * (v3end - v3start) + v3start
+    theta = _fraction * (theta_end - theta_start) + theta_start
+
+    return v2ref,v3ref,theta
+
+def readWCSCoeffs(header):
+   
+    #Read distortion coeffs from WCS header keywords and
+    #populate distortion coeffs arrays.
+    
+    # Read in order for polynomials
+    _xorder = header['a_order']
+    _yorder = header['b_order']
+    order = max(max(_xorder,_yorder),3)
+
+    fx = N.zeros(shape=(order+1,order+1),dtype=N.float64)
+    fy = N.zeros(shape=(order+1,order+1),dtype=N.float64)
+
+    # Read in CD matrix
+    _cd11 = header['cd1_1']
+    _cd12 = header['cd1_2']
+    _cd21 = header['cd2_1']
+    _cd22 = header['cd2_2']
+    _cdmat = N.array([[_cd11,_cd12],[_cd21,_cd22]])
+    _theta = N.arctan2(-_cd12,_cd22)
+    _rotmat = N.array([[N.cos(_theta),N.sin(_theta)],
+                      [-N.sin(_theta),N.cos(_theta)]])
+    _rCD = N.dot(_rotmat,_cdmat)
+    _skew = N.arcsin(-_rCD[1][0] / _rCD[0][0])
+    _scale = _rCD[0][0] * N.cos(_skew) * 3600.
+    _scale2 = _rCD[1][1] * 3600.
+
+    # Set up refpix
+    refpix = {}
+    refpix['XREF'] = header['crpix1']
+    refpix['YREF'] = header['crpix2']
+    refpix['XSIZE'] = header['naxis1']
+    refpix['YSIZE'] = header['naxis2']
+    refpix['PSCALE'] = _scale
+    refpix['V2REF'] = 0.
+    refpix['V3REF'] = 0.
+    refpix['THETA'] = RADTODEG(_theta)
+    refpix['XDELTA'] = 0.0
+    refpix['YDELTA'] = 0.0
+    refpix['DEFAULT_SCALE'] = yes
+    refpix['centered'] = yes
+
+
+    # Set up template for coeffs keyword names
+    cxstr = 'A_'
+    cystr = 'B_'
+    # Read coeffs into their own matrix
+    for i in xrange(_xorder+1):
+        for j in xrange(i+1):
+            xcname = cxstr+str(j)+'_'+str(i-j)
+            if header.has_key(xcname):
+                fx[i,j] = header[xcname]
+
+    # Extract Y coeffs separately as a different order may
+    # have been used to fit it.
+    for i in xrange(_yorder+1):
+        for j in xrange(i+1):
+            ycname = cystr+str(j)+'_'+str(i-j)
+            if header.has_key(ycname):
+                fy[i,j] = header[ycname]
+
+    # Now set the linear terms
+    fx[0][0] = 1.0
+    fy[0][0] = 1.0
+
+    return fx,fy,refpix,order
+
+
+def readTraugerTable(idcfile,wavelength):
+
+    # Return a default geometry model if no coefficients filename
+    # is given.  This model will not distort the data in any way.
+    if idcfile == None:
+        return fileutil.defaultModel()
+
+    # Trauger coefficients only result in a cubic file...
+    order = 3
+    numco = 10
+    a_coeffs = [0] * numco
+    b_coeffs = [0] * numco
+    indx = _MgF2(wavelength)
+
+    ifile = open(idcfile,'r')
+    # Search for the first line of the coefficients
+    _line = fileutil.rAsciiLine(ifile)
+    while string.lower(_line[:7]) != 'trauger':
+        _line = fileutil.rAsciiLine(ifile)
+    # Read in each row of coefficients,split them into their values,
+    # and convert them into cubic coefficients based on
+    # index of refraction value for the given wavelength
+    # Build X coefficients from first 10 rows of Trauger coefficients
+    j = 0
+    while j < 20:
+        _line = fileutil.rAsciiLine(ifile)
+        if _line == '': continue
+        _lc = string.split(_line)
+        if j < 10:
+            a_coeffs[j] = float(_lc[0])+float(_lc[1])*(indx-1.5)+float(_lc[2])*(indx-1.5)**2
+        else:
+            b_coeffs[j-10] = float(_lc[0])+float(_lc[1])*(indx-1.5)+float(_lc[2])*(indx-1.5)**2
+        j = j + 1
+
+    ifile.close()
+    del ifile
+
+    # Now, convert the coefficients into a Numeric array
+    # with the right coefficients in the right place.
+    # Populate output values now...
+    fx = N.zeros(shape=(order+1,order+1),dtype=N.float64)
+    fy = N.zeros(shape=(order+1,order+1),dtype=N.float64)
+    # Assign the coefficients to their array positions
+    fx[0,0] = 0.
+    fx[1] = N.array([a_coeffs[2],a_coeffs[1],0.,0.],dtype=N.float64)
+    fx[2] = N.array([a_coeffs[5],a_coeffs[4],a_coeffs[3],0.],dtype=N.float64)
+    fx[3] = N.array([a_coeffs[9],a_coeffs[8],a_coeffs[7],a_coeffs[6]],dtype=N.float64)
+    fy[0,0] = 0.
+    fy[1] = N.array([b_coeffs[2],b_coeffs[1],0.,0.],dtype=N.float64)
+    fy[2] = N.array([b_coeffs[5],b_coeffs[4],b_coeffs[3],0.],dtype=N.float64)
+    fy[3] = N.array([b_coeffs[9],b_coeffs[8],b_coeffs[7],b_coeffs[6]],dtype=N.float64)
+
+    # Used in Pattern.computeOffsets()
+    refpix = {}
+    refpix['XREF'] = None
+    refpix['YREF'] = None
+    refpix['V2REF'] = None
+    refpix['V3REF'] = None
+    refpix['XDELTA'] = 0.
+    refpix['YDELTA'] = 0.
+    refpix['PSCALE'] = None
+    refpix['DEFAULT_SCALE'] = no
+    refpix['centered'] = yes
+
+    return fx,fy,refpix,order
+
+
+def readCubicTable(idcfile):
+    # Assumption: this will only be used for cubic file...
+    order = 3
+    # Also, this function does NOT perform any scaling on
+    # the coefficients, it simply passes along what is found
+    # in the file as is...
+
+    # Return a default geometry model if no coefficients filename
+    # is given.  This model will not distort the data in any way.
+    if idcfile == None:
+        return fileutil.defaultModel()
+
+    ifile = open(idcfile,'r')
+    # Search for the first line of the coefficients
+    _line = fileutil.rAsciiLine(ifile)
+
+    _found = no
+    while _found == no:
+        if _line[:7] in  ['cubic','quartic','quintic'] or _line[:4] == 'poly':
+            found = yes
+            break
+        _line = fileutil.rAsciiLine(ifile)
+
+    # Read in each row of coefficients, without line breaks or newlines
+    # split them into their values, and create a list for A coefficients
+    # and another list for the B coefficients
+    _line = fileutil.rAsciiLine(ifile)
+    a_coeffs = string.split(_line)
+
+    x0 = float(a_coeffs[0])
+    _line = fileutil.rAsciiLine(ifile)
+    a_coeffs[len(a_coeffs):] = string.split(_line)
+    # Scale coefficients for use within PyDrizzle
+    for i in range(len(a_coeffs)):
+        a_coeffs[i] = float(a_coeffs[i])
+
+    _line = fileutil.rAsciiLine(ifile)
+    b_coeffs = string.split(_line)
+    y0 = float(b_coeffs[0])
+    _line = fileutil.rAsciiLine(ifile)
+    b_coeffs[len(b_coeffs):] = string.split(_line)
+    # Scale coefficients for use within PyDrizzle
+    for i in range(len(b_coeffs)):
+        b_coeffs[i] = float(b_coeffs[i])
+
+    ifile.close()
+    del ifile
+    # Now, convert the coefficients into a Numeric array
+    # with the right coefficients in the right place.
+    # Populate output values now...
+    fx = N.zeros(shape=(order+1,order+1),dtype=N.float64)
+    fy = N.zeros(shape=(order+1,order+1),dtype=N.float64)
+    # Assign the coefficients to their array positions
+    fx[0,0] = 0.
+    fx[1] = N.array([a_coeffs[2],a_coeffs[1],0.,0.],dtype=N.float64)
+    fx[2] = N.array([a_coeffs[5],a_coeffs[4],a_coeffs[3],0.],dtype=N.float64)
+    fx[3] = N.array([a_coeffs[9],a_coeffs[8],a_coeffs[7],a_coeffs[6]],dtype=N.float64)
+    fy[0,0] = 0.
+    fy[1] = N.array([b_coeffs[2],b_coeffs[1],0.,0.],dtype=N.float64)
+    fy[2] = N.array([b_coeffs[5],b_coeffs[4],b_coeffs[3],0.],dtype=N.float64)
+    fy[3] = N.array([b_coeffs[9],b_coeffs[8],b_coeffs[7],b_coeffs[6]],dtype=N.float64)
+
+    # Used in Pattern.computeOffsets()
+    refpix = {}
+    refpix['XREF'] = None
+    refpix['YREF'] = None
+    refpix['V2REF'] = x0
+    refpix['V3REF'] = y0
+    refpix['XDELTA'] = 0.
+    refpix['YDELTA'] = 0.
+    refpix['PSCALE'] = None
+    refpix['DEFAULT_SCALE'] = no
+    refpix['centered'] = yes
+
+    return fx,fy,refpix,order
+
+def factorial(n):
+    """ Compute a factorial for integer n. """
+    m = 1
+    for i in range(int(n)):
+        m = m * (i+1)
+    return m
+
+def combin(j,n):
+    """ Return the combinatorial factor for j in n."""
+    return (factorial(j) / (factorial(n) * factorial( (j-n) ) ) )
+
+
+def defaultModel():
+    """ This function returns a default, non-distorting model
+        that can be used with the data.
+    """
+    order = 3
+
+    fx = N.zeros(shape=(order+1,order+1),dtype=N.float64)
+    fy = N.zeros(shape=(order+1,order+1),dtype=N.float64)
+
+    fx[1,1] = 1.
+    fy[1,0] = 1.
+
+    # Used in Pattern.computeOffsets()
+    refpix = {}
+    refpix['empty_model'] = yes
+    refpix['XREF'] = None
+    refpix['YREF'] = None
+    refpix['V2REF'] = 0.
+    refpix['XSIZE'] = 0.
+    refpix['YSIZE'] = 0.
+    refpix['V3REF'] = 0.
+    refpix['XDELTA'] = 0.
+    refpix['YDELTA'] = 0.
+    refpix['PSCALE'] = None
+    refpix['DEFAULT_SCALE'] = no
+    refpix['THETA'] = 0.
+    refpix['centered'] = yes
+    return fx,fy,refpix,order
+
+# Function to compute the index of refraction for MgF2 at
+# the specified wavelength for use with Trauger coefficients
+def _MgF2(lam):
+    _sig = pow((1.0e7/lam),2)
+    return N.sqrt(1.0 + 2.590355e10/(5.312993e10-_sig) +
+        4.4543708e9/(11.17083e9-_sig) + 4.0838897e5/(1.766361e5-_sig))
+
+
+def convertDate(date):
+    """ Converts the DATE-OBS date string into an integer of the
+        number of seconds since 1970.0 using calendar.timegm().
+
+        INPUT: DATE-OBS in format of 'YYYY-MM-DD'.
+        OUTPUT: Date (integer) in seconds.
+    """
+
+    _dates = date.split('-')
+    _val = 0
+    _date_tuple = (int(_dates[0]), int(_dates[1]), int(_dates[2]), 0, 0, 0, 0, 0, 0)
+
+    return calendar.timegm(_date_tuple)