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

4 files changed, 892 insertions, 0 deletions

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
+