Index: /trunk/archive/modules/src/pmAstrom.c
===================================================================
--- /trunk/archive/modules/src/pmAstrom.c	(revision 5494)
+++ /trunk/archive/modules/src/pmAstrom.c	(revision 5495)
@@ -228,5 +228,4 @@
 }
 
-
 // rotate the focal-plane coordinates about the center coordinate
 // angle specified in radians
Index: /trunk/psastro/src/pmAstrom.c
===================================================================
--- /trunk/psastro/src/pmAstrom.c	(revision 5494)
+++ /trunk/psastro/src/pmAstrom.c	(revision 5495)
@@ -1,203 +1,11 @@
 # include "pmAstrom.h"
-
-// input: st1, st2 are psArray of psPlane? pmAstromObj
-
-static double maxOffset;  // maximum allowed offset between lists, in raw pixels
-static double Scale;      // grid pixel scale
-static double Offset;     // deltas to pixels
-
-bool p_pmAstromGridBin (int *pOut, int *qOut, double dP, int dQ) {
-
-    if (fabs(dP) > maxOffset) return false;
-    if (fabs(dQ) > maxOffset) return false;
-
-    *pOut = dP / Scale + Offset;
-    *qOut = dQ / Scale + Offset;
-    return true;
-}
-
-/* Illustration of the grid bins
-   dP        bin
-   -35:-25 -> 0     bin = dP / Scale + Offset
-   -25:-15 -> 1     Scale = 10
-   -15:-05 -> 2     Offset = 3.5
-   -05:+05 -> 3     nPix = 3 (maxOffset = 35 = (nPix + 0.5)*dPix
-   +05:+15 -> 4     dPix = 10
-   +15:+25 -> 5
-   +25:+35 -> 6
-
-   maxOffsetRequest = 30
-   nPix = (int) (maxOffset / dPix + 0.5);
-   maxOffset = (nPix + 0.5)*Scale;
-*/
-
-// apply the measured FPA offset and rotation (stat) to the fpa astrom structures
-psFPA *pmAstromApplyGridMatch (psFPA *fpa, pmAstromGridMatchStat stat) {
-
-    // stat.angle modifies fpa.toTangentPlane (effective angle)
-    // stat.dP, stat.dQ modifies fpa.projection (Ro, Do)
-
-    return (fpa);
-}
-
-pmAstromGridMatchStat pmAstromGridMatch (psArray *st1, psArray *st2, psMetadata *config) {
-
-    pmAstromGridMatchStat minStat, newStat;
-
-    // find center of the st2 field (focal-plane coords)
-    pMin = 1e10; pMax = -1e10;
-    qMin = 1e10; qMax = -1e10;
-    for (int i = 0; i < st2->n; i++) {
-	ob2 = (pmAstromObj *)st2->data[i];
-	pMin = PS_MIN (ob2->P, pMin);
-	pMax = PS_MAX (ob2->P, pMax);
-	qMin = PS_MIN (ob2->Q, qMin);
-	qMax = PS_MAX (ob2->Q, qMax);
-    }
-    pCenter = 0.5*(pMin + pMax);
-    qCenter = 0.5*(qMin + qMax);
-
-    double minAngle = psMetadataLookupF32 (&status, config, "MIN_ANGLE");
-    double maxAngle = psMetadataLookupF32 (&status, config, "MAX_ANGLE");
-    double delAngle = psMetadataLookupF32 (&status, config, "DEL_ANGLE");
-
-    minStat.minMetric = 1e10;
-    for (angle = minAngle; angle <= maxAngle; angle += delAngle) {
-	st2r = pmAstromRotateObj (st2, angle, pCenter, qCenter);
-	newStat = pmAstromGridMatchAngle (st1, st2r, config);
-	newStat.angle = angle;
-	if (newStat.minMetric < minStat.minMetric) {
-	    minStat = newStat;
-	}
-    }
-    return (minStat);
-}
-
-// rotate the focal-plane coordinates about the center coordinate
-// angle specified in radians
-psArray *pmAstromRotateObj (psArray *old, double angle, double pCenter, double qCenter) {
-
-    pmAstromObj *ob;
-
-    psArray *new = psArrayAlloc (old->n);
-
-    cs = cos(angle);
-    sn = sin(angle);
-    
-    for (int i = 0; i < old->n; i++) {
-
-	oldObj = (pmAstromObj *)old->data[i];
-	newObj = pmAstromObjCopy (oldObj);
-
-	P = oldObj->P - pCenter;
-	Q = oldObj->Q - qCenter;
-	
-	newObj->P = P*cs + Q*sn;
-	newObj->Q = Q*cs - P*sn;
-	
-	new->data[i] = newObj;
-    }
-    return (new);
-}
-
-// match the two lists using the binned delta-delta max
-pmAstromGridMatchStats pmAstromGridMatchAngle (psArray *st1, psArray *st2, psMetadata *config) {
-
-    // input lists must be projected to common focal plane
-    pmAstromGridMatchStats matchStats;
-    
-    double gridOffset = psMetadataLookupF32 (&status, config, "GRID_OFFSET");
-    double gridScale  = psMetadataLookupF32 (&status, config, "GRID_SCALE");
-
-    // set the static scaling factors
-    nPixHalf = (int)(gridOffset / gridScale + 0.5);
-    nPix     = 2*nPixHalf + 1;
-
-    maxOffpix = gridScale * (nPix + 0.5);
-    Offset    = maxOffpix / gridScale; 
-    Scale     = gridScale;
-
-    // ** can we assume the allocated image is init-ed?
-    gridNP = psImageAlloc (nPix, nPix, PS_TYPE_S32);
-    gridDP = psImageAlloc (nPix, nPix, PS_TYPE_F32);
-    gridDQ = psImageAlloc (nPix, nPix, PS_TYPE_F32);
-    gridD2 = psImageAlloc (nPix, nPix, PS_TYPE_F32);
-    psImageInit (gridNP);
-    psImageInit (gridDP);
-    psImageInit (gridDQ);
-    psImageInit (gridD2);
-
-    // short-cut names for grid images
-    NP = gridNP->data.F32;
-    DP = gridDP->data.F32;
-    DQ = gridDQ->data.F32;
-    D2 = gridD2->data.F32;
-
-    // accumulate grids for focal plane (L,M) matches
-    for (int i = 0; i < st1->n; i++) {
-	ob1 = (pmAstromObj *)st1->data[i];
-	for (int j = 0; j < st2->n; j++) {
-	    ob2 = (pmAstromObj *)st2->data[i];
-	    dP = ob1->P - ob2->P;
-	    if (fabs(dP) > maxOffset) continue;
-
-	    dQ = ob1->Q - ob2->Q;
-	    if (fabs(dQ) > maxOffset) continue;
-
-	    // find bin coordinates for this delta-delta
-	    p_pmAstromGridBin (&iP, &iQ, dP, dQ);
-		
-	    // accumulate dP2, dQ2? 
-	    NP[iQ][iP] ++;
-	    DP[iQ][iP] += dP;
-	    DQ[iQ][iP] += dQ;
-	    D2[iQ][iP] += PS_SQR(dP) + PS_SQR(dQ);
-	}
-    }
-
-    // find the max pixel
-    stats = psStatsAlloc (PS_STAT_MAX);
-    stats = psImageStats (stats, gridNP);
-
-    // only check bins with at least 1/2 of max bin
-    minNpt = 0.5*stats->max;
-
-    // find the 'best' bin
-    minMetric = minVar = 1e10;
-    minP = minQ = -1;
-    for (int j = 0; j < gridNP->nY; j++) {
-	for (int i = 0; i < gridNP->nX; i++) {
-
-	    if (NP[j][i] < minNpts) continue;
-
-	    // this metric emphasize a narrow peak with lots of sources over one with few.
-	    var = fabs((D2[j][i]/NP[j][i]) - PS_SQR(DP[j][i]/NP[j][i]) - PS_SQR(DQ[j][i]/NP[j][i]));
-	    metric = var / PS_SQR(PS_SQR(NP[j][i]));
-	    
-	    if (metric < minMetric) {
-		minMetric = metric;
-		minVar    = var;
-		minP      = i;
-		minQ      = j;
-	    }
-	}
-    }
-
-    // convert the bin to delta-delta
-    matchStats.dP        = DP[minQ][minP] / NP[minQ][minP];
-    matchStats.dQ        = DQ[minQ][minP] / NP[minQ][minP];
-    matchStats.minMetric = minMetric;
-    matchStats.minVar    = minVar;
-    matchStats.nMatch    = NP[minQ][minP];
-
-    return (matchStats);
-}
 
 psArray *pmAstromRadiusMatch (psArray *st1, psArray *st2, psMetadata *config) {
 
     // assume the lists are sorted, or sort in place?
-
     // sort st1 by P
     // sort st2 by P
+
+    double dX, dY;
 
     int i = 0;
@@ -212,10 +20,10 @@
     while ((i < st1->n) && (j < st2->n)) {
 
-	dP = ((pmAstromObj *)st1->data[i])->P - ((pmAstromObj *)st2->data[i])->P;
-	if (dP < -RADIUS) {
+	dX = ((pmAstromObj *)st1->data[i])->FP.x - ((pmAstromObj *)st2->data[j])->FP.x;
+	if (dX < -RADIUS) {
 	    i++;
 	    continue;
 	}
-	if (dP > +RADIUS) {
+	if (dX > +RADIUS) {
 	    j++;
 	    continue;
@@ -225,7 +33,7 @@
 	while ((dX > -RADIUS) && (j < st2->n)) {
 	    
-	    dP = ((pmAstromObj *)st1->data[i])->P - ((pmAstromObj *)st2->data[i])->P;
-	    dQ = ((pmAstromObj *)st1->data[i])->Q - ((pmAstromObj *)st2->data[i])->Q;
-	    dR = dP*dP + dQ*dQ;
+	    dX = ((pmAstromObj *)st1->data[i])->FP.x - ((pmAstromObj *)st2->data[j])->FP.x;
+	    dQ = ((pmAstromObj *)st1->data[i])->FP.y - ((pmAstromObj *)st2->data[j])->FP.y;
+	    dR = dX*dX + dQ*dQ;
 
 	    if (dR > RADIUS_SQR) {
@@ -276,2 +84,85 @@
     // apply fitted polynomials to fpa
 }
+
+// rotate the focal-plane coordinates about the center coordinate
+// angle specified in radians
+psArray *pmAstromRotateObj (psArray *old, psPlane center, double angle) {
+
+    pmAstromObj *newObj;
+
+    psArray *new = psArrayAlloc (old->n);
+
+    cs = cos(angle);
+    sn = sin(angle);
+    xCenter = center.x;
+    yCenter = center.y;
+    
+    for (int i = 0; i < old->n; i++) {
+
+	oldObj = (pmAstromObj *)old->data[i];
+	newObj = pmAstromObjCopy (oldObj);
+
+	P = oldObj->FP.x - xCenter;
+	Q = oldObj->FP.y - yCenter;
+	
+	newObj->FP.x = P*cs + Q*sn;
+	newObj->FP.y = Q*cs - P*sn;
+	
+	new->data[i] = newObj;
+    }
+    return (new);
+}
+
+static void pmAstromObjFree (pmAstromObj *obj) {
+
+  if (obj == NULL) return;
+  return;
+}
+
+pmAstromObj *pmAstromObjAlloc (void) {
+
+  pmAstromObj *obj = psAlloc (sizeof(pmAstromObj));
+  psMemSetDeallocator (obj, (psFreeFunc) pmAstromObjFree);
+
+  obj->pix.x = 0;
+  obj->pix.y = 0;
+
+  obj->FP.x = 0;
+  obj->FP.y = 0;
+
+  obj->TP.x = 0;
+  obj->TP.y = 0;
+
+  obj->sky.x = 0;
+  obj->sky.y = 0;
+
+  obj->mag = 0;
+
+  return (obj);
+}
+
+pmAstromObj *pmAstromObjCopy (pmAstromObj *old) {
+
+  pmAstromObj *obj = pmAstromObjAlloc ();
+
+  obj[0] = old[0];
+
+  return (obj);
+}
+
+static void pmAstromMatchFree (pmAstromMatch *match) {
+
+  if (match == NULL) return;
+  return;
+}
+
+pmAstromMatch *pmAstromMatchAlloc (int i1, int i2) {
+
+  pmAstromMatch *match = psAlloc (sizeof(pmAstromMatch));
+  psMemSetDeallocator(match, (psFreeFunc) pmAstromMatchFree);
+
+  match->i1 = i1;
+  match->i2 = i2;
+
+  return (match);
+}
Index: /trunk/psastro/src/pmAstromGrid.c
===================================================================
--- /trunk/psastro/src/pmAstromGrid.c	(revision 5495)
+++ /trunk/psastro/src/pmAstromGrid.c	(revision 5495)
@@ -0,0 +1,188 @@
+# include "pmAstrom.h"
+
+static double maxOffset;  // maximum allowed offset between lists, in raw pixels
+static double Scale;      // grid pixel scale
+static double Offset;     // deltas to pixels
+
+// local function to convert x,y coords to grid bins
+// it requires the globals defined above
+static bool AstromGridBin (int *pOut, int *qOut, double dP, int dQ) {
+
+    if (fabs(dP) > maxOffset) return false;
+    if (fabs(dQ) > maxOffset) return false;
+
+    *pOut = dP / Scale + Offset;
+    *qOut = dQ / Scale + Offset;
+    return true;
+}
+
+// match two star lists
+pmAstromGridMatchStats pmAstromGridMatch (psArray *st1, psArray *st2, psMetadata *config) {
+
+    double xMin, xMax, yMin, yMax;
+    pmAstromObj *ob1, *ob2;
+
+    pmAstromGridMatchStat minStat, newStat;
+    psPlane center;
+
+    // find center of the st2 field (focal-plane coords)
+    xMin = yMin = +1e10; 
+    xMax = yMax = -1e10;
+    for (int i = 0; i < st2->n; i++) {
+	ob2 = (pmAstromObj *)st2->data[i];
+	xMin = PS_MIN (ob2->FP.x, xMin);
+	xMax = PS_MAX (ob2->FP.x, xMax);
+	yMin = PS_MIN (ob2->FP.y, yMin);
+	yMax = PS_MAX (ob2->FP.y, yMax);
+    }
+    center.x = 0.5*(xMin + xMax);
+    center.y = 0.5*(yMin + yMax);
+
+    double minAngle = psMetadataLookupF32 (&status, config, "GRID.MIN.ANGLE");
+    double maxAngle = psMetadataLookupF32 (&status, config, "GRID.MAX.ANGLE");
+    double delAngle = psMetadataLookupF32 (&status, config, "GRID.DEL.ANGLE");
+
+    minStat.minMetric = 1e10;
+    for (angle = minAngle; angle <= maxAngle; angle += delAngle) {
+	st2r = pmAstromRotateObj (st2, center, angle);
+	newStat = pmAstromGridMatchAngle (st1, st2r, config);
+	newStat.angle = angle;
+	if (newStat.minMetric < minStat.minMetric) {
+	    minStat = newStat;
+	}
+	psFree (st2r);
+    }
+    return (minStat);
+}
+
+// match the two lists using the binned delta-delta max
+pmAstromGridMatchStats pmAstromGridMatchAngle (psArray *st1, psArray *st2, psMetadata *config) {
+
+    int nPix;       // size of matching grid
+    double dX, dY;  // offset between a possible matched pair
+    int iX, iY;     // corresponding grid bin
+
+    pmAstromObj *ob1, *ob2 // short-cut pointers to the objects
+    pmAstromGridMatchStats matchStats; // output match statistics
+
+    // max allowed offset in either X or Y directions
+    double gridOffset = psMetadataLookupF32 (&status, config, "GRID.OFFSET"); 
+
+    // sampling scale of the grid
+    double gridScale  = psMetadataLookupF32 (&status, config, "GRID.SCALE");
+
+    // set the static scaling factors
+    nPix = (int)(gridOffset / gridScale + 0.5);  // half-grid
+    nPix = 2*nPix + 1;
+
+    // these are globals used by p_pmAstromGridBin
+    maxOffpix = gridScale * (nPix + 0.5);
+    Offset    = maxOffpix / gridScale; 
+    Scale     = gridScale;
+
+    // XXX EAM : can we assume the allocated image is init-ed?
+    psImage *gridNP = psImageAlloc (nPix, nPix, PS_TYPE_S32);
+    psImage *gridDX = psImageAlloc (nPix, nPix, PS_TYPE_F32);
+    psImage *gridDY = psImageAlloc (nPix, nPix, PS_TYPE_F32);
+    psImage *gridD2 = psImageAlloc (nPix, nPix, PS_TYPE_F32);
+    psImageInit (gridNP);
+    psImageInit (gridDX);
+    psImageInit (gridDY);
+    psImageInit (gridD2);
+
+    // short-cut names for grid images
+    psS32 *NP = gridNP->data.S32;
+    psF32 *DX = gridDP->data.F32;
+    psF32 *DY = gridDQ->data.F32;
+    psF32 *D2 = gridD2->data.F32;
+
+    // accumulate grids for focal plane (L,M) matches
+    for (int i = 0; i < st1->n; i++) {
+	ob1 = (pmAstromObj *)st1->data[i];
+	for (int j = 0; j < st2->n; j++) {
+	    ob2 = (pmAstromObj *)st2->data[i];
+	    dX = ob1->FP.x - ob2->FP.x;
+	    dY = ob1->FP.y - ob2->FP.y;
+
+	    // find bin coordinates for this delta-delta
+	    if (!p_pmAstromGridBin (&iX, &iY, dX, dY)) {
+		continue; // matched pair is too far offset
+
+	    }
+
+	    // accumulate bin stats
+	    NP[iY][iX] ++;
+	    DP[iY][iX] += dX;
+	    DQ[iY][iX] += dY;
+	    D2[iY][iX] += PS_SQR(dX) + PS_SQR(dY);
+	}
+    }
+
+    // now assess the grid images
+    { 
+	double minMetric = 1e10;
+	double minVar = 1e10;
+	double minX = -1;
+	double minY = -1;
+	double metric, var;
+
+	// find the max pixel
+	psStats *stats = psStatsAlloc (PS_STAT_MAX);
+	stats = psImageStats (stats, gridNP);
+
+	// only check bins with at least 1/2 of max bin
+	int minNpt = 0.5*stats->max;
+
+	// find the 'best' bin
+	for (int j = 0; j < gridNP->nY; j++) {
+	    for (int i = 0; i < gridNP->nX; i++) {
+
+		if (NP[j][i] < minNpts) continue;
+
+		// this metric emphasize a narrow peak with lots of sources over one with few.
+		var = fabs((D2[j][i]/NP[j][i]) - PS_SQR(DX[j][i]/NP[j][i]) - PS_SQR(DY[j][i]/NP[j][i]));
+		metric = var / PS_SQR(PS_SQR(NP[j][i]));
+	    
+		if (metric < minMetric) {
+		    minMetric = metric;
+		    minVar    = var;
+		    minP      = i;
+		    minQ      = j;
+		}
+	    }
+	}
+
+	// convert the bin to delta-delta
+	matchStats.dX        = DP[minY][minX] / NP[minY][minX];
+	matchStats.dY        = DQ[minY][minX] / NP[minY][minX];
+	matchStats.minMetric = minMetric;
+	matchStats.minVar    = minVar;
+	matchStats.nMatch    = NP[minY][minX];
+    }
+    return (matchStats);
+}
+
+// apply the measured FPA offset and rotation (stat) to the fpa astrom structures
+psFPA *pmAstromApplyGridMatch (psFPA *fpa, pmAstromGridMatchStat stat) {
+
+    // stat.angle modifies fpa.toTangentPlane (effective angle)
+    // stat.dP, stat.dQ modifies fpa.projection (Ro, Do)
+
+    return (fpa);
+}
+
+/* Illustration of the grid bins
+   dX        bin
+   -35:-25 -> 0     bin = dX / Scale + Offset
+   -25:-15 -> 1     Scale = 10
+   -15:-05 -> 2     Offset = 3.5
+   -05:+05 -> 3     nPix = 3 (maxOffset = 35 = (nPix + 0.5)*dXix
+   +05:+15 -> 4     dPix = 10
+   +15:+25 -> 5
+   +25:+35 -> 6
+
+   maxOffsetRequest = 30
+   nPix = (int) (maxOffset / dPix + 0.5);
+   maxOffset = (nPix + 0.5)*Scale;
+*/
+