first round of pep8 changes

13 files changed, 474 insertions, 497 deletions

diff --git a/stwcs/distortion/coeff_converter.py b/stwcs/distortion/coeff_converter.py
index 415b512..99cf74d 100644
--- a/stwcs/distortion/coeff_converter.py
+++ b/stwcs/distortion/coeff_converter.py
@@ -1,9 +1,10 @@
-from __future__ import division, print_function # confidence high
+from __future__ import (absolute_import, unicode_literals, division, print_function)
 
 import numpy as np
 from astropy.io import fits
 from astropy import wcs as pywcs
 
+
 def sip2idc(wcs):
     """
     Converts SIP style coefficients to IDCTAB coefficients.
@@ -49,28 +50,28 @@ def sip2idc(wcs):
         print('Input to sip2idc must be a PyFITS header or a wcsutil.HSTWCS object\n')
         return
 
-
     if None in [ocx10, ocx11, ocy10, ocy11]:
         print('First order IDC coefficients not found, exiting ...\n')
         return None, None
     idc_coeff = np.array([[ocx11, ocx10], [ocy11, ocy10]])
-    cx = np.zeros((order+1,order+1), dtype=np.double)
-    cy = np.zeros((order+1,order+1), dtype=np.double)
-    for n in range(order+1):
-        for m in range(order+1):
-            if n >= m and n>=2:
-                sipval = np.array([[sipa[m,n-m]],[sipb[m,n-m]]])
+    cx = np.zeros((order + 1, order + 1), dtype=np.double)
+    cy = np.zeros((order + 1, order + 1), dtype=np.double)
+    for n in range(order + 1):
+        for m in range(order + 1):
+            if n >= m and n >= 2:
+                sipval = np.array([[sipa[m, n - m]], [sipb[m, n - m]]])
                 idcval = np.dot(idc_coeff, sipval)
-                cx[n,m] = idcval[0]
-                cy[n,m] = idcval[1]
+                cx[n, m] = idcval[0]
+                cy[n, m] = idcval[1]
 
-    cx[1,0] = ocx10
-    cx[1,1] = ocx11
-    cy[1,0] = ocy10
-    cy[1,1] = ocy11
+    cx[1, 0] = ocx10
+    cx[1, 1] = ocx11
+    cy[1, 0] = ocy10
+    cy[1, 1] = ocy11
 
     return cx, cy
 
+
 def _read_sip_kw(header):
     """
     Reads SIP header keywords and returns an array of coefficients.
@@ -84,15 +85,15 @@ def _read_sip_kw(header):
                 "keyword for SIP distortion")
 
         m = int(header["A_ORDER"])
-        a = np.zeros((m+1, m+1), np.double)
-        for i in range(m+1):
-            for j in range(m-i+1):
+        a = np.zeros((m + 1, m + 1), np.double)
+        for i in range(m + 1):
+            for j in range(m - i + 1):
                 a[i, j] = header.get("A_%d_%d" % (i, j), 0.0)
 
         m = int(header["B_ORDER"])
-        b = np.zeros((m+1, m+1), np.double)
-        for i in range(m+1):
-            for j in range(m-i+1):
+        b = np.zeros((m + 1, m + 1), np.double)
+        for i in range(m + 1):
+            for j in range(m - i + 1):
                 b[i, j] = header.get("B_%d_%d" % (i, j), 0.0)
     elif "B_ORDER" in header:
         raise ValueError(
@@ -102,7 +103,7 @@ def _read_sip_kw(header):
         a = None
         b = None
 
-    return a , b
+    return a, b
 
 
 """
diff --git a/stwcs/distortion/models.py b/stwcs/distortion/models.py
index 231a9f1..136c6cf 100644
--- a/stwcs/distortion/models.py
+++ b/stwcs/distortion/models.py
@@ -1,14 +1,9 @@
-from __future__ import absolute_import, division, print_function # confidence high
-
+from __future__ import absolute_import, division, print_function
 
 import numpy as np
-
-# Import PyDrizzle utility modules
 from . import mutil
 from .mutil import combin
 
-yes = True
-no = False
 
 #################
 #
@@ -72,8 +67,8 @@ class GeometryModel:
         to the reference position of the chip.
         """
 
-        _cxs = np.zeros(shape=self.cx.shape,dtype=self.cx.dtype)
-        _cys = np.zeros(shape=self.cy.shape,dtype=self.cy.dtype)
+        _cxs = np.zeros(shape=self.cx.shape, dtype=self.cx.dtype)
+        _cys = np.zeros(shape=self.cy.shape, dtype=self.cy.dtype)
         _k = self.norder + 1
         # loop over each input coefficient
         for m in range(_k):
@@ -83,15 +78,17 @@ class GeometryModel:
                     _ilist = list(range(m, _k))
                     # sum from m to k
                     for i in _ilist:
-                        _jlist = list(range(n, i - (m-n)+1))
+                        _jlist = list(range(n, i - (m - n) + 1))
                         # sum from n to i-(m-n)
                         for j in _jlist:
-                            _cxs[m,n] += self.cx[i,j]*combin(j,n)*combin((i-j),(m-n))*pow(xs,(j-n))*pow(ys,((i-j)-(m-n)))
-                            _cys[m,n] += self.cy[i,j]*combin(j,n)*combin((i-j),(m-n))*pow(xs,(j-n))*pow(ys,((i-j)-(m-n)))
+                            _cxs[m, n] += self.cx[i, j] * combin(j, n) * combin((i - j), (m - n)) * \
+                                pow(xs, (j - n)) * pow(ys, ((i - j) - (m - n)))
+                            _cys[m, n] += self.cy[i, j] * combin(j, n) * combin((i - j), (m - n)) * \
+                                pow(xs, (j - n)) * pow(ys, ((i - j) - (m - n)))
         self.cx = _cxs.copy()
         self.cy = _cys.copy()
 
-    def convert(self, tmpname, xref=None,yref=None,delta=yes):
+    def convert(self, tmpname, xref=None, yref=None, delta=True):
         """
          Open up an ASCII file, output coefficients in drizzle
           format after converting them as necessary.
@@ -104,89 +101,62 @@ class GeometryModel:
             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']
+        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['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
+        self.x0 = x0
+        self.y0 = y0
 
         # 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('# 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:
+        elif self.norder == 4:
             lines.append('quartic\n')
-        elif self.norder==5:
+        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])
+        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])
+        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])
+        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])
+        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])
+        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])
+        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])
+        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])
+        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 = open(tmpname, 'w')
         output.writelines(lines)
         output.close()
 
-
-    def apply(self, pixpos,scale=1.0,order=None):
+    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.
@@ -195,13 +165,13 @@ class GeometryModel:
 
         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 self.cx is None:
+            return pixpos[:, 0], pixpos[:, 1]
 
         if order is None:
             order = self.norder
@@ -209,29 +179,29 @@ class GeometryModel:
         # Apply in the same way that 'drizzle' would...
         _cx = self.cx / (self.pscale * scale)
         _cy = self.cy / (self.pscale * scale)
-        _convert = no
+        _convert = False
         _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.
+        _cx[0, 0] = 0.
+        _cy[0, 0] = 0.
 
         if isinstance(_p, list) or isinstance(_p, tuple):
-            _p = np.array(_p,dtype=np.float64)
-            _convert = yes
+            _p = np.array(_p, dtype=np.float64)
+            _convert = True
 
-        dxy = _p - (self.refpix['XREF'],self.refpix['YREF'])
+        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]
+        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:
@@ -242,11 +212,11 @@ class GeometryModel:
                 xc = xc[0]
                 yc = yc[0]
 
-        return xc,yc
+        return xc, yc
 
-    def setPScaleCoeffs(self,pscale):
-        self.cx[1,1] = pscale
-        self.cy[1,0] = pscale
+    def setPScaleCoeffs(self, pscale):
+        self.cx[1, 1] = pscale
+        self.cy[1, 0] = pscale
 
         self.refpix['PSCALE'] = pscale
         self.pscale = pscale
@@ -259,7 +229,7 @@ class IDCModel(GeometryModel):
     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):
+                 filter1='CLEAR1', filter2='CLEAR2', offtab=None, binned=1):
         GeometryModel.__init__(self)
         #
         # Norder must be derived from the coeffs file itself,
@@ -269,9 +239,9 @@ class IDCModel(GeometryModel):
         # 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)
+        self.cx, self.cy, self.refpix, self.norder = mutil.readIDCtab(idcfile,
+                                                                      chip=chip, direction=direction, filter1=filter1, filter2=filter2,
+                                                                      date=date, offtab=offtab)
 
         if 'empty_model' in self.refpix and self.refpix['empty_model']:
             pass
@@ -285,61 +255,54 @@ class IDCModel(GeometryModel):
             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):
+    def __init__(self, header, rootname):
         GeometryModel.__init__(self)
 
-
         if 'rootname' in header:
             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()
+        # 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.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):
+    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)
+        self.cx, self.cy, self.refpix, self.norder = mutil.readCubicTable(idcfile)
 
-        # scale is the ratio wcs.pscale/model.pscale. 
+        # 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:
+
+        if scale is not 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):
     """
@@ -349,14 +312,8 @@ class TraugerModel(GeometryModel):
     """
     NORDER = 3
 
-    def __init__(self, idcfile,lam):
+    def __init__(self, idcfile, lam):
         GeometryModel.__init__(self)
         self.name = idcfile
-        self.cx,self.cy,self.refpix,self.norder = mutil.readTraugerTable(idcfile,lam)
+        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/stwcs/distortion/mutil.py b/stwcs/distortion/mutil.py
index ed6a1ea..0e193c8 100644
--- a/stwcs/distortion/mutil.py
+++ b/stwcs/distortion/mutil.py
@@ -1,25 +1,22 @@
-from __future__ import division, print_function # confidence high
+from __future__ import division, print_function
 
 from stsci.tools import fileutil
 import numpy as np
 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.
+#    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):
+
+def readIDCtab(tabname, chip=1, date=None, direction='forward',
+               filter1=None, filter2=None, offtab=None):
 
     """
         Read IDCTAB, and optional OFFTAB if sepcified, and generate
@@ -30,17 +27,17 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
 
     """
 
- # Return a default geometry model if no IDCTAB filename
+    # Return a default geometry model if no IDCTAB filename
     # is given.  This model will not distort the data in any way.
-    if tabname == None:
+    if tabname is 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 or filter1.strip() == '':
+    if filter1 is None or filter1.find('CLEAR') == 0 or filter1.strip() == '':
         filter1 = 'CLEAR'
-    if filter2 == None or filter2.find('CLEAR') == 0 or filter2.strip() == '':
+    if filter2 is None or filter2.find('CLEAR') == 0 or filter2.strip() == '':
         filter2 = 'CLEAR'
 
     # Insure that tabname is full filename with fully expanded
@@ -51,7 +48,7 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
     try:
         ftab = fileutil.openImage(tabname)
     except:
-        err_str =  "------------------------------------------------------------------------ \n"
+        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"
@@ -63,7 +60,7 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
         err_str += "------------------------------------------------------------------------ \n"
         raise IOError(err_str)
 
-    #First thing we need, is to read in the coefficients from the IDC
+    # First thing we need, is to read in the coefficients from the IDC
     # table and populate the Fx and Fy matrices.
 
     if 'DETECTOR' in ftab['PRIMARY'].header:
@@ -96,10 +93,10 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
     else:
         order = norder
 
-    fx = np.zeros(shape=(order+1,order+1),dtype=np.float64)
-    fy = np.zeros(shape=(order+1,order+1),dtype=np.float64)
+    fx = np.zeros(shape=(order + 1, order + 1), dtype=np.float64)
+    fy = np.zeros(shape=(order + 1, order + 1), dtype=np.float64)
 
-    #Determine row from which to get the coefficients.
+    # 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
@@ -111,14 +108,14 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
 
         try:
             # Match FILTER combo to appropriate row,
-            #if there is a filter column in the IDCTAB...
+            # 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]
+                if filt1.find('CLEAR') > -1: filt1 = filt1[: 5]
 
                 filt2 = ftab[1].data.field('FILTER2')[i]
-                if filt2.find('CLEAR') > -1: filt2 = filt2[:5]
+                if filt2.find('CLEAR') > -1: filt2 = filt2[: 5]
             else:
                 if 'OPT_ELEM' in colnames:
                     filt1 = ftab[1].data.field('OPT_ELEM')
@@ -170,16 +167,17 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
         joinstr = ''
     else:
         f2str = filter2.strip()
-    filtstr = (joinstr.join([f1str,f2str])).strip()
+    filtstr = (joinstr.join([f1str, f2str])).strip()
     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: '+filtstr+'\n'
+        err_str += '        CHIP: ' + str(detchip) + '\n'
+        err_str += '     FILTERS: ' + filtstr + '\n'
         ftab.close()
         del ftab
         raise LookupError(err_str)
     else:
-        print('- IDCTAB: Distortion model from row',str(row+1),'for chip',detchip,':',filtstr)
+        print('- IDCTAB: Distortion model from row', str(row + 1), 'for chip',
+              detchip, ':', filtstr)
 
     # Read in V2REF and V3REF: this can either come from current table,
     # or from an OFFTAB if time-dependent (i.e., for WFPC2)
@@ -190,12 +188,12 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
     else:
         # Read V2REF/V3REF from offset table (OFFTAB)
         if offtab:
-            v2ref,v3ref,theta = readOfftab(offtab, date, chip=detchip)
+            v2ref, v3ref, theta = readOfftab(offtab, date, chip=detchip)
         else:
             v2ref = 0.0
             v3ref = 0.0
 
-    if theta == None:
+    if theta is None:
         if 'THETA' in colnames:
             theta = ftab[1].data.field('THETA')[row]
         else:
@@ -206,14 +204,14 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
     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['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
+    refpix['DEFAULT_SCALE'] = True
+    refpix['centered'] = False
     refpix['skew_coeffs'] = skew_coeffs
     # Now that we know which row to look at, read coefficients into the
     #   numeric arrays we have set up...
@@ -226,13 +224,13 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
         cxstr = 'A'
         cystr = 'B'
 
-    for i in range(norder+1):
+    for i in range(norder + 1):
         if i > 0:
-            for j in range(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]
+            for j in range(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
@@ -240,21 +238,21 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
     # 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']:
+    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
-#
-#
-# Time-dependent skew correction coefficients (only ACS/WFC)
-#
-#
+    return fx, fy, refpix, order
+
+
 def read_tdd_coeffs(phdr, chip=1):
-    ''' Read in the TDD related keywords from the PRIMARY header of the IDCTAB
-    '''
+    """
+    Time-dependent skew correction coefficients (only ACS/WFC).
+
+    Read in the TDD related keywords from the PRIMARY header of the IDCTAB
+    """
     # Insure we have an integer form of chip
     ic = int(chip)
 
@@ -269,7 +267,7 @@ def read_tdd_coeffs(phdr, chip=1):
     skew_coeffs['TDD_CX_ALPHA'] = None
 
     # Skew-based TDD coefficients
-    skew_terms = ['TDD_CTB','TDD_CTA','TDD_CYA','TDD_CYB','TDD_CXA','TDD_CXB']
+    skew_terms = ['TDD_CTB', 'TDD_CTA', 'TDD_CYA', 'TDD_CYB', 'TDD_CXA', 'TDD_CXB']
     for s in skew_terms:
         skew_coeffs[s] = None
 
@@ -280,7 +278,7 @@ def read_tdd_coeffs(phdr, chip=1):
         print("Using 2015-calibrated VAFACTOR-corrected TDD correction...")
         skew_coeffs['TDD_DATE'] = phdr['TDD_DATE']
         for s in skew_terms:
-            skew_coeffs[s] = phdr.get('{0}{1}'.format(s,ic),None)
+            skew_coeffs[s] = phdr.get('{0}{1}'.format(s, ic), None)
 
     elif "TDD_CYB1" in phdr:
         # We have 2014-calibrated TDD correction to apply, not J.A.-derived values
@@ -288,17 +286,17 @@ def read_tdd_coeffs(phdr, chip=1):
         skew_coeffs['TDD_DATE'] = phdr['TDD_DATE']
         # Read coefficients for TDD Y coefficient
         cyb_kw = 'TDD_CYB{0}'.format(int(chip))
-        skew_coeffs['TDD_CY_BETA'] = phdr.get(cyb_kw,None)
+        skew_coeffs['TDD_CY_BETA'] = phdr.get(cyb_kw, None)
         cya_kw = 'TDD_CYA{0}'.format(int(chip))
-        tdd_cya = phdr.get(cya_kw,None)
+        tdd_cya = phdr.get(cya_kw, None)
         if tdd_cya == 0 or tdd_cya == 'N/A': tdd_cya = None
         skew_coeffs['TDD_CY_ALPHA'] = tdd_cya
 
         # Read coefficients for TDD X coefficient
         cxb_kw = 'TDD_CXB{0}'.format(int(chip))
-        skew_coeffs['TDD_CX_BETA'] = phdr.get(cxb_kw,None)
+        skew_coeffs['TDD_CX_BETA'] = phdr.get(cxb_kw, None)
         cxa_kw = 'TDD_CXA{0}'.format(int(chip))
-        tdd_cxa = phdr.get(cxa_kw,None)
+        tdd_cxa = phdr.get(cxa_kw, None)
         if tdd_cxa == 0 or tdd_cxa == 'N/A': tdd_cxa = None
         skew_coeffs['TDD_CX_ALPHA'] = tdd_cxa
 
@@ -309,12 +307,12 @@ def read_tdd_coeffs(phdr, chip=1):
             print('TDDORDER kw not present, using default TDD correction')
             return None
 
-        a = np.zeros((n+1,), np.float64)
-        b = np.zeros((n+1,), np.float64)
-        for i in range(n+1):
+        a = np.zeros((n + 1,), np.float64)
+        b = np.zeros((n + 1,), np.float64)
+        for i in range(n + 1):
             a[i] = phdr.get(("TDD_A%d" % i), 0.0)
             b[i] = phdr.get(("TDD_B%d" % i), 0.0)
-        if (a==0).all() and (b==0).all():
+        if (a == 0).all() and (b == 0).all():
             print('Warning: TDD_A and TDD_B coeffiecients have values of 0, \n \
                    but TDDORDER is %d.' % TDDORDER)
 
@@ -325,15 +323,15 @@ def read_tdd_coeffs(phdr, chip=1):
 
     return skew_coeffs
 
-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.
+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 is None:
+        return 0., 0.
 
     # Provide a default value for chip
     if chip:
@@ -347,7 +345,7 @@ def readOfftab(offtab, date, chip=None):
     except:
         raise IOError("Offset table '%s' not valid as specified!" % offtab)
 
-    #Determine row from which to get the coefficients.
+    # 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
@@ -374,7 +372,7 @@ def readOfftab(offtab, date, chip=None):
         obsdate = convertDate(ftab[1].data.field('OBSDATE')[i])
 
         # If the row is appropriate for the chip...
-            # Interpolate between dates
+        # Interpolate between dates
         if int(detchip) == int(chip) or int(detchip) == -999:
             if num_date <= obsdate:
                 date_end = obsdate
@@ -384,7 +382,7 @@ def readOfftab(offtab, date, chip=None):
                 row_end = i
                 continue
 
-            if row_end == None and (num_date > obsdate):
+            if row_end is None and (num_date > obsdate):
                 date_end = obsdate
                 v2end = ftab[1].data.field('V2REF')[i]
                 v3end = ftab[1].data.field('V3REF')[i]
@@ -403,16 +401,17 @@ def readOfftab(offtab, date, chip=None):
     ftab.close()
     del ftab
 
-    if row_start == None and row_end == None:
-        print('Row corresponding to DETCHIP of ',detchip,' was not found!')
+    if row_start is None and row_end is 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))
+    elif row_start is 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))
+        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:
+    if row_start is None or row_end is row_start:
         # We are processing an observation taken after the last calibration
         date_start = date_end
         v2start = v2end
@@ -426,31 +425,32 @@ def readOfftab(offtab, date, chip=None):
     v3ref = _fraction * (v3end - v3start) + v3start
     theta = _fraction * (theta_end - theta_start) + theta_start
 
-    return v2ref,v3ref,theta
+    return v2ref, v3ref, theta
 
-def readWCSCoeffs(header):
-
-    #Read distortion coeffs from WCS header keywords and
-    #populate distortion coeffs arrays.
 
+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)
+    order = max(max(_xorder, _yorder), 3)
 
-    fx = np.zeros(shape=(order+1,order+1),dtype=np.float64)
-    fy = np.zeros(shape=(order+1,order+1),dtype=np.float64)
+    fx = np.zeros(shape=(order + 1, order + 1), dtype=np.float64)
+    fy = np.zeros(shape=(order + 1, order + 1), dtype=np.float64)
 
     # Read in CD matrix
     _cd11 = header['cd1_1']
     _cd12 = header['cd1_2']
     _cd21 = header['cd2_1']
     _cd22 = header['cd2_2']
-    _cdmat = np.array([[_cd11,_cd12],[_cd21,_cd22]])
-    _theta = np.arctan2(-_cd12,_cd22)
-    _rotmat = np.array([[np.cos(_theta),np.sin(_theta)],
-                      [-np.sin(_theta),np.cos(_theta)]])
-    _rCD = np.dot(_rotmat,_cdmat)
+    _cdmat = np.array([[_cd11, _cd12], [_cd21, _cd22]])
+    _theta = np.arctan2(-_cd12, _cd22)
+    _rotmat = np.array([[np.cos(_theta), np.sin(_theta)],
+                        [-np.sin(_theta), np.cos(_theta)]])
+    _rCD = np.dot(_rotmat, _cdmat)
     _skew = np.arcsin(-_rCD[1][0] / _rCD[0][0])
     _scale = _rCD[0][0] * np.cos(_skew) * 3600.
     _scale2 = _rCD[1][1] * 3600.
@@ -467,40 +467,40 @@ def readWCSCoeffs(header):
     refpix['THETA'] = np.rad2deg(_theta)
     refpix['XDELTA'] = 0.0
     refpix['YDELTA'] = 0.0
-    refpix['DEFAULT_SCALE'] = yes
-    refpix['centered'] = yes
-
+    refpix['DEFAULT_SCALE'] = True
+    refpix['centered'] = True
 
     # Set up template for coeffs keyword names
     cxstr = 'A_'
     cystr = 'B_'
     # Read coeffs into their own matrix
-    for i in range(_xorder+1):
-        for j in range(i+1):
-            xcname = cxstr+str(j)+'_'+str(i-j)
+    for i in range(_xorder + 1):
+        for j in range(i + 1):
+            xcname = cxstr + str(j) + '_' + str(i - j)
             if xcname in header:
-                fx[i,j] = header[xcname]
+                fx[i, j] = header[xcname]
 
     # Extract Y coeffs separately as a different order may
     # have been used to fit it.
-    for i in range(_yorder+1):
-        for j in range(i+1):
-            ycname = cystr+str(j)+'_'+str(i-j)
+    for i in range(_yorder + 1):
+        for j in range(i + 1):
+            ycname = cystr + str(j) + '_' + str(i - j)
             if ycname in header:
-                fy[i,j] = header[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
+    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:
+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 is None:
         return fileutil.defaultModel()
 
     # Trauger coefficients only result in a cubic file...
@@ -510,10 +510,10 @@ def readTraugerTable(idcfile,wavelength):
     b_coeffs = [0] * numco
     indx = _MgF2(wavelength)
 
-    ifile = open(idcfile,'r')
+    ifile = open(idcfile, 'r')
     # Search for the first line of the coefficients
     _line = fileutil.rAsciiLine(ifile)
-    while _line[:7].lower() != 'trauger':
+    while _line[: 7].lower() != 'trauger':
         _line = fileutil.rAsciiLine(ifile)
     # Read in each row of coefficients,split them into their values,
     # and convert them into cubic coefficients based on
@@ -525,9 +525,11 @@ def readTraugerTable(idcfile,wavelength):
         if _line == '': continue
         _lc = _line.split()
         if j < 10:
-            a_coeffs[j] = float(_lc[0])+float(_lc[1])*(indx-1.5)+float(_lc[2])*(indx-1.5)**2
+            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
+            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()
@@ -536,17 +538,17 @@ def readTraugerTable(idcfile,wavelength):
     # Now, convert the coefficients into a Numeric array
     # with the right coefficients in the right place.
     # Populate output values now...
-    fx = np.zeros(shape=(order+1,order+1),dtype=np.float64)
-    fy = np.zeros(shape=(order+1,order+1),dtype=np.float64)
+    fx = np.zeros(shape=(order + 1, order + 1), dtype=np.float64)
+    fy = np.zeros(shape=(order + 1, order + 1), dtype=np.float64)
     # Assign the coefficients to their array positions
-    fx[0,0] = 0.
-    fx[1] = np.array([a_coeffs[2],a_coeffs[1],0.,0.],dtype=np.float64)
-    fx[2] = np.array([a_coeffs[5],a_coeffs[4],a_coeffs[3],0.],dtype=np.float64)
-    fx[3] = np.array([a_coeffs[9],a_coeffs[8],a_coeffs[7],a_coeffs[6]],dtype=np.float64)
-    fy[0,0] = 0.
-    fy[1] = np.array([b_coeffs[2],b_coeffs[1],0.,0.],dtype=np.float64)
-    fy[2] = np.array([b_coeffs[5],b_coeffs[4],b_coeffs[3],0.],dtype=np.float64)
-    fy[3] = np.array([b_coeffs[9],b_coeffs[8],b_coeffs[7],b_coeffs[6]],dtype=np.float64)
+    fx[0, 0] = 0.
+    fx[1] = np.array([a_coeffs[2], a_coeffs[1], 0., 0.], dtype=np.float64)
+    fx[2] = np.array([a_coeffs[5], a_coeffs[4], a_coeffs[3], 0.], dtype=np.float64)
+    fx[3] = np.array([a_coeffs[9], a_coeffs[8], a_coeffs[7], a_coeffs[6]], dtype=np.float64)
+    fy[0, 0] = 0.
+    fy[1] = np.array([b_coeffs[2], b_coeffs[1], 0., 0.], dtype=np.float64)
+    fy[2] = np.array([b_coeffs[5], b_coeffs[4], b_coeffs[3], 0.], dtype=np.float64)
+    fy[3] = np.array([b_coeffs[9], b_coeffs[8], b_coeffs[7], b_coeffs[6]], dtype=np.float64)
 
     # Used in Pattern.computeOffsets()
     refpix = {}
@@ -557,10 +559,10 @@ def readTraugerTable(idcfile,wavelength):
     refpix['XDELTA'] = 0.
     refpix['YDELTA'] = 0.
     refpix['PSCALE'] = None
-    refpix['DEFAULT_SCALE'] = no
-    refpix['centered'] = yes
+    refpix['DEFAULT_SCALE'] = False
+    refpix['centered'] = True
 
-    return fx,fy,refpix,order
+    return fx, fy, refpix, order
 
 
 def readCubicTable(idcfile):
@@ -572,17 +574,17 @@ def readCubicTable(idcfile):
 
     # Return a default geometry model if no coefficients filename
     # is given.  This model will not distort the data in any way.
-    if idcfile == None:
+    if idcfile is None:
         return fileutil.defaultModel()
 
-    ifile = open(idcfile,'r')
+    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
+    _found = False
+    while not _found:
+        if _line[:7] in ['cubic', 'quartic', 'quintic'] or _line[: 4] == 'poly':
+            found = True
             break
         _line = fileutil.rAsciiLine(ifile)
 
@@ -613,17 +615,17 @@ def readCubicTable(idcfile):
     # Now, convert the coefficients into a Numeric array
     # with the right coefficients in the right place.
     # Populate output values now...
-    fx = np.zeros(shape=(order+1,order+1),dtype=np.float64)
-    fy = np.zeros(shape=(order+1,order+1),dtype=np.float64)
+    fx = np.zeros(shape=(order + 1, order + 1), dtype=np.float64)
+    fy = np.zeros(shape=(order + 1, order + 1), dtype=np.float64)
     # Assign the coefficients to their array positions
-    fx[0,0] = 0.
-    fx[1] = np.array([a_coeffs[2],a_coeffs[1],0.,0.],dtype=np.float64)
-    fx[2] = np.array([a_coeffs[5],a_coeffs[4],a_coeffs[3],0.],dtype=np.float64)
-    fx[3] = np.array([a_coeffs[9],a_coeffs[8],a_coeffs[7],a_coeffs[6]],dtype=np.float64)
-    fy[0,0] = 0.
-    fy[1] = np.array([b_coeffs[2],b_coeffs[1],0.,0.],dtype=np.float64)
-    fy[2] = np.array([b_coeffs[5],b_coeffs[4],b_coeffs[3],0.],dtype=np.float64)
-    fy[3] = np.array([b_coeffs[9],b_coeffs[8],b_coeffs[7],b_coeffs[6]],dtype=np.float64)
+    fx[0, 0] = 0.
+    fx[1] = np.array([a_coeffs[2], a_coeffs[1], 0., 0.], dtype=np.float64)
+    fx[2] = np.array([a_coeffs[5], a_coeffs[4], a_coeffs[3], 0.], dtype=np.float64)
+    fx[3] = np.array([a_coeffs[9], a_coeffs[8], a_coeffs[7], a_coeffs[6]], dtype=np.float64)
+    fy[0, 0] = 0.
+    fy[1] = np.array([b_coeffs[2], b_coeffs[1], 0., 0.], dtype=np.float64)
+    fy[2] = np.array([b_coeffs[5], b_coeffs[4], b_coeffs[3], 0.], dtype=np.float64)
+    fy[3] = np.array([b_coeffs[9], b_coeffs[8], b_coeffs[7], b_coeffs[6]], dtype=np.float64)
 
     # Used in Pattern.computeOffsets()
     refpix = {}
@@ -634,21 +636,23 @@ def readCubicTable(idcfile):
     refpix['XDELTA'] = 0.
     refpix['YDELTA'] = 0.
     refpix['PSCALE'] = None
-    refpix['DEFAULT_SCALE'] = no
-    refpix['centered'] = yes
+    refpix['DEFAULT_SCALE'] = False
+    refpix['centered'] = True
+
+    return fx, fy, refpix, order
 
-    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)
+        m = m * (i + 1)
     return m
 
-def combin(j,n):
+
+def combin(j, n):
     """ Return the combinatorial factor for j in n."""
-    return (factorial(j) / (factorial(n) * factorial( (j-n) ) ) )
+    return (factorial(j) / (factorial(n) * factorial((j - n))))
 
 
 def defaultModel():
@@ -657,15 +661,15 @@ def defaultModel():
     """
     order = 3
 
-    fx = np.zeros(shape=(order+1,order+1),dtype=np.float64)
-    fy = np.zeros(shape=(order+1,order+1),dtype=np.float64)
+    fx = np.zeros(shape=(order + 1, order + 1), dtype=np.float64)
+    fy = np.zeros(shape=(order + 1, order + 1), dtype=np.float64)
 
-    fx[1,1] = 1.
-    fy[1,0] = 1.
+    fx[1, 1] = 1.
+    fy[1, 0] = 1.
 
     # Used in Pattern.computeOffsets()
     refpix = {}
-    refpix['empty_model'] = yes
+    refpix['empty_model'] = True
     refpix['XREF'] = None
     refpix['YREF'] = None
     refpix['V2REF'] = 0.
@@ -675,17 +679,20 @@ def defaultModel():
     refpix['XDELTA'] = 0.
     refpix['YDELTA'] = 0.
     refpix['PSCALE'] = None
-    refpix['DEFAULT_SCALE'] = no
+    refpix['DEFAULT_SCALE'] = False
     refpix['THETA'] = 0.
-    refpix['centered'] = yes
-    return fx,fy,refpix,order
+    refpix['centered'] = True
+    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 np.sqrt(1.0 + 2.590355e10/(5.312993e10-_sig) +
-        4.4543708e9/(11.17083e9-_sig) + 4.0838897e5/(1.766361e5-_sig))
+    """
+    Function to compute the index of refraction for MgF2 at
+    the specified wavelength for use with Trauger coefficients
+    """
+    _sig = pow((1.0e7 / lam), 2)
+    return np.sqrt(1.0 + 2.590355e10 / (5.312993e10 - _sig) +
+                   4.4543708e9 / (11.17083e9 - _sig) + 4.0838897e5 / (1.766361e5 - _sig))
 
 
 def convertDate(date):
diff --git a/stwcs/distortion/utils.py b/stwcs/distortion/utils.py
index 4228f62..1baad3f 100644
--- a/stwcs/distortion/utils.py
+++ b/stwcs/distortion/utils.py
@@ -1,4 +1,4 @@
-from __future__ import division, print_function # confidence high
+from __future__ import absolute_import, division, print_function
 
 import os
 
@@ -6,11 +6,11 @@ import numpy as np
 from numpy import linalg
 from astropy import wcs as pywcs
 
-from stwcs import wcsutil
-from stwcs import updatewcs
+from .. import updatewcs
 from numpy import sqrt, arctan2
 from stsci.tools import fileutil
 
+
 def output_wcs(list_of_wcsobj, ref_wcs=None, owcs=None, undistort=True):
     """
     Create an output WCS.
@@ -33,63 +33,65 @@ def output_wcs(list_of_wcsobj, ref_wcs=None, owcs=None, undistort=True):
 
     # This new algorithm may not be strictly necessary, but it may be more
     # robust in handling regions near the poles or at 0h RA.
-    crval1,crval2 = computeFootprintCenter(fra_dec)
+    crval1, crval2 = computeFootprintCenter(fra_dec)
 
-    crval = np.array([crval1,crval2], dtype=np.float64) # this value is now zero-based
+    crval = np.array([crval1, crval2], dtype=np.float64)  # this value is now zero-based
     if owcs is None:
         if ref_wcs is None:
             ref_wcs = list_of_wcsobj[0].deepcopy()
         if undistort:
-            #outwcs = undistortWCS(ref_wcs)
+            # outwcs = undistortWCS(ref_wcs)
             outwcs = make_orthogonal_cd(ref_wcs)
         else:
             outwcs = ref_wcs.deepcopy()
         outwcs.wcs.crval = crval
         outwcs.wcs.set()
-        outwcs.pscale = sqrt(outwcs.wcs.cd[0,0]**2 + outwcs.wcs.cd[1,0]**2)*3600.
-        outwcs.orientat = arctan2(outwcs.wcs.cd[0,1],outwcs.wcs.cd[1,1]) * 180./np.pi
+        outwcs.pscale = sqrt(outwcs.wcs.cd[0, 0] ** 2 + outwcs.wcs.cd[1, 0] ** 2) * 3600.
+        outwcs.orientat = arctan2(outwcs.wcs.cd[0, 1], outwcs.wcs.cd[1, 1]) * 180. / np.pi
     else:
         outwcs = owcs.deepcopy()
-        outwcs.pscale = sqrt(outwcs.wcs.cd[0,0]**2 + outwcs.wcs.cd[1,0]**2)*3600.
-        outwcs.orientat = arctan2(outwcs.wcs.cd[0,1],outwcs.wcs.cd[1,1]) * 180./np.pi
+        outwcs.pscale = sqrt(outwcs.wcs.cd[0, 0] ** 2 + outwcs.wcs.cd[1, 0] ** 2) * 3600.
+        outwcs.orientat = arctan2(outwcs.wcs.cd[0, 1], outwcs.wcs.cd[1, 1]) * 180. / np.pi
 
     tanpix = outwcs.wcs.s2p(fra_dec, 0)['pixcrd']
 
-    outwcs._naxis1 = int(np.ceil(tanpix[:,0].max() - tanpix[:,0].min()))
-    outwcs._naxis2 = int(np.ceil(tanpix[:,1].max() - tanpix[:,1].min()))
-    crpix = np.array([outwcs._naxis1/2., outwcs._naxis2/2.], dtype=np.float64)
+    outwcs._naxis1 = int(np.ceil(tanpix[:, 0].max() - tanpix[:, 0].min()))
+    outwcs._naxis2 = int(np.ceil(tanpix[:, 1].max() - tanpix[:, 1].min()))
+    crpix = np.array([outwcs._naxis1 / 2., outwcs._naxis2 / 2.], dtype=np.float64)
     outwcs.wcs.crpix = crpix
     outwcs.wcs.set()
     tanpix = outwcs.wcs.s2p(fra_dec, 0)['pixcrd']
 
     # shift crpix to take into account (floating-point value of) position of
     # corner pixel relative to output frame size: no rounding necessary...
-    newcrpix = np.array([crpix[0]+tanpix[:,0].min(), crpix[1]+
-                         tanpix[:,1].min()])
+    newcrpix = np.array([crpix[0] + tanpix[:, 0].min(), crpix[1] +
+                         tanpix[:, 1].min()])
 
     newcrval = outwcs.wcs.p2s([newcrpix], 1)['world'][0]
     outwcs.wcs.crval = newcrval
     outwcs.wcs.set()
-    outwcs.wcs.name = wcsname # keep track of label for this solution
+    outwcs.wcs.name = wcsname  # keep track of label for this solution
     return outwcs
 
+
 def computeFootprintCenter(edges):
     """ Geographic midpoint in spherical coords for points defined by footprints.
         Algorithm derived from: http://www.geomidpoint.com/calculation.html
 
         This algorithm should be more robust against discontinuities at the poles.
     """
-    alpha = np.deg2rad(edges[:,0])
-    dec = np.deg2rad(edges[:,1])
+    alpha = np.deg2rad(edges[:, 0])
+    dec = np.deg2rad(edges[:, 1])
 
-    xmean = np.mean(np.cos(dec)*np.cos(alpha))
-    ymean = np.mean(np.cos(dec)*np.sin(alpha))
+    xmean = np.mean(np.cos(dec) * np.cos(alpha))
+    ymean = np.mean(np.cos(dec) * np.sin(alpha))
     zmean = np.mean(np.sin(dec))
 
-    crval1 = np.rad2deg(np.arctan2(ymean,xmean))%360.0
-    crval2 = np.rad2deg(np.arctan2(zmean,np.sqrt(xmean*xmean+ymean*ymean)))
+    crval1 = np.rad2deg(np.arctan2(ymean, xmean)) % 360.0
+    crval2 = np.rad2deg(np.arctan2(zmean, np.sqrt(xmean * xmean + ymean * ymean)))
+
+    return crval1, crval2
 
-    return crval1,crval2
 
 def make_orthogonal_cd(wcs):
     """ Create a perfect (square, orthogonal, undistorted) CD matrix from the
@@ -98,21 +100,20 @@ def make_orthogonal_cd(wcs):
     # get determinant of the CD matrix:
     cd = wcs.celestial.pixel_scale_matrix
 
-
     if hasattr(wcs, 'idcv2ref') and wcs.idcv2ref is not None:
         # 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 = updatewcs.makewcs.troll(wcs.pav3,wcs.wcs.crval[1],wcs.idcv2ref,wcs.idcv3ref)
+        pv = updatewcs.makewcs.troll(wcs.pav3, wcs.wcs.crval[1], wcs.idcv2ref, wcs.idcv3ref)
         # Add the chip rotation angle
         if wcs.idctheta:
             pv += wcs.idctheta
         cs = np.cos(np.deg2rad(pv))
         sn = np.sin(np.deg2rad(pv))
-        pvmat = np.dot(np.array([[cs,sn],[-sn,cs]]),wcs.parity)
-        rot = np.arctan2(pvmat[0,1],pvmat[1,1])
-        scale = wcs.idcscale/3600.
+        pvmat = np.dot(np.array([[cs, sn], [-sn, cs]]), wcs.parity)
+        rot = np.arctan2(pvmat[0, 1], pvmat[1, 1])
+        scale = wcs.idcscale / 3600.
 
         det = linalg.det(wcs.parity)
 
@@ -122,36 +123,37 @@ def make_orthogonal_cd(wcs):
 
         # find pixel scale:
         if hasattr(wcs, 'idcscale'):
-            scale = (wcs.idcscale) / 3600. # HST pixel scale provided
+            scale = (wcs.idcscale) / 3600.  # HST pixel scale provided
         else:
-            scale = np.sqrt(np.abs(det)) # find as sqrt(pixel area)
+            scale = np.sqrt(np.abs(det))  # find as sqrt(pixel area)
 
         # find Y-axis orientation:
         if hasattr(wcs, 'orientat') and not ignoreHST:
-            rot = np.deg2rad(wcs.orientat) # use HST ORIENTAT
+            rot = np.deg2rad(wcs.orientat)  # use HST ORIENTAT
         else:
-            rot = np.arctan2(wcs.wcs.cd[0,1], wcs.wcs.cd[1,1]) # angle of the Y-axis
+            rot = np.arctan2(wcs.wcs.cd[0, 1], wcs.wcs.cd[1, 1])  # angle of the Y-axis
 
     par = -1 if det < 0.0 else 1
 
     # create a perfectly square, orthogonal WCS
     sn = np.sin(rot)
     cs = np.cos(rot)
-    orthogonal_cd = scale * np.array([[par*cs, sn], [-par*sn, cs]])
+    orthogonal_cd = scale * np.array([[par * cs, sn], [-par * sn, cs]])
 
     lin_wcsobj = pywcs.WCS()
     lin_wcsobj.wcs.cd = orthogonal_cd
     lin_wcsobj.wcs.set()
-    lin_wcsobj.orientat = arctan2(lin_wcsobj.wcs.cd[0,1],lin_wcsobj.wcs.cd[1,1]) * 180./np.pi
-    lin_wcsobj.pscale = sqrt(lin_wcsobj.wcs.cd[0,0]**2 + lin_wcsobj.wcs.cd[1,0]**2)*3600.
-    lin_wcsobj.wcs.crval = np.array([0.,0.])
-    lin_wcsobj.wcs.crpix = np.array([0.,0.])
+    lin_wcsobj.orientat = arctan2(lin_wcsobj.wcs.cd[0, 1], lin_wcsobj.wcs.cd[1, 1]) * 180. / np.pi
+    lin_wcsobj.pscale = sqrt(lin_wcsobj.wcs.cd[0, 0] ** 2 + lin_wcsobj.wcs.cd[1, 0] ** 2) * 3600.
+    lin_wcsobj.wcs.crval = np.array([0., 0.])
+    lin_wcsobj.wcs.crpix = np.array([0., 0.])
     lin_wcsobj.wcs.ctype = ['RA---TAN', 'DEC--TAN']
     lin_wcsobj.wcs.set()
 
     return lin_wcsobj
 
-def  undistortWCS(wcsobj):
+
+def undistortWCS(wcsobj):
     """
     Creates an undistorted linear WCS by applying the IDCTAB distortion model
     to a 3-point square. The new ORIENTAT angle is calculated as well as the
@@ -175,25 +177,26 @@ def  undistortWCS(wcsobj):
             return wcsobj.copy()
     crpix1 = wcsobj.wcs.crpix[0]
     crpix2 = wcsobj.wcs.crpix[1]
-    xy = np.array([(crpix1,crpix2),(crpix1+1.,crpix2),(crpix1,crpix2+1.)],dtype=np.double)
+    xy = np.array([(crpix1, crpix2), (crpix1 + 1., crpix2),
+                   (crpix1, crpix2 + 1.)], dtype=np.double)
     offsets = np.array([wcsobj.ltv1, wcsobj.ltv2])
     px = xy + offsets
-    #order = wcsobj.sip.a_order
+    # order = wcsobj.sip.a_order
     pscale = wcsobj.idcscale
-    #pixref = np.array([wcsobj.sip.SIPREF1, wcsobj.sip.SIPREF2])
+    # pixref = np.array([wcsobj.sip.SIPREF1, wcsobj.sip.SIPREF2])
 
     tan_pix = apply_idc(px, cx, cy, wcsobj.wcs.crpix, pscale, order=1)
-    xc = tan_pix[:,0]
-    yc = tan_pix[:,1]
+    xc = tan_pix[:, 0]
+    yc = tan_pix[:, 1]
     am = xc[1] - xc[0]
     bm = xc[2] - xc[0]
     cm = yc[1] - yc[0]
     dm = yc[2] - yc[0]
-    cd_mat = np.array([[am,bm],[cm,dm]],dtype=np.double)
+    cd_mat = np.array([[am, bm], [cm, dm]], dtype=np.double)
 
     # Check the determinant for singularity
     _det = (am * dm) - (bm * cm)
-    if ( _det == 0.0):
+    if (_det == 0.0):
         print('Singular matrix in updateWCS, aborting ...')
         return
 
@@ -202,15 +205,16 @@ def  undistortWCS(wcsobj):
     cd = np.dot(wcsobj.wcs.cd, cd_inv).astype(np.float64)
     lin_wcsobj.wcs.cd = cd
     lin_wcsobj.wcs.set()
-    lin_wcsobj.orientat = arctan2(lin_wcsobj.wcs.cd[0,1],lin_wcsobj.wcs.cd[1,1]) * 180./np.pi
-    lin_wcsobj.pscale = sqrt(lin_wcsobj.wcs.cd[0,0]**2 + lin_wcsobj.wcs.cd[1,0]**2)*3600.
-    lin_wcsobj.wcs.crval = np.array([0.,0.])
-    lin_wcsobj.wcs.crpix = np.array([0.,0.])
+    lin_wcsobj.orientat = arctan2(lin_wcsobj.wcs.cd[0, 1], lin_wcsobj.wcs.cd[1, 1]) * 180. / np.pi
+    lin_wcsobj.pscale = sqrt(lin_wcsobj.wcs.cd[0, 0] ** 2 + lin_wcsobj.wcs.cd[1, 0] ** 2) * 3600.
+    lin_wcsobj.wcs.crval = np.array([0., 0.])
+    lin_wcsobj.wcs.crpix = np.array([0., 0.])
     lin_wcsobj.wcs.ctype = ['RA---TAN', 'DEC--TAN']
     lin_wcsobj.wcs.set()
     return lin_wcsobj
 
-def apply_idc(pixpos, cx, cy, pixref, pscale= None, order=None):
+
+def apply_idc(pixpos, cx, cy, pixref, pscale=None, order=None):
     """
     Apply the IDCTAB polynomial distortion model to pixel positions.
     pixpos must be already corrected for ltv1/2.
@@ -233,27 +237,28 @@ def apply_idc(pixpos, cx, cy, pixref, pscale= None, order=None):
         return pixpos
 
     # Apply in the same way that 'drizzle' would...
-    _cx = cx/pscale
-    _cy = cy/ pscale
+    _cx = cx / pscale
+    _cy = cy / pscale
     _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.
+    _cx[0, 0] = 0.
+    _cy[0, 0] = 0.
 
     dxy = _p - pixref
     # 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))
+    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))
+
+    return c
 
-    return  c
 
 def foundIDCTAB(idctab):
     idctab_found = True
diff --git a/stwcs/gui/apply_headerlet.py b/stwcs/gui/apply_headerlet.py
index d517e9f..0cef612 100644
--- a/stwcs/gui/apply_headerlet.py
+++ b/stwcs/gui/apply_headerlet.py
@@ -1,31 +1,34 @@
+from __future__ import absolute_import, division, print_function
 import os
 from stsci.tools import teal, parseinput
 
-import stwcs
-from stwcs.wcsutil import headerlet
+from .. import __version__
+from ..wcsutil import headerlet
 
-__taskname__ = __name__.split('.')[-1] # needed for help string
+__taskname__ = __name__.split('.')[-1]  # needed for help string
 __package__ = headerlet.__name__
-__version__ = stwcs.__version__
-#
-#### Interfaces used by TEAL
-#
+# __version__ = stwcs.__version__
+
+
+############### Interfaces used by TEAL ###############
+
 def getHelpAsString(docstring=False):
     """
     return useful help from a file in the script directory called __taskname__.help
     """
     install_dir = os.path.dirname(__file__)
-    htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html')
-    helpfile = os.path.join(install_dir,__taskname__+'.help')
+    htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html')
+    helpfile = os.path.join(install_dir, __taskname__ + '.help')
     if docstring or (not docstring and not os.path.exists(htmlfile)):
-        helpString = __taskname__+' Version '+__version__+'\n\n'
+        helpString = __taskname__ + ' Version ' + __version__ + '\n\n'
         if os.path.exists(helpfile):
-            helpString += teal.getHelpFileAsString(__taskname__,__file__)
+            helpString += teal.getHelpFileAsString(__taskname__, __file__)
     else:
-        helpString = 'file://'+htmlfile
+        helpString = 'file://' + htmlfile
 
     return helpString
 
+
 def run(configObj=None):
 
     # start by interpreting filename and hdrlet inputs
@@ -37,18 +40,18 @@ def run(configObj=None):
         # Syntax: apply_headerlet_as_primary(filename, hdrlet, attach=True,
         #            archive=True, force=False, verbose=False)
         headerlet.apply_headerlet_as_primary(filename,
-                        hdrlet,attach=configObj['attach'],
-                        archive=configObj['archive'],force=configObj['force'],
-                        logging=configObj['logging'])
+                                             hdrlet, attach=configObj['attach'],
+                                             archive=configObj['archive'], force=configObj['force'],
+                                             logging=configObj['logging'])
     else:
         wcsname = configObj['wcsname']
-        if wcsname in ['',' ','INDEF']: wcsname = None
+        if wcsname in ['', ' ', 'INDEF']: wcsname = None
         wcskey = configObj['wcskey']
         if wcskey == '': wcskey = None
         # Call function with properly interpreted input parameters
         #         apply_headerlet_as_alternate(filename, hdrlet, attach=True,
         #                        wcskey=None, wcsname=None, verbose=False)
         headerlet.apply_headerlet_as_alternate(filename,
-                        hdrlet, attach=configObj['attach'],
-                        wcsname=wcsname, wcskey=wcskey,
-                        logging=configObj['logging'])
+                                               hdrlet, attach=configObj['attach'],
+                                               wcsname=wcsname, wcskey=wcskey,
+                                               logging=configObj['logging'])
diff --git a/stwcs/gui/archive_headerlet.py b/stwcs/gui/archive_headerlet.py
index 7ad3d4d..577a348 100644
--- a/stwcs/gui/archive_headerlet.py
+++ b/stwcs/gui/archive_headerlet.py
@@ -1,49 +1,49 @@
-from __future__ import print_function
+from __future__ import absolute_import, division, print_function
 import os
 
-from astropy.io import fits
 from stsci.tools import teal
+from .. import __version__
+from ..wcsutil import headerlet
 
-import stwcs
-from stwcs.wcsutil import headerlet
-
-__taskname__ = __name__.split('.')[-1] # needed for help string
+__taskname__ = __name__.split('.')[-1]  # needed for help string
 __package__ = headerlet.__name__
-__version__ = stwcs.__version__
+# __version__ = stwcs.__version__
+
 #
 #### Interfaces used by TEAL
 #
 def getHelpAsString(docstring=False):
     """
-    return useful help from a file in the script directory called __taskname__.help
+    Return useful help from a file in the script directory called __taskname__.help
     """
     install_dir = os.path.dirname(__file__)
-    htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html')
-    helpfile = os.path.join(install_dir,__taskname__+'.help')
+    htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html')
+    helpfile = os.path.join(install_dir, __taskname__ + '.help')
     if docstring or (not docstring and not os.path.exists(htmlfile)):
-        helpString = __taskname__+' Version '+__version__+'\n\n'
+        helpString = __taskname__ + ' Version ' + __version__ + '\n\n'
         if os.path.exists(helpfile):
-            helpString += teal.getHelpFileAsString(__taskname__,__file__)
+            helpString += teal.getHelpFileAsString(__taskname__, __file__)
         else:
             helpString += headerlet.archive_as_headerlet.__doc__
 
     else:
-        helpString = 'file://'+htmlfile
+        helpString = 'file://' + htmlfile
 
     return helpString
 
+
 def run(configObj=None):
 
-    if configObj['hdrname'] in ['',' ','INDEF']:
-        print('='*60)
+    if configObj['hdrname'] in ['', ' ', 'INDEF']:
+        print('=' * 60)
         print('ERROR:')
         print('    No valid "hdrname" parameter value provided!')
         print('    Please restart this task and provide a value for this parameter.')
-        print('='*60)
+        print('=' * 60)
         return
 
-    str_kw = ['wcsname','destim','sipname','npolfile','d2imfile',
-            'descrip','history','author']
+    str_kw = ['wcsname', 'destim', 'sipname', 'npolfile', 'd2imfile',
+              'descrip', 'history', 'author']
 
     # create dictionary of remaining parameters, deleting extraneous ones
     # such as those above
@@ -66,4 +66,4 @@ def run(configObj=None):
     #                    author=None, descrip=None, history=None,
     #                    hdrlet=None, clobber=False)
     headerlet.archive_as_headerlet(configObj['filename'], configObj['hdrname'],
-                                **cdict)
+                                   **cdict)
diff --git a/stwcs/gui/attach_headerlet.py b/stwcs/gui/attach_headerlet.py
index 873c549..f3ab690 100644
--- a/stwcs/gui/attach_headerlet.py
+++ b/stwcs/gui/attach_headerlet.py
@@ -1,13 +1,13 @@
+from __future__ import absolute_import, division, print_function
 import os
 
 from stsci.tools import teal
+from .. import __version
+from ..wcsutil import headerlet
 
-import stwcs
-from stwcs.wcsutil import headerlet
-
-__taskname__ = __name__.split('.')[-1] # needed for help string
+__taskname__ = __name__.split('.')[-1]  # needed for help string
 __package__ = headerlet.__name__
-__version__ = stwcs.__version__
+#__version__ = stwcs.__version__
 #
 #### Interfaces used by TEAL
 #
@@ -16,21 +16,21 @@ def getHelpAsString(docstring=False):
     return useful help from a file in the script directory called __taskname__.help
     """
     install_dir = os.path.dirname(__file__)
-    htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html')
-    helpfile = os.path.join(install_dir,__taskname__+'.help')
+    htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html')
+    helpfile = os.path.join(install_dir, __taskname__ + '.help')
     if docstring or (not docstring and not os.path.exists(htmlfile)):
-        helpString = __taskname__+' Version '+__version__+'\n\n'
+        helpString = __taskname__ + ' Version ' + __version__ + '\n\n'
         if os.path.exists(helpfile):
-            helpString += teal.getHelpFileAsString(__taskname__,__file__)
+            helpString += teal.getHelpFileAsString(__taskname__, __file__)
         else:
-            helpString += eval('.'.join([__package__,__taskname__,'__doc__']))
+            helpString += eval('.'.join([__package__, __taskname__, '__doc__']))
     else:
-        helpString = 'file://'+htmlfile
+        helpString = 'file://' + htmlfile
 
     return helpString
 
+
 def run(configObj=None):
 
-    headerlet.attach_headerlet(configObj['filename'],configObj['hdrlet'],
+    headerlet.attach_headerlet(configObj['filename'], configObj['hdrlet'],
                                configObj['logging'])
-
diff --git a/stwcs/gui/delete_headerlet.py b/stwcs/gui/delete_headerlet.py
index b3df5a7..29937f7 100644
--- a/stwcs/gui/delete_headerlet.py
+++ b/stwcs/gui/delete_headerlet.py
@@ -1,15 +1,15 @@
-from __future__ import print_function
-import os
+from __future__ import absolute_import, division, print_function
 
+import os
 from stsci.tools import teal
 from stsci.tools import parseinput
 
-import stwcs
-from stwcs.wcsutil import headerlet
+from .. import __version__
+from ..wcsutil import headerlet
 
-__taskname__ = __name__.split('.')[-1] # needed for help string
+__taskname__ = __name__.split('.')[-1]  # needed for help string
 __package__ = headerlet.__name__
-__version__ = stwcs.__version__
+# __version__ = stwcs.__version__
 #
 #### Interfaces used by TEAL
 #
@@ -18,35 +18,34 @@ def getHelpAsString(docstring=False):
     return useful help from a file in the script directory called __taskname__.help
     """
     install_dir = os.path.dirname(__file__)
-    htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html')
-    helpfile = os.path.join(install_dir,__taskname__+'.help')
+    htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html')
+    helpfile = os.path.join(install_dir, __taskname__ + '.help')
     if docstring or (not docstring and not os.path.exists(htmlfile)):
-        helpString = __taskname__+' Version '+__version__+'\n\n'
+        helpString = __taskname__ + ' Version ' + __version__ + '\n\n'
         if os.path.exists(helpfile):
-            helpString += teal.getHelpFileAsString(__taskname__,__file__)
+            helpString += teal.getHelpFileAsString(__taskname__, __file__)
         else:
-            helpString += eval('.'.join([__package__,__taskname__,'__doc__']))
+            helpString += eval('.'.join([__package__, __taskname__, '__doc__']))
     else:
-        helpString = 'file://'+htmlfile
+        helpString = 'file://' + htmlfile
 
     return helpString
 
 def run(configObj=None):
 
     if configObj['hdrname'] == '' and configObj['hdrext'] is None and \
-        configObj['distname'] == '':
-        print('='*60)
+            configObj['distname'] == '':
+        print('=' * 60)
         print('ERROR:')
         print('    No valid "hdrname", "hdrext" or "distname" parameter value provided!')
         print('    Please restart this task and provide a value for one of these parameters.')
-        print('='*60)
+        print('=' * 60)
         return
     filename = parseinput.parseinput(configObj['filename'])[0]
     # Call function with properly interpreted input parameters
     # Syntax: delete_headerlet(filename, hdrname=None, hdrext=None, distname=None)
     headerlet.delete_headerlet(filename,
-                                hdrname = configObj['hdrname'],
-                                hdrext = configObj['hdrext'],
-                                distname = configObj['distname'],
-                                logging = configObj['logging'])
-
+                               hdrname=configObj['hdrname'],
+                               hdrext=configObj['hdrext'],
+                               distname=configObj['distname'],
+                               logging=configObj['logging'])
diff --git a/stwcs/gui/extract_headerlet.py b/stwcs/gui/extract_headerlet.py
index 02ecd7a..1e88221 100644
--- a/stwcs/gui/extract_headerlet.py
+++ b/stwcs/gui/extract_headerlet.py
@@ -1,14 +1,14 @@
-from __future__ import print_function
+from __future__ import absolute_import, division, print_function
+
 import os
 
 from stsci.tools import teal
-
-import stwcs
-from stwcs.wcsutil import headerlet
+from .. import __version__
+from ..wcsutil import headerlet
 
 __taskname__ = __name__.split('.')[-1] # needed for help string
 __package__ = headerlet.__name__
-__version__ = stwcs.__version__
+# __version__ = stwcs.__version__
 #
 #### Interfaces used by TEAL
 #
@@ -17,28 +17,29 @@ def getHelpAsString(docstring=False):
     return useful help from a file in the script directory called __taskname__.help
     """
     install_dir = os.path.dirname(__file__)
-    htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html')
-    helpfile = os.path.join(install_dir,__taskname__+'.help')
+    htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html')
+    helpfile = os.path.join(install_dir, __taskname__ + '.help')
     if docstring or (not docstring and not os.path.exists(htmlfile)):
-        helpString = __taskname__+' Version '+__version__+'\n\n'
+        helpString = __taskname__ + ' Version ' + __version__ + '\n\n'
         if os.path.exists(helpfile):
-            helpString += teal.getHelpFileAsString(__taskname__,__file__)
+            helpString += teal.getHelpFileAsString(__taskname__, __file__)
         else:
-            helpString += eval('.'.join([__package__,__taskname__,'__doc__']))
+            helpString += eval('.'.join([__package__, __taskname__, '__doc__']))
 
     else:
-        helpString = 'file://'+htmlfile
+        helpString = 'file://' + htmlfile
 
     return helpString
 
+
 def run(configObj=None):
 
-    if configObj['output'] in ['',' ',None]:
-        print('='*60)
+    if configObj['output'] in ['', ' ', None]:
+        print('=' * 60)
         print('ERROR:')
         print('    No valid "output" parameter value provided!')
         print('    Please restart this task and provide a value for this parameter.')
-        print('='*60)
+        print('=' * 60)
         return
 
     # create dictionary of remaining parameters, deleting extraneous ones
@@ -55,4 +56,3 @@ def run(configObj=None):
     #                clobber=False, verbose=100)
     headerlet.extract_headerlet(configObj['filename'], configObj['output'],
                                 **cdict)
-
diff --git a/stwcs/gui/headerlet_summary.py b/stwcs/gui/headerlet_summary.py
index 82a3e0c..eac59ee 100644
--- a/stwcs/gui/headerlet_summary.py
+++ b/stwcs/gui/headerlet_summary.py
@@ -1,12 +1,13 @@
+from __future__ import absolute_import, division, print_function
 import os
 from stsci.tools import teal
 
-import stwcs
-from stwcs.wcsutil import headerlet
+from .. import __version__
+from ..wcsutil import headerlet
 
 __taskname__ = __name__.split('.')[-1] # needed for help string
 __package__ = headerlet.__name__
-__version__ = stwcs.__version__
+# __version__ = stwcs.__version__
 #
 #### Interfaces used by TEAL
 #
@@ -15,19 +16,20 @@ def getHelpAsString(docstring=False):
     return useful help from a file in the script directory called __taskname__.help
     """
     install_dir = os.path.dirname(__file__)
-    htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html')
-    helpfile = os.path.join(install_dir,__taskname__+'.help')
+    htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html')
+    helpfile = os.path.join(install_dir, __taskname__ + '.help')
     if docstring or (not docstring and not os.path.exists(htmlfile)):
-        helpString = __taskname__+' Version '+__version__+'\n\n'
+        helpString = __taskname__ + ' Version ' + __version__ + '\n\n'
         if os.path.exists(helpfile):
-            helpString += teal.getHelpFileAsString(__taskname__,__file__)
+            helpString += teal.getHelpFileAsString(__taskname__, __file__)
         else:
-            helpString += eval('.'.join([__package__,__taskname__,'__doc__']))
+            helpString += eval('.'.join([__package__, __taskname__, '__doc__']))
     else:
-        helpString = 'file://'+htmlfile
+        helpString = 'file://' + htmlfile
 
     return helpString
 
+
 def run(configObj=None):
 
     # create dictionary of remaining parameters, deleting extraneous ones
@@ -43,5 +45,4 @@ def run(configObj=None):
 
     # Syntax: headerlet_summary(filename,columns=None,pad=2,maxwidth=None,
     #                    output=None,clobber=True,quiet=False)
-    headerlet.headerlet_summary(configObj['filename'],**cdict)
-
+    headerlet.headerlet_summary(configObj['filename'], **cdict)
diff --git a/stwcs/gui/restore_headerlet.py b/stwcs/gui/restore_headerlet.py
index 7570d76..790c239 100644
--- a/stwcs/gui/restore_headerlet.py
+++ b/stwcs/gui/restore_headerlet.py
@@ -1,13 +1,14 @@
+from __future__ import absolute_import, division, print_function
 import os
 
 from stsci.tools import teal
 
-import stwcs
-from stwcs.wcsutil import headerlet
+from .. import __version__
+from ..wcsutil import headerlet
 
 __taskname__ = __name__.split('.')[-1] # needed for help string
 __package__ = headerlet.__name__
-__version__ = stwcs.__version__
+#__version__ = stwcs.__version__
 #
 #### Interfaces used by TEAL
 #
@@ -16,33 +17,33 @@ def getHelpAsString(docstring=False):
     return useful help from a file in the script directory called __taskname__.help
     """
     install_dir = os.path.dirname(__file__)
-    htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html')
-    helpfile = os.path.join(install_dir,__taskname__+'.help')
+    htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html')
+    helpfile = os.path.join(install_dir, __taskname__ + '.help')
     if docstring or (not docstring and not os.path.exists(htmlfile)):
-        helpString = __taskname__+' Version '+__version__+'\n\n'
+        helpString = __taskname__ + ' Version ' + __version__ + '\n\n'
         if os.path.exists(helpfile):
-            helpString += teal.getHelpFileAsString(__taskname__,__file__)
+            helpString += teal.getHelpFileAsString(__taskname__, __file__)
     else:
-        helpString = 'file://'+htmlfile
+        helpString = 'file://' + htmlfile
 
     return helpString
 
+
 def run(configObj=None):
 
-    if configObj['distname'] not in ['',' ','INDEF']:
+    if configObj['distname'] not in ['', ' ', 'INDEF']:
         # Call function with properly interpreted input parameters
         # Syntax: restore_all_with_distname(filename, distname, primary,
         #           archive=True, sciext='SCI', verbose=False)
         headerlet.restore_all_with_distname(configObj['filename'],
-                        configObj['distname'],configObj['primary'],
-                        archive=configObj['archive'],sciext=configObj['sciext'],
-                        logging=configObj['logging'])
+                                            configObj['distname'], configObj['primary'],
+                                            archive=configObj['archive'], sciext=configObj['sciext'],
+                                            logging=configObj['logging'])
     else:
         # Call function with properly interpreted input parameters
         #         restore_from_headerlet(filename, hdrname=None, hdrext=None,
         #           archive=True, force=False)
         headerlet.restore_from_headerlet(configObj['filename'],
-                        hdrname=configObj['hdrname'],hdrext=configObj['hdrext'],
-                        archive=configObj['archive'], force=configObj['force'],
-                        logging=configObj['logging'])
-
+                                         hdrname=configObj['hdrname'], hdrext=configObj['hdrext'],
+                                         archive=configObj['archive'], force=configObj['force'],
+                                         logging=configObj['logging'])
diff --git a/stwcs/gui/updatewcs.py b/stwcs/gui/updatewcs.py
index 3dacb67..b5dc700 100644
--- a/stwcs/gui/updatewcs.py
+++ b/stwcs/gui/updatewcs.py
@@ -1,19 +1,22 @@
-from __future__ import print_function
+from __future__ import absolute_import, division, print_function
+
 import os
 
 from astropy.io import fits
 from stsci.tools import parseinput
 from stsci.tools import fileutil
 from stsci.tools import teal
-import stwcs
-from stwcs import updatewcs
-from stwcs.wcsutil import convertwcs
+from .. import __version__
+from .. import updatewcs
+
+
+allowed_corr_dict = {'vacorr': 'VACorr', 'tddcorr': 'TDDCorr', 'npolcorr': 'NPOLCorr',
+                     'd2imcorr': 'DET2IMCorr'}
 
-allowed_corr_dict = {'vacorr':'VACorr','tddcorr':'TDDCorr','npolcorr':'NPOLCorr','d2imcorr':'DET2IMCorr'}
 
 __taskname__ = __name__.split('.')[-1] # needed for help string
 __package__ = updatewcs.__name__
-__version__ = stwcs.__version__
+#__version__ = stwcs.__version__
 
 #
 #### Interfaces used by TEAL
@@ -23,24 +26,24 @@ def getHelpAsString(docstring=False):
     return useful help from a file in the script directory called __taskname__.help
     """
     install_dir = os.path.dirname(__file__)
-    htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html')
-    helpfile = os.path.join(install_dir,__taskname__+'.help')
+    htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html')
+    helpfile = os.path.join(install_dir, __taskname__ + '.help')
     if docstring or (not docstring and not os.path.exists(htmlfile)):
-        helpString = __taskname__+' Version '+__version__+'\n\n'
+        helpString = __taskname__ + ' Version ' + __version__ + '\n\n'
         if os.path.exists(helpfile):
-            helpString += teal.getHelpFileAsString(__taskname__,__file__)
+            helpString += teal.getHelpFileAsString(__taskname__, __file__)
         else:
-            helpString += eval('.'.join([__package__,__taskname__,'__doc__']))
+            helpString += eval('.'.join([__package__, __taskname__, '__doc__']))
 
     else:
-        helpString = 'file://'+htmlfile
+        helpString = 'file://' + htmlfile
 
     return helpString
 
+
 def run(configObj=None):
 
     # Interpret primary parameters from configObj instance
-    extname = configObj['extname']
     input = configObj['input']
 
     # create dictionary of remaining parameters, deleting extraneous ones
@@ -53,7 +56,7 @@ def run(configObj=None):
     del cdict['extname']
 
     # parse input
-    input,altfiles = parseinput.parseinput(configObj['input'])
+    input, altfiles = parseinput.parseinput(configObj['input'])
 
     # Insure that all input files have a correctly archived
     #    set of OPUS WCS keywords
@@ -69,7 +72,7 @@ def run(configObj=None):
         # Check to insure that there is a valid reference file to be used
         idctab = fits.getval(file, 'idctab')
         if not os.path.exists(fileutil.osfn(idctab)):
-            print('No valid distortion reference file ',idctab,' found in ',file,'!')
+            print('No valid distortion reference file ', idctab, ' found in ', file, '!')
             raise ValueError
 
     # Re-define 'cdict' to only have switches for steps supported by that instrument
@@ -78,7 +81,7 @@ def run(configObj=None):
     #
     for file in input:
         # get instrument name from input file
-        instr = fits.getval(file,'INSTRUME')
+        instr = fits.getval(file, 'INSTRUME')
         # make copy of input parameters dict for this file
         fdict = cdict.copy()
         # Remove any parameter that is not part of this instrument's allowed corrections
@@ -86,5 +89,4 @@ def run(configObj=None):
             if allowed_corr_dict[step] not in updatewcs.apply_corrections.allowed_corrections[instr]:
                 fdict[step]
         # Call 'updatewcs' on correctly archived file
-        updatewcs.updatewcs(file,**fdict)
-
+        updatewcs.updatewcs(file, **fdict)
diff --git a/stwcs/gui/write_headerlet.py b/stwcs/gui/write_headerlet.py
index e18bed8..56520b3 100644
--- a/stwcs/gui/write_headerlet.py
+++ b/stwcs/gui/write_headerlet.py
@@ -1,15 +1,15 @@
-from __future__ import print_function
+from __future__ import absolute_import, division, print_function
 import os
 
 from stsci.tools import teal
 from stsci.tools import parseinput
 
-import stwcs
-from stwcs.wcsutil import headerlet
+from .. import __version__
+from ..wcsutil import headerlet
 
 __taskname__ = __name__.split('.')[-1] # needed for help string
 __package__ = headerlet.__name__
-__version__ = stwcs.__version__
+#__version__ = stwcs.__version__
 #
 #### Interfaces used by TEAL
 #
@@ -18,43 +18,44 @@ def getHelpAsString(docstring=False):
     return useful help from a file in the script directory called __taskname__.help
     """
     install_dir = os.path.dirname(__file__)
-    htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html')
-    helpfile = os.path.join(install_dir,__taskname__+'.help')
+    htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html')
+    helpfile = os.path.join(install_dir, __taskname__ + '.help')
     if docstring or (not docstring and not os.path.exists(htmlfile)):
-        helpString = __taskname__+' Version '+__version__+'\n\n'
+        helpString = __taskname__ + ' Version ' + __version__ + '\n\n'
         if os.path.exists(helpfile):
-            helpString += teal.getHelpFileAsString(__taskname__,__file__)
+            helpString += teal.getHelpFileAsString(__taskname__, __file__)
         else:
             helpString += headerlet.write_headerlet.__doc__
 
     else:
-        helpString = 'file://'+htmlfile
+        helpString = 'file://' + htmlfile
 
     return helpString
 
+
 def run(configObj=None):
-    flist,oname = parseinput.parseinput(configObj['filename'])
+    flist, oname = parseinput.parseinput(configObj['filename'])
     if len(flist) == 0:
-        print('='*60)
+        print('=' * 60)
         print('ERROR:')
         print('    No valid "filename" parameter value provided!')
         print('    Please check the working directory and restart this task.')
-        print('='*60)
+        print('=' * 60)
         return
 
-    if configObj['hdrname'] in ['',' ','INDEF']:
-        print('='*60)
+    if configObj['hdrname'] in ['', ' ', 'INDEF']:
+        print('=' * 60)
         print('ERROR:')
         print('    No valid "hdrname" parameter value provided!')
         print('    Please restart this task and provide a value for this parameter.')
-        print('='*60)
+        print('=' * 60)
         return
 
-    if configObj['output'] in ['',' ','INDEF']:
+    if configObj['output'] in ['', ' ', 'INDEF']:
         configObj['output'] = None
 
-    str_kw = ['wcsname','destim','sipname','npolfile','d2imfile',
-            'descrip','history','author','output','catalog']
+    str_kw = ['wcsname', 'destim', 'sipname', 'npolfile', 'd2imfile',
+              'descrip', 'history', 'author', 'output', 'catalog']
 
     # create dictionary of remaining parameters, deleting extraneous ones
     # such as those above