Index: trunk/psModules/src/astrom/pmAstrometryWCS.c
===================================================================
--- trunk/psModules/src/astrom/pmAstrometryWCS.c	(revision 15254)
+++ trunk/psModules/src/astrom/pmAstrometryWCS.c	(revision 15896)
@@ -7,6 +7,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-10-09 19:27:04 $
+ *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-12-22 17:52:57 $
  *
  *  Copyright 2006 Institute for Astronomy, University of Hawaii
@@ -56,14 +56,7 @@
 bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip, double tol)
 {
-    // XXX require chip->toFPA->x->nX == chip->toFPA->x->nY
-    // XXX require chip->toFPA->y->nX == chip->toFPA->y->nY
-    // XXX require chip->toFPA->x->nX == chip->toFPA->y->nX
-    // XXX require chip->toFPA->nX == 1,2,3
-
-    // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
-    if (!psPlaneTransformIsDiagonal (fpa->toTPA))
-        psAbort("invalid TPA transformation");
-
     pmAstromWCS *wcs = pmAstromWCSfromFPA(fpa, chip, tol);
+    if (!wcs) return false;
+
     pmAstromWCStoHeader (header, wcs);
 
@@ -118,5 +111,4 @@
 {
     pmAstromWCS *wcs = pmAstromWCSBilevelChipFromFPA (chip, tol);
-
     pmAstromWCStoHeader (header, wcs);
 
@@ -337,4 +329,6 @@
     char *type;
 
+    if (!wcs) return false;
+
     type = psProjectTypeToString (wcs->toSky->type, "RA--");
     psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", type);
@@ -385,5 +379,5 @@
     }
 
-    return (true);
+    return true;
 }
 
@@ -550,10 +544,4 @@
 pmAstromWCS *pmAstromWCSfromFPA (const pmFPA *fpa, const pmChip *chip, double tol)
 {
-
-    // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
-    // XXX not really: toTPA needs to have unity scale for distortion fitting function
-    if (!psPlaneTransformIsDiagonal (fpa->toTPA))
-        psAbort("invalid TPA transformation");
-
     // XXX require chip->toFPA->x->nX == chip->toFPA->x->nY
     // XXX require chip->toFPA->y->nX == chip->toFPA->y->nY
@@ -561,5 +549,36 @@
     // XXX require chip->toFPA->nX == 1,2,3
 
-    pmAstromWCS *wcs = pmAstromWCSAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY);
+    // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
+    // XXX not really: toTPA needs to have unity scale for distortion fitting function
+    // if (!psPlaneTransformIsDiagonal (fpa->toTPA))
+    // psAbort("invalid TPA transformation");
+
+    // create a temporary transform which combines toTPA and toFPA.  allow toTPA to have 0th
+    // and 1st order terms
+    
+    // XXX require fpa->toTPA->x->nX == fpa->toTPA->x->nY
+    // XXX require fpa->toTPA->y->nX == fpa->toTPA->y->nY
+    // XXX require fpa->toTPA->x->nX == fpa->toTPA->y->nX
+    // XXX require fpa->toTPA->x->coeffMask[1][1]
+    // XXX require fpa->toTPA->y->coeffMask[1][1]
+    // XXX require fpa->toTPA->nX == 1
+
+    psPlaneTransform *toTPA = psPlaneTransformAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY);
+
+    for (int i = 0; i <= toTPA->x->nX; i++) {
+	for (int j = 0; j <= toTPA->x->nY; j++) {
+	    double f1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[1][0]*chip->toFPA->x->coeff[i][j];
+	    double f2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[0][1]*chip->toFPA->y->coeff[i][j];
+	    toTPA->x->coeff[i][j] = f1 + f2;
+
+	    double g1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[1][0]*chip->toFPA->x->coeff[i][j];
+	    double g2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[0][1]*chip->toFPA->y->coeff[i][j];
+	    toTPA->y->coeff[i][j] = g1 + g2;
+	}
+    }
+    toTPA->x->coeff[0][0] += fpa->toTPA->x->coeff[0][0];
+    toTPA->y->coeff[0][0] += fpa->toTPA->y->coeff[0][0];
+
+    pmAstromWCS *wcs = pmAstromWCSAlloc(toTPA->x->nX, toTPA->x->nY);
 
     // convert projection from FPA to SKY into wcs projection (degrees to radians)
@@ -569,9 +588,9 @@
 
     // given transformation, solve for coordinates which yields output coordinates of 0,0
-    psPlane *center = psPlaneTransformGetCenter (chip->toFPA, tol);
+    psPlane *center = psPlaneTransformGetCenter (toTPA, tol);
 
     // create wcs transform from toFPA, resulting transformation has units of microns/pixel
     // adjust wcs transform to use center as reference coordinate
-    psPlaneTransformSetCenter (wcs->trans, chip->toFPA, center->x, center->y);
+    psPlaneTransformSetCenter (wcs->trans, toTPA, center->x, center->y);
 
     // calculated center is crpix1,2
@@ -595,4 +614,6 @@
     wcs->cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
     wcs->cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
+
+    psFree (toTPA);
 
     return wcs;