Added TDD correction to WCS

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

3 files changed, 73 insertions, 37 deletions

diff --git a/lib/utils.py b/lib/utils.py
index 87ff6d9..b29c67c 100644
--- a/lib/utils.py
+++ b/lib/utils.py
@@ -35,6 +35,11 @@ def restoreWCS(fnames):
                             nkw = ('O'+kw)[:7]
                             if backup.has_key(kw):
                                 hdr.update(kw, hdr[nkw])
+                    tddalpha = hdr.get('TDDALPHA', None)
+                    tddbeta = hdr.get('TDDBETA', None)
+                    if tddalpha or tddbeta:
+                        hdr.update('TDDALPHA', 0.0)
+                        hdr.update('TDDBETA', 0.0)
             fobj.close()
 
 def get_archive(header):
@@ -74,7 +79,7 @@ def diff_angles(a,b):
     """
     
     diff = a - b
-
+    
     if diff > 180.0:
        diff -= 360.0
 
diff --git a/updatewcs/corrections.py b/updatewcs/corrections.py
index 4699b1b..c76e918 100644
--- a/updatewcs/corrections.py
+++ b/updatewcs/corrections.py
@@ -29,19 +29,22 @@ class TDDCorr(object):
         """
         - Calculates alpha and beta for ACS/WFC data.
         - Writes 2 new kw to the extension header: TDDALPHA and TDDBETA
-        - sets TDDCORR to COMPLETE
         """
         
         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}
         
         return newkw
     updateWCS = classmethod(updateWCS)
     
-    def apply_tdd2idc(cls,ext_wcs, alpha, beta):
+    def apply_tdd2idc(cls, hwcs, alpha, beta):
         """
         Applies TDD to the idctab coefficients of a ACS/WFC observation.
         This should be always the first correction.
@@ -52,13 +55,12 @@ class TDDCorr(object):
         tdd_mat = numpy.array([[1+(beta/2048.), alpha/2048.],[alpha/2048.,1-(beta/2048.)]],numpy.float64)
         abmat1 = numpy.dot(tdd_mat, mrotn)
         abmat2 = numpy.dot(mrotp,abmat1)
-        xshape, yshape = ext_wcs.idcmodel.cx.shape, ext_wcs.idcmodel.cy.shape
-        icxy = numpy.dot(abmat2,[ext_wcs.idcmodel.cx.ravel(),ext_wcs.idcmodel.cy.ravel()])
-        ext_wcs.idcmodel.cx = icxy[0]
-        ext_wcs.idcmodel.cy = icxy[1]
-        ext_wcs.idcmodel.cx.shape = xshape
-        ext_wcs.idcmodel.cy.shape = yshape
-        #ext_wcs.wcs.set()
+        xshape, yshape = hwcs.idcmodel.cx.shape, hwcs.idcmodel.cy.shape
+        icxy = numpy.dot(abmat2,[hwcs.idcmodel.cx.ravel(), hwcs.idcmodel.cy.ravel()])
+        hwcs.idcmodel.cx = icxy[0]
+        hwcs.idcmodel.cy = icxy[1]
+        hwcs.idcmodel.cx.shape = xshape
+        hwcs.idcmodel.cy.shape = yshape
         
     apply_tdd2idc = classmethod(apply_tdd2idc)
         
@@ -93,11 +95,11 @@ class VACorr(object):
             ext_wcs.wcs.cd = ext_wcs.wcs.cd * ext_wcs.vafactor
             ext_wcs.wcs.crval[0] = ref_wcs.wcs.crval[0] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[0], ref_wcs.wcs.crval[0]) 
             ext_wcs.wcs.crval[1] = ref_wcs.wcs.crval[1] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[1], ref_wcs.wcs.crval[1]) 
+            ext_wcs.wcs.set()
             
             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]}
-            ext_wcs.wcs.set()
+                    'CRVAL1':ext_wcs.wcs.crval[0], 'CRVAL2':ext_wcs.wcs.crval[1]}            
         else:
             pass
         return kw2update
@@ -137,6 +139,8 @@ class CompSIP(object):
                     Bkey="B_%d_%d" % (m,n-m)
                     kw2update[Akey] = sipval[0,0]
                     kw2update[Bkey] = sipval[1,0]
+        #kw2update['CTYPE1'] = 'RA---TAN-SIP'
+        #kw2update['CTYPE2'] = 'DEC--TAN-SIP'
         return kw2update
     
     updateWCS = classmethod(updateWCS)
diff --git a/updatewcs/makewcs.py b/updatewcs/makewcs.py
index 1f3e640..d2a3cb6 100644
--- a/updatewcs/makewcs.py
+++ b/updatewcs/makewcs.py
@@ -4,6 +4,7 @@ import numpy
 from numpy import math 
 from math import sin, sqrt, pow, cos, asin, atan2,pi
 from hstwcs import utils
+from pytools import fileutil
 
 class MakeWCS(object):
     """
@@ -18,16 +19,18 @@ class MakeWCS(object):
     - update extension header
     
     """
- 
+    tdd_xyref = {1: [2048, 3072], 2:[2048, 1024]}
     def updateWCS(cls, ext_wcs, ref_wcs):    
         """
         recomputes the basic WCS kw 
         """
         ltvoff, offshift = cls.getOffsets(ext_wcs)
-            
-        cls.uprefwcs(ext_wcs, ref_wcs, ltvoff, offshift)
-        cls.upextwcs(ext_wcs, ref_wcs, ltvoff, offshift)
+        v23_corr = cls.zero_point_corr(ext_wcs)
+        rv23_corr = cls.zero_point_corr(ref_wcs)
+        cls.uprefwcs(ext_wcs, ref_wcs, rv23_corr, ltvoff, offshift)
+        cls.upextwcs(ext_wcs, ref_wcs, v23_corr, rv23_corr, ltvoff, offshift)
         
+
         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],
@@ -41,7 +44,7 @@ class MakeWCS(object):
     
     updateWCS = classmethod(updateWCS)
     
-    def upextwcs(self, ext_wcs, ref_wcs, loff, offsh):
+    def upextwcs(cls, ext_wcs, ref_wcs, v23_corr, rv23_corr, loff, offsh):
         """
         updates an extension wcs
         """  
@@ -58,13 +61,15 @@ class MakeWCS(object):
             fx, fy = ext_wcs.idcmodel.cx, ext_wcs.idcmodel.cy
         
         ridcmodel = ref_wcs.idcmodel
-        v2 = ext_wcs.idcmodel.refpix['V2REF']
-        v3 = ext_wcs.idcmodel.refpix['V3REF']
-        v2ref = ref_wcs.idcmodel.refpix['V2REF']
-
-        v3ref = ref_wcs.idcmodel.refpix['V3REF']
+        
+        v2 = ext_wcs.idcmodel.refpix['V2REF'] + v23_corr[0,0]
+        v3 = ext_wcs.idcmodel.refpix['V3REF'] - v23_corr[1,0]
+        v2ref = ref_wcs.idcmodel.refpix['V2REF'] + rv23_corr[0,0]
+        v3ref = ref_wcs.idcmodel.refpix['V3REF'] - rv23_corr[1,0]
+        
         R_scale = ref_wcs.idcmodel.refpix['PSCALE']/3600.0
         off = sqrt((v2-v2ref)**2 + (v3-v3ref)**2)/(R_scale*3600.0)
+        
         if v3 == v3ref:
            theta=0.0
         else:
@@ -74,13 +79,14 @@ class MakeWCS(object):
 
         dX=(off*sin(theta)) + offshiftx
         dY=(off*cos(theta)) + offshifty
+        
         px = numpy.array([[dX,dY]])
-        newcrval = ref_wcs.all_pix2sky_fits(px)[0]
+        newcrval = ref_wcs.wcs.p2s_fits(px)['world'][0]
         newcrpix = numpy.array([ext_wcs.idcmodel.refpix['XREF'] + ltvoffx, 
                                 ext_wcs.idcmodel.refpix['YREF'] + ltvoffy])
         ext_wcs.wcs.crval = newcrval
         ext_wcs.wcs.crpix = newcrpix
-        
+        ext_wcs.wcs.set()
         # Account for subarray offset
         # Angle of chip relative to chip
         if ext_wcs.idcmodel.refpix['THETA']:
@@ -96,7 +102,7 @@ class MakeWCS(object):
         delYX=fy[1,1]/R_scale/3600.
         delXY=fx[1,0]/R_scale/3600.
         delYY=fy[1,0]/R_scale/3600.
-
+        
         # Convert to radians
         rr=dtheta*pi/180.0
         
@@ -106,54 +112,59 @@ class MakeWCS(object):
 
         dXY = cos(rr)*delXY - sin(rr)*delYY
         dYY = sin(rr)*delXY + cos(rr)*delYY
+        
         px = numpy.array([[dX+dXX,dY+dYX]])
         
         # Transform to sky coordinates
-        wc = ref_wcs.all_pix2sky_fits(px)
+        wc = ref_wcs.wcs.p2s_fits(px)['world']
         
         a = wc[0,0]
         b = wc[0,1]
         px = numpy.array([[dX+dXY,dY+dYY]])
         
-        wc = ref_wcs.all_pix2sky_fits(px)
+        wc = ref_wcs.wcs.p2s_fits(px)['world']
+        
         c = wc[0,0]
         d = wc[0,1]
-
+        
         # Calculate the new CDs and convert to degrees
         cd11 = utils.diff_angles(a,newcrval[0])*cos(newcrval[1]*pi/180.0)
         cd12 = utils.diff_angles(c,newcrval[0])*cos(newcrval[1]*pi/180.0)
         cd21 = utils.diff_angles(b,newcrval[1])
         cd22 = utils.diff_angles(d,newcrval[1])
         cd = numpy.array([[cd11, cd12], [cd21, cd22]])
-        ext_wcs.wcs.cd = cd  #???
-        
+        ext_wcs.wcs.cd = cd  
+        ext_wcs.wcs.set()
     upextwcs = classmethod(upextwcs)
         
-    def uprefwcs(self, ext_wcs, ref_wcs, loff, offsh):
+    def uprefwcs(cls, ext_wcs, ref_wcs, rv23_corr_tdd, loff, offsh):
         """
         Updates the reference chip
         """
         ltvoffx, ltvoffy = loff[0], loff[1]
         offshift = offsh
         dec = ref_wcs.wcs.crval[1]
+        
+        rv23 = [ref_wcs.idcmodel.refpix['V2REF'] + rv23_corr_tdd[0,0], 
+            ref_wcs.idcmodel.refpix['V3REF'] - rv23_corr_tdd[1,0]]
         # Get an approximate reference position on the sky
-        rref = numpy.array([[ref_wcs.idcmodel.refpix['XREF']+ltvoffx, 
+        rref = numpy.array([[ref_wcs.idcmodel.refpix['XREF']+ltvoffx , 
                             ref_wcs.idcmodel.refpix['YREF']+ltvoffy]])
         
-        crval = ref_wcs.all_pix2sky_fits(rref)
+        crval = ref_wcs.wcs.p2s_fits(rref)['world'][0]
         # Convert the PA_V3 orientation to the orientation at the aperture
         # This is for the reference chip only - we use this for the
         # reference tangent plane definition
         # It has the same orientation as the reference chip
-        pv = troll(ext_wcs.pav3,dec,ref_wcs.idcmodel.refpix['V2REF'],
-                    ref_wcs.idcmodel.refpix['V3REF'])
+        pv = troll(ext_wcs.pav3,dec,rv23[0], rv23[1])
+        
         # Add the chip rotation angle
         if ref_wcs.idcmodel.refpix['THETA']:
             pv += ref_wcs.idcmodel.refpix['THETA']
             
         
         # Set values for the rest of the reference WCS
-        ref_wcs.wcs.crval = crval[0]
+        ref_wcs.wcs.crval = crval
         ref_wcs.wcs.crpix = numpy.array([0.0,0.0])+offsh
         parity = ref_wcs.parity
         R_scale = ref_wcs.idcmodel.refpix['PSCALE']/3600.0
@@ -168,7 +179,23 @@ 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:
+            return numpy.array([[0., 0.],[0.,0.]])
+        
+        tdd = numpy.array([[beta, alpha], [alpha, -beta]])
+        mrotp = fileutil.buildRotMatrix(2.234529)/2048.
+        xy0 = numpy.array([[cls.tdd_xyref[hwcs.chip][0]-2048.], [cls.tdd_xyref[hwcs.chip][1]-2048.]])
+        
+        v23_corr = numpy.dot(mrotp,numpy.dot(tdd,xy0)) * 0.05
+        
+        return v23_corr 
+    
+    zero_point_corr = classmethod(zero_point_corr)
+        
     def getOffsets(cls, ext_wcs):
         ltv1 = ext_wcs.ltv1
         ltv2 = ext_wcs.ltv2