Index: trunk/psModules/src/astrom/pmAstrometryDistortion.c
===================================================================
--- trunk/psModules/src/astrom/pmAstrometryDistortion.c	(revision 10851)
+++ trunk/psModules/src/astrom/pmAstrometryDistortion.c	(revision 10859)
@@ -7,6 +7,6 @@
 *  @author EAM, IfA
 *
-*  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-12-29 18:29:59 $
+*  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-12-30 03:27:19 $
 *
 *  Copyright 2006 Institute for Astronomy, University of Hawaii
@@ -128,6 +128,6 @@
             fscanf (stdin, "%c", &key);
 
-            // fit the collection of positions and offsets with a local 1st order gradient
-            // 3hi/3lo sigma clipping on the rflux vs metric fit
+            // stats structure and mask for use in measuring the clipping statistic
+            // this analysis has too few data points to use the robust median method
             psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
             psVector *mask = psVectorAlloc (Npts, PS_TYPE_MASK);
@@ -136,11 +136,7 @@
             pmAstromGradient *grad = pmAstromGradientAlloc ();
 
-            psVectorStats (stats, L, NULL, NULL, 0);
-            grad->FP.x = stats->sampleMedian;
-
-            psVectorStats (stats, M, NULL, NULL, 0);
-            grad->FP.y = stats->sampleMedian;
-
-            // test how many survive and include / exclude points on that basis
+            // fit the collection of positions and offsets with a local 1st order gradient
+            // apply 3hi/3lo sigma clipping to the fitted data values
+            // the mask is used to mark the points which pass / fail the fit
             if (!psVectorClipFitPolynomial2D (local, stats, mask, 0xff, dP, NULL, L, M)) {
                 psFree (grad);
@@ -151,4 +147,7 @@
             grad->dTPdM.x = local->coeff[0][1];
 
+            // fit the collection of positions and offsets with a local 1st order gradient
+            // apply 3hi/3lo sigma clipping to the fitted data values
+            // the mask is used to mark the points which pass / fail the fit
             if (!psVectorClipFitPolynomial2D (local, stats, mask, 0xff, dQ, NULL, L, M)) {
                 psFree (grad);
@@ -159,6 +158,15 @@
             grad->dTPdM.y = local->coeff[0][1];
 
+            // also measure the L and M median positions as a representative coordinate
+            psVectorStats (stats, L, NULL, NULL, 0);
+            grad->FP.x = stats->sampleMedian;
+
+            psVectorStats (stats, M, NULL, NULL, 0);
+            grad->FP.y = stats->sampleMedian;
+
             psArrayAdd (gradients, 100, grad);
             psFree (grad);
+            psFree (stats);
+            psFree (mask);
         }
     }
@@ -210,4 +218,6 @@
     fscanf (stdin, "%c", &key);
 
+    // mask and stats structure for measuring the clipping statistic
+    // this analysis has too few data points to use the robust median method
     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     psVector *mask = psVectorAlloc (gradients->n, PS_TYPE_MASK);
@@ -217,8 +227,7 @@
     // determine the gradient order(s) from the fpa->toTP structure
     localX = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTPA->x->nX-1, fpa->toTPA->x->nY);
-
     localY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTPA->x->nX,   fpa->toTPA->x->nY-1);
 
-    // set masks based on fpa->toTP
+    // set masks based on fpa->toTPA
     for (int i = 0; i <= fpa->toTPA->x->nX; i++) {
         for (int j = 0; j <= fpa->toTPA->x->nY; j++) {
@@ -246,6 +255,6 @@
 
     // XXX test: generate fit
-    psVector *dPdLf = psPolynomial2DEvalVector (localX, L, M);
-    psVector *dPdMf = psPolynomial2DEvalVector (localY, L, M);
+    psVector *dPdLfit = psPolynomial2DEvalVector (localX, L, M);
+    psVector *dPdMfit = psPolynomial2DEvalVector (localY, L, M);
 
     // update fpa->toTP distortion terms
@@ -291,17 +300,21 @@
 
     // XXX test: generate fit
-    psVector *dQdLf = psPolynomial2DEvalVector (localX, L, M);
-    psVector *dQdMf = psPolynomial2DEvalVector (localY, L, M);
+    psVector *dQdLfit = psPolynomial2DEvalVector (localX, L, M);
+    psVector *dQdMfit = psPolynomial2DEvalVector (localY, L, M);
 
     FILE *f = fopen ("gradients.dat", "w");
-    for (int i = 0; i < dPdLf->n; i++) {
+    for (int i = 0; i < dPdLfit->n; i++) {
         fprintf (f, "%f %f  | %f %f  %f %f | %f %f  %f %f\n",
                  L->data.F32[i], M->data.F32[i],
                  dPdL->data.F32[i], dPdM->data.F32[i],
                  dQdL->data.F32[i], dQdM->data.F32[i],
-                 dPdLf->data.F32[i], dPdMf->data.F32[i],
-                 dQdLf->data.F32[i], dQdMf->data.F32[i]);
+                 dPdLfit->data.F32[i], dPdMfit->data.F32[i],
+                 dQdLfit->data.F32[i], dQdMfit->data.F32[i]);
     }
     fclose (f);
+    psFree (dPdLfit);
+    psFree (dPdMfit);
+    psFree (dQdLfit);
+    psFree (dQdMfit);
 
     // update fpa->toTP distortion terms
@@ -337,6 +350,9 @@
     psFree (mask);
 
-    // XXX need to reset the fromFPA terms here
-
+    // XXX need to reset the fromTPA terms here
+    // XXX choose an appropriate region based on the range of values
+    // in L and M?
+    psRegion region = psRegionSet (0, 4000, 0, 4000);
+    psPlaneTransformInvert(fpa->fromTPA, fpa->toTPA, region, 50);
 
     return true;