Index: /trunk/psastro/doc/notes.txt
===================================================================
--- /trunk/psastro/doc/notes.txt	(revision 5510)
+++ /trunk/psastro/doc/notes.txt	(revision 5510)
@@ -0,0 +1,27 @@
+
+we have a few different astrometry circumstances for our complete
+image/chip hierarchy structure:
+
+- ChipAstrom with one readout, one cell per chip:
+
+  whether or not there is more than one chip, we need to have two
+  stages
+  
+    - in the first stage, the per-chip results are applied to the
+      chip.toFPA parameter set.  
+    - in the second stage, the chip results are collectively used to
+      modify the fpa.toSky terms.
+
+      - the average offsets of the chip positions relative to starting
+	coordinates are used to calculate a single average offset to
+	the boresite. 
+
+      - the average rotations of the chips relative to their starting
+        rotations are used to calculate a single average rotation of
+        the boresite
+
+
+- ChipAstrom with more than one readout and/or cell:
+
+  - match all stars up relative to first readout?
+   
Index: /trunk/psastro/src/pmAstrom.c
===================================================================
--- /trunk/psastro/src/pmAstrom.c	(revision 5509)
+++ /trunk/psastro/src/pmAstrom.c	(revision 5510)
@@ -54,5 +54,17 @@
 }
 
-pmAstromMatchedListFit (pmFPA *fpa, psArray *st1, psArray *st2, psArray *match, psMetadata *config) {
+// take two matched star lists and fit a psPlaneTransform between them
+// 
+psPlaneTransform *pmAstromMatchedListFit (psPlaneTransform *map, psArray *st1, psArray *st2, psArray *match, psMetadata *config) {
+
+    if (map == NULL) {
+	// int nX = psMetadataLookupS32 (&status, myHeader, "CHIP.NX");
+	// int nY = psMetadataLookupS32 (&status, myHeader, "CHIP.NY");
+	map = psPlaneTransform (2, 2);
+    }
+
+    psStats *stats = psStatsAlloc (CLIPPED_MEAN);
+    stats->nSigma = psMetadataLookupS32 (&status, myHeader, "CHIP.NITER");
+    stats->nSigma = psMetadataLookupS32 (&status, myHeader, "CHIP.NSIGMA");
 
     psVector *X = psVectorAlloc (match->n);
@@ -62,5 +74,5 @@
   
     // take the matched stars, first fit 
-    for (i = 0; i < match->n; i++) {
+    for (int i = 0; i < match->n; i++) {
 
 	pair = match->data[i];
@@ -68,19 +80,21 @@
 	ob2 = st2->data[pair->i2];
 
-	X->data.F64[i] = ob1->P;
-	Y->data.F64[i] = ob1->Q;
+	X->data.F64[i] = ob1->chip.x;
+	Y->data.F64[i] = ob1->chip.y;
 
-	x->data.F64[i] = ob2->P;
-	y->data.F64[i] = ob2->Q;
-
-	// use one or the other...
-	psPolynomial2D *xt = psVectorFitPolynomial2D (NULL, stats, NULL, 0, X, NULL, x, y);
-	psPolynomial2D *xt = psVectorClipFitPolynomial2D (NULL, stats, NULL, 0, X, NULL, x, y);
-
-	psPolynomial2D *yt = psVectorFitPolynomial2D (NULL, stats, NULL, 0, Y, NULL, x, y);
-	psPolynomial2D *yt = psVectorClipFitPolynomial2D (NULL, stats, NULL, 0, Y, NULL, x, y);
+	x->data.F64[i] = ob2->FP.x;
+	y->data.F64[i] = ob2->FP.y;
     }
 
-    // apply fitted polynomials to fpa
+    // no masking, no errors
+    psVectorClipFitPolynomial2D (map->x, stats, NULL, 0, X, NULL, x, y);
+    psVectorClipFitPolynomial2D (map->y, stats, NULL, 0, Y, NULL, x, y);
+    psFree (x);
+    psFree (y);
+    psFree (X);
+    psFree (Y);
+    psFree (stats);
+
+    return (map);
 }
 
@@ -103,9 +117,9 @@
 	newObj = pmAstromObjCopy (oldObj);
 
-	P = oldObj->FP.x - xCenter;
-	Q = oldObj->FP.y - yCenter;
+	X = oldObj->FP.x - xCenter;
+	Y = oldObj->FP.y - yCenter;
 	
-	newObj->FP.x = P*cs + Q*sn;
-	newObj->FP.y = Q*cs - P*sn;
+	newObj->FP.x = X*cs + Y*sn;
+	newObj->FP.y = Y*cs - X*sn;
 	
 	new->data[i] = newObj;
Index: /trunk/psastro/src/pmAstrom.h
===================================================================
--- /trunk/psastro/src/pmAstrom.h	(revision 5509)
+++ /trunk/psastro/src/pmAstrom.h	(revision 5510)
@@ -7,8 +7,10 @@
 
 typedef struct {
-    double X, Y;
-    double P, Q;
-    double L, M;
-    double R, D;
+    psPlane pix;
+    psPlane cell;
+    psPlane chip;
+    psPlane FP;
+    psPlane TP;
+    psSphere sky;
     double Mag, dMag;
 } pmAstromObj;
@@ -26,7 +28,7 @@
 
 typedef struct {
+    psPlane center;
+    psPlane offset;
     double angle;
-    double dP;
-    double dQ;
     double minMetric;
     double minVar;
Index: /trunk/psastro/src/pmAstromGrid.c
===================================================================
--- /trunk/psastro/src/pmAstromGrid.c	(revision 5509)
+++ /trunk/psastro/src/pmAstromGrid.c	(revision 5510)
@@ -47,5 +47,6 @@
 	st2r = pmAstromRotateObj (st2, center, angle);
 	newStat = pmAstromGridMatchAngle (st1, st2r, config);
-	newStat.angle = angle;
+	newStat.angle  = angle;
+	newStat.center = center;
 	if (newStat.minMetric < minStat.minMetric) {
 	    minStat = newStat;
@@ -147,6 +148,6 @@
 		    minMetric = metric;
 		    minVar    = var;
-		    minP      = i;
-		    minQ      = j;
+		    minX      = i;
+		    minY      = j;
 		}
 	    }
@@ -154,6 +155,6 @@
 
 	// convert the bin to delta-delta
-	matchStats.dX        = DP[minY][minX] / NP[minY][minX];
-	matchStats.dY        = DQ[minY][minX] / NP[minY][minX];
+	matchStats.offset.x  = DP[minY][minX] / NP[minY][minX];
+	matchStats.offset.y  = DQ[minY][minX] / NP[minY][minX];
 	matchStats.minMetric = minMetric;
 	matchStats.minVar    = minVar;
@@ -164,8 +165,8 @@
 
 // apply the measured FPA offset and rotation (stat) to the fpa astrom structures
-psFPA *pmAstromApplyGridMatch (psFPA *fpa, pmAstromGridMatchStat stat) {
+psFPA *pmAstromGridApply (psPlaneTransform *map, pmAstromGridMatchStat stat) {
 
-    // stat.angle modifies fpa.toTangentPlane (effective angle)
-    // stat.dP, stat.dQ modifies fpa.projection (Ro, Do)
+    // stat.angle, stat.center, stat.offse
+    // I think i need to know the center reference....
 
     return (fpa);
Index: /trunk/psastro/src/psastro.c
===================================================================
--- /trunk/psastro/src/psastro.c	(revision 5509)
+++ /trunk/psastro/src/psastro.c	(revision 5510)
@@ -22,12 +22,6 @@
     psastroProjectRawstars (subset);
 
-    // load the corresponding reference data (DVO command)
-    psArray *refstars = psastroLoadReference (argv[3]);
-
-    // use the header & config info to project refstars on the focal plane
-    psastroProjectRefstars (refstars, subset);
-
     // fpa and subset point to the same astrometry terms
-    psastroChipAstrom (refstars, subset);
+    psastroChipAstrom (subset, config);
 
     // write out data (cmp file)
Index: /trunk/psastro/src/psastroBuildFPA.c
===================================================================
--- /trunk/psastro/src/psastroBuildFPA.c	(revision 5509)
+++ /trunk/psastro/src/psastroBuildFPA.c	(revision 5510)
@@ -136,4 +136,7 @@
     // set toChip to identity as default
     // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581)
+    // XXX EAM : define these in the config?
+    // int nX = psMetadataLookupS32 (&status, myHeader, "CHIP.NX");
+    // int nY = psMetadataLookupS32 (&status, myHeader, "CHIP.NY");
     cell->toChip   = psPlaneTransformAlloc (2, 2);
     cell->toChip->x->coeff[1][0] = 1;
Index: /trunk/psastro/src/psastroUtils.c
===================================================================
--- /trunk/psastro/src/psastroUtils.c	(revision 5509)
+++ /trunk/psastro/src/psastroUtils.c	(revision 5510)
@@ -197,5 +197,5 @@
 }
 
-bool psastroProjectRawstars (pmFPA *fpa) {
+bool psastroProjectFPA (pmFPA *fpa, bool toSky) {
 
     for (int i = 0; i < fpa->chips->n; i++) {
@@ -205,37 +205,141 @@
 	    for (int k = 0; k < cell->readouts->n; k++) {
 		pmReadout *readout = cell->readouts->data[k];
-		for (int m = 0; m < readout->stars->n; m++) {
-		    pmAstromObj *star = readout->data[m];
-		    // apply readout offsets
-		    // apply cell PlaneTransform
-		    // apply chip PlaneTransform
-		    psPlaneTransformApply (star-> 
-
-	    }
-	    chip->cells->data[j] = cell;
-	}
-	fpa->chips->data[i] = chip;
-    }
-    return (fpa);
+		if (toSky) {
+		    psastroProjectRawstars (readout->stars, readout);
+		} else {
+		    psastroProjectRefstars (readout->stars, readout);
+		}
+	    }
+	}
+    }
+    return true;
 }
  
-bool psastroChipAstrom (psArray *refstars, pmFPA *subset) {
-
-    // do this loop over the readouts 
-
-    // find initial offset / rotation
-    stat = pmAstromGridMatch (subset, refstars, config);
-
-    pmAstromModifyFPA (fpa, stat);
-
-    // use fit result to re-project rawstars
-    psastroProjectRawstars (fpa, subset);
-    psastroProjectRefstars (fpa, refstars);
+bool psastroProjectRawstars (psArray *stars, pmReadout *readout) {
+
+    pmCell *cell = readout->cell;
+    pmChip *chip = cell->chip;
+    pmFPA  *fpa  = chip->cell;
+
+    for (int i = 0; i < readout->stars->n; i++) {
+	pmAstromObj *star = readout->stars->data[i];
+	psCoordReadoutToCell (&star->cell, &star->pix, readout);
+	psCoordCellToChip (&star->chip, &star->cell, cell);
+	psCoordChipToFP (&star->FP, &star->chip, chip);
+	psCoordFPtoTP (&star->TP, &star->FP, fpa);
+	psCoordTPtoSky (&star->sky, &star->TP, fpa);
+    }
+    return true;
+}
+ 
+bool psastroProjectRefstars (psArray *stars, pmReadout *readout) {
+
+    pmCell *cell = readout->cell;
+    pmChip *chip = cell->chip;
+    pmFPA  *fpa  = chip->cell;
+
+    for (int i = 0; i < stars->n; i++) {
+	pmAstromObj *star = stars->data[i];
+
+	psCoordSkytoTP (&star->TP, &star->sky, fpa);
+	psCoordTPtoFP (&star->FP, &star->TP, fpa);
+	psCoordFPtoChip (&star->chip, &star->FP, chip);
+	psCoordChipToCell (&star->cell, &star->chip, cell);
+	psCoordReadoutToCell (&star->pix, &star->cell, readout);
+    }
+    return true;
+}
+ 
+// measure per-chip astrometry terms 
+bool psastroChipAstrom (pmFPA *fpa, psMetadata *config) {
+
+    pmFPA *raw = pmFPACopy (fpa);
+
+    // first pass: measure the per-chip solutions, modify the chip.toFPA terms
+    for (int i = 0; i < fpa->chips->n; i++) {
+	pmChip *chip = fpa->chips->data[i];
+
+	// cells->n > 1 is not yet well-defined
+	if (chip->cells->n > 1) {
+	    psLogMsg ("pmSourcesReadCMP", 3, "undefined behavior for nCells > 1");
+	    return false;
+	}
+
+	for (int j = 0; j < chip->cells->n; j++) {
+	    pmCell *cell = chip->cells->data[j];
+
+	    // readouts->n > 1 is not yet well-defined
+	    if (cell->readouts->n > 1) {
+		psLogMsg ("pmSourcesReadCMP", 3, "undefined behavior for nReadouts > 1");
+		return false;
+	    }
+
+	    // load the corresponding reference data (DVO command)
+	    psArray *refstars = psastroLoadReference (cell->header, config);
+
+	    for (int k = 0; k < cell->readouts->n; k++) {
+
+		pmReadout *readout = cell->readouts->data[k];
+
+		// use the header & config info to project refstars on the focal plane
+		psastroProjectRefstars (refstars, readout);
+
+		// find initial offset / rotation
+		stat = pmAstromGridMatch (readout->stars, refstars, config);
+
+		// adjust the chip.toFPA terms only
+		pmAstromGridApply (chip->toFPA, stat);
+
+		// use fit result to re-project rawstars
+		psastroProjectRawstars (readout->stars, readout);
+		psastroProjectRefstars (refstars,       readout);
     
-    // use small radius to match stars
-    match = pmAstromRadiusMatch (rawstars, refstars, options);
-
-    // fit astrometric terms
-    output = pmAstromMatchedListFit (fpa, subset, refstars, match, options);
-}
-
+		// use small radius to match stars
+		match = pmAstromRadiusMatch (rawstars, refstars, options);
+
+		// fit astrometric terms
+		pmAstromMatchedListFit (chip->toFPA, readout->stars, refstars, match, options);
+	    }
+	}
+    }
+
+    // second stage: re-normalize the chip terms, passing the 
+    // average rotation and offset values to the fpa.toSky
+	    
+# if (0)
+    // code disabled for now 
+    // this is not as trivial as it seems at first:
+    // applying a new rotation to the FPA implies changes to the 
+    // (x,y) reference coordinates for the chips as well as a rotation
+    // I need to calculate both the x,y offset effect and the rotation
+    // in one step.  what about the effect of higher order terms in either
+    // the chip.toFPA or the fpa.toTPA? 
+    // this calculation also needs to account for scale effects between the 
+    // fpa and the chips.
+
+    // calculate the average rotation and boresite offset relative to raw
+    for (int i = 0; i < fpa->chips->n; i++) {
+	pmChip *iChip = raw->chips->data[i];
+	pmChip *oChip = fpa->chips->data[i];
+
+	T1 = psPlaneTransformGetRotation (iChip->toFPA);
+	T2 = psPlaneTransformGetRotation (oChip->toFPA);
+	dT += T1 - T2;
+	dN ++;
+
+	
+    }
+
+    dT /= dN;
+
+    psPlaneTransformSetRotation (fpa->toTPA, dT);
+	
+    // apply the new modifcations in rotation and boresite
+    for (int i = 0; i < fpa->chips->n; i++) {
+	pmChip *iChip = raw->chips->data[i];
+	pmChip *oChip = fpa->chips->data[i];
+
+	psPlaneTransformSetRotation (oChip->toFPA, -dT);
+    }
+# endif 
+}