Index: trunk/psastro/src/Makefile.am
===================================================================
--- trunk/psastro/src/Makefile.am	(revision 19313)
+++ trunk/psastro/src/Makefile.am	(revision 19314)
@@ -62,5 +62,8 @@
 	psastroUseModel.c           \
 	psastroMosaicAstrom.c       \
+	psastroMosaicDistortion.c   \
+	psastroMosaicFPtoTP.c       \
 	psastroMosaicGradients.c    \
+	psastroMosaicCorrectDistortion.c   \
 	psastroMosaicChipAstrom.c   \
 	psastroMosaicOneChip.c      \
Index: trunk/psastro/src/psastro.h
===================================================================
--- trunk/psastro/src/psastro.h	(revision 19313)
+++ trunk/psastro/src/psastro.h	(revision 19314)
@@ -55,5 +55,9 @@
 
 // mosaic fitting functions
-bool 		  psastroMosaicGradients (pmFPA *fpa, psMetadata *recipe);
+bool 		  psastroMosaicDistortion (pmFPA *fpa, psMetadata *recipe);
+psPlaneTransform *psastroMosaicFitRotAndScale (pmFPA *fpa);
+bool              psastroMosaicApplyRotAndScale (pmFPA *fpa, psPlaneTransform *TPtoFP);
+bool 		  psastroMosaicDistortionFromGradients (pmFPA *fpa, psMetadata *recipe);
+bool              psastroMosaicCorrectDistortion (pmFPA *fpa, psPlaneTransform *TPtoFP);
 bool 		  psastroMosaicCommonScale (pmFPA *fpa, psMetadata *recipe);
 bool 		  psastroMosaicAstrom (pmConfig *config);
Index: trunk/psastro/src/psastroMosaicCorrectDistortion.c
===================================================================
--- trunk/psastro/src/psastroMosaicCorrectDistortion.c	(revision 19314)
+++ trunk/psastro/src/psastroMosaicCorrectDistortion.c	(revision 19314)
@@ -0,0 +1,81 @@
+# include "psastroInternal.h"
+
+bool psastroMosaicCorrectDistortion (pmFPA *fpa, psPlaneTransform *TPtoFP) {
+
+    // invert the linear TPtoFP transform
+    psPlaneTransform *FPtoTP = p_psPlaneTransformLinearInvert (TPtoFP);
+    if (FPtoTP == NULL) {
+        psError (PS_ERR_UNKNOWN, false, "failed to invert TPtoFP\n");
+        return false;
+    }
+
+    // store the new coeffs in a new structure
+    psPlaneTransform *toTPAnew = psPlaneTransformAlloc(fpa->toTPA->x->nX, fpa->toTPA->x->nY);
+
+    // set the new coeffs, or set the mask
+    for (int i = 0; i <= toTPAnew->y->nX; i++) {
+        for (int j = 0; j <= toTPAnew->y->nY; j++) {
+
+	    // init the coeffs
+	    toTPAnew->x->coeff[i][j] = 0.0;
+	    toTPAnew->y->coeff[i][j] = 0.0;
+
+	    // if both are masked, mask the outpu
+            if ((fpa->toTPA->x->coeffMask[i][j] & PS_POLY_MASK_SET) && (fpa->toTPA->x->coeffMask[i][j] & PS_POLY_MASK_SET)) {
+		toTPAnew->x->coeffMask[i][j] = PS_POLY_MASK_SET;
+		toTPAnew->y->coeffMask[i][j] = PS_POLY_MASK_SET;
+		continue;
+	    } 
+
+	    // set the X terms
+	    toTPAnew->x->coeff[i][j] += (fpa->toTPA->x->coeffMask[i][j] & PS_POLY_MASK_SET) ? 0.0 : FPtoTP->x->coeff[1][0] * fpa->toTPA->x->coeff[i][j];
+	    toTPAnew->x->coeff[i][j] += (fpa->toTPA->y->coeffMask[i][j] & PS_POLY_MASK_SET) ? 0.0 : FPtoTP->x->coeff[0][1] * fpa->toTPA->y->coeff[i][j];
+
+	    // set the Y terms
+	    toTPAnew->y->coeff[i][j] += (fpa->toTPA->x->coeffMask[i][j] & PS_POLY_MASK_SET) ? 0.0 : FPtoTP->y->coeff[1][0] * fpa->toTPA->x->coeff[i][j];
+	    toTPAnew->y->coeff[i][j] += (fpa->toTPA->y->coeffMask[i][j] & PS_POLY_MASK_SET) ? 0.0 : FPtoTP->y->coeff[0][1] * fpa->toTPA->y->coeff[i][j];
+        }
+    }
+
+    // adjust the 0,0 terms:
+    toTPAnew->x->coeff[0][0] += FPtoTP->x->coeff[0][0];
+    toTPAnew->y->coeff[0][0] += FPtoTP->y->coeff[0][0];
+
+    psFree (fpa->toTPA);
+    fpa->toTPA = toTPAnew;
+
+    // invert toTPA to determine fromTPA. choose an appropriate region based on the dimensions
+    // of the complete FPA
+    psRegion *region = pmAstromFPAExtent (fpa);
+
+    psFree (fpa->fromTPA);
+    fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, *region, 50);
+    psFree (region);
+
+    if (fpa->fromTPA == NULL) {
+        psError (PS_ERR_UNKNOWN, false, "failed to invert fpa->toTPA\n");
+	psFree (FPtoTP);
+        return false;
+    }
+
+    psFree (FPtoTP);
+    return true;
+}
+
+// we have three coordinate systems and two polynomial transformations:
+// TP is the raw tangent plane coordinate system (aligned with RA,DEC)
+// FP is the raw focal plane coordinate system (aligned with camera X,Y)
+// dFP is the distorted focal plane coordinate system
+
+// fpa->toTPA is a polynomial transformation from FP to dFP
+//   L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j and 
+//   M(x,y) = \sum_i \sum_j B_{i,j} x^i y^j 
+//   where (x,y) are the FP coords and L,M  are the dFP coords
+
+// FPtoTP is a linear transformation from dFP to TP (we recover this from TPtoFP by inversion)
+// P(L,M) = r_xo + r_xx L + r_xy M
+// Q(L,M) = r_yo + r_yx L + r_yy M
+
+// we need to merge them into a single transformation:
+// P(x,y) = r_xo + r_xx * \sum_i \sum_j A_{i,j} x^i y^j + r_xy * \sum_i \sum_j B_{i,j} x^i y^j 
+// Q(x,y) = r_yo + r_yx * \sum_i \sum_j A_{i,j} x^i y^j + r_yy * \sum_i \sum_j B_{i,j} x^i y^j 
Index: trunk/psastro/src/psastroMosaicDistortion.c
===================================================================
--- trunk/psastro/src/psastroMosaicDistortion.c	(revision 19314)
+++ trunk/psastro/src/psastroMosaicDistortion.c	(revision 19314)
@@ -0,0 +1,44 @@
+# include "psastroInternal.h"
+
+bool psastroMosaicDistortion (pmFPA *fpa, psMetadata *recipe) {
+
+    // fit a linear transformation from reference TP to observed FP coords
+    psPlaneTransform *TPtoFP = psastroMosaicFitRotAndScale (fpa);
+    if (!TPtoFP) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to fit a linear TP correction\n");
+        return false;
+    }
+
+    // Correct the current reference star TP coordinates to the nearly-FP
+    // system.  We note two points here: 1) the corrected TP coordinates are
+    // NOT consistent with the sky coordinates, 2) the remaining differnce
+    // between the new TP reference coords and the observerd FP coordinates is
+    // only distortion
+    if (!psastroMosaicApplyRotAndScale (fpa, TPtoFP)) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to apply the linear TP correction\n");
+	psFree (TPtoFP);
+        return false;
+    }
+
+    if (!psastroMosaicDistortionFromGradients (fpa, recipe)) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to fit the distortion field\n");
+	psFree (TPtoFP);
+        return false;
+    }
+
+    if (!psastroMosaicCorrectDistortion (fpa, TPtoFP)) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to correct the distortion for the linear fit\n");
+	psFree (TPtoFP);
+        return false;
+    }
+	
+    if (!psastroMosaicSetAstrom (fpa)) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to apply mosaic distortion terms\n");
+	psFree (TPtoFP);
+        return false;
+    }
+
+    psFree (TPtoFP);
+    return true;
+}
+
Index: trunk/psastro/src/psastroMosaicFPtoTP.c
===================================================================
--- trunk/psastro/src/psastroMosaicFPtoTP.c	(revision 19314)
+++ trunk/psastro/src/psastroMosaicFPtoTP.c	(revision 19314)
@@ -0,0 +1,164 @@
+# include "psastroInternal.h"
+
+/************************************************/
+psPlaneTransform *psastroMosaicFitRotAndScale (pmFPA *fpa) {
+
+    pmChip *chip = NULL;
+    pmCell *cell = NULL;
+    pmReadout *readout = NULL;
+
+    // XXX first pass: fit and remove any linear fp->tp transformation
+    // accumulate FP(x,y) & TP(x,y) in a single set of vectors
+
+    psVector *X  = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *Y  = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *x  = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *y  = psVectorAllocEmpty(100, PS_TYPE_F32);
+
+    pmFPAview *view = pmFPAviewAlloc (0);
+
+    // int nPts = 0;
+    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
+        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+        if (!chip->process || !chip->file_exists) { continue; }
+	if (!chip->toFPA) { continue; }
+	
+	while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
+            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+            if (!cell->process || !cell->file_exists) { continue; }
+
+	    // process each of the readouts
+	    // XXX there can only be one readout per chip, right?
+	    while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
+		if (! readout->data_exists) { continue; }
+
+		// select the raw objects for this readout
+		psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
+		if (rawstars == NULL) { continue; }
+
+		// select the raw objects for this readout
+		psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
+		if (refstars == NULL) { continue; }
+
+		psArray *match = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH");
+		if (match == NULL) { continue; }
+
+		// take the matched stars, first fit
+		for (int i = 0; i < match->n; i++) {
+		    pmAstromMatch *pair = match->data[i];
+		    pmAstromObj *rawStar = rawstars->data[pair->raw];
+		    pmAstromObj *refStar = refstars->data[pair->ref];
+
+		    // independent variables
+		    X->data.F32[X->n] = refStar->TP->x;
+		    Y->data.F32[Y->n] = refStar->TP->y;
+
+		    // fitted values
+		    x->data.F32[x->n] = rawStar->FP->x;
+		    y->data.F32[y->n] = rawStar->FP->y;
+
+		    psVectorExtend (X, 100, 1);
+		    psVectorExtend (Y, 100, 1);
+		    psVectorExtend (x, 100, 1);
+		    psVectorExtend (y, 100, 1);
+		}
+	    }
+	}
+    }
+    // x->n = y->n = X->n = Y->n = nPts;
+
+    // linear fit without xy cross term
+    psPlaneTransform *map = psPlaneTransformAlloc (1, 1);
+    map->x->coeffMask[1][1] = PS_POLY_MASK_SET;
+    map->y->coeffMask[1][1] = PS_POLY_MASK_SET;
+
+    // constant errors
+    psVector *mask = psVectorAlloc (X->n, PS_TYPE_U8);
+    psVectorInit (mask, 0);
+
+    // the stats options supplied are used to perform the clip fitting
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    stats->clipIter = 1;
+
+    // fit TP-to-FP transformation
+
+    // we run 3 cycles clipping in each of x and y, with only one iteration each.
+    // XXX use the stats lookups functions to get the width and center
+    for (int i = 0; i < 3; i++) {
+	if (!psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, x, NULL, X, Y)) {
+            psError(PS_ERR_UNKNOWN, false, "failure in clip-fitting for x\n");
+	    psFree (map);
+	    map = NULL;
+	    goto escape;
+	}
+        psTrace ("psastro", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, x->n);
+
+        if (!psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, y, NULL, X, Y)) {
+            psError(PS_ERR_UNKNOWN, false, "failure in clip-fitting for y\n");
+	    psFree (map);
+	    map = NULL;
+	    goto escape;
+	}
+        psTrace ("psastro", 3, "y resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, y->n);
+    }
+
+escape:
+    psFree (x);
+    psFree (y);
+    psFree (X);
+    psFree (Y);
+    psFree (mask);
+    psFree (view);
+    psFree (stats);
+
+    return (map);
+}
+
+/*****************************/
+bool psastroMosaicApplyRotAndScale (pmFPA *fpa, psPlaneTransform *TPtoFP) {
+
+    pmChip *chip = NULL;
+    pmCell *cell = NULL;
+    pmReadout *readout = NULL;
+
+    pmFPAview *view = pmFPAviewAlloc (0);
+    psPlane newTP;
+
+    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
+        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+        if (!chip->process || !chip->file_exists) { continue; }
+	if (!chip->toFPA) { continue; }
+	
+	while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
+            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+            if (!cell->process || !cell->file_exists) { continue; }
+
+	    // process each of the readouts
+	    // XXX there can only be one readout per chip, right?
+	    while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
+		if (! readout->data_exists) { continue; }
+
+		// select the raw objects for this readout
+		psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
+		if (refstars == NULL) { continue; }
+
+		// take the matched stars, first fit
+		for (int i = 0; i < refstars->n; i++) {
+		    pmAstromObj *refStar = refstars->data[i];
+
+		    // Correct the current reference star TP coordinates to the nearly-FP
+		    // system.  We note two points here: 1) the corrected TP coordinates are
+		    // NOT consistent with the sky coordinates, 2) the remaining differnce
+		    // between the new TP reference coords and the observerd FP coordinates is
+		    // only distortion
+
+		    psPlaneTransformApply (&newTP, TPtoFP, refStar->TP);
+		    refStar->TP->x = newTP.x;
+		    refStar->TP->y = newTP.y;
+		}
+	    }
+	}
+    }
+    psFree (view);
+    return true;
+}