Index: trunk/psastro/src/pmAstrom.c
===================================================================
--- trunk/psastro/src/pmAstrom.c	(revision 5510)
+++ trunk/psastro/src/pmAstrom.c	(revision 5560)
@@ -1,21 +1,132 @@
 # include "pmAstrom.h"
 
+// sort by mag (descending)
+int pmAstromObjSortByFPX (const void **a, const void **b)
+{
+    pmAstromObj *A = *(pmAstromObj **)a;
+    pmAstromObj *B = *(pmAstromObj **)b;
+
+    psF32 diff = A->FP.x - B->FP.x;
+    if (diff > FLT_EPSILON) return (+1);
+    if (diff < FLT_EPSILON) return (-1);
+    return (0);
+}
+
+psPlane *psCoordReadoutToCell (psPlane *cellpix, psPlane *readpix, pmReadout *readout) {
+
+    if (cellpix == NULL) {
+	cellpix = psPlaneAlloc ();
+    }
+
+    cellpix->x = readpix->x*readout->colBins + readout->col0;
+    cellpix->y = readpix->y*readout->rowBins + readout->row0;
+
+    return (cellpix);
+}
+
+psPlane *psCoordCellToReadout (psPlane *readpix, psPlane *cellpix, pmReadout *readout) {
+
+    if (readpix == NULL) {
+	readpix = psPlaneAlloc ();
+    }
+
+    readpix->x = (cellpix->x - readout->col0) / readout->colBins;
+    readpix->y = (cellpix->y - readout->row0) / readout->rowBins;
+
+    return (readpix);
+}
+
+psPlane* psCoordChipToCell_EAM(psPlane* cellCoord,
+                           const psPlane* chipCoord,
+                           const pmCell* cell)
+{
+    PS_ASSERT_PTR_NON_NULL(chipCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
+
+    // XXX EAM : why was this being done?
+    // pmCell *tmpCell = pmCellInChip(chipCoord, cell->parent);
+    PS_ASSERT_PTR_NON_NULL(cell->toChip, NULL);
+    psPlaneTransform *tmpChipToCell = p_psPlaneTransformLinearInvert(cell->toChip);
+    PS_ASSERT_PTR_NON_NULL(tmpChipToCell, NULL);
+    cellCoord = psPlaneTransformApply(cellCoord, tmpChipToCell, chipCoord);
+    psFree(tmpChipToCell);
+    return(cellCoord);
+}
+
+bool psastroProjectFPA (pmFPA *fpa, char *starlist, bool toSky) {
+
+    bool status;
+
+    for (int i = 0; i < fpa->chips->n; i++) {
+	pmChip *chip = fpa->chips->data[i];
+	for (int j = 0; j < chip->cells->n; j++) {
+	    pmCell *cell = chip->cells->data[j];
+	    for (int k = 0; k < cell->readouts->n; k++) {
+		pmReadout *readout = cell->readouts->data[k];
+		psArray *stars = psMetadataLookupPtr (&status, readout->analysis, starlist);
+		if (toSky) {
+		    psastroProjectRawstars (stars, readout);
+		} else {
+		    psastroProjectRefstars (stars, readout);
+		}
+	    }
+	}
+    }
+    return true;
+}
+ 
+bool psastroProjectRawstars (psArray *stars, pmReadout *readout) {
+
+    pmCell *cell = readout->parent;
+    pmChip *chip = cell->parent;
+    pmFPA  *fpa  = chip->parent;
+
+    for (int i = 0; i < stars->n; i++) {
+	pmAstromObj *star = stars->data[i];
+	psCoordReadoutToCell (&star->cell, &star->pix, readout);
+	psCoordCellToChip (&star->chip, &star->cell, cell);
+	psCoordChipToFPA (&star->FP, &star->chip, chip);
+	psCoordFPAToTP (&star->TP, &star->FP, 0.0, 0.0, fpa);
+	psCoordTPToSky (&star->sky, &star->TP, fpa->projection);
+    }
+    return true;
+}
+ 
+bool psastroProjectRefstars (psArray *stars, pmReadout *readout) {
+
+    pmCell *cell = readout->parent;
+    pmChip *chip = cell->parent;
+    pmFPA  *fpa  = chip->parent;
+
+    for (int i = 0; i < stars->n; i++) {
+	pmAstromObj *star = stars->data[i];
+	psCoordSkyToTP (&star->TP, &star->sky, fpa->projection);
+	psCoordTPToFPA (&star->FP, &star->TP, 0.0, 0.0, fpa);
+	psCoordFPAToChip (&star->chip, &star->FP, chip);
+	psCoordChipToCell_EAM (&star->cell, &star->chip, cell);
+	psCoordCellToReadout (&star->pix, &star->cell, readout);
+    }
+    return true;
+}
+ 
 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;
+    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, pmAstromObjSortByFPX);
+    st2 = psArraySort (st2, pmAstromObjSortByFPX);
+
+    psArray *matches = psArrayAlloc (100);
+    matches->n = 0;
 
     int i = 0;
     int j = 0;
-
-    double RADIUS = psMetadataLookupR32 (&status, myHeader, "ASTROM_MATCH_RADIUS");
-    double RADIUS_SQR = PS_SQR (RADIUS);
-
-    psArray *matches = psArrayAlloc (100);
-    matches->n = 0;
-
     while ((i < st1->n) && (j < st2->n)) {
 
@@ -34,6 +145,6 @@
 	    
 	    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;
+	    dY = ((pmAstromObj *)st1->data[i])->FP.y - ((pmAstromObj *)st2->data[j])->FP.y;
+	    dR = dX*dX + dY*dY;
 
 	    if (dR > RADIUS_SQR) {
@@ -44,5 +155,5 @@
 	    // got a match; add to output list
 	    pmAstromMatch *match = pmAstromMatchAlloc (i, j);
-	    psArrayAdd (matches, match, 100);
+	    psArrayAdd (matches, 100, match);
 
 	    j++;
@@ -51,47 +162,60 @@
 	i++;
     }
-    return (match);
+    return (matches);
 }
 
 // take two matched star lists and fit a psPlaneTransform between them
 // 
-psPlaneTransform *pmAstromMatchedListFit (psPlaneTransform *map, psArray *st1, psArray *st2, psArray *match, psMetadata *config) {
+psPlaneTransform *pmAstromMatchFit (psPlaneTransform *map, psArray *raw, psArray *ref, psArray *match, psMetadata *config) {
+
+    bool status;
+    pmAstromObj *rawStar, *refStar;
+    pmAstromMatch *pair;
 
     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);
-    psVector *Y = psVectorAlloc (match->n);
-    psVector *x = psVectorAlloc (match->n);
-    psVector *y = psVectorAlloc (match->n);
+	int nX = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NX");
+	if (!status) nX = 1;
+	int nY = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NY");
+	if (!status) nY = 1;
+	map = psPlaneTransformAlloc (nX, nY);
+	// XXX EAM : not clear that we are allowed to use orders > 1
+    }
+
+    psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+    stats->clipSigma = psMetadataLookupF32 (&status, config, "PSASTRO.CHIP.NSIGMA");
+    stats->clipIter  = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NITER");
+
+    // XXX EAM : clip fit seems to only work for F32!
+    // XXX EAM : clip fit fails to handle NULL error
+    psVector *X = psVectorAlloc (match->n, PS_TYPE_F32);
+    psVector *Y = psVectorAlloc (match->n, PS_TYPE_F32);
+    psVector *x = psVectorAlloc (match->n, PS_TYPE_F32);
+    psVector *y = psVectorAlloc (match->n, PS_TYPE_F32);
+    psVector *wt = psVectorAlloc (match->n, PS_TYPE_F32);
   
     // take the matched stars, first fit 
     for (int i = 0; i < match->n; i++) {
 
-	pair = match->data[i];
-	ob1 = st1->data[pair->i1];
-	ob2 = st2->data[pair->i2];
-
-	X->data.F64[i] = ob1->chip.x;
-	Y->data.F64[i] = ob1->chip.y;
-
-	x->data.F64[i] = ob2->FP.x;
-	y->data.F64[i] = ob2->FP.y;
+	pair    = match->data[i];
+	rawStar = raw->data[pair->i1];
+	refStar = ref->data[pair->i2];
+
+	X->data.F32[i] = rawStar->chip.x;
+	Y->data.F32[i] = rawStar->chip.y;
+
+	x->data.F32[i] = refStar->FP.x;
+	y->data.F32[i] = refStar->FP.y;
+
+	wt->data.F32[i] = 1.0;
     }
 
     // no masking, no errors
-    psVectorClipFitPolynomial2D (map->x, stats, NULL, 0, X, NULL, x, y);
-    psVectorClipFitPolynomial2D (map->y, stats, NULL, 0, Y, NULL, x, y);
+    psVectorClipFitPolynomial2D (map->x, stats, NULL, 0, X, wt, x, y);
+    psVectorClipFitPolynomial2D (map->y, stats, NULL, 0, Y, wt, x, y);
     psFree (x);
     psFree (y);
     psFree (X);
     psFree (Y);
+    psFree (wt);
     psFree (stats);
 
@@ -103,12 +227,14 @@
 psArray *pmAstromRotateObj (psArray *old, psPlane center, double angle) {
 
+    double X, Y;
     pmAstromObj *newObj;
+    pmAstromObj *oldObj;
 
     psArray *new = psArrayAlloc (old->n);
 
-    cs = cos(angle);
-    sn = sin(angle);
-    xCenter = center.x;
-    yCenter = center.y;
+    double cs = cos(angle);
+    double sn = sin(angle);
+    double xCenter = center.x;
+    double yCenter = center.y;
     
     for (int i = 0; i < old->n; i++) {
@@ -120,6 +246,6 @@
 	Y = oldObj->FP.y - yCenter;
 	
-	newObj->FP.x = X*cs + Y*sn;
-	newObj->FP.y = Y*cs - X*sn;
+	newObj->FP.x = X*cs + Y*sn + xCenter;
+	newObj->FP.y = Y*cs - X*sn + yCenter;
 	
 	new->data[i] = newObj;
@@ -148,8 +274,9 @@
   obj->TP.y = 0;
 
-  obj->sky.x = 0;
-  obj->sky.y = 0;
-
-  obj->mag = 0;
+  obj->sky.r = 0;
+  obj->sky.d = 0;
+
+  obj->Mag = 0;
+  obj->dMag = 0;
 
   return (obj);