Addresses ticket #563

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

2 files changed, 92 insertions, 10 deletions

diff --git a/distortion/mutil.py b/distortion/mutil.py
index 44b46cd..780de0b 100644
--- a/distortion/mutil.py
+++ b/distortion/mutil.py
@@ -74,7 +74,23 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
             detector = str(ftab['PRIMARY'].header['CAMERA'])
         else:
             detector = 1
-
+    # First, read in TDD coeffs if present
+    phdr = ftab['PRIMARY'].header
+    skew_coeffs = read_tdd_coeffs(phdr)
+    """
+    # Read in ACS TDD correction coefficients
+    if ftab['PRIMARY'].header.has_key('INSTRUME'):
+        instrument = ftab['PRIMARY'].header['INSTRUME']
+    if instrument == 'ACS' and detector == 'WFC':
+        skew_coeffs = None    
+        if phdr.has_key('TDD_DATE'):
+            print 'Reading TDD coefficients from ',tabname
+            skew_coeffs = read_tdd_coeffs(phdr)
+            # Using coefficients read in from IDCTAB Primary header
+            #skew_coeffs = {'TDD_A0':phdr['TDD_A0'],'TDD_A1':phdr["TDD_A1"],
+            #            'TDD_B0':phdr['TDD_B0'],'TDD_B1':phdr['TDD_B1'],
+            #            'TDD_D0':phdr['TDD_D0'],'TDD_DATE':phdr['TDD_DATE']}
+    """
     # Set default filters for SBC
     if detector == 'SBC':
         if filter1 == 'CLEAR':
@@ -195,6 +211,7 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
     refpix['YDELTA'] = 0.0
     refpix['DEFAULT_SCALE'] = yes
     refpix['centered'] = no
+    refpix['skew_coeffs'] = skew_coeffs
 
     # Now that we know which row to look at, read coefficients into the
     #   numeric arrays we have set up...
@@ -228,6 +245,41 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
     # 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)
+#
+#
+def read_tdd_coeffs(phdr):
+    ''' Read in the TDD related keywords from the PRIMARY header of the IDCTAB
+    '''
+    skew_coeffs = {}
+    skew_coeffs['TDDORDER'] = 0
+    skew_coeffs['TDD_DATE'] = ""
+    skew_coeffs['TDD_A'] = None
+    skew_coeffs['TDD_B'] = None
+
+    if phdr.has_key("TDDORDER"):
+        n = int(phdr["TDDORDER"])
+    else:
+        print 'TDDORDER kw not present, TDD correction will not be applied.'
+        return skew_coeffs
+    
+    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():
+        print 'Warning: TDD_A and TDD_B coeffiecients have values of 0, \n \
+               but TDDORDER is %d.' % TDDORDER
+
+    skew_coeffs['TDDORDER'] = n
+    skew_coeffs['TDD_DATE'] = phdr['TDD_DATE']
+    skew_coeffs['TDD_A'] = a
+    skew_coeffs['TDD_B'] = b
+    
+    return skew_coeffs
 
 def readOfftab(offtab, date, chip=None):
 
diff --git a/updatewcs/corrections.py b/updatewcs/corrections.py
index f296c10..eac4b6e 100644
--- a/updatewcs/corrections.py
+++ b/updatewcs/corrections.py
@@ -70,18 +70,48 @@ class TDDCorr(object):
         
     def compute_alpha_beta(cls, ext_wcs):
         """
-        Compute the time dependent distortion skew terms
-        default date of 2004.5 = 2004-7-1
-        """
-        dday = 2004.5
-        year,month,day = ext_wcs.date_obs.split('-')
-        rdate = datetime.datetime(int(year),int(month),int(day))
-        rday = float(rdate.strftime("%j"))/365.25 + rdate.year
+        Compute the ACS time dependent distortion skew terms
+        as described in ACS ISR 07-08 by J. Anderson.
+        
+        Jay's code only computes the alpha/beta values based on a decimal year
+        with only 3 digits, so this line reproduces that when needed for comparison
+        with his results.
+        rday = float(('%0.3f')%rday)
+        
+        The zero-point terms account for the skew accumulated between
+        2002.0 and 2004.5, when the latest IDCTAB was delivered.
         alpha = 0.095 + 0.090*(rday-dday)/2.5
         beta = -0.029 - 0.030*(rday-dday)/2.5
+        """
+        skew_coeffs = ext_wcs.idcmodel.refpix['skew_coeffs']
+        if skew_coeffs is None:
+            err_str =  "------------------------------------------------------------------------  \n"
+            err_str += "WARNING: the IDCTAB geometric distortion file specified in the image      \n"
+            err_str += "         header did not have the time-dependent distortion coefficients.  \n"
+            err_str += "         The pre-SM4 time-dependent skew solution will be used by default.\n"
+            err_str += "         Please update IDCTAB with new reference file from HST archive.   \n"
+            err_str +=  "------------------------------------------------------------------------  \n"
+            print err_str
+            # Using default pre-SM4 coefficients
+            skew_coeffs = {'TDD_A':[0.095,0.090/2.5],
+                        'TDD_B':[-0.029,-0.030/2.5],
+                        'TDD_DATE':2004.5}
         
-        return alpha, beta
-
+        if not isinstance(ext_wcs.date_obs,float):
+            year,month,day = ext_wcs.date_obs.split('-')
+            rdate = datetime.datetime(int(year),int(month),int(day))
+            rday = float(rdate.strftime("%j"))/365.25 + rdate.year
+        else:
+            rday = ext_wcs.date_obs
+        
+        alpha = 0
+        beta = 0
+        # Compute skew terms, allowing for non-linear coefficients as well
+        for c in range(skew_coeffs['TDDORDER']+1):
+            alpha += skew_coeffs['TDD_A'][c]* np.power((rday-skew_coeffs['TDD_DATE']),c)
+            beta += skew_coeffs['TDD_B'][c]*np.power((rday-skew_coeffs['TDD_DATE']),c)
+            
+        return alpha,beta
     compute_alpha_beta = classmethod(compute_alpha_beta)