Initial check-in to add support for new ACS TDD correction described in Ticket #1108

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

5 files changed, 138 insertions, 99 deletions

diff --git a/lib/stwcs/distortion/mutil.py b/lib/stwcs/distortion/mutil.py
index 02fc4c2..ea131f9 100644
--- a/lib/stwcs/distortion/mutil.py
+++ b/lib/stwcs/distortion/mutil.py
@@ -78,7 +78,7 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
     phdr = ftab['PRIMARY'].header
     instrument = phdr['INSTRUME']
     if instrument == 'ACS' and detector == 'WFC':
-        skew_coeffs = read_tdd_coeffs(phdr)
+        skew_coeffs = read_tdd_coeffs(phdr, chip=chip)
     else:
         skew_coeffs = None
 
@@ -216,7 +216,6 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
     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...
     # Setup which column name convention the IDCTAB follows
@@ -254,7 +253,7 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward',
 # Time-dependent skew correction coefficients (only ACS/WFC)
 #
 #
-def read_tdd_coeffs(phdr):
+def read_tdd_coeffs(phdr, chip=1):
     ''' Read in the TDD related keywords from the PRIMARY header of the IDCTAB
     '''
     skew_coeffs = {}
@@ -262,26 +261,34 @@ def read_tdd_coeffs(phdr):
     skew_coeffs['TDD_DATE'] = ""
     skew_coeffs['TDD_A'] = None
     skew_coeffs['TDD_B'] = None
-
-    if "TDDORDER" in phdr:
-        n = int(phdr["TDDORDER"])
+    skew_coeffs['TDD_CY_BETA'] = None
+
+    if "TDD_CYB1" in phdr:
+        # We have 2014-calibrated TDD correction to apply, not J.A.-derived values
+        print "Using 2014-calibrated TDD correction..."
+        skew_coeffs['TDD_DATE'] = phdr['TDD_DATE']
+        cyb_kw = 'TDD_CYB{0}'.format(int(chip))
+        skew_coeffs['TDD_CY_BETA'] = phdr[cyb_kw]
     else:
-        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[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
+        if "TDDORDER" in phdr:
+            n = int(phdr["TDDORDER"])
+        else:
+            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[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
 
diff --git a/lib/stwcs/updatewcs/__init__.py b/lib/stwcs/updatewcs/__init__.py
index a0e39f6..8778065 100644
--- a/lib/stwcs/updatewcs/__init__.py
+++ b/lib/stwcs/updatewcs/__init__.py
@@ -111,6 +111,7 @@ def makecorr(fname, allowed_corr):
     `acorr`: list
              list of corrections to be applied
     """
+    logger.info("Allowed corrections: {0}".format(allowed_corr))
     f = pyfits.open(fname, mode='update')
     #Determine the reference chip and create the reference HSTWCS object
     nrefchip, nrefext = getNrefchip(f)
diff --git a/lib/stwcs/updatewcs/apply_corrections.py b/lib/stwcs/updatewcs/apply_corrections.py
index 3640c56..740d424 100644
--- a/lib/stwcs/updatewcs/apply_corrections.py
+++ b/lib/stwcs/updatewcs/apply_corrections.py
@@ -110,7 +110,6 @@ def applyTDDCorr(fname, utddcorr):
             #print "***IDCTAB file not found - not applying TDD correction***\n"
     else:
         tddcorr = False
-
     return tddcorr
 
 def applyNpolCorr(fname, unpolcorr):
@@ -224,4 +223,3 @@ def applyD2ImCorr(fname, d2imcorr):
         print 'D2IMFILE keyword not found in primary header'
         applyD2IMCorr = False
         return applyD2IMCorr
-
diff --git a/lib/stwcs/updatewcs/corrections.py b/lib/stwcs/updatewcs/corrections.py
index 70bca3b..6dc2ecd 100644
--- a/lib/stwcs/updatewcs/corrections.py
+++ b/lib/stwcs/updatewcs/corrections.py
@@ -17,46 +17,46 @@ class TDDCorr(object):
     """
     Apply time dependent distortion correction to distortion coefficients and basic
     WCS keywords. This correction **must** be done before any other WCS correction.
-        
+
     Parameters
     ----------
     ext_wcs: HSTWCS object
              An HSTWCS object to be modified
     ref_wcs: HSTWCS object
              A reference HSTWCS object
-             
+
     Notes
     -----
-    Compute the ACS/WFC time dependent distortion terms as described 
+    Compute the ACS/WFC time dependent distortion terms as described
     in [1]_ and apply the correction to the WCS of the observation.
-    
+
     The model coefficients are stored in the primary header of the IDCTAB.
-    :math:`D_{ref}` is the reference date. The computed corrections are saved 
+    :math:`D_{ref}` is the reference date. The computed corrections are saved
     in the science extension header as TDDALPHA and TDDBETA keywords.
-    
+
     .. math:: TDDALPHA = A_{0} + {A_{1}*(obsdate - D_{ref})}
-    
+
     .. math:: TDDBETA =  B_{0} + B_{1}*(obsdate - D_{ref})
-    
-       
+
+
     The time dependent distortion affects the IDCTAB coefficients, and the
-    relative location of the two chips. Because the linear order IDCTAB 
+    relative location of the two chips. Because the linear order IDCTAB
     coefficients ar eused in the computatuion of the NPOL extensions,
     the TDD correction affects all components of the distortion model.
-    
+
     Application of TDD to the IDCTAB polynomial coefficients:
-    The TDD model is computed in Jay's frame, while the IDCTAB coefficients 
-    are in the HST V2/V3 frame. The coefficients are transformed to Jay's frame, 
-    TDD is applied and they are transformed back to the V2/V3 frame. This 
-    correction is performed in this class. 
-    
-    Application of TDD to the relative location of the two chips is 
+    The TDD model is computed in Jay's frame, while the IDCTAB coefficients
+    are in the HST V2/V3 frame. The coefficients are transformed to Jay's frame,
+    TDD is applied and they are transformed back to the V2/V3 frame. This
+    correction is performed in this class.
+
+    Application of TDD to the relative location of the two chips is
     done in `makewcs`.
-    
+
     References
     ----------
     .. [1] Jay Anderson, "Variation of the Distortion Solution of the WFC", ACS ISR 2007-08.
-    
+
     """
     def updateWCS(cls, ext_wcs, ref_wcs):
         """
@@ -64,21 +64,51 @@ class TDDCorr(object):
         - Writes 2 new kw to the extension header: TDDALPHA and TDDBETA
         """
         logger.info("\n\tStarting TDDCorr: %s" % time.asctime())
-        alpha, beta = cls.compute_alpha_beta(ext_wcs)
-        cls.apply_tdd2idc(ref_wcs, alpha, beta)
-        cls.apply_tdd2idc(ext_wcs, alpha, beta)
-        ext_wcs.idcmodel.refpix['TDDALPHA'] = alpha
-        ext_wcs.idcmodel.refpix['TDDBETA'] = beta
-        ref_wcs.idcmodel.refpix['TDDALPHA'] = alpha
-        ref_wcs.idcmodel.refpix['TDDBETA'] = beta
-        
-        newkw = {'TDDALPHA': alpha, 'TDDBETA':beta, 'OCX10':ext_wcs.idcmodel.cx[1,0],
-                'OCX11':ext_wcs.idcmodel.cx[1,1],'OCY10':ext_wcs.idcmodel.cy[1,0],
-                'OCY11':ext_wcs.idcmodel.cy[1,1],}
-        
+
+        if ext_wcs.idcmodel.refpix['skew_coeffs']is not None and \
+            ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CY_BETA'] is not None:
+            logger.info("\n\t Applying 2014-calibrated TDD: %s" % time.asctime())
+            # We have 2014-calibrated TDD, not J.A.-style TDD
+            cls.apply_tdd2idc2(ref_wcs)
+            cls.apply_tdd2idc2(ext_wcs)
+            newkw = {'TDDALPHA': None, 'TDDBETA':None, 'OCX10':ext_wcs.idcmodel.cx[1,0],
+                    'OCX11':ext_wcs.idcmodel.cx[1,1],'OCY10':ext_wcs.idcmodel.cy[1,0],
+                    'OCY11':ext_wcs.idcmodel.cy[1,1],}
+
+        else:
+            alpha, beta = cls.compute_alpha_beta(ext_wcs)
+            cls.apply_tdd2idc(ref_wcs, alpha, beta)
+            cls.apply_tdd2idc(ext_wcs, alpha, beta)
+            ext_wcs.idcmodel.refpix['TDDALPHA'] = alpha
+            ext_wcs.idcmodel.refpix['TDDBETA'] = beta
+            ref_wcs.idcmodel.refpix['TDDALPHA'] = alpha
+            ref_wcs.idcmodel.refpix['TDDBETA'] = beta
+
+            newkw = {'TDDALPHA': alpha, 'TDDBETA':beta, 'OCX10':ext_wcs.idcmodel.cx[1,0],
+                    'OCX11':ext_wcs.idcmodel.cx[1,1],'OCY10':ext_wcs.idcmodel.cy[1,0],
+                    'OCY11':ext_wcs.idcmodel.cy[1,1],}
+
         return newkw
     updateWCS = classmethod(updateWCS)
-    
+
+    def apply_tdd2idc2(cls, hwcs):
+        """ Applies 2014-calibrated TDD correction to single IDCTAB coefficient
+            of an ACS/WFC observation.
+        """
+        if not isinstance(hwcs.date_obs,float):
+            year,month,day = hwcs.date_obs.split('-')
+            rdate = datetime.datetime(int(year),int(month),int(day))
+            rday = float(rdate.strftime("%j"))/365.25 + rdate.year
+        else:
+            rday = hwcs.date_obs
+
+        skew_coeffs = hwcs.idcmodel.refpix['skew_coeffs']
+        cy_beta = skew_coeffs['TDD_CY_BETA']
+
+        hwcs.idcmodel.cy[1,1] += cy_beta*(rday - skew_coeffs['TDD_DATE'])
+
+    apply_tdd2idc2 = classmethod(apply_tdd2idc2)
+
     def apply_tdd2idc(cls, hwcs, alpha, beta):
         """
         Applies TDD to the idctab coefficients of a ACS/WFC observation.
@@ -96,19 +126,19 @@ class TDDCorr(object):
         hwcs.idcmodel.cy = icxy[1]
         hwcs.idcmodel.cx.shape = xshape
         hwcs.idcmodel.cy.shape = yshape
-        
+
     apply_tdd2idc = classmethod(apply_tdd2idc)
-        
+
     def compute_alpha_beta(cls, ext_wcs):
         """
         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
@@ -136,64 +166,64 @@ class TDDCorr(object):
             skew_coeffs = {'TDD_A':[0.095,0.090/2.5],
                         'TDD_B':[-0.029,-0.030/2.5],
                         'TDD_DATE':2004.5,'TDDORDER':1}
-                    
+
         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)
-    
-        
+
+
 class VACorr(object):
     """
     Apply velocity aberation correction to WCS keywords.
-    
+
     Notes
     -----
     Velocity Aberration is stored in the extension header keyword 'VAFACTOR'.
     The correction is applied to the CD matrix and CRVALs.
-    
+
     """
     def updateWCS(cls, ext_wcs, ref_wcs):
         logger.info("\n\tStarting VACorr: %s" % time.asctime())
         if ext_wcs.vafactor != 1:
             ext_wcs.wcs.cd = ext_wcs.wcs.cd * ext_wcs.vafactor
-            crval0 = ref_wcs.wcs.crval[0] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[0], 
-                                                                         ref_wcs.wcs.crval[0]) 
-            crval1 = ref_wcs.wcs.crval[1] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[1], 
-                                                                         ref_wcs.wcs.crval[1]) 
+            crval0 = ref_wcs.wcs.crval[0] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[0],
+                                                                         ref_wcs.wcs.crval[0])
+            crval1 = ref_wcs.wcs.crval[1] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[1],
+                                                                         ref_wcs.wcs.crval[1])
             crval = np.array([crval0,crval1])
-            ext_wcs.wcs.crval = crval            
+            ext_wcs.wcs.crval = crval
             ext_wcs.wcs.set()
         else:
             pass
-        kw2update={'CD1_1': ext_wcs.wcs.cd[0,0], 'CD1_2':ext_wcs.wcs.cd[0,1], 
-                    'CD2_1':ext_wcs.wcs.cd[1,0], 'CD2_2':ext_wcs.wcs.cd[1,1], 
-                    'CRVAL1':ext_wcs.wcs.crval[0], 'CRVAL2':ext_wcs.wcs.crval[1]}    
+        kw2update={'CD1_1': ext_wcs.wcs.cd[0,0], 'CD1_2':ext_wcs.wcs.cd[0,1],
+                    'CD2_1':ext_wcs.wcs.cd[1,0], 'CD2_2':ext_wcs.wcs.cd[1,1],
+                    'CRVAL1':ext_wcs.wcs.crval[0], 'CRVAL2':ext_wcs.wcs.crval[1]}
         return kw2update
-    
+
     updateWCS = classmethod(updateWCS)
 
-        
+
 class CompSIP(object):
     """
     Compute Simple Imaging Polynomial (SIP) coefficients as defined in [2]_
     from IDC table coefficients.
-    
+
     This class transforms the TDD corrected IDCTAB coefficients into SIP format.
-    It also applies a binning factor to the coefficients if the observation was 
+    It also applies a binning factor to the coefficients if the observation was
     binned.
-    
+
     References
     ----------
     .. [2] David Shupe, et al, "The SIP Convention of representing Distortion
        in FITS Image headers",  Astronomical Data Analysis Software And Systems, ASP
        Conference Series, Vol. 347, 2005
-    
+
     """
     def updateWCS(cls, ext_wcs, ref_wcs):
         logger.info("\n\tStarting CompSIP: %s" %time.asctime())
@@ -205,10 +235,10 @@ class CompSIP(object):
         kw2update['A_ORDER'] = order
         kw2update['B_ORDER'] = order
         pscale = ext_wcs.idcmodel.refpix['PSCALE']
-        
+
         cx = ext_wcs.idcmodel.cx
         cy = ext_wcs.idcmodel.cy
-        
+
         matr = np.array([[cx[1,1],cx[1,0]], [cy[1,1],cy[1,0]]], dtype=np.float64)
         imatr = linalg.inv(matr)
         akeys1 = np.zeros((order+1,order+1), dtype=np.float64)
@@ -218,8 +248,8 @@ class CompSIP(object):
                 if n >= m and n>=2:
                     idcval = np.array([[cx[n,m]],[cy[n,m]]])
                     sipval = np.dot(imatr, idcval)
-                    akeys1[m,n-m] = sipval[0] 
-                    bkeys1[m,n-m] = sipval[1] 
+                    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] * ext_wcs.binned
@@ -227,7 +257,5 @@ class CompSIP(object):
         kw2update['CTYPE1'] = 'RA---TAN-SIP'
         kw2update['CTYPE2'] = 'DEC--TAN-SIP'
         return kw2update
-    
-    updateWCS = classmethod(updateWCS)
-
 
+    updateWCS = classmethod(updateWCS)
diff --git a/lib/stwcs/updatewcs/makewcs.py b/lib/stwcs/updatewcs/makewcs.py
index 9f95a13..68ceded 100644
--- a/lib/stwcs/updatewcs/makewcs.py
+++ b/lib/stwcs/updatewcs/makewcs.py
@@ -1,5 +1,7 @@
 from __future__ import division # confidence high
 
+import datetime
+
 import numpy as np
 from math import sin, sqrt, pow, cos, asin, atan2,pi
 import utils
@@ -80,9 +82,9 @@ class MakeWCS(object):
             offsetx = ext_wcs.wcs.crpix[0] - ltv1 - ext_wcs.idcmodel.refpix['XREF']
             offsety = ext_wcs.wcs.crpix[1] - ltv2 - ext_wcs.idcmodel.refpix['YREF']
             ext_wcs.idcmodel.shift(offsetx,offsety)
-        
+
         fx, fy = ext_wcs.idcmodel.cx, ext_wcs.idcmodel.cy
-        
+
         ext_wcs.setPscale()
         tddscale = (ref_wcs.pscale/fx[1,1])
         v2 = ext_wcs.idcmodel.refpix['V2REF'] + v23_corr[0,0] * tddscale
@@ -149,7 +151,7 @@ class MakeWCS(object):
             offsetx = ref_wcs.wcs.crpix[0] - ltv1 - ref_wcs.idcmodel.refpix['XREF']
             offsety = ref_wcs.wcs.crpix[1] - ltv2 - ref_wcs.idcmodel.refpix['YREF']
             ref_wcs.idcmodel.shift(offsetx,offsety)
-        
+
         rfx, rfy = ref_wcs.idcmodel.cx, ref_wcs.idcmodel.cy
 
         offshift = offsh
@@ -189,16 +191,20 @@ class MakeWCS(object):
     uprefwcs = classmethod(uprefwcs)
 
     def zero_point_corr(cls,hwcs):
-        try:
-            alpha = hwcs.idcmodel.refpix['TDDALPHA']
-            beta = hwcs.idcmodel.refpix['TDDBETA']
-        except KeyError:
-            alpha = 0.0
-            beta = 0.0
+        if hwcs.idcmodel.refpix['skew_coeffs'] is not None and 'TDD_CY_BETA' in hwcs.idcmodel.refpix['skew_coeffs']:
             v23_corr = np.array([[0.],[0.]])
-            logger.debug("\n\tTDD Zero point correction for chip %s defaulted to: %s" % (hwcs.chip, v23_corr))
             return v23_corr
- 
+        else:
+            try:
+                alpha = hwcs.idcmodel.refpix['TDDALPHA']
+                beta = hwcs.idcmodel.refpix['TDDBETA']
+            except KeyError:
+                alpha = 0.0
+                beta = 0.0
+                v23_corr = np.array([[0.],[0.]])
+                logger.debug("\n\tTDD Zero point correction for chip %s defaulted to: %s" % (hwcs.chip, v23_corr))
+                return v23_corr
+
         tdd = np.array([[beta, alpha], [alpha, -beta]])
         mrotp = fileutil.buildRotMatrix(2.234529)/2048.
         xy0 = np.array([[cls.tdd_xyref[hwcs.chip][0]-2048.], [cls.tdd_xyref[hwcs.chip][1]-2048.]])
@@ -264,4 +270,3 @@ def troll(roll, dec, v2, v3):
     troll = np.rad2deg(pi - (gamma+B))
 
     return troll
-