Index: trunk/psastro/src/psastroWCS.c
===================================================================
--- trunk/psastro/src/psastroWCS.c	(revision 7014)
+++ trunk/psastro/src/psastroWCS.c	(revision 7332)
@@ -3,9 +3,6 @@
 // interpret header WCS (only handles traditional WCS for the moment)
 // plateScale is nominal physical scale on tangent plane (microns / arcsecond)
-bool pmAstromReadWCS (psProjection **toSkyOut, psPlaneDistort **toTPAout, psPlaneTransform **toFPAout, psMetadata *header, double plateScale) { 
-
-    psPlaneTransform *toFPA;
-    psPlaneDistort *toTPA;
-    psProjection *toSky;
+bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic) { 
+
     psProjectionType type;
     bool status;
@@ -13,5 +10,5 @@
     float pc1_1, pc1_2, pc2_1, pc2_2;
 
-    // *** interpret header data, convert to crval(i), etc
+    // interpret header data, convert to crval(i), etc
     char *ctype = psMetadataLookupPtr (&status, header, "CTYPE2");
     if (!status) {
@@ -88,51 +85,109 @@
 got_matrix:
 
-    // XXX EAM : if fpa->toSky and fpa->toTPA are already defined, then the
-    //           toFPA must be modified to match the crval(i), scale(i) and crpix(i)
-    //           scale = scale(i)/scale(0) (i == chip #)
-    //           project crval1(0),crval2(0 using projection
-    toFPA = psPlaneTransformAlloc (1, 1);
+    /*****
+
+    For mosaic astrometry, we need to have a starting set of projection terms in which the
+    chip-to-FPA terms result in a fixed physical unit on the focal plane (eg, pixels or
+    microns).  This set of projections, coupled with an identity toTPA (ie, no distortion) will
+    result in substantial errors between the observed and predicted star positions on the focal
+    plane: this is the measurement of the optical distortion in the camera.  At the same time,
+    we need to carry around the transformations which allow us to make an accurate calculation
+    of the position of the stars based on the input (per-chip) astrometry.  These
+    transformations will allow us to match the raw and ref stars robustly.  To convert the
+    per-chip astrometry (which may have been calculated with a different plate scale for each
+    chip) to a collection of astrometry terms for chips in a single mosaic, we need to adjust
+    the chip-to-FPA scaling (eg, pc11) to match the variations in the effective plate scale for
+    each chip (eg, cdelt1).  Thus, we need to carry around both the
+
+    *****/
+
+    {
+	double rX = 1.0;
+	double rY = 1.0;
+
+	// XXX free an existing toFPA?
+	psPlaneTransform *toFPA = psPlaneTransformAlloc (1, 1);
     
-    double cdelt = hypot (cdelt1, cdelt2) / plateScale;  // degrees / micron (eg, in fact, whatever unit user chooses for focal plane)
-    cdelt1 /= cdelt;
-    cdelt2 /= cdelt;
-
-    toFPA->x->coeff[0][0] = -(pc1_1*cdelt1*crpix1 + pc1_2*cdelt2*crpix2);
-    toFPA->x->coeff[1][0] = +(pc1_1*cdelt1);
-    toFPA->x->coeff[0][1] = +(pc1_2*cdelt2);
-    toFPA->x->mask[1][1]  = 1;
-
-    toFPA->y->coeff[0][0] = -(pc2_1*cdelt1*crpix1 + pc2_2*cdelt2*crpix2);
-    toFPA->y->coeff[1][0] = +(pc2_1*cdelt1);
-    toFPA->y->coeff[0][1] = +(pc2_2*cdelt2);
-    toFPA->y->mask[1][1]  = 1;
-    *toFPAout = toFPA;
-
-    if (*toSkyOut != NULL) {
-	if (*toTPAout == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined");
-
-	psSphere sky;
-	psPlane tpa;
-
-	sky.r = crval1*RAD_DEG;
-	sky.d = crval2*RAD_DEG;
-	p_psProject (&tpa, &sky, *toSkyOut);
+	// basic transformation from chip to FPA
+	toFPA->x->coeff[0][0] = -(pc1_1*crpix1 + pc1_2*crpix2);
+	toFPA->x->coeff[1][0] = pc1_1;
+	toFPA->x->coeff[0][1] = pc1_2;
+	toFPA->x->mask[1][1]  = 1;
+
+	toFPA->y->coeff[0][0] = -(pc2_1*crpix1 + pc2_2*crpix2);
+	toFPA->y->coeff[1][0] = pc2_1;
+	toFPA->y->coeff[0][1] = pc2_2;
+	toFPA->y->mask[1][1]  = 1;
+
+	// projection from TPA to SKY
+	psProjection *toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type);
+
+	if (fpa->toSky == NULL) {
+	    // XXX for now, use the identity for TPA <--> FPA
+	    fpa->toTangentPlane = psPlaneDistortIdentity (1);
+	    fpa->fromTangentPlane = psPlaneDistortIdentity (1);
+	    fpa->toSky = toSky;
+	} else {
+	    if (fpa->toTangentPlane == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined");
+	    if (fpa->fromTangentPlane == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined");
+
+	    // adjust for common toSky for mosaic:
+	    // ignore the TPA since toTPA is identity
+	    // find the FPA coordinate of 0,0 for this chip.
+	    psPlane *fp = psPlaneAlloc();
+	    psPlane *chip = psPlaneAlloc();
+	    psSphere *sky = psSphereAlloc();
+	    chip->x = chip->y = 0;
 	
-	// XXX for the moment, assume toTPA is the identity transformation
-	toFPA->x->coeff[0][0] = tpa.x;
-	toFPA->y->coeff[0][0] = tpa.y;
-    } else {
-	// XXX for now, use the identity for TPA <--> FPA
-	toTPA = psPlaneDistortIdentity ();
-
-	// center of projection is (0,0) coordinate of TPA
-	toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, RAD_DEG*cdelt, RAD_DEG*cdelt, type);
-
-	*toTPAout = toTPA;
-	*toSkyOut = toSky;
+	    psPlaneTransformApply (fp, toFPA, chip); // find the focal-plane coordinate of this chip's 0,0 coordinate
+	    p_psDeproject (sky, fp, toSky); // find the RA,DEC coord of the focal-plane coordinate
+	    p_psProject (fp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
+
+	    toFPA->x->coeff[0][0] = fp->x;
+	    toFPA->y->coeff[0][0] = fp->y;
+
+	    rX = toSky->Xs/fpa->toSky->Xs;
+	    rY = toSky->Ys/fpa->toSky->Ys;
+
+	    // correct to common plate scale
+	    toFPA->x->coeff[1][0] *= rX;
+	    toFPA->x->coeff[0][1] *= rX;
+	    toFPA->y->coeff[1][0] *= rY;
+	    toFPA->y->coeff[0][1] *= rY;
+
+	    psFree (fp);
+	    psFree (sky);
+	    psFree (chip);
+	    psFree (toSky);
+	}
+
+	chip->toFPA = toFPA;
+	chip->fromFPA = p_psPlaneTransformLinearInvert(toFPA);
+
+	// remove the correction to the common plate scale
+	if (isMosaic) {
+	    chip->toFPA->x->coeff[1][0] /= rX;
+	    chip->toFPA->x->coeff[0][1] /= rX;
+	    chip->toFPA->y->coeff[1][0] /= rY;
+	    chip->toFPA->y->coeff[0][1] /= rY;
+	}
+
+	fprintf (stderr, "toFPA: %f %f  (%f,%f),(%f,%f)\n", 
+		 chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0], 
+		 chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1], 
+		 chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]);
+
+	fprintf (stderr, "frFPA: %f %f  (%f,%f),(%f,%f)\n", 
+		 chip->fromFPA->x->coeff[0][0], chip->fromFPA->y->coeff[0][0], 
+		 chip->fromFPA->x->coeff[1][0], chip->fromFPA->x->coeff[0][1], 
+		 chip->fromFPA->y->coeff[1][0], chip->fromFPA->y->coeff[0][1]);
+
+	fprintf (stderr, "field: %f,%f  %f,%f\n", 
+		 DEG_RAD*fpa->toSky->R, DEG_RAD*fpa->toSky->D, 
+		 3600*DEG_RAD*fpa->toSky->Xs, 3600*DEG_RAD*fpa->toSky->Ys);
+
     }
     return true;
 }
-
 
 // convert toFPA / toSky components to traditional WCS
@@ -169,11 +224,11 @@
     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", 	PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", 	PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
-    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", 	PS_META_REPLACE, "", toSky->Xs*DEG_RAD);
-    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", 	PS_META_REPLACE, "", toSky->Ys*DEG_RAD);
-
-    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]);
-    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]);
-    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]);
-    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", 	PS_META_REPLACE, "", toSky->Xs*DEG_RAD*plateScale);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", 	PS_META_REPLACE, "", toSky->Ys*DEG_RAD*plateScale);
+
+    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]/plateScale);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]/plateScale);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]/plateScale);
 
     // alternative representations use
@@ -184,18 +239,133 @@
 }
 
-psPlaneDistort *psPlaneDistortIdentity () {
+// convert toFPA / toSky components to traditional WCS
+// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
+bool pmAstromWriteBilevelChip (psPlaneTransform *toFPA, psMetadata *header, double plateScale) { 
+
+    psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---WRP");
+    psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--WRP");
+
+    // XXX not really right: needs to deal with non-identity toTP coeffs & plateScale
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", 	PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", 	PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", 	PS_META_REPLACE, "", 0.0);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", 	PS_META_REPLACE, "", 0.0);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", 	PS_META_REPLACE, "", 1.0);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", 	PS_META_REPLACE, "", 1.0);
+
+    if (toFPA->x->nX != toFPA->x->nY) psAbort ("psastro", "mis-matched tangent plane orders (1)");
+    if (toFPA->x->nX != toFPA->y->nX) psAbort ("psastro", "mis-matched tangent plane orders (2)");
+    if (toFPA->x->nX != toFPA->y->nY) psAbort ("psastro", "mis-matched tangent plane orders (3)");
+
+    switch (toFPA->x->nX) {
+      case 3:
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X3Y0", PS_META_REPLACE, "", toFPA->x->coeff[3][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y1", PS_META_REPLACE, "", toFPA->x->coeff[2][1]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y2", PS_META_REPLACE, "", toFPA->x->coeff[1][2]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y3", PS_META_REPLACE, "", toFPA->x->coeff[0][3]);
+
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X3Y0", PS_META_REPLACE, "", toFPA->y->coeff[3][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y1", PS_META_REPLACE, "", toFPA->y->coeff[2][1]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y2", PS_META_REPLACE, "", toFPA->y->coeff[1][2]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y3", PS_META_REPLACE, "", toFPA->y->coeff[0][3]);
+
+      case 2:
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y0", PS_META_REPLACE, "", toFPA->x->coeff[2][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y1", PS_META_REPLACE, "", toFPA->x->coeff[1][1]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y2", PS_META_REPLACE, "", toFPA->x->coeff[0][2]);
+
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y0", PS_META_REPLACE, "", toFPA->y->coeff[2][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y1", PS_META_REPLACE, "", toFPA->y->coeff[1][1]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y2", PS_META_REPLACE, "", toFPA->y->coeff[0][2]);
+
+      case 1:
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]/plateScale);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]/plateScale);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]/plateScale);
+	break;
+
+      case 0:
+	psAbort ("psastro", "invalid tangent plane order");
+    }
+    return true;
+}
+
+// convert toFPA / toSky components to traditional WCS
+// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
+psMetadata *pmAstromWriteBilevelMosaic (psProjection *toSky, psPlaneDistort *toTP, double plateScale) { 
+
+    psMetadata *header = psMetadataAlloc ();
+
+    psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---DIS");
+    psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--DIS");
+
+    // XXX need to handle the plateScale??
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", 	PS_META_REPLACE, "", toSky->R*DEG_RAD);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", 	PS_META_REPLACE, "", toSky->D*DEG_RAD);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", 	PS_META_REPLACE, "", 0.0);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", 	PS_META_REPLACE, "", 0.0);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", 	PS_META_REPLACE, "", toSky->Xs*DEG_RAD*plateScale);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", 	PS_META_REPLACE, "", toSky->Ys*DEG_RAD*plateScale);
+
+    if (toTP->x->nX != toTP->x->nY) psAbort ("psastro", "mis-matched tangent plane orders (1)");
+    if (toTP->x->nX != toTP->y->nX) psAbort ("psastro", "mis-matched tangent plane orders (2)");
+    if (toTP->x->nX != toTP->y->nY) psAbort ("psastro", "mis-matched tangent plane orders (3)");
+
+    switch (toTP->x->nX) {
+      case 3:
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X3Y0", PS_META_REPLACE, "", toTP->x->coeff[3][0][0][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y1", PS_META_REPLACE, "", toTP->x->coeff[2][1][0][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y2", PS_META_REPLACE, "", toTP->x->coeff[1][2][0][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y3", PS_META_REPLACE, "", toTP->x->coeff[0][3][0][0]);
+
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X3Y0", PS_META_REPLACE, "", toTP->y->coeff[3][0][0][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y1", PS_META_REPLACE, "", toTP->y->coeff[2][1][0][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y2", PS_META_REPLACE, "", toTP->y->coeff[1][2][0][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y3", PS_META_REPLACE, "", toTP->y->coeff[0][3][0][0]);
+
+      case 2:
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y0", PS_META_REPLACE, "", toTP->x->coeff[2][0][0][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y1", PS_META_REPLACE, "", toTP->x->coeff[1][1][0][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y2", PS_META_REPLACE, "", toTP->x->coeff[0][2][0][0]);
+
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y0", PS_META_REPLACE, "", toTP->y->coeff[2][0][0][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y1", PS_META_REPLACE, "", toTP->y->coeff[1][1][0][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y2", PS_META_REPLACE, "", toTP->y->coeff[0][2][0][0]);
+
+      case 1:
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toTP->x->coeff[1][0][0][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toTP->x->coeff[0][1][0][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toTP->y->coeff[1][0][0][0]);
+	psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toTP->y->coeff[0][1][0][0]);
+	break;
+
+      case 0:
+	psAbort ("psastro", "invalid tangent plane order");
+    }
+    return header;
+}
+
+// XXX for the moment, we mimic the limited Elixir mosastro orders
+psPlaneDistort *psPlaneDistortIdentity (int order) {
 
     psPlaneDistort *distort;
 
-    distort = psPlaneDistortAlloc (1, 1, 0, 0);
-    distort->x->coeff[0][0][0][0] = 0;
+    if (order < 1) psAbort ("psastro", "invalid order");
+    if (order > 3) psAbort ("psastro", "invalid order");
+    
+    // all coeffs and masks initially set to 0
+    distort = psPlaneDistortAlloc (order, order, 0, 0);
+
+    for (int i = 0; i <= order; i++) {
+	for (int j = 0; j <= order; j++) {
+	    if (i + j > order) {
+		distort->x->mask [1][1][0][0] = 1;
+		distort->y->mask [1][1][0][0] = 1;
+	    } 
+	}
+    }
     distort->x->coeff[1][0][0][0] = 1;
-    distort->x->coeff[0][1][0][0] = 0;
-    distort->x->mask [1][1][0][0] = 1;
-
-    distort->y->coeff[0][0][0][0] = 0;
-    distort->y->coeff[1][0][0][0] = 0;
     distort->y->coeff[0][1][0][0] = 1;
-    distort->y->mask [1][1][0][0] = 1;
 
     return distort;