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

13 files changed, 2337 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)
diff --git a/lib/__init__.py b/lib/__init__.py
new file mode 100644
index 0000000..08e421b
--- /dev/null
+++ b/lib/__init__.py
@@ -0,0 +1,47 @@
+import wcsutil
+from  wcsutil import HSTWCS
+from pytools import fileutil, parseinput
+import pyfits
+
+#to avoid relative imports
+import mappings
+from mappings import basic_wcs
+import distortion
+import pywcs
+
+def restoreWCS(fnames):
+    """
+    Given a list of fits file names restore the original basic WCS kw
+    and write out the files. This overwrites the original files.
+    """
+    files = parseinput.parseinput(fnames)[0]
+    for f in files:
+        isfits, ftype = fileutil.isFits(f)
+        if not isfits or (isfits and ftype == 'waiver'):
+            print "RestoreWCS works only on multiextension fits files."
+            return
+        else:
+            fobj = pyfits.open(f, mode='update')
+            hdr0 = fobj[0].header
+            for ext in fobj:
+                try:
+                    extname = ext.header['EXTNAME'].lower()
+                except KeyError:
+                    continue
+                if extname in ['sci', 'err', 'sdq']:
+                    hdr = ext.header
+                    owcs = HSTWCS(hdr0, hdr)
+                    #Changed restore to update the attributes ???
+                    # this may need to be changed
+                    backup = owcs.restore()
+                    if not backup:
+                        print '\No archived keywords found.\n'
+                        continue
+                    else:
+                        for kw in basic_wcs:
+                            nkw = ('O'+kw)[:7]
+                            if backup.has_key(kw):
+                                hdr.update(kw, hdr[nkw])
+            fobj.close()
+            
+                    
diff --git a/lib/mappings.py b/lib/mappings.py
new file mode 100644
index 0000000..fed9b60
--- /dev/null
+++ b/lib/mappings.py
@@ -0,0 +1,38 @@
+from pytools import fileutil
+
+# This dictionary maps an instrument into an instrument class
+# The instrument class handles instrument specific keywords
+
+inst_mappings={'WFPC2': 'WFPC2WCS',
+                'ACS': 'ACSWCS'
+                }
+
+DEGTORAD = fileutil.DEGTORAD
+RADTODEG = fileutil.RADTODEG
+
+# A dictionary which lists the allowed corrections for each instrument.
+# These are the default corrections applied also in the pipeline.
+#Dgeo correction is applied separately.
+allowed_corrections={'WFPC2': ['MakeWCS','CompSIP'],
+                    'ACS': ['TDDCorr', 'MakeWCS','CompSIP', 'VACorr']
+                    }
+
+# A list of instrument specific keywords
+# Every instrument class must have methods which define each of these
+# as class attributes.
+ins_spec_kw = ['ltv1', 'ltv2', 'parity', 'binned','vafactor', 'chip', 
+            'naxis1', 'naxis2', 'filter1', 'filter2']
+
+# A list of keywords defined in the primary header.
+# The HSTWCS class sets this as attributes 
+prim_hdr_kw = ['detector', 'offtab', 'idctab', 'date-obs', 
+              'pa_v3', 'ra_targ', 'dec_targ']
+
+# These are the keywords which are archived before MakeWCS is run
+basic_wcs = ['CD1_1', 'CD1_2', 'CD2_1', 'CD2_2', 'CRVAL1','CRVAL2','CTYPE1', 'CTYPE2',
+            'CRPIX1', 'CRPIX2', 'CTYPE1', 'CTYPE2', 'ORIENTAT']
+            
+dgeo_vals = {'ACS':
+                {'naxis1':65, 'naxis2':33, 'extver':1, 'crpix1':33.5, 'crpix2':16.5, 
+                'cdelt1':1., 'cdelt2':1., 'crval1':2048., 'crval2':1024.}
+            }
\ No newline at end of file
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..cc150f1
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,23 @@
+from distutils.core import setup
+from os.path import join
+# PyFITS
+try:
+    import pyfits
+except ImportError:
+    print "WARNING: PyFITS must be installed to use hstwcs."
+    print "         Since this is not a build-time dependency, the build will proceed."
+
+# PyWCS
+try:
+    import pywcs
+except ImportError:
+    print "WARNING: PyWCS must be installed to use hstwcs."
+    print "         Since this is not a build-time dependency, the build will proceed."
+
+setup(name="hstwcs",
+      version="0.1",
+      description="HST WCS Corrections",
+      packages=['hstwcs', 'hstwcs/updatewcs', 'hstwcs/wcsutil', 'hstwcs/distortion'],
+      package_dir={'hstwcs':'lib', 'hstwcs/updatewcs': 'updatewcs',
+                    'hstwcs/wcsutil': 'wcsutil', 'hstwcs/distortion': 'distortion'}
+      )
diff --git a/setup.py.old b/setup.py.old
new file mode 100644
index 0000000..a4a57a2
--- /dev/null
+++ b/setup.py.old
@@ -0,0 +1,89 @@
+from distutils.core import setup, Extension
+from os.path import join
+
+######################################################################
+# NUMPY
+try:
+    import numpy
+except ImportError:
+    print "numpy must be installed to build pywcs."
+    print "ABORTING."
+    sys.exit(1)
+
+try:
+    numpy_include = numpy.get_include()
+except AttributeError:
+    numpy_include = numpy.get_numpy_include()
+
+######################################################################
+# PyFITS
+try:
+    import pyfits
+except ImportError:
+    print "WARNING: PyFITS must be installed to use pywcs."
+    print "         Since this is not a build-time dependency, the build will proceed."
+
+######################################################################
+######################################################################
+# WCSLIB
+WCSVERSION = "4.3"
+WCSLIB = "wcslib-%s" % WCSVERSION # Path to wcslib
+WCSLIBC = join('pywcs', WCSLIB, "C") # Path to wcslib source files
+WCSFILES = [ # List of wcslib files to compile
+    'flexed/wcsulex.c',
+    'flexed/wcsutrn.c',
+    'cel.c',
+    'lin.c',
+    'log.c',
+    'prj.c',
+    'spc.c',
+    'sph.c',
+    'spx.c',
+    'tab.c',
+    'wcs.c',
+    'wcsfix.c',
+    'wcshdr.c',
+    'wcspih.c',
+    'wcstrig.c',
+    'wcsunits.c',
+    'wcsutil.c']
+WCSFILES = [join(WCSLIBC, x) for x in WCSFILES]
+######################################################################
+
+######################################################################
+# GENERATE DOCSTRINGS IN C
+docstrings = {}
+execfile("pywcs/doc/docstrings.py", docstrings)
+keys = docstrings.keys()
+keys.sort()
+fd = open("pywcs/src/docstrings.h", "w")
+fd.write('/* This file is autogenerated by setup.py.  To edit its contents\n')
+fd.write('   edit doc/docstrings.py\n')
+fd.write('*/\n')
+for key in keys:
+    if key.startswith('__'):
+        continue
+    val = docstrings[key].lstrip().encode("string_escape").replace('"', '\\"')
+    fd.write("static const char doc_%s[] = \"%s\";\n\n" % (key, val))
+fd.close()
+
+######################################################################
+
+setup(name="hstwcs",
+      version="0.1",
+      description="HST WCS Corrections",
+      packages=['hstwcs', 'hstwcs/pywcs', 'hstwcs/updatewcs', 'hstwcs/wcsutil', 'hstwcs/distortion'],
+      package_dir={'hstwcs':'lib', 'hstwcs/pywcs': 'pywcs/pywcs', 'hstwcs/updatewcs': 'updatewcs',
+                    'hstwcs/wcsutil': 'wcsutil', 'hstwcs/distortion': 'distortion'},
+      ext_modules=[
+        Extension('hstwcs.pywcs._pywcs',
+                  ['pywcs/src/pywcs.c'] +
+                  WCSFILES,
+                  include_dirs=[
+                    numpy_include,
+                    WCSLIBC,
+                    "src"
+                    ]
+                  )
+     ]
+      )
diff --git a/updatewcs/__init__.py b/updatewcs/__init__.py
new file mode 100644
index 0000000..db89d13
--- /dev/null
+++ b/updatewcs/__init__.py
@@ -0,0 +1,262 @@
+import os
+import pyfits
+#from .. wcsutil import HSTWCS
+from hstwcs.wcsutil import HSTWCS
+from hstwcs.mappings import allowed_corrections
+#from .. mappings import allowed_corrections
+import corrections, makewcs
+import dgeo
+import time
+from pytools import parseinput, fileutil
+
+#NB! the order of corrections matters
+
+__docformat__ = 'restructuredtext'
+
+
+def updatewcs(input, vacorr=True, tddcorr=True, checkfiles=True):
+    """
+    Purpose
+    =======
+    Applies corrections to the WCS keywords.
+    
+    Example
+    =======
+    >>>from hstwcs import updatewcs
+    >>>updatewcs.updatewcs(filename)
+    
+    Dependencies 
+    ============
+    `pytools`
+    `pyfits`
+    `pywcs`
+    `numpy`
+
+    :Parameters:
+    `input`: a python list of file names or a string (wild card characters allowed)
+             input files may be in fits, geis or waiver fits format
+    `vacorr`: boolean
+              If True, vecocity aberration correction will be applied
+    `tddcorr`: boolean
+              If True, time dependen correction will be applied to the distortion model
+    `checkfiles`: boolean
+              If True, the format of the input files will be checked,
+              geis and waiver fits files will be converted to MEF format.
+              Default value is True for standalone mode.
+    """
+    
+    files = parseinput.parseinput(input)[0]
+    if checkfiles:
+        files = checkFiles(files)
+        if not files:
+            print 'No valid input, quitting ...\n'
+            return
+    for f in files:
+        instr = pyfits.getval(f, 'INSTRUME')
+        try:
+            acorr = setCorrections(instr,vacorr=vacorr, tddcorr=tddcorr)
+        except KeyError:
+            print 'Unsupported instrument %s ' %instr
+            print 'Removing %s from list of processed files\n' % f
+            files.remove(f)
+            continue
+        
+        makecorr(f, acorr)
+    return files
+
+def makecorr(fname, acorr):
+    """
+    Purpose
+    =======
+    Applies corrections to a single file
+    
+    :Parameters:
+    `fname`: string
+             file name
+    `acorr`: list
+             list of corrections to be applied
+            
+    """
+    f = pyfits.open(fname, mode='update')
+    nrefchip, nrefext = getNrefchip(f)
+    primhdr = f[0].header
+    
+    
+    for extn in f:
+        refwcs = HSTWCS(primhdr, f[nrefext].header)
+        refwcs.archive_kw()
+        refwcs.readModel()
+        if extn.header.has_key('extname') and extn.header['extname'].lower() == 'sci':
+            hdr = extn.header
+            owcs = HSTWCS(primhdr, hdr)
+            owcs.archive_kw()
+            owcs.readModel()
+            for c in acorr:
+                owcs.__setattr__('DO'+c, 'PERFORM')
+                corr_klass = corrections.__getattribute__(c)
+                corr_klass(owcs, refwcs)
+            
+            
+    
+    #always do dgeo correction
+    if applyDgeoCorr(fname):
+        dgeo.DGEO(f)
+    f.close()
+
+        
+def setCorrections(instrument, vacorr=True, tddcorr=True):
+    """
+    Purpose
+    =======
+    Creates a list of corrections to be applied to a file.
+    based on user input paramters and allowed corrections
+    for the instrument, which are defined in mappings.py.
+    """
+    acorr = allowed_corrections[instrument]
+    if 'VACorr' in acorr and not vacorr:  acorr.remove('VACorr')
+    if 'TDDCorr' in acorr and not tddcorr: acorr.remove('TDDCorr')
+    if 'DGEOCorr' in acorr and not dgeocorr: acorr.remove('DGEOCorr')
+
+    return acorr
+
+
+
+def applyDgeoCorr(fname):
+    """
+    Purpose
+    =======
+    Adds dgeo extensions to files based on the DGEOFILE keyword in the primary 
+    header. This is a default correction and will always run in the pipeline.
+    The file used to generate the extensions is 
+    recorded in the DGEOFILE keyword in each science extension.
+    If 'DGEOFILE' in the primary header is different from 'DGEOFILE' in the 
+    extension header and the file exists on disk and is a 'new type' dgeofile, 
+    then the dgeo extensions will be updated.
+    """
+    applyDGEOCorr = True
+    try:
+        # get DGEOFILE kw from primary header
+        fdgeo0 = pyfits.getval(fname, 'DGEOFILE')
+        fdgeo0 = fileutil.osfn(fdgeo0)
+        if not fileutil.findFile(fdgeo0):
+            print 'Kw DGEOFILE exists in primary header but file %s not found\n' % fdgeo0
+            print 'DGEO correction will not be applied\n'
+            applyDGEOCorr = False
+            return applyDGEOCorr 
+        try:
+            # get DGEOFILE kw from first extension header
+            fdgeo1 = pyfits.getval(fname, 'DGEOFILE', ext=1)
+            fdgeo1 = fileutil.osfn(fdgeo1)
+            if fdgeo1 and fileutil.findFile(fdgeo1):
+                if fdgeo0 != fdgeo1:
+                    applyDGEOCorr = True
+                else:
+                    applyDGEOCorr = False
+            else: 
+                # dgeo file defined in first extension may not be found
+                # but if a valid kw exists in the primary header, dgeo should be applied.
+                applyDGEOCorr = True
+        except KeyError:
+            # the case of DGEOFILE kw present in primary header but missing 
+            # in first extension header
+            applyDGEOCorr = True
+    except KeyError:
+
+        print 'DGEOFILE keyword not found in primary header'
+        applyDGEOCorr = False
+    
+    if isOldStyleDGEO(fname, fdgeo0):
+        applyDGEOCorr = False
+        
+    return applyDGEOCorr
+
+def isOldStyleDGEO(fname, dgname):
+    # checks if the file defined in a DGEOFILE kw is a full size 
+    # (old style) image
+    
+    sci_naxis1 = pyfits.getval(fname, 'NAXIS1', ext=1)
+    sci_naxis2 = pyfits.getval(fname, 'NAXIS2', ext=1)
+    dg_naxis1 = pyfits.getval(dgname, 'NAXIS1', ext=1)
+    dg_naxis2 = pyfits.getval(dgname, 'NAXIS2', ext=1)
+    if sci_naxis1 <= dg_naxis1 or sci_naxis2 <= dg_naxis2:
+        print 'Only full size (old style) XY file was found.'
+        print 'DGEO correction will not be applied.\n'
+        return True
+    else:
+        return False
+    
+def getNrefchip(fobj):
+    """
+    This handles the fact that WFPC2 subarray observations
+    may not include chip 3 which is the default reference chip for
+    full observations. Also for subarrays chip 3  may not be the third
+    extension in a MEF file. 
+    """
+    Nrefext = 1
+    instrument = fobj[0].header['INSTRUME']
+    if instrument == 'WFPC2':
+        detectors = [img.header['DETECTOR'] for img in fobj[1:]]
+
+        if 3 not in detectors:
+            Nrefchip=detectors[0]
+            Nrefext = 1
+        else:
+            Nrefchip = 3
+            Nrefext = detectors.index(3) + 1
+    elif instrument == 'ACS':
+        detector = fobj[0].header['DETECTOR']
+        if detector == 'WCS':
+            Nrefchip =2
+        else:
+            Nrefchip = 1
+    elif instrument == 'NICMOS':
+        Nrefchip = fobj[0].header['CAMERA']
+    return Nrefchip, Nrefext
+
+def checkFiles(input):
+    """
+    Purpose
+    =======
+    Checks that input files are in the correct format.
+    Converts geis and waiver fits files to multietension fits.
+    """
+    from pytools.check_files import geis2mef, waiver2mef
+    removed_files = []
+    newfiles = []
+    for file in input:
+        try:
+                imgfits,imgtype = fileutil.isFits(file)
+        except IOError:
+            print "Warning:  File %s could not be found\n" %file
+            print "Removing file %s from input list" %file
+            removed_files.append(file)
+            continue
+        # Check for existence of waiver FITS input, and quit if found.
+        # Or should we print a warning and continue but not use that file
+        if imgfits: 
+            if imgtype == 'waiver':
+                newfilename = waiver2mef(file, convert_dq=True)
+                if newfilename == None:
+                    print "Removing file %s from input list - could not convert waiver to mef" %file
+                    removed_files.append(file)
+                else:
+                    newfiles.append(newfilename)
+            else:
+                newfiles.append(file)
+
+        # If a GEIS image is provided as input, create a new MEF file with 
+        # a name generated using 'buildFITSName()'
+        # Convert the corresponding data quality file if present    
+        if not imgfits:
+            newfilename = geis2mef(file, convert_dq=True)
+            if newfilename == None:
+                print "Removing file %s from input list - could not convert geis to mef" %file
+                removed_files.append(file)
+            else:
+                newfiles.append(newfilename)
+    if removed_files:
+        print 'The following files will be removed from the list of files to be processed :\n'
+        for f in removed_files:
+            print f
+    return newfiles
+
diff --git a/updatewcs/corrections.py b/updatewcs/corrections.py
new file mode 100644
index 0000000..800c1c6
--- /dev/null
+++ b/updatewcs/corrections.py
@@ -0,0 +1,164 @@
+import datetime
+import numpy
+from numpy import linalg
+#from pytools import wcsutil
+
+from makewcs import MakeWCS
+from dgeo import DGEO
+
+class TDDCorr(object):
+    """
+    Purpose
+    =======
+    Apply time dependent distortion correction to SIP coefficients and basic
+    WCS keywords. Atthi stime it is applicable only to ACS/WFC data.
+    
+    Ref: Jay Anderson, ACS ISR 2007-08, Variation of the Distortion 
+    Solution of the WFC 
+    
+    """
+    def __init__(self, owcs=None, refwcs=None):
+        """
+        :Parameters:
+        `owcs`: HSTWCS object
+                An extension HSTWCS object to be modified
+        `refwcs`: HSTWCS object
+                 A reference HSTWCS object
+        """
+        self.wcsobj = owcs
+        if self.wcsobj.DOTDDCorr == 'PERFORM':
+            self.updateWCS()
+            self.wcsobj.DOTDDCorr = 'COMPLETE'
+        else:
+            pass
+
+    def updateWCS(self):
+        """
+        - Calculates alpha and beta for ACS/WFC data.
+        - Writes 2 new kw to the extension header: TDDALPHA and TDDBETA
+        """
+        if self.wcsobj.instrument == 'ACS' and self.wcsobj.detector == 'WFC':
+            newkw = ['TDDALPHA', 'TDDBETA']
+            self.set_alpha_beta()
+            self.updatehdr(newkeywords=newkw)
+            
+        else:
+            pass
+        
+    def set_alpha_beta(self):
+        # Compute the time dependent distortion skew terms
+        # default date of 2004.5 = 2004-7-1
+        
+        datedefault = datetime.datetime(2004,7,1)
+        year,month,day = self.wcsobj.date_obs.split('-')
+        rdate = datetime.datetime(int(year),int(month),int(day))
+        self.TDDALPHA = 0.095+0.090*((rdate-datedefault).days/365.25)/2.5
+        self.TDDBETA = -0.029-0.030*((rdate-datedefault).days/365.25)/2.5
+        self.TDDALPHA = 0.0
+        self.TDDBETA = 0.0
+        
+    def updatehdr(self, newkeywords=None):
+        
+        for kw in newkeywords:
+            self.wcsobj.hdr.update(kw, self.__getattribute__(kw))
+        
+        
+class VACorr(object):
+    """
+    Purpose
+    =======
+    Apply velocity aberation correction to WCS keywords.
+    Modifies the CD matrix and CRVAL1/2
+    """
+    def __init__(self, owcs=None, refwcs=None):
+        self.wcsobj = owcs
+        self.refwcs = refwcs
+        if self.wcsobj.DOVACorr == 'PERFORM':
+            self.updateWCS()
+            self.wcsobj.DOVACorr == 'COMPLETE'
+        else:
+            pass
+        
+    def updateWCS(self):
+        if self.wcsobj.vafactor != 1:
+            self.wcsobj.cd = self.wcsobj.cd * self.wcsobj.vafactor
+            
+            #self.wcsobj.crval[0] = self.wcsobj.ra_targ + self.wcsobj.vafactor*(self.wcsobj.crval[0] - self.wcsobj.ra_targ) 
+            #self.wcsobj.crval[1] = self.wcsobj.dec_targ + self.wcsobj.vafactor*(self.wcsobj.crval[1] - self.wcsobj.dec_targ) 
+            ref_backup = self.refwcs.restore()
+            crval1 = ref_backup['CRVAL1']
+            crval2 = ref_backup['CRVAL2']
+            self.wcsobj.crval[0] = crval1 + self.wcsobj.vafactor*diff_angles(self.wcsobj.crval[0], crval1) 
+            self.wcsobj.crval[1] = crval2 + self.wcsobj.vafactor*diff_angles(self.wcsobj.crval[1], crval2) 
+
+            kw2update={'CD1_1': self.wcsobj.cd[0,0], 'CD1_2':self.wcsobj.cd[0,1], 
+                    'CD2_1':self.wcsobj.cd[1,0], 'CD2_2':self.wcsobj.cd[1,1], 
+                    'CRVAL1':self.wcsobj.crval[0], 'CRVAL2':self.wcsobj.crval[1]}
+            self.updatehdr(newkeywords=kw2update)
+        else:
+            pass
+            
+    def updatehdr(self, newkeywords=None):
+        for kw in newkeywords:
+            self.wcsobj.hdr.update(kw, newkeywords[kw])
+        
+class CompSIP(object):
+    """
+    Purpose
+    =======
+    Compute SIP coefficients from idc table coefficients.
+        
+    """
+    def __init__(self, owcs, refwcs):
+        self.wcsobj = owcs
+        self.updateWCS()
+        self.DOCOMPSIP = 'COMPLETE'
+        
+    def updateWCS(self):
+        kw2update = {}
+        order = self.wcsobj.idcmodel.norder
+        kw2update['A_ORDER'] = order
+        kw2update['B_ORDER'] = order
+        pscale = self.wcsobj.idcmodel.refpix['PSCALE']
+        
+        cx = self.wcsobj.idcmodel.cx
+        cy = self.wcsobj.idcmodel.cy
+        
+        matr = numpy.array([[cx[1,1],cx[1,0]], [cy[1,1],cy[1,0]]], dtype=numpy.float)
+        imatr = linalg.inv(matr)
+        akeys1 = numpy.zeros((order+1,order+1), dtype=numpy.float)
+        bkeys1 = numpy.zeros((order+1,order+1), dtype=numpy.float)
+        for n in range(order+1):
+            for m in range(order+1):
+                if n >= m and n>=2:
+                    idcval = numpy.array([[cx[n,m]],[cy[n,m]]])
+                    sipval = numpy.dot(imatr, idcval)
+                    akeys1[m,n-m] = sipval[0]
+                    bkeys1[m,n-m] = sipval[1]
+                    Akey="A_%d_%d" % (m,n-m)
+                    Bkey="B_%d_%d" % (m,n-m)
+                    kw2update[Akey] = sipval[0,0]
+                    kw2update[Bkey] = sipval[1,0]
+        self.updatehdr(newkw=kw2update)
+                    
+    def updatehdr(self, newkw=None):
+        if not newkw: return
+        for kw in newkw.keys():
+            self.wcsobj.hdr.update(kw, newkw[kw])
+            
+            
+def diff_angles(a,b):
+    """ 
+    Perform angle subtraction a-b taking into account
+    small-angle differences across 360degree line. 
+    """
+    
+    diff = a - b
+
+    if diff > 180.0:
+       diff -= 360.0
+
+    if diff < -180.0:
+       diff += 360.0
+    
+    return diff
diff --git a/updatewcs/dgeo.py b/updatewcs/dgeo.py
new file mode 100644
index 0000000..32e1a38
--- /dev/null
+++ b/updatewcs/dgeo.py
@@ -0,0 +1,206 @@
+import pyfits
+from pytools import fileutil
+from hstwcs.mappings import dgeo_vals
+
+class DGEO(object):
+    """
+    Purpose
+    =======
+    Defines a Lookup table prior distortion correction as per WCS paper IV.
+    It uses a reference file defined by the DGEOFILE keyword in the primary header.
+    
+    Algorithm
+    =========
+    - Using extensions in the reference file create a WCSDVARR extension 
+      and add it to the file object.
+    - Add record-valued keywords which describe the looikup tables to the 
+      science extension header
+    - Add a keyword 'DGEOFILE' to the science extension header, whose
+      value is the reference file used to create the WCSVARR extension
+    
+    """
+    def __init__(self, fobj):
+        """
+        :Parameters:
+        `fobj`: pyfits object
+                Science file, for which a distortion correc tion in a DGEOFILE is available
+                
+        """
+        assert isinstance(fobj, pyfits.NP_pyfits.HDUList)
+        self.fobj = fobj
+        self.applyDgeoCorr()
+        
+        
+    def applyDgeoCorr(self):
+        """
+        For each science extension in a pyfits file object:
+            - create a WCSDVARR extension
+            - update science header
+            - add/update DGEOFILE keyword
+        """
+        
+        dgeover = 0
+        dgfile = fileutil.osfn(self.fobj[0].header['DGEOFILE'])
+        wcsdvarr_ind = self.getwcsindex()
+        for ext in self.fobj:
+            try:
+                extname = ext.header['EXTNAME'].lower()
+            except KeyError:
+                continue
+            if extname == 'sci':
+                extversion = ext.header['EXTVER']
+                for ename in ['DX', 'DY']:
+                    dgeover +=1
+                              
+                    self.addSciExtKw(ext.header, extver=dgeover, ename=ename)
+                    hdu = self.createDgeoHDU(dgeofile=dgfile, dgeover=dgeover,ename=ename, extver=extversion)
+                    if wcsdvarr_ind:
+                        self.fobj[wcsdvarr_ind[dgeover]] = hdu
+                    else:
+                        self.fobj.append(hdu)
+                self.updateDGEOkw(ext.header)
+                
+    def getwcsindex(self):
+        """
+        Returns the index of a WCSDVARR extension in a pyfits file object if it exists.
+        If it exists subsequent updates will overwrite it. If not, it will be
+        added to the file object.
+        """
+        wcsd = {}
+        for e in range(len(self.fobj)):
+            try:
+                ename = self.fobj[e].header['EXTNAME']
+            except KeyError:
+                continue
+            if ename == 'WCSDVARR':
+                wcsd[self.fobj[e].header['EXTVER']] = e
+        return wcsd
+        
+        return self.fobj.index_of(('WCSDVARR', dgeover))
+    
+    def addSciExtKw(self, hdr, extver=None, ename=None):
+        """
+        Adds kw to sci extension to define dgeo correction extensions
+        kw to be added to sci ext:
+        CPERRORj
+        CPDISj
+        DPj.EXTVER
+        DPj.NAXES
+        DPj.AXIS.i (i=DP1.NAXES)
+        """
+        if ename =='DX':
+            j=1
+        else:
+            j=2
+        
+        cperror = 'CPERROR%s' %j
+        cpdis = 'CPDIS%s' %j
+        dpext = 'DP%s.' %j + 'EXTVER'
+        dpnaxes = 'DP%s.' %j +'NAXES'
+        dpaxis1 = 'DP%s.' %j+'AXIS.1'
+        dpaxis2 = 'DP%s.' %j+'AXIS.2'
+        keys = [cperror, cpdis, dpext, dpnaxes, dpaxis1, dpaxis2]
+        values = {cperror: 0.0, cpdis: 'Lookup',  dpext: extver, dpnaxes: 2,
+                dpaxis1: 1, dpaxis2: 2}
+                
+        comments = {cperror: 'Maximum error of dgeo correction for axis %s' % (extver/2), 
+                    cpdis: 'Prior distortion funcion type',  
+                    dpext: 'Version number of WCSDVARR extension containing lookup distortion table', 
+                    dpnaxes: 'Number of independent variables in distortion function',
+                    dpaxis1: 'Axis number of the jth independent variable in a distortion function', 
+                    dpaxis2: 'Axis number of the jth independent variable in a distortion function'
+                    }
+        
+        for key in keys:
+            hdr.update(key=key, value=values[key], comment=comments[key], before='HISTORY')
+        
+        dgfile = self.fobj[0].header['DGEOFILE']
+        hdr.update(key='DGEOFILE', value=dgfile, comment='DGEOFILE used to create WCSDVARR extensions', before='HISTORY')
+        
+        
+    def getDgeoData(self, dgfile=None, ename=None, extver=1):
+        """
+        Given a dgeo file name, creates an array to be used 
+        as a data array in the dgeo extension.
+        """ 
+        return pyfits.getdata(dgfile, ext=(ename,extver))
+            
+    def createDgeoHDU(self, dgeofile=None, dgeover=1, ename=None,extver=1):
+        """
+        Creates an HDU to be added to the file object.
+        """
+        dgeokw = {'naxis1':None, 'naxis2':None, 'extver':dgeover, 'crpix1':None, 
+                'crpix2':None, 'cdelt1':None, 'cdelt2':None, 'crval1':None, 'crval2':None}
+        hdr = self.createDgeoHdr(**dgeokw)
+        data = self.getDgeoData(dgfile=dgeofile, ename=ename, extver=extver)
+        hdu=pyfits.ImageHDU(header=hdr, data=data)
+        return hdu
+    
+    def createDgeoHdr(self, **kw):
+        """
+        Creates a header for the dgeo extension based on dgeo file 
+        and sci extension header.
+        **kw = {'naxis1':None, 'naxis2':None, 'extver':None, 'crpix1':None, 
+                    'crpix2':None, 'cdelt1':None, 'cdelt2':None, 'crval1':None, 'crval2':None}
+        """
+        instr = self.fobj[0].header['INSTRUME']
+        instr_vals = dgeo_vals[instr]
+        naxis1 = kw['naxis1'] or instr_vals['naxis1']
+        naxis2 = kw['naxis2'] or instr_vals['naxis2']
+        extver = kw['extver'] or instr_vals['extver']
+        crpix1 = kw['crpix1'] or instr_vals['crpix1']
+        crpix2 = kw['crpix2'] or instr_vals['crpix2']
+        cdelt1 = kw['cdelt1'] or instr_vals['cdelt1']
+        cdelt2 = kw['cdelt2'] or instr_vals['cdelt2']
+        crval1 = kw['crval1'] or instr_vals['crval1']
+        crval2 = kw['crval2'] or instr_vals['crval2']
+            
+        keys = ['XTENSION','BITPIX','NAXIS','NAXIS1','NAXIS2',
+              'EXTNAME','EXTVER','PCOUNT','GCOUNT','CRPIX1',
+                        'CDELT1','CRVAL1','CRPIX2','CDELT2','CRVAL2']
+                        
+        comments = {'XTENSION': 'Image extension',
+                    'BITPIX': 'IEEE floating point',
+                    'NAXIS': 'Number of image columns',
+                    'NAXIS1': 'Number of image columns',
+                    'NAXIS2': 'Number of image rows',
+                    'EXTNAME': 'WCS distortion array',
+                    'EXTVER': 'Distortion array version number',
+                    'PCOUNT': 'Special data area of size 0',
+                    'GCOUNT': 'One data group',
+                    'CRPIX1': 'Distortion array reference pixel',
+                    'CDELT1': 'Grid step size in first coordinate',
+                    'CRVAL1': 'Image array pixel coordinate',
+                    'CRPIX2': 'Distortion array reference pixel',
+                    'CDELT2': 'Grid step size in second coordinate',
+                    'CRVAL2': 'Image array pixel coordinate'}
+        
+        values = {'XTENSION': 'IMAGE',
+                'BITPIX': -32,
+                'NAXIS': 2,
+                'NAXIS1': naxis1,
+                'NAXIS2': naxis2,
+                'EXTNAME': 'WCSDVARR',
+                'EXTVER':  extver,
+                'PCOUNT': 0,
+                'GCOUNT': 1,
+                'CRPIX1': crpix1,
+                'CDELT1': cdelt1,
+                'CRVAL1': crval1,
+                'CRPIX2': crpix1,
+                'CDELT2': cdelt2,
+                'CRVAL2': crval2
+                }
+                    
+        
+        cdl = pyfits.CardList()
+        for c in keys:
+            cdl.append(pyfits.Card(key=c, value=values[c], comment=comments[c]))
+
+        hdr = pyfits.Header(cards=cdl)
+        return hdr
+        
+    def updateDGEOkw(self, hdr):
+        dgfile = self.fobj[0].header['DGEOFILE']
+        hdr.update(key='DGEOFILE', value=dgfile, comment='DGEOFILE used to create WCSDVARR extensions', before='HISTORY')
+        
\ No newline at end of file
diff --git a/updatewcs/makewcs.py b/updatewcs/makewcs.py
new file mode 100644
index 0000000..451da13
--- /dev/null
+++ b/updatewcs/makewcs.py
@@ -0,0 +1,260 @@
+#from .. mappings import DEGTORAD, RADTODEG
+from hstwcs.mappings import DEGTORAD, RADTODEG
+import numpy
+from numpy import math 
+from math import sin, sqrt, pow, cos, asin, atan2,pi
+
+class MakeWCS(object):
+    """
+    Purpose
+    =======
+    Recompute WCS keywords based on PA_V3 and distortion model.
+    
+    Algorithm
+    =========
+    - update reference chip wcs
+    - update extension wcs
+    - update extension header
+    
+    """
+    
+    def __init__(self, owcs, refwcs):
+        self.wcsobj = owcs
+        self.refwcs = refwcs
+        self.updateWCS()
+        self.DOMakeWCS = 'COMPLETE'
+        
+    def updateWCS(self):    
+        """
+        recomputes the basic WCS kw 
+        """
+        backup_wcs = self.wcsobj.restore()
+        ltvoff, offshift = self.getOffsets(backup_wcs)
+            
+        self.uprefwcs(ltvoff, offshift)
+        self.upextwcs(ltvoff, offshift)
+        self.updatehdr()
+        
+    def upextwcs(self, loff, offsh):
+        """
+        updates an extension wcs
+        """  
+        ltvoffx, ltvoffy = loff[0], loff[1]
+        offshiftx, offshifty = offsh[0], offsh[1]
+        backup_wcs = self.wcsobj.restore()
+        ltv1 = self.wcsobj.ltv1
+        ltv2 = self.wcsobj.ltv2  
+        
+        
+        
+        if ltv1 != 0. or ltv2 != 0.:
+            offsetx = backup_wcs['CRPIX1'] - ltv1 - self.wcsobj.idcmodel.refpix['XREF']
+            offsety = backup_wcs['CRPIX2'] - ltv2 - self.wcsobj.idcmodel.refpix['YREF']
+            fx,fy = self.wcsobj.idcmodel.shift(self.wcsobj.idcmodel.cx,self.wcsobj.idcmodel.cy,offsetx,offsety)
+        else:
+            fx, fy = self.wcsobj.idcmodel.cx, self.wcsobj.idcmodel.cy
+        
+        ridcmodel = self.refwcs.idcmodel
+        v2 = self.wcsobj.idcmodel.refpix['V2REF']
+        v3 = self.wcsobj.idcmodel.refpix['V3REF']
+        v2ref = self.refwcs.idcmodel.refpix['V2REF']
+
+        v3ref = self.refwcs.idcmodel.refpix['V3REF']
+        R_scale = self.refwcs.idcmodel.refpix['PSCALE']/3600.0
+        off = sqrt((v2-v2ref)**2 + (v3-v3ref)**2)/(R_scale*3600.0)
+
+        # Here we must include the PARITY
+        if v3 == v3ref:
+           theta=0.0
+        else:
+           theta = atan2(self.wcsobj.parity[0][0]*(v2-v2ref),self.wcsobj.parity[1][1]*(v3-v3ref))
+        
+        if self.refwcs.idcmodel.refpix['THETA']: theta += self.refwcs.idcmodel.refpix['THETA']*pi/180.0
+
+        dX=(off*sin(theta)) + offshiftx
+        dY=(off*cos(theta)) + offshifty
+        px = numpy.array([[dX,dY]])
+        newcrval = self.refwcs.pixel2world_fits(px)[0]
+        newcrpix = numpy.array([self.wcsobj.idcmodel.refpix['XREF'] + ltvoffx, 
+                                self.wcsobj.idcmodel.refpix['YREF'] + ltvoffy])
+        
+        self.wcsobj.crval = newcrval
+        self.wcsobj.crpix = newcrpix
+        
+        # Account for subarray offset
+        # Angle of chip relative to chip
+        if self.wcsobj.idcmodel.refpix['THETA']:
+           dtheta = self.wcsobj.idcmodel.refpix['THETA'] - self.refwcs.idcmodel.refpix['THETA']
+        else:
+           dtheta = 0.0
+        
+        
+        # Create a small vector, in reference image pixel scale
+        # There is no parity effect here ???
+        
+        delXX=fx[1,1]/R_scale/3600.
+        delYX=fy[1,1]/R_scale/3600.
+        delXY=fx[1,0]/R_scale/3600.
+        delYY=fy[1,0]/R_scale/3600.
+        
+        
+        # Convert to radians
+        rr=dtheta*pi/180.0
+        
+        # Rotate the vectors
+        dXX = cos(rr)*delXX - sin(rr)*delYX
+        dYX = sin(rr)*delXX + cos(rr)*delYX
+
+        dXY = cos(rr)*delXY - sin(rr)*delYY
+        dYY = sin(rr)*delXY + cos(rr)*delYY
+        px = numpy.array([[dX+dXX,dY+dYX]])
+        
+        # Transform to sky coordinates
+        
+        wc = self.refwcs.pixel2world_fits(px)
+        
+        a = wc[0,0]
+        b = wc[0,1]
+        px = numpy.array([[dX+dXY,dY+dYY]])
+        
+        wc = self.refwcs.pixel2world_fits(px)
+        c = wc[0,0]
+        d = wc[0,1]
+        
+        # Calculate the new CDs and convert to degrees
+        cd11 = diff_angles(a,newcrval[0])*cos(newcrval[1]*pi/180.0)
+        cd12 = diff_angles(c,newcrval[0])*cos(newcrval[1]*pi/180.0)
+        cd21 = diff_angles(b,newcrval[1])
+        cd22 = diff_angles(d,newcrval[1])
+        
+        cd = numpy.array([[cd11, cd12], [cd21, cd22]])
+        self.wcsobj.cd = cd
+        
+    def uprefwcs(self, loff, offsh):
+        """
+        Updates the reference chip
+        """
+        ltvoffx, ltvoffy = loff[0], loff[1]
+        offshift = offsh
+        backup_refwcs = self.refwcs.restore()
+        dec = backup_refwcs['CRVAL2']
+    
+        # Get an approximate reference position on the sky
+        rref = numpy.array([[self.refwcs.idcmodel.refpix['XREF']+ltvoffx, 
+                            self.refwcs.idcmodel.refpix['YREF']+ltvoffy]])
+        
+        
+        crval = self.refwcs.pixel2world_fits(rref)
+        
+        # Convert the PA_V3 orientation to the orientation at the aperture
+        # This is for the reference chip only - we use this for the
+        # reference tangent plane definition
+        # It has the same orientation as the reference chip
+        pv = troll(self.wcsobj.pav3,dec,self.refwcs.idcmodel.refpix['V2REF'],
+                    self.refwcs.idcmodel.refpix['V3REF'])
+        # Add the chip rotation angle
+        if self.refwcs.idcmodel.refpix['THETA']:
+            pv += self.refwcs.idcmodel.refpix['THETA']
+            
+        
+        # Set values for the rest of the reference WCS
+        self.refwcs.crval = crval[0]
+        self.refwcs.crpix = numpy.array([0.0,0.0])+offsh
+        parity = self.refwcs.parity
+        R_scale = self.refwcs.idcmodel.refpix['PSCALE']/3600.0
+        cd11 = parity[0][0] *  cos(pv*pi/180.0)*R_scale
+        cd12 = parity[0][0] * -sin(pv*pi/180.0)*R_scale
+        cd21 = parity[1][1] *  sin(pv*pi/180.0)*R_scale
+        cd22 = parity[1][1] *  cos(pv*pi/180.0)*R_scale
+        rcd = numpy.array([[cd11, cd12], [cd21, cd22]])
+        self.refwcs.cd = rcd
+        self.refwcs.set()
+    
+    
+        
+    def updatehdr(self):
+        """
+        Keywords to update:
+        CD1_1, CD1_2, CD2_1, CD2_2, CRPIX1, CRPIX2, CRVAL1, CRVAL2
+        """
+        
+        self.wcsobj.hdr.update('CD1_1', self.wcsobj.cd[0,0])
+        self.wcsobj.hdr.update('CD1_2', self.wcsobj.cd[0,1])
+        self.wcsobj.hdr.update('CD2_1', self.wcsobj.cd[1,0])
+        self.wcsobj.hdr.update('CD2_2', self.wcsobj.cd[1,1])
+        self.wcsobj.hdr.update('CRVAL1', self.wcsobj.crval[0])
+        self.wcsobj.hdr.update('CRVAL2', self.wcsobj.crval[1])
+        self.wcsobj.hdr.update('CRPIX1', self.wcsobj.crpix[0])
+        self.wcsobj.hdr.update('CRPIX2', self.wcsobj.crpix[1])
+        self.wcsobj.hdr.update('ORIENTAT', self.wcsobj.orientat)
+        
+    def getOffsets(self, backup_wcs):
+        ltv1 = self.wcsobj.ltv1
+        ltv2 = self.wcsobj.ltv2
+        
+        offsetx = backup_wcs['CRPIX1'] - ltv1 - self.wcsobj.idcmodel.refpix['XREF']
+        offsety = backup_wcs['CRPIX2'] - ltv2 - self.wcsobj.idcmodel.refpix['YREF']
+        
+        shiftx = self.wcsobj.idcmodel.refpix['XREF'] + ltv1
+        shifty = self.wcsobj.idcmodel.refpix['YREF'] + ltv2
+        if ltv1 != 0. or ltv2 != 0.:
+            ltvoffx = ltv1 + offsetx
+            ltvoffy = ltv2 + offsety
+            offshiftx = offsetx + shiftx
+            offshifty = offsety + shifty
+        else:
+            ltvoffx = 0.
+            ltvoffy = 0.
+            offshiftx = 0.
+            offshifty = 0.
+            
+        ltvoff = numpy.array([ltvoffx, ltvoffy])
+        offshift = numpy.array([offshiftx, offshifty])
+        return ltvoff, offshift
+    
+def diff_angles(a,b):
+    """ 
+    Perform angle subtraction a-b taking into account
+    small-angle differences across 360degree line. 
+    """
+    
+    diff = a - b
+
+    if diff > 180.0:
+       diff -= 360.0
+
+    if diff < -180.0:
+       diff += 360.0
+    
+    return diff
+
+def troll(roll, dec, v2, v3):
+    """ Computes the roll angle at the target position based on:
+            the roll angle at the V1 axis(roll),
+            the dec of the target(dec), and
+            the V2/V3 position of the aperture (v2,v3) in arcseconds.
+
+        Based on the algorithm provided by Colin Cox that is used in
+        Generic Conversion at STScI.
+    """
+    # Convert all angles to radians
+    _roll = DEGTORAD(roll)
+    _dec = DEGTORAD(dec)
+    _v2 = DEGTORAD(v2 / 3600.)
+    _v3 = DEGTORAD(v3 / 3600.)
+
+    # compute components
+    sin_rho = sqrt((pow(sin(_v2),2)+pow(sin(_v3),2)) - (pow(sin(_v2),2)*pow(sin(_v3),2)))
+    rho = asin(sin_rho)
+    beta = asin(sin(_v3)/sin_rho)
+    if _v2 < 0: beta = pi - beta
+    gamma = asin(sin(_v2)/sin_rho)
+    if _v3 < 0: gamma = pi - gamma
+    A = pi/2. + _roll - beta
+    B = atan2( sin(A)*cos(_dec), (sin(_dec)*sin_rho - cos(_dec)*cos(rho)*cos(A)))
+
+    # compute final value
+    troll = RADTODEG(pi - (gamma+B))
+
+    return troll
+
diff --git a/wcsutil/__init__.py b/wcsutil/__init__.py
new file mode 100644
index 0000000..46ad42e
--- /dev/null
+++ b/wcsutil/__init__.py
@@ -0,0 +1,166 @@
+#from .. pywcs.sip import SIP
+from pywcs.sip import SIP
+import pyfits
+import instruments
+#from .. distortion import models
+from hstwcs.distortion import models
+import numpy as N
+
+#from .. mappings import inst_mappings, ins_spec_kw, DEGTORAD, RADTODEG, basic_wcs
+from hstwcs.mappings import inst_mappings, ins_spec_kw, DEGTORAD, RADTODEG
+from hstwcs.mappings import basic_wcs, prim_hdr_kw
+
+__docformat__ = 'restructuredtext'
+
+class HSTWCS(SIP):
+    """
+    Purpose
+    =======
+    Create a wcs object based on the instrument.
+    It has all basic WCS kw as attribbutes (set by pywcs).
+    It also uses the primary and extension header to define 
+    instrument specific attributes needed by the correction classes.
+    """
+    
+    def __init__(self, hdr0, hdr):
+        """
+        :Parameters:
+        `hdr0`: Pyfits Header
+                primary header
+        `hdr`: Pyfits Header
+               extension header
+        """
+        self.hdr = hdr
+        self.hdr0 = hdr0
+        SIP.__init__(self, self.hdr)
+        self.setHDR0kw(hdr0)
+        self.detector = self.setDetector()
+        self.inst_kw = ins_spec_kw
+        self.setInstrSpecKw()
+        self.pscale = self.setPscale()
+        self.orientat = self.setOrient()
+        
+        
+       
+    def setHDR0kw(self, primhdr):
+        #sets attributes from kw defined in the primary header
+        self.instrument = primhdr.get('INSTRUME', None)
+        self.offtab = primhdr.get('OFFTAB', None) 
+        self.idctab = primhdr.get('IDCTAB', None)
+        self.date_obs = primhdr.get('DATE-OBS', None)
+        self.pav3 = primhdr.get('PA_V3', None)
+        self.ra_targ = primhdr.get('RA_TARG', None)
+        self.dec_targ = primhdr.get('DEC_TARG', None)
+    
+    
+    def setDetector(self):
+        #sets detector attribute for instuments which have more than one detector
+        if self.instrument in ['ACS', 'WFC3']:
+            return self.hdr0.get('DETECTOR', None)
+    
+    
+    def setInstrSpecKw(self):
+        #Based on the instrument kw creates an instalnce of an instrument WCS class
+        #and sets attributes from instrument specific kw 
+        if self.instrument in inst_mappings.keys():
+            inst_kl = inst_mappings[self.instrument]
+            inst_kl = instruments.__dict__[inst_kl]
+            insobj = inst_kl(self.hdr0, self.hdr)
+            for key in self.inst_kw:
+                self.__setattr__(key, insobj.__getattribute__(key))
+        else:
+            raise KeyError, "Unsupported instrument - %s" %self.instrument
+    
+    def setPscale(self):
+        #Calculates the plate scale from the cd matrix
+        #this may not be needed now and shoufl probably be done after all 
+        # corrections
+        cd11 = self.cd[0][0]
+        cd21 = self.cd[1][0]
+        return N.sqrt(N.power(cd11,2)+N.power(cd21,2)) * 3600.
+    
+    def setOrient(self):
+        #same as setPscale
+        cd12 = self.cd[0][1]
+        cd22 = self.cd[1][1]
+        return RADTODEG(N.arctan2(cd12,cd22))
+    
+    def verifyPa_V3(self):
+        """make sure PA_V3 is populated"""    
+        
+    def verifyKw(self):
+        """verify that all required kw have meaningful values"""
+            
+    def readModel(self):
+        """
+        Purpose
+        =======
+        Read distortion model from idc table.
+        Save some of the information as kw needed later by pydrizzle
+        
+        list of kw - to be revised later
+        """
+        self.idcmodel = models.IDCModel(self.idctab,
+                    chip=self.chip, direction='forward', date=self.date_obs,
+                    filter1=self.filter1, filter2=self.filter2, 
+                    offtab=self.offtab, binned=self.binned)
+        
+        
+        self.updatehdr()
+                   
+    
+    def updatehdr(self, newkeywords=None):
+        #kw2add : OCX10, OCX11, OCY10, OCY11 
+        # record the model in the header for use by pydrizzle
+        self.hdr.update('OCX10', self.idcmodel.cx[1,0])
+        self.hdr.update('OCX11', self.idcmodel.cx[1,1])
+        self.hdr.update('OCY10', self.idcmodel.cy[1,0])
+        self.hdr.update('OCY11', self.idcmodel.cy[1,1])
+        self.hdr.update('IDCSCALE', self.idcmodel.refpix['PSCALE'])
+        self.hdr.update('IDCTHETA', self.idcmodel.refpix['THETA'])
+        self.hdr.update('IDCXREF', self.idcmodel.refpix['XREF'])
+        self.hdr.update('IDCYREF', self.idcmodel.refpix['YREF'])
+        self.hdr.update('IDCV2REF', self.idcmodel.refpix['V2REF'])
+        self.hdr.update('IDCV3REF', self.idcmodel.refpix['V3REF'])
+        self.hdr.update('IDCXSIZE', self.idcmodel.refpix['XSIZE'])
+        self.hdr.update('IDCYSIZE', self.idcmodel.refpix['YSIZE'])
+        self.hdr.update('IDCXDELT', self.idcmodel.refpix['XDELTA'])
+        self.hdr.update('IDCYDELT', self.idcmodel.refpix['YDELTA'])
+        self.hdr.update('CENTERED', self.idcmodel.refpix['centered'])
+        
+    
+    def restore(self):
+        """
+        restores basic wcs keywords from archive
+        this should be modified to always return a populated 
+        dictionary, although the values may be None.
+        """    
+
+        backup = {}
+        for k in basic_wcs:
+            try:
+                nkw = ('O'+k)[:7]
+                #backup[k] = self.hdr.__getitem__('O'+k)
+                backup[k] = self.hdr[nkw]
+                #self.__setattr__(k, self.hdr.__getitem__('O'+k))
+                self.__setattr__(k, self.hdr[nkw])
+            except KeyError:
+                print 'Keyword %s not found' % k
+                
+        return backup
+    
+    def archive_kw(self):
+        """
+        archive original WCS kw before recalculating them.
+        """
+        backup_kw = self.restore()
+        if backup_kw != {}:
+            return
+        else:
+            # keep the archived kw 8 characters long
+            cmt = 'archived value'
+            for kw in basic_wcs:
+                nkw = 'O'+kw
+                self.hdr.update(nkw[:7], self.hdr[kw], comment=cmt)
+        
+
diff --git a/wcsutil/instruments.py b/wcsutil/instruments.py
new file mode 100644
index 0000000..2790d8a
--- /dev/null
+++ b/wcsutil/instruments.py
@@ -0,0 +1,117 @@
+import pyfits
+import numpy as N
+#from .. mappings import ins_spec_kw
+from hstwcs.mappings import ins_spec_kw, prim_hdr_kw
+    
+class InstrWCS(object):
+    """
+    A base class for instrument specific keyword definition.
+    It prvides a default implementation (modeled by ACS) for
+    all set_kw methods.
+    """
+    def __init__(self, hdr0, hdr):
+        self.exthdr = hdr   
+        self.primhdr = hdr0
+        
+    def set_ins_spec_kw(self):
+        """
+        This method MUST call all set_kw methods.
+        There should be a set_kw method for all kw listed in 
+        mappings.ins_spec_kw
+        """
+        self.set_filter1()
+        self.set_filter2()
+        self.set_vafactor()
+        self.set_naxis1()
+        self.set_naxis2()
+        self.set_ltv1()
+        self.set_ltv2()
+        self.set_binned()
+        self.set_chip()
+        self.set_parity()
+        
+    def set_filter1(self):
+        self.filter1 = self.primhdr.get('FILTER1', None)
+        
+    def set_filter2(self):
+        self.filter2 = self.primhdr.get('FILTER2', None) 
+               
+    def set_vafactor(self):
+        self.vafactor = self.exthdr.get('vafactor', 1) 
+    
+    def set_naxis1(self):
+        self.naxis1 = self.exthdr.get('naxis1', None)
+        
+    def set_naxis2(self):
+        self.naxis2 = self.exthdr.get('naxis2', None)
+        
+    def set_ltv1(self):
+        self.ltv1 = self.exthdr.get('ltv1', 0.0)
+        
+    def set_ltv2(self):
+        self.ltv2 = self.exthdr.get('ltv2', 0.0)
+        
+    def set_binned(self):
+        self.binned = self.exthdr.get('BINAXIS1', 1)
+
+    def set_chip(self):
+        self.chip = self.exthdr.get('CCDCHIP', 1)
+    
+    def set_parity(self):
+        self.parity = [[1.0,0.0],[0.0,-1.0]]
+        
+class ACSWCS(InstrWCS):
+    """
+    get instrument specific kw    
+    """
+    
+    def __init__(self, hdr0, hdr):
+        self.primhdr = hdr0
+        self.exthdr = hdr
+        InstrWCS.__init__(self,hdr0, hdr)
+        self.set_ins_spec_kw()
+    
+    
+    def set_parity(self):
+        parity = {'WFC':[[1.0,0.0],[0.0,-1.0]],
+                'HRC':[[-1.0,0.0],[0.0,1.0]],
+                'SBC':[[-1.0,0.0],[0.0,1.0]]}
+        detector = self.primhdr.get('detector', None)  
+        print 'detector', detector   
+        self.parity = parity[detector]
+    
+        
+class WFPC2WCS(InstrWCS):   
+
+
+    def __init__(self, hdr0, hdr):
+        self.primhdr = hdr0
+        self.exthdr = hdr
+        InstrWCS.__init__(self,hdr0, hdr)
+        self.set_ins_spec_kw()
+    
+    def set_filter1(self):
+        self.filter1 = self.primhdr.get('FILTNAM1', None)
+        if self.filter1 == " " or self.filter1 == None:
+            self.filter1 = 'CLEAR1'
+
+    def set_filter2(self):
+        self.filter2 = self.primhdr.get('FILTNAM2', None)
+        if self.filter2 == " " or self.filter2 == None:
+            self.filter2 = 'CLEAR1'
+            
+    
+    def set_binned(self):
+        mode = self.primhdr.get('MODE', 1)
+        if mode == 'FULL':
+            self.binned = 1
+        elif mode == 'AREA':
+            self.binned = 2
+
+    def set_chip(self):
+        self.chip = self.exthdr.get('DETECTOR', 1)
+    
+    def set_parity(self):
+        self.parity = [[-1.0,0.],[0.,1.0]]
+        
+