Index: /trunk/psModules/src/astrom/pmAstrometryDistortion.c
===================================================================
--- /trunk/psModules/src/astrom/pmAstrometryDistortion.c	(revision 7281)
+++ /trunk/psModules/src/astrom/pmAstrometryDistortion.c	(revision 7282)
@@ -7,6 +7,6 @@
 *  @author EAM, IfA
 *
-*  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-05-19 15:49:37 $
+*  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-06-02 02:15:35 $
 *
 *  Copyright 2006 Institute for Astronomy, University of Hawaii
@@ -57,5 +57,5 @@
     // XXX we need to get the chip dimensions from the metadata somewhere
     int Nx = 2;
-    int Ny = 5;
+    int Ny = 2;
     int DX = 2080 / Nx;
     int DY = 4600 / Ny;
@@ -108,5 +108,8 @@
                 psVectorExtend (dP, 100, 1);
                 psVectorExtend (dQ, 100, 1);
-            }
+                Npts++;
+            }
+            if (Npts < 5)
+                continue;
 
             // fit the collection of positions and offsets with a local 1st order gradient
@@ -114,4 +117,6 @@
             psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
             psVector *mask = psVectorAlloc (Npts, PS_TYPE_MASK);
+            mask->n = Npts;
+            psVectorInit (mask, 0);
 
             pmAstromGradient *grad = pmAstromGradientAlloc ();
@@ -141,4 +146,7 @@
 bool pmAstromFitDistortion(pmFPA *fpa, psArray *grads, psMetadata *config)
 {
+
+    psPolynomial2D *localX = NULL;
+    psPolynomial2D *localY = NULL;
 
     // assign the gradient elements to psVectors for fitting
@@ -168,21 +176,117 @@
     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     psVector *mask = psVectorAlloc (grads->n, PS_TYPE_MASK);
-
-    // XXX the order of the gradient fits need to be 1 less than the distortion term
-    // XXX determine the distortion order(s) from the fpa->toTP structure
-
-    psPolynomial2D *localX = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
-    psPolynomial2D *localY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
-    // XXX set masks based on fpa->toTP
-
+    mask->n = grads->n;
+    psVectorInit (mask, 0);
+
+    // the order of the gradient fits need to be 1 less than the distortion term
+    // determine the gradient order(s) from the fpa->toTP structure
+    localX = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTangentPlane->x->nX-1, fpa->toTangentPlane->x->nY);
+    localY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTangentPlane->x->nX,   fpa->toTangentPlane->x->nY-1);
+
+    // set masks based on fpa->toTP
+    for (int i = 0; i < fpa->toTangentPlane->x->nX; i++) {
+        for (int j = 0; j < fpa->toTangentPlane->x->nY; j++) {
+            if (fpa->toTangentPlane->x->mask[i][j][0][0]) {
+                localX->mask[i-1][j] = 1;
+                localY->mask[i][j-1] = 1;
+            }
+        }
+    }
+
+    // fit the local gradients in each direction
     psVectorClipFitPolynomial2D (localX, stats, mask, 0xff, dPdL, NULL, L, M);
     psVectorClipFitPolynomial2D (localY, stats, mask, 0xff, dPdM, NULL, L, M);
 
-    // XXX update fpa->toTP distortion terms
-
+    // XXX test: generate fit
+    psVector *dPdLf = psPolynomial2DEvalVector (localX, L, M);
+    psVector *dPdMf = psPolynomial2DEvalVector (localY, L, M);
+
+    // update fpa->toTP distortion terms
+    fpa->toTangentPlane->x->coeff[0][0][0][0] = 0;
+    for (int i = 1; i < fpa->toTangentPlane->x->nX; i++) {
+        if (!fpa->toTangentPlane->x->mask[i][0][0][0]) {
+            fpa->toTangentPlane->x->coeff[i][0][0][0] = localX->coeff[i-1][0] / i;
+        }
+    }
+    for (int j = 1; j < fpa->toTangentPlane->x->nY; j++) {
+        if (!fpa->toTangentPlane->x->mask[0][j][0][0]) {
+            fpa->toTangentPlane->x->coeff[0][j][0][0] = localY->coeff[0][j-1] / j;
+        }
+    }
+    for (int i = 1; i < fpa->toTangentPlane->x->nX; i++) {
+        for (int j = 1; j < fpa->toTangentPlane->x->nY; j++) {
+            if (!fpa->toTangentPlane->x->mask[i][j][0][0]) {
+                fpa->toTangentPlane->x->coeff[i][j][0][0] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
+            }
+        }
+    }
+    psFree (localX);
+    psFree (localY);
+
+    // the order of the gradient fits need to be 1 less than the distortion term
+    // determine the gradient order(s) from the fpa->toTP structure
+    localX = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTangentPlane->y->nX-1, fpa->toTangentPlane->y->nY);
+    localY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, fpa->toTangentPlane->y->nX,   fpa->toTangentPlane->y->nY-1);
+
+    // set masks based on fpa->toTP
+    for (int i = 0; i < fpa->toTangentPlane->y->nX; i++) {
+        for (int j = 0; j < fpa->toTangentPlane->y->nY; j++) {
+            if (fpa->toTangentPlane->y->mask[i][j][0][0]) {
+                localX->mask[i-1][j] = 1;
+                localY->mask[i][j-1] = 1;
+            }
+        }
+    }
+
+    // fit the local gradients in each direction
     psVectorClipFitPolynomial2D (localX, stats, mask, 0xff, dQdL, NULL, L, M);
     psVectorClipFitPolynomial2D (localY, stats, mask, 0xff, dQdM, NULL, L, M);
 
+    // XXX test: generate fit
+    psVector *dQdLf = psPolynomial2DEvalVector (localX, L, M);
+    psVector *dQdMf = psPolynomial2DEvalVector (localY, L, M);
+
+    FILE *f = fopen ("grads.dat", "w");
+    for (int i = 0; i < dPdLf->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]);
+    }
+    fclose (f);
+
+    // update fpa->toTP distortion terms
+    fpa->toTangentPlane->y->coeff[0][0][0][0] = 0;
+    for (int i = 1; i < fpa->toTangentPlane->y->nX; i++) {
+        if (!fpa->toTangentPlane->y->mask[i][0][0][0]) {
+            fpa->toTangentPlane->y->coeff[i][0][0][0] = localX->coeff[i-1][0] / i;
+        }
+    }
+    for (int j = 1; j < fpa->toTangentPlane->y->nY; j++) {
+        if (!fpa->toTangentPlane->y->mask[0][j][0][0]) {
+            fpa->toTangentPlane->y->coeff[0][j][0][0] = localY->coeff[0][j-1] / j;
+        }
+    }
+    for (int i = 1; i < fpa->toTangentPlane->y->nX; i++) {
+        for (int j = 1; j < fpa->toTangentPlane->y->nY; j++) {
+            if (!fpa->toTangentPlane->y->mask[i][j][0][0]) {
+                fpa->toTangentPlane->y->coeff[i][j][0][0] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
+            }
+        }
+    }
+    psFree (localX);
+    psFree (localY);
+
     // XXX free unneeded structures
+    psFree (dPdL);
+    psFree (dPdM);
+    psFree (dQdL);
+    psFree (dQdM);
+    psFree (L);
+    psFree (M);
+    psFree (stats);
+    psFree (mask);
 
     return true;
Index: /trunk/psModules/src/astrom/pmAstrometryObjects.c
===================================================================
--- /trunk/psModules/src/astrom/pmAstrometryObjects.c	(revision 7281)
+++ /trunk/psModules/src/astrom/pmAstrometryObjects.c	(revision 7282)
@@ -8,6 +8,6 @@
 *  @author EAM, IfA
 *
-*  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-05-19 15:49:37 $
+*  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-06-02 02:15:35 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -34,4 +34,52 @@
  ******************************************************************************/
 
+// XXX an example (might be faster)
+void pmAstromObjSortByChipX2 (psArray *objects)
+{
+
+    int l,j,ir,i;
+    float v1, v2;
+    pmAstromObj *tmp;
+    pmAstromObj *obj;
+
+    int N = objects->n;
+
+    if (N < 2)
+        return;
+    l = N >> 1;
+    ir = N - 1;
+    for (;;) {
+        if (l > 0) {
+            tmp = objects->data[--l];
+        } else {
+            tmp = objects->data[ir];
+            objects->data[ir] = objects->data[0];
+            if (--ir == 0) {
+                objects->data[0] = tmp;
+                return;
+            }
+        }
+        i = l;
+        j = (l << 1) + 1;
+        while (j <= ir) {
+            obj = objects->data[j];
+            v1 = obj->chip->x;
+            obj = objects->data[j+1];
+            v2 = obj->chip->x;
+            if (j < ir && v1 < v2)
+                ++j;
+            obj = objects->data[j];
+            v1 = obj->chip->x;
+            v2 = tmp->chip->x;
+            if (v2 < v1) {
+                objects->data[i] = objects->data[j];
+                j += (i=j) + 1;
+            } else
+                j = ir + 1;
+        }
+        objects->data[i] = tmp;
+    }
+}
+
 int pmAstromObjSortByFPX(
     const void **a,
@@ -59,5 +107,28 @@
 }
 
-
+int pmAstromObjSortByChipX(
+    const void **a,
+    const void **b)
+{
+    if (PM_ASTROMETRYOBJECTS_DEBUG) {
+        PS_ASSERT_PTR_NON_NULL(a, 0);
+        PS_ASSERT_PTR_NON_NULL(*a, 0);
+        PS_ASSERT_PTR_NON_NULL(b, 0);
+        PS_ASSERT_PTR_NON_NULL(*b, 0);
+    }
+
+    pmAstromObj *A = *(pmAstromObj **)a;
+    pmAstromObj *B = *(pmAstromObj **)b;
+
+    psF32 diff = A->chip->x - B->chip->x;
+    if (diff > FLT_EPSILON) {
+        return (+1);
+    }
+
+    if (diff < FLT_EPSILON) {
+        return (-1);
+    }
+    return (0);
+}
 
 /******************************************************************************
@@ -116,4 +187,7 @@
     st2 = psArraySort (st2, pmAstromObjSortByFPX);
 
+    // pmAstromObjSortByChipX2 (st1);
+    // pmAstromObjSortByChipX2 (st2);
+
     psArray *matches = psArrayAlloc (100);
     matches->n = 0;
@@ -124,9 +198,9 @@
 
         dX = ((pmAstromObj *)st1->data[i])->FP->x - ((pmAstromObj *)st2->data[j])->FP->x;
-        if (dX < -RADIUS) {
+        if (dX <= -RADIUS) {
             i++;
             continue;
         }
-        if (dX > +RADIUS) {
+        if (dX >= +RADIUS) {
             j++;
             continue;
@@ -134,5 +208,5 @@
 
         jStart = j;
-        while ((dX > -RADIUS) && (j < st2->n)) {
+        while ((fabs(dX) < RADIUS) && (j < st2->n)) {
 
             dX = ((pmAstromObj *)st1->data[i])->FP->x - ((pmAstromObj *)st2->data[j])->FP->x;
@@ -155,5 +229,72 @@
         i++;
     }
-    psLogMsg (__func__, 3, "radius match: %d pairs\n", matches->n);
+    psLogMsg (__func__, 3, "radius match: %d pairs (radius: %f)\n", matches->n, RADIUS);
+    return (matches);
+}
+
+
+
+/******************************************************************************
+pmAstromRadiusMatchChip(st1, st2, config) : match objects on the chip
+ ******************************************************************************/
+psArray *pmAstromRadiusMatchChip(
+    psArray *st1,
+    psArray *st2,
+    psMetadata *config)
+{
+    PS_ASSERT_PTR_NON_NULL(st1, NULL);
+    PS_ASSERT_PTR_NON_NULL(st2, NULL);
+    PS_ASSERT_PTR_NON_NULL(config, NULL);
+    bool status;
+    int jStart;
+    double dX, dY, dR;
+
+    double RADIUS = psMetadataLookupF32 (&status, config, "PSASTRO.MATCH.RADIUS");
+    double RADIUS_SQR = PS_SQR (RADIUS);
+
+    // sort both lists by Focal Plane X coord
+    st1 = psArraySort (st1, pmAstromObjSortByChipX);
+    st2 = psArraySort (st2, pmAstromObjSortByChipX);
+
+    psArray *matches = psArrayAlloc (100);
+    matches->n = 0;
+
+    int i = 0;
+    int j = 0;
+    while ((i < st1->n) && (j < st2->n)) {
+
+        dX = ((pmAstromObj *)st1->data[i])->chip->x - ((pmAstromObj *)st2->data[j])->chip->x;
+        if (dX <= -RADIUS) {
+            i++;
+            continue;
+        }
+        if (dX >= +RADIUS) {
+            j++;
+            continue;
+        }
+
+        jStart = j;
+        while ((dX > -RADIUS) && (j < st2->n)) {
+
+            dX = ((pmAstromObj *)st1->data[i])->chip->x - ((pmAstromObj *)st2->data[j])->chip->x;
+            dY = ((pmAstromObj *)st1->data[i])->chip->y - ((pmAstromObj *)st2->data[j])->chip->y;
+            dR = dX*dX + dY*dY;
+
+            if (dR > RADIUS_SQR) {
+                j++;
+                continue;
+            }
+
+            // got a match; add to output list
+            pmAstromMatch *match = pmAstromMatchAlloc (i, j);
+            psArrayAdd (matches, 100, match);
+            psFree (match);
+
+            j++;
+        }
+        j = jStart;
+        i++;
+    }
+    psLogMsg (__func__, 3, "radius match: %d pairs (radius: %f)\n", matches->n, RADIUS);
     return (matches);
 }
@@ -235,5 +376,8 @@
     // fit chip-to-FPA transformation
     psVectorClipFitPolynomial2D (map->x, stats, xMask, 0, x, wt, X, Y);
+    psLogMsg ("psModules.astrom", 3, "x resid: %f +/- %f (%d of %d)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, x->n);
+
     psVectorClipFitPolynomial2D (map->y, stats, yMask, 0, y, wt, X, Y);
+    psLogMsg ("psModules.astrom", 3, "y resid: %f +/- %f (%d of %d)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, y->n);
 
     psFree (x);
Index: /trunk/psModules/src/astrom/pmAstrometryObjects.h
===================================================================
--- /trunk/psModules/src/astrom/pmAstrometryObjects.h	(revision 7281)
+++ /trunk/psModules/src/astrom/pmAstrometryObjects.h	(revision 7282)
@@ -8,6 +8,6 @@
 *  @author EAM, IfA
 *
-*  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-05-19 15:49:37 $
+*  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-06-02 02:15:35 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -92,4 +92,10 @@
  */
 psArray *pmAstromRadiusMatch(
+    psArray *st1,
+    psArray *st2,
+    psMetadata *config
+);
+
+psArray *pmAstromRadiusMatchChip(
     psArray *st1,
     psArray *st2,