Index: /trunk/psastro/Makefile
===================================================================
--- /trunk/psastro/Makefile	(revision 5559)
+++ /trunk/psastro/Makefile	(revision 5560)
@@ -1,5 +1,5 @@
-default: psphot
+default: psastro
 help:
-	@echo "USAGE: make psphot"
+	@echo "USAGE: make psastro"
 
 CC      =       gcc
@@ -12,69 +12,40 @@
 IPSLIB  :=      $(shell pslib-config --cflags)
 
-INCS	= 	$(IPSLIB)
-# LIBS	= 	-lpsmodule $(LPSLIB)
-LIBS	= 	$(LPSLIB)
+LPSMOD  :=      $(shell psmodule-config --libs)
+IPSMOD  :=      $(shell psmodule-config --cflags)
+
+INCS	= 	$(IPSLIB) $(IPSMOD)
+LIBS	= 	$(LPSLIB) $(LPSMOD)
 CFLAGS	=	$(INCS) -std=c99 -Wall -Werror -g
 # CFLAGS	=	$(INCS) -std=c99 -g
 LFLAGS	=	$(LIBS) 
 
-PSPHOT = \
-$(SRC)/psphot.$(ARCH).o            \
-$(SRC)/psphotArguments.$(ARCH).o   \
-$(SRC)/psphotSetup.$(ARCH).o	   \
-$(SRC)/psphotImageStats.$(ARCH).o  \
-$(SRC)/psphotSourceStats.$(ARCH).o \
-$(SRC)/psphotChoosePSF.$(ARCH).o   \
-$(SRC)/psphotApplyPSF.$(ARCH).o	   \
-$(SRC)/psphotFitGalaxies.$(ARCH).o \
-$(SRC)/psphotOutput.$(ARCH).o      \
-$(SRC)/psphotMarkPSF.$(ARCH).o     \
-$(SRC)/psphotSubtractPSF.$(ARCH).o \
-$(SRC)/psphotSortBySN.$(ARCH).o    \
-$(SRC)/psphotDefinePixels.$(ARCH).o \
-$(SRC)/psphotMagnitudes.$(ARCH).o  \
+PSASTRO = \
+$(SRC)/psastro.$(ARCH).o           \
+$(SRC)/psastroIO.$(ARCH).o	   \
+$(SRC)/psastroBuildFPA.$(ARCH).o   \
+$(SRC)/psastroArguments.$(ARCH).o  \
+$(SRC)/psastroUtils.$(ARCH).o	   \
 $(SRC)/psLibUtils.$(ARCH).o	   \
-$(SRC)/psLine.$(ARCH).o		   \
-$(SRC)/psMinimize.$(ARCH).o	   \
-$(SRC)/psPolynomials.$(ARCH).o	   \
-$(SRC)/psEllipse.$(ARCH).o         \
-$(SRC)/psModulesUtils.$(ARCH).o	   \
-$(SRC)/pmPeaksSigmaLimit.$(ARCH).o \
-$(SRC)/pmPSF.$(ARCH).o             \
-$(SRC)/pmPSFtry.$(ARCH).o          \
-$(SRC)/psImageData.$(ARCH).o	   \
-$(SRC)/pmModelGroup.$(ARCH).o       \
-$(SRC)/pmObjects_EAM.$(ARCH).o
+$(SRC)/pmAstrom.$(ARCH).o	   \
+$(SRC)/pmAstromGrid.$(ARCH).o
 
-MODELS = \
-$(SRC)/models/pmModel_GAUSS.c      \
-$(SRC)/models/pmModel_PGAUSS.c     \
-$(SRC)/models/pmModel_QGAUSS.c     \
-$(SRC)/models/pmModel_SGAUSS.c     \
-$(SRC)/models/pmModel_TGAUSS.c      
+PSASTRO-MKTEST = \
+$(SRC)/psastro-mktest.$(ARCH).o    \
+$(SRC)/psastroArguments.$(ARCH).o  \
+$(SRC)/psLibUtils.$(ARCH).o        \
+$(SRC)/pmAstrom.$(ARCH).o
 
-$(SRC)/pmModelGroup.$(ARCH).o: $(MODELS)
+# $(SRC)/psModUtils.$(ARCH).o	   \
 
-# deprecated
+psastro: $(BIN)/psastro.$(ARCH)
+$(BIN)/psastro.$(ARCH) : $(PSASTRO)
+$(PSASTRO): $(SRC)/psastro.h $(SRC)/pmAstrom.h
 
-FITSOURCE = \
-$(SRC)/fitsource.$(ARCH).o \
-$(SRC)/onesource.$(ARCH).o \
-$(SRC)/fs_args.$(ARCH).o \
-$(SRC)/psphot-utils.$(ARCH).o \
-$(SRC)/psPolynomials.$(ARCH).o \
-$(SRC)/psUtils.$(ARCH).o \
-$(SRC)/setup.$(ARCH).o \
-$(SRC)/LocalSky.$(ARCH).o
+psastro-mktest: $(BIN)/psastro-mktest.$(ARCH)
+$(BIN)/psastro-mktest.$(ARCH) : $(PSASTRO-MKTEST)
+$(PSASTRO-MKTEST): $(SRC)/psastro.h $(SRC)/pmAstrom.h
 
-psphot: $(BIN)/psphot.$(ARCH)
-$(BIN)/psphot.$(ARCH) : $(PSPHOT)
-$(PSPHOT): $(SRC)/psphot.h
-
-fitsource: $(BIN)/fitsource.$(ARCH)
-$(BIN)/fitsource.$(ARCH) : $(FITSOURCE)
-$(FITSOURCE): $(SRC)/psphot.h
-
-INSTALL = psphot
+INSTALL = psastro psastro-mktest
 
 # dependancy rules for binary code #########################
Index: /trunk/psastro/src/pmAstrom.c
===================================================================
--- /trunk/psastro/src/pmAstrom.c	(revision 5559)
+++ /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);
Index: /trunk/psastro/src/pmAstrom.h
===================================================================
--- /trunk/psastro/src/pmAstrom.h	(revision 5559)
+++ /trunk/psastro/src/pmAstrom.h	(revision 5560)
@@ -3,6 +3,5 @@
 # include <unistd.h>   // for unlink
 # include <pslib.h>
-
-psArray *pmAstromGridMatch (psArray *s1, psArray *s2, pmAstromGridMatchOpt *opt);
+# include <pmAstrometry.h>
 
 typedef struct {
@@ -22,16 +21,30 @@
 
 typedef struct {
-    double minAngle;
-    double maxAngle;
-    double delAngle;
-} pmAstromOpt;
-
-typedef struct {
     psPlane center;
     psPlane offset;
-    double angle;
-    double minMetric;
-    double minVar;
-    int    nMatch;
-} pmAstromGridMatchStats;
+    double  angle;
+    double  minMetric;
+    double  minVar;
+    int     nMatch;
+} pmAstromStats;
 
+pmAstromObj 	 *pmAstromObjAlloc (void);
+pmAstromObj 	 *pmAstromObjCopy (pmAstromObj *old);
+pmAstromMatch    *pmAstromMatchAlloc (int i1, int i2);
+psArray          *pmAstromRadiusMatch (psArray *st1, psArray *st2, psMetadata *config);
+psArray          *pmAstromRotateObj (psArray *old, psPlane center, double angle);
+psPlaneTransform *pmAstromMatchFit (psPlaneTransform *map, psArray *st1, psArray *st2, psArray *match, psMetadata *config);
+
+int               pmAstromObjSortByFPX (const void **a, const void **b);
+
+bool 		  psastroProjectFPA (pmFPA *fpa, char *starlist, bool toSky);
+bool 		  psastroProjectRawstars (psArray *stars, pmReadout *readout);
+bool 		  psastroProjectRefstars (psArray *stars, pmReadout *readout);
+
+psPlane* psCoordChipToCell_EAM(
+    psPlane* out,                      ///< a plane struct to recycle. If NULL, a new struct is created
+    const psPlane* in,                 ///< the Chip coordinate
+    const pmCell* cell                 ///< the cell of interest
+);
+
+psPlaneTransform *p_psPlaneTransformLinearInvert_EAM(psPlaneTransform *transform);
Index: /trunk/psastro/src/pmAstromGrid.c
===================================================================
--- /trunk/psastro/src/pmAstromGrid.c	(revision 5559)
+++ /trunk/psastro/src/pmAstromGrid.c	(revision 5560)
@@ -1,5 +1,5 @@
 # include "pmAstrom.h"
 
-static double maxOffset;  // maximum allowed offset between lists, in raw pixels
+static double maxOffpix;  // maximum allowed offset between lists, in raw pixels
 static double Scale;      // grid pixel scale
 static double Offset;     // deltas to pixels
@@ -7,113 +7,75 @@
 // 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;
+static bool AstromGridBin (int *dx, int *dy, double dX, double dY) {
+
+    if (fabs(dX) > maxOffpix) return false;
+    if (fabs(dY) > maxOffpix) return false;
+
+    *dx = dX / Scale + Offset;
+    *dy = dY / 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;
-	newStat.center = center;
-	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) {
-
+pmAstromStats pmAstromGridAngle (psArray *raw, psArray *ref, psMetadata *config) {
+
+    bool status;
     int nPix;       // size of matching grid
+    int nPixHalf;   // half-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
+    pmAstromObj *ob1, *ob2; // short-cut pointers to the objects
+    pmAstromStats stats;    // output match statistics
 
     // max allowed offset in either X or Y directions
-    double gridOffset = psMetadataLookupF32 (&status, config, "GRID.OFFSET"); 
+    double gridOffset = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.OFFSET"); 
 
     // sampling scale of the grid
-    double gridScale  = psMetadataLookupF32 (&status, config, "GRID.SCALE");
+    double gridScale  = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.SCALE");
 
     // set the static scaling factors
-    nPix = (int)(gridOffset / gridScale + 0.5);  // half-grid
-    nPix = 2*nPix + 1;
+    nPixHalf = (int)(gridOffset / gridScale + 0.5);  // half-grid
+    nPix = 2*nPixHalf + 1;                           // full grid width
 
     // these are globals used by p_pmAstromGridBin
-    maxOffpix = gridScale * (nPix + 0.5);
+    maxOffpix = gridScale * (nPixHalf + 0.5);            // max offset from true center
     Offset    = maxOffpix / gridScale; 
     Scale     = gridScale;
 
-    // XXX EAM : can we assume the allocated image is init-ed?
-    psImage *gridNP = psImageAlloc (nPix, nPix, PS_TYPE_S32);
+    // images used as accumulators for the loop below
+    psImage *gridNP = psImageAlloc (nPix, nPix, PS_TYPE_U32);
     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);
+    psImageInit (gridNP, 0);
+    psImageInit (gridDX, 0);
+    psImageInit (gridDY, 0);
+    psImageInit (gridD2, 0);
 
     // 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;
+    psS32 **NP = gridNP->data.S32;
+    psF32 **DX = gridDX->data.F32;
+    psF32 **DY = gridDY->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];
+    for (int i = 0; i < raw->n; i++) {
+	ob1 = (pmAstromObj *)raw->data[i];
+	for (int j = 0; j < ref->n; j++) {
+	    ob2 = (pmAstromObj *)ref->data[j];
 	    dX = ob1->FP.x - ob2->FP.x;
 	    dY = ob1->FP.y - ob2->FP.y;
 
+	    // fprintf (stderr, "dX,dY: %8.2f %8.2f : %8.2f %8.2f : %8.2f %8.2f\n", dX, dY, ob1->FP.x, ob2->FP.x, ob1->FP.y, ob2->FP.y);
 	    // find bin coordinates for this delta-delta
-	    if (!p_pmAstromGridBin (&iX, &iY, dX, dY)) {
+	    if (!AstromGridBin (&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;
+	    DX[iY][iX] += dX;
+	    DY[iY][iX] += dY;
 	    D2[iY][iX] += PS_SQR(dX) + PS_SQR(dY);
 	}
@@ -124,18 +86,19 @@
 	double minMetric = 1e10;
 	double minVar = 1e10;
-	double minX = -1;
-	double minY = -1;
+	int minX = -1;
+	int minY = -1;
 	double metric, var;
 
 	// find the max pixel
-	psStats *stats = psStatsAlloc (PS_STAT_MAX);
-	stats = psImageStats (stats, gridNP);
+	psStats *imStats = psStatsAlloc (PS_STAT_MAX);
+	imStats = psImageStats (imStats, gridNP, NULL, 0);
 
 	// only check bins with at least 1/2 of max bin
-	int minNpt = 0.5*stats->max;
+	int minNpts = 0.5*imStats->max;
+	fprintf (stderr, "minNpts: %d, max: %d\n", minNpts, (int)(imStats->max));
 
 	// find the 'best' bin
-	for (int j = 0; j < gridNP->nY; j++) {
-	    for (int i = 0; i < gridNP->nX; i++) {
+	for (int j = 0; j < gridNP->numRows; j++) {
+	    for (int i = 0; i < gridNP->numCols; i++) {
 
 		if (NP[j][i] < minNpts) continue;
@@ -145,4 +108,6 @@
 		metric = var / PS_SQR(PS_SQR(NP[j][i]));
 	    
+		fprintf (stderr, "try : %f %f (%d pts, %f var, %f met)\n", DX[j][i]/NP[j][i], DY[j][i]/NP[j][i], NP[j][i], var, metric);
+
 		if (metric < minMetric) {
 		    minMetric = metric;
@@ -155,20 +120,84 @@
 
 	// convert the bin to delta-delta
-	matchStats.offset.x  = DP[minY][minX] / NP[minY][minX];
-	matchStats.offset.y  = DQ[minY][minX] / NP[minY][minX];
-	matchStats.minMetric = minMetric;
-	matchStats.minVar    = minVar;
-	matchStats.nMatch    = NP[minY][minX];
-    }
-    return (matchStats);
+	stats.offset.x  = DX[minY][minX] / NP[minY][minX];
+	stats.offset.y  = DY[minY][minX] / NP[minY][minX];
+	stats.minMetric = minMetric;
+	stats.minVar    = minVar;
+	stats.nMatch    = NP[minY][minX];
+
+	// XXX EAM : This routine, and pmAstromGridMatch, need to handle failure cases better
+    }
+    return (stats);
+}
+
+// match two star lists
+pmAstromStats pmAstromGridMatch (psArray *raw, psArray *ref, psMetadata *config) {
+
+    bool status;
+    double xMin, xMax, yMin, yMax;
+    pmAstromObj *obj;
+    psArray *rot;
+
+    pmAstromStats minStat, newStat;
+    psPlane center;
+
+    // find center of the raw field (focal-plane coords)
+    xMin = yMin = +1e10; 
+    xMax = yMax = -1e10;
+    for (int i = 0; i < raw->n; i++) {
+	obj = (pmAstromObj *)raw->data[i];
+	xMin = PS_MIN (obj->FP.x, xMin);
+	xMax = PS_MAX (obj->FP.x, xMax);
+	yMin = PS_MIN (obj->FP.y, yMin);
+	yMax = PS_MAX (obj->FP.y, yMax);
+    }
+    center.x = 0.5*(xMin + xMax);
+    center.y = 0.5*(yMin + yMax);
+
+    double minAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MIN.ANGLE");
+    double maxAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MAX.ANGLE");
+    double delAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.DEL.ANGLE");
+
+    minStat.minMetric = 1e10;
+    for (double angle = minAngle; angle <= maxAngle; angle += delAngle) {
+	rot = pmAstromRotateObj (raw, center, angle);
+	newStat = pmAstromGridAngle (rot, ref, config);
+	newStat.angle  = angle;
+	newStat.center = center;
+	if (newStat.minMetric < minStat.minMetric) {
+	    minStat = newStat;
+	}
+	psFree (rot);
+    }
+    fprintf (stderr, "best: %f %f (%d pts, %f var, %f met)\n", newStat.offset.x, newStat.offset.y, newStat.nMatch, newStat.minVar, newStat.minMetric);
+    return (minStat);
 }
 
 // apply the measured FPA offset and rotation (stat) to the fpa astrom structures
-psFPA *pmAstromGridApply (psPlaneTransform *map, pmAstromGridMatchStat stat) {
-
-    // stat.angle, stat.center, stat.offse
-    // I think i need to know the center reference....
-
-    return (fpa);
+psPlaneTransform *pmAstromGridApply (psPlaneTransform *map, pmAstromStats stat) {
+
+    double cs = cos (stat.angle);
+    double sn = sin (stat.angle);
+    
+    double dx = (map->x->coeff[0][0] - stat.center.x);
+    double dy = (map->y->coeff[0][0] - stat.center.y);
+
+    // new offset 
+    map->x->coeff[0][0] =  cs*dx + sn*dy - stat.offset.x + stat.center.x;
+    map->y->coeff[0][0] = -sn*dx + cs*dy - stat.offset.y + stat.center.y;
+
+    // original rotation matrix
+    double pc1_1 = map->x->coeff[1][0];
+    double pc1_2 = map->x->coeff[0][1];
+    double pc2_1 = map->y->coeff[1][0];
+    double pc2_2 = map->y->coeff[0][1];
+
+    // new rotation matrix
+    map->x->coeff[1][0] = +cs*pc1_1 + sn*pc2_1;
+    map->x->coeff[0][1] = +cs*pc1_2 + sn*pc2_2;
+    map->y->coeff[1][0] = -sn*pc1_1 + cs*pc2_1;
+    map->y->coeff[0][1] = -sn*pc1_2 + cs*pc2_2;
+
+    return (map);
 }
 
Index: /trunk/psastro/src/psModUtils.c
===================================================================
--- /trunk/psastro/src/psModUtils.c	(revision 5559)
+++ /trunk/psastro/src/psModUtils.c	(revision 5560)
@@ -4,4 +4,5 @@
 
 // template sample alloc/free pair:
+# if (0)
 static void psFooFree (psFoo *foo) {
 
@@ -9,5 +10,4 @@
   return;
 }
-
 psFoo *psFooAlloc (int i1, int i2) {
 
@@ -17,4 +17,5 @@
   return (foo);
 }
+# endif
 
 // pmFPA
@@ -31,5 +32,5 @@
 }
 
-pmFPA *pmFPAAlloc (int i1, int i2) {
+pmFPA *pmFPAAlloc (void) {
 
   pmFPA *fpa = psAlloc (sizeof(pmFPA));
@@ -56,5 +57,5 @@
 }
 
-pmChip *pmChipAlloc (int i1, int i2) {
+pmChip *pmChipAlloc (void) {
 
   pmChip *chip = psAlloc (sizeof(pmChip));
@@ -79,5 +80,5 @@
 }
 
-pmCell *pmCellAlloc (int i1, int i2) {
+pmCell *pmCellAlloc (void) {
 
   pmCell *cell = psAlloc (sizeof(pmCell));
@@ -95,10 +96,10 @@
   if (readout == NULL) return;
 
-  psFree (readout->objects);
+  psFree (readout->stars);
 
   return;
 }
 
-pmReadout *pmReadoutAlloc (int i1, int i2) {
+pmReadout *pmReadoutAlloc (void) {
 
   pmReadout *readout = psAlloc (sizeof(pmReadout));
@@ -109,5 +110,5 @@
   readout->colBins = 1;
   readout->rowBins = 1;
-  readout->objects = NULL;
+  readout->stars = NULL;
 
   return (readout);
Index: /trunk/psastro/src/psastro.c
===================================================================
--- /trunk/psastro/src/psastro.c	(revision 5559)
+++ /trunk/psastro/src/psastro.c	(revision 5560)
@@ -4,6 +4,4 @@
 
     psMetadata *header = NULL;
-    pmAstromGridMatchStat stat;
-    pmAstromMatch *match;
 
     // load configuration information
@@ -11,20 +9,20 @@
 
     // load the input data (cmp file)
-    psArray *rawstars = pmSourcesReadCMP (&header, argv[1]);
+    psArray *rawstars = psastroReadCMP (&header, argv[2]);
 
     // simple layout interpretation, basic WCS conversion
-    psFPA *fpa = psastroBuildFPA (header, rawstars);
+    pmFPA *fpa = psastroBuildFPA (header, rawstars);
 
     // limit fit to bright stars only 
-    psFPA *subset = psastroSelectBrightStars (fpa, config);
-
-    // use the header & config info to project rawstars on the focal plane
-    psastroProjectRawstars (subset);
+    psastroSelectBrightStars (fpa, config);
 
     // fpa and subset point to the same astrometry terms
-    psastroChipAstrom (subset, config);
+    psastroChipAstrom (fpa, config);
 
     // write out data (cmp file)
-    psastroWriteCMP (fpa, argv[2]);
+    psastroWriteCMP (fpa, argv[3]);
+
+    psFree (config);
+    psFree (fpa);
 
     exit (0);
Index: /trunk/psastro/src/psastro.h
===================================================================
--- /trunk/psastro/src/psastro.h	(revision 5559)
+++ /trunk/psastro/src/psastro.h	(revision 5560)
@@ -3,5 +3,13 @@
 # include <unistd.h>   // for unlink
 # include <pslib.h>
+# include <pmAstrometry.h>
+# include "psLibUtils.h"
+# include "pmAstrom.h"
 
+# define fromTPA fromTangentPlane
+# define toTPA toTangentPlane
+# define toSky projection
+
+# if (0)
 typedef struct
 {
@@ -19,4 +27,5 @@
     psPlaneTransform *fromFPA;          ///< Transformation from FPA to chip coordinates
     psArray *cells;                     ///< The cells
+    pmFPA *fpa;
 }
 pmChip;
@@ -28,4 +37,5 @@
     psArray *readouts;                  ///< The readouts (referred to by number)
     psMetadata *header;                 ///< header always corresponds to a cell
+    pmChip *chip;
 }
 pmCell;
@@ -39,11 +49,43 @@
     int rowBins;                        ///< Amount of binning in y-dimension
     psArray *stars;			///< sources detected / measured on readout
+    pmCell *cell;
 }
 pmReadout;
+# endif
 
-// dummy structure for template code (psFooAlloc, etc)
-typedef struct
-{
-    int i;
-} 
-psFoo;
+pmAstromStats     pmAstromGridAngle (psArray *st1, psArray *st2, psMetadata *config);
+psPlaneTransform *pmAstromGridApply (psPlaneTransform *map, pmAstromStats stat);
+pmAstromStats     pmAstromGridMatch (psArray *st1, psArray *st2, psMetadata *config);
+int               pmAstromObjSortByMag (const void **a, const void **b);
+bool              pmCellInterpretWCS (pmCell *cell, psMetadata *header) ;
+pmFPA            *pmFPACopyAstrom (pmFPA *inFPA);
+psMetadata       *psastroArguments (int *argc, char **argv);
+pmFPA            *psastroBuildFPA (psMetadata *header, psArray *stars);
+bool 		  psastroChipAstrom (pmFPA *fpa, psMetadata *config);
+psArray          *psastroLoadReference (psMetadata *header, psMetadata *config);
+psArray          *psastroReadCMP (psMetadata **header, char *filename);
+bool              psastroSelectBrightStars (pmFPA *fpa, psMetadata *config);
+bool              psastroWriteCMP (pmFPA *fpa, char *filename);
+
+psMetadata *testArguments (int *argc, char **argv);
+
+# if (0)
+pmFPA            *pmFPAAlloc (void);
+pmChip           *pmChipAlloc (void);
+pmCell           *pmCellAlloc (void);
+pmReadout        *pmReadoutAlloc (void);
+# endif
+
+# define SIGN(X)  (((X) == 0) ? 0 : ((fabs((double)(X))) / (X)))
+
+psPolynomial2D *psPolynomial2DCopy (psPolynomial2D *input);
+psPolynomial4D *psPolynomial4DCopy (psPolynomial4D *input);
+psPlaneDistort *psPlaneDistortCopy (psPlaneDistort *input);
+psPlaneTransform *psPlaneTransformCopy (psPlaneTransform *input);
+psProjection *psProjectionCopy (psProjection *input);
+psPlane *psCoordReadoutToCell (psPlane *cellpix, psPlane *readpix, pmReadout *readout);
+psPlane *psCoordCellToReadout (psPlane *readpix, psPlane *cellpix, pmReadout *readout);
+double psPlaneTransformGetRotation (psPlaneTransform *map);
+
+psPlaneDistort *psPlaneDistortInvert(psPlaneDistort *distort);
+psPlaneTransform *psPlaneTransformLinearInvert(psPlaneTransform *transform);
Index: /trunk/psastro/src/psastroArguments.c
===================================================================
--- /trunk/psastro/src/psastroArguments.c	(revision 5559)
+++ /trunk/psastro/src/psastroArguments.c	(revision 5560)
@@ -1,28 +1,50 @@
-# include "psphot.h"
+# include "psastro.h"
 
-static void usage (void) ;
+static void usage (void);
+static void usage_test (void);
 
-psMetadata *psphotArguments (int *argc, char **argv) {
+psMetadata *psastroArguments (int *argc, char **argv) {
 
-  int N, Nfail;
-  int mode = PS_DATA_STRING | PS_META_REPLACE;
+    unsigned int Nfail;
 
-  // basic pslib options
-  fprintf (stderr, "starting... %s\n", psLibVersion());
-  psLogSetFormat ("M");
-  psLogArguments (argc, argv);
-  psTraceArguments (argc, argv);
+    // basic pslib options
+    fprintf (stderr, "starting... %s\n", psLibVersion());
+    psLogSetFormat ("M");
+    psLogArguments (argc, argv);
+    psTraceArguments (argc, argv);
 
-  if (*argc != 2) usage ();
+    if (*argc != 4) usage ();
 
-  // load config information
-  psMetadata *config = psMetadataAlloc ();
-  config = psMetadataConfigParse (config, &Nfail, argv[3], FALSE);
-  return (config);
+    // load config information
+    psMetadata *config = psMetadataAlloc ();
+    config = psMetadataConfigParse (config, &Nfail, argv[1], FALSE);
+    return (config);
 }
 
 static void usage (void) {
-
-    fprintf (stderr, "USAGE: psastro (cmpfile)\n");
+    fprintf (stderr, "USAGE: psastro (config) (input) (output)\n");
     exit (2);
 }
+
+psMetadata *testArguments (int *argc, char **argv) {
+
+    unsigned int Nfail;
+
+    // basic pslib options
+    fprintf (stderr, "starting... %s\n", psLibVersion());
+    psLogSetFormat ("M");
+    psLogArguments (argc, argv);
+    psTraceArguments (argc, argv);
+
+    if (*argc != 6) usage_test ();
+
+    // load config information
+    psMetadata *config = psMetadataAlloc ();
+    config = psMetadataConfigParse (config, &Nfail, argv[1], FALSE);
+    return (config);
+}
+
+static void usage_test (void) {
+    fprintf (stderr, "USAGE: psastro-mktest (config) (cmpfile) (reflist) (rawfile) (rawfile)\n");
+    exit (2);
+}
Index: /trunk/psastro/src/psastroBuildFPA.c
===================================================================
--- /trunk/psastro/src/psastroBuildFPA.c	(revision 5559)
+++ /trunk/psastro/src/psastroBuildFPA.c	(revision 5560)
@@ -18,30 +18,32 @@
     // allocate the structures
 
-    pmFPA *fpa = pmFPAAlloc ();
+    pmFPA *fpa = pmFPAAlloc (NULL, NULL);
 
-    fpa->chips = psArrayAlloc (Nchips);
+    psArrayRealloc (fpa->chips, Nchips);
     for (int i = 0; i < Nchips; i++) {
+	pmChip *chip = pmChipAlloc (fpa);
 
-	pmChip *chip = pmChipAlloc ();
-	chip->fpa = fpa; // assign parent fpa (view only; don't free)
+	psArrayRealloc (chip->cells, Ncells);
+	for (int j = 0; j < Ncells; j++) {
+	    pmCell *cell = pmCellAlloc (chip);
 
-	chip->cells = psArrayAlloc (Ncells);
-	for (int j = 0; j < Ncells; j++) {
-
-	    pmCell *cell = pmCellAlloc ();
-	    cell->chip = chip;      // assign parent chip (view only; don't free)
-	    cell->header = header;  // XXX EAM : extend this function to take more than header
-
+	    // XXX EAM : extend this function to take more than one header
+	    cell->header = header;	
 	    pmCellInterpretWCS (cell, header);
 
-	    cell->readouts = psArrayAlloc (Nreadouts);
+	    psArrayRealloc (cell->readouts, Nreadouts);
 	    for (int k = 0; k < Nreadouts; k++) {
+		pmReadout *readout = pmReadoutAlloc (cell);
 
-		pmReadout *readout = pmReadoutAlloc ();
-		readout->stars  = stars;   // XXX EAM : need more than one stars...
+		// XXX EAM : extend this function to take more than one stars...
+		psMetadataAdd (readout->analysis, PS_LIST_TAIL, "STARS.FULLSET", PS_DATA_ARRAY, "stars from analysis", stars);
 
-		cells->readouts->data[j] = readout;
+		// XXX EAM : this information should be in the header...
+		readout->col0 = readout->row0 = 0;
+		readout->colBins = readout->rowBins = 1;
+
+		cell->readouts->data[j] = readout;
 	    }
-	    chips->cells->data[j] = cell;
+	    chip->cells->data[j] = cell;
 	}
 	fpa->chips->data[i] = chip;
@@ -53,9 +55,13 @@
 bool pmCellInterpretWCS (pmCell *cell, psMetadata *header) { 
 
+    pmChip *chip;
+    pmFPA  *fpa;
+    bool status;
+    psProjectionType type;
     float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2;
     float pc1_1, pc1_2, pc2_1, pc2_2;
 
     // *** interpret header data, convert to crval(i), etc
-    char *ctype = pmMetadataLookupPtr (&status, header, "CTYPE2");
+    char *ctype = psMetadataLookupPtr (&status, header, "CTYPE2");
     if (!status) {
 	psLogMsg ("psastro", 2, "warning: no WCS metadata in header\n");
@@ -85,7 +91,7 @@
 
 	// test the CROTAi varient:
-	float rotate = psMetadataLookupF32 (&status, header, "CROTA2");
+	double rotate = psMetadataLookupF32 (&status, header, "CROTA2");
 	if (status) {
-	    Lambda = cdelt2 / cdelt1;
+	    double Lambda = cdelt2 / cdelt1;
 	    pc1_1 =  cos(rotate*RAD_DEG);
 	    pc1_2 = -sin(rotate*RAD_DEG) * Lambda;
@@ -131,13 +137,15 @@
 
     // XXX EAM : these must already be set
-    chip = cell->chip;
-    fpa = chip->fpa;
+    chip = cell->parent;
+    fpa  = chip->parent;
 
     // set toChip to identity as default
-    // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581)
-    // XXX EAM : define these in the config?
-    // int nX = psMetadataLookupS32 (&status, myHeader, "CHIP.NX");
-    // int nY = psMetadataLookupS32 (&status, myHeader, "CHIP.NY");
-    cell->toChip   = psPlaneTransformAlloc (2, 2);
+    // XXX EAM : can we actually use higher order?
+    int nX = psMetadataLookupS32 (&status, header, "PSASTRO.CHIP.NX");
+    if (!status) nX = 1;
+    int nY = psMetadataLookupS32 (&status, header, "PSASTRO.CHIP.NY");
+    if (!status) nY = 1;
+
+    cell->toChip   = psPlaneTransformAlloc (1, 1);
     cell->toChip->x->coeff[1][0] = 1;
     cell->toChip->x->mask[1][1]  = 1;
@@ -152,5 +160,6 @@
 
     // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581)
-    chip->toFPA   = psPlaneTransformAlloc (2, 2);
+    // XXX EAM : I've fixed this in pslib eam_rel8_b2
+    chip->toFPA   = psPlaneTransformAlloc (1, 1);
     
     chip->toFPA->x->coeff[0][0] = crpix1;
@@ -164,8 +173,10 @@
     chip->toFPA->y->mask[1][1]  = 1;
 
+    chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA);
+
     // set toTPA to identity as default
     // XXX EAM : psPlaneDistortAlloc uses nTerm not nOrder (bug 581)
     if (fpa->toTPA == NULL) {
-	fpa->toTPA   = psPlaneDistortAlloc (2, 2, 1, 1);
+	fpa->toTPA   = psPlaneDistortAlloc (1, 1, 0, 0);
 	fpa->toTPA->x->coeff[1][0][0][0] = 1;
 	fpa->toTPA->x->mask[1][1][0][0]  = 1;
@@ -173,4 +184,5 @@
 	fpa->toTPA->y->coeff[0][1][0][0] = 1;
 	fpa->toTPA->y->mask[1][1][0][0]  = 1;
+	fpa->fromTangentPlane = psPlaneDistortInvert (fpa->toTangentPlane);
     } else {
 	psLogMsg ("psastro", 2, "warning: fpa distortion already defined\n");
@@ -178,5 +190,5 @@
 
     // center of projection is (0,0) coordinate of TPA
-    fpa->toSky = psProjectionAlloc (crval1, crval2, cdelt1, cdelt2, type);
+    fpa->toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type);
     return true;
 }
Index: /trunk/psastro/src/psastroUtils.c
===================================================================
--- /trunk/psastro/src/psastroUtils.c	(revision 5559)
+++ /trunk/psastro/src/psastroUtils.c	(revision 5560)
@@ -1,67 +1,30 @@
 # include "psastro.h"
-
-// crude function to load CMP data (header + ascii)
-// XXX EAM : assumes fixed line, expects NSTARS entries
-psArrray *pmSourcesReadCMP (psMetadata **header, char *filename) {
-
-    psMetadata *myHeader;
-    char line[80];
-
-    psFits *fits = psFitsAlloc (filename);
-    myHeader = psFitsReadHeader (fits);
-    psFree (fits);
-
-    // how many lines in the header?
-    // XXX EAM : is this calculation robust?
-    nLines = myHeader->list->n;
-    nBytes = nLines * 80;
-    if (nBytes % 2880) {
-	nBlock = 1 + (int)(nBytes / 2880);
-	nBytes = nBlock * 2880;
-    }
-
-    // re-open, seek to end of header
-    FILE *f = fopen (filename, "r");
-    if (f == NULL) {
-	psLogMsg ("pmSourcesReadCMP", 3, "can't open output file for input %s\n", filename);
-	return NULL;
-    }
-    fseek (f, nBytes, SEEK_SET);
-
-    // prepare array to store data
-    nStars = psMetadataLookupS32 (&status, myHeader, "NSTARS");
-    psArray *stars = psArrayAlloc (nStars);
-    stars->n = 0;
-
-    // we have fixed bytes / line : use that info
-    for (i = 0; i < nStars; i++) {
-	fread (line, 1, 67, f);
-	sscanf (line, "%lf %lf %lf %lf", &X, &Y, &Mag, &dMag);
-	
-	pmAstromObj *obj = pmAstromObjAlloc ();
-	obj->pix.X    = X;
-	obj->pix.Y    = Y;
-	obj->pix.Mag  = Mag;
-	obj->pix.dMag = dMag;
-	psArrayAdd (stars, obj, 100);
-    }
-    fclose (f);
-
-    *header = myHeader;
-    return (stars);
+# define RENORM 0
+
+bool testWriteRaw (char *filename, psArray *sources);
+
+void psastroDumpStars (psArray *sources) {
+
+    for (int i = 0; i < sources->n; i++) {
+	pmAstromObj *star = sources->data[i];
+
+	fprintf (stderr, "%8.2f %8.2f   %8.2f %8.2f   %8.2f %8.2f   %8.2f %8.2f   %8.2f %8.2f   %10.6f %10.6f   %8.2f %8.2f\n", 
+		 star->pix.x, star->pix.y, 
+		 star->cell.x, star->cell.y, 
+		 star->chip.x, star->chip.y, 
+		 star->FP.x, star->FP.y, 
+		 star->TP.x, star->TP.y, 
+		 star->sky.r*DEG_RAD, star->sky.d*DEG_RAD, 
+		 star->Mag, star->dMag);
+    }
 }
 
 // sort by mag (descending)
-int psastroSortByMag (const void **a, const void **b)
+int pmAstromObjSortByMag (const void **a, const void **b)
 {
-    pmSource *A = *(pmSource **)a;
-    pmSource *B = *(pmSource **)b;
-
-    psF32 fA = (A->moments == NULL) ? 0 : A->moments->SN;
-    psF32 fB = (B->moments == NULL) ? 0 : B->moments->SN;
-    if (isnan (fA)) fA = 0;
-    if (isnan (fB)) fB = 0;
-
-    psF32 diff = fA - fB;
+    pmAstromObj *A = *(pmAstromObj **)a;
+    pmAstromObj *B = *(pmAstromObj **)b;
+
+    psF32 diff = A->Mag - B->Mag;
     if (diff > FLT_EPSILON) return (-1);
     if (diff < FLT_EPSILON) return (+1);
@@ -69,133 +32,10 @@
 }
 
-psArray *psastroSelectBrightStars (pmFPA *inFPA, psMetadata *config) {
-
-    pmFPA *sbFPA = pmFPACopy (fpa);
+bool psastroSelectBrightStars (pmFPA *fpa, psMetadata *config) {
+
+    bool status;
 
     // add exclusions for objects on some basis?
-    int MAX_NSTARS = pmMetadataLookupS32 (&status, config, "MAX_NSTARS");
-
-    for (int i = 0; i < inFPA->chips->n; i++) {
-	pmChip *inChip = inFPA->chips->data[i];
-	pmChip *sbChip = sbFPA->chips->data[i];
-	for (int j = 0; j < Ncells; j++) {
-	    pmCell *inCell = inChip->cells->data[j];
-o	    pmCell *sbCell = sbChip->cells->data[j];
-	    for (int k = 0; k < Nreadouts; k++) {
-		pmReadout *inReadout = inCell->readouts->data[k];
-		pmReadout *sbReadout = sbCell->readouts->data[k];
-
-		inReadout->stars = psArraySort (inReadout->stars, psastroSortByMag);
-
-		nSubset = PS_MIN (MAX_NSTARS, stars->n);
-
-		psArray *subset = psArrayAlloc (nSubset);
-
-		for (int i = 0; i < nSubset; i++) {
-		    subset->data[i] = inReadout->stars->data[i];
-		}
-		sbReadout->stars = subset;
-	    }
-	}
-    }
-    return (subset);
-}
-
-// fake this by loading from a text file
-psArray *psastroLoadReference (char *filename) {
-
-    psMetadata *myHeader;
-    char line[80];
-
-    psFits *fits = psFitsAlloc (filename);
-    myHeader = psFitsReadHeader (fits);
-    psFree (fits);
-
-    // how many lines in the header?
-    // XXX EAM : is this calculation robust?
-    nLines = header->list->n;
-    nBytes = nLines * 80;
-    if (nBytes % 2880) {
-	nBlock = 1 + (int)(nBytes / 2880);
-	nBytes = nBlock * 2880;
-    }
-
-    // re-open, seek to end of header
-    FILE *f = fopen (filename, "r");
-    if (f == NULL) {
-	psLogMsg ("pmSourcesReadCMP", 3, "can't open output file for input %s\n", filename);
-	return NULL;
-    }
-    fseek (f, nBytes, SEEK_SET);
-
-    // prepare array to store data
-    nStars = psMetadataLookupS32 (&status, myHeader, "NSTARS");
-    psArray *stars = psArrayAlloc (nStars);
-    stars->n = 0;
-
-    // we have fixed bytes / line : use that info
-    for (i = 0; i < nStars; i++) {
-	fread (line, 1, 67, f);
-	sscanf (line, "%lf %lf %lf %lf", &X, &Y, &Mag, &dMag);
-	
-	pmAstromObj *obj = pmAstromObjAlloc (X, Y, Mag, dMag);
-	psArrayAdd (stars, obj, 100);
-    }
-    fclose (f);
-
-    *header = myHeader;
-    return (stars);
-}
-
-pmFPA *pmFPACopy (pmFPA *inFPA) {
-
-    pmFPA *fpa = pmFPAAlloc ();
-
-    fpa->toSky   = psMemCopy (inFPA->toSky);
-    fpa->toTPA   = psMemCopy (inFPA->toTPA);
-    fpa->fromTPA = psMemCopy (inFPA->fromTPA);
-
-    fpa->chips = psArrayAlloc (inFPA->chips->n);
-    for (int i = 0; i < inFPA->chips->n; i++) {
-
-	pmChip *inChip = inFPA->chips->data[i];
-
-	pmChip *chip = pmChipAlloc ();
-	chip->fpa = fpa; // assign parent fpa (view only; don't free)
-
-	chip->toFPA = psMemCopy (inChip->toFPA);
-	chip->fromFPA = psMemCopy (inChip->fromFPA);
-
-	chip->cells = psArrayAlloc (inChip->cells->n);
-	for (int j = 0; j < Ncells; j++) {
-
-	    pmCell *inCell = inChip->cells->data[j];
-
-	    pmCell *cell = pmCellAlloc ();
-	    cell->chip = chip;      // assign parent chip (view only; don't free)
-
-	    cell->header = psMemCopy (inCell->header);
-	    cell->toChip = psMemCopy (inCell->toChip);
-
-	    cell->readouts = psArrayAlloc (inCell);
-	    for (int k = 0; k < Nreadouts; k++) {
-
-		pmReadout *inReadout = inCell->readouts->data[k];
-
-		pmReadout *readout = pmReadoutAlloc ();
-
-		*readout = *inReadout;
-		readout->stars = NULL;
-		
-		cell->readouts->data[k] = readout;
-	    }
-	    chip->cells->data[j] = cell;
-	}
-	fpa->chips->data[i] = chip;
-    }
-    return (fpa);
-}
-
-bool psastroProjectFPA (pmFPA *fpa, bool toSky) {
+    int MAX_NSTARS = psMetadataLookupS32 (&status, config, "PSASTRO.STARS.MAX");
 
     for (int i = 0; i < fpa->chips->n; i++) {
@@ -205,9 +45,15 @@
 	    for (int k = 0; k < cell->readouts->n; k++) {
 		pmReadout *readout = cell->readouts->data[k];
-		if (toSky) {
-		    psastroProjectRawstars (readout->stars, readout);
-		} else {
-		    psastroProjectRefstars (readout->stars, readout);
+
+		psArray *stars = psMetadataLookupPtr (&status, readout->analysis, "STARS.FULLSET");
+		stars = psArraySort (stars, pmAstromObjSortByMag);
+
+		int nSubset = PS_MIN (MAX_NSTARS, stars->n);
+		psArray *subset = psArrayAlloc (nSubset);
+
+		for (int i = 0; i < nSubset; i++) {
+		    subset->data[i] = stars->data[i];
 		}
+		psMetadataAdd (readout->analysis, PS_LIST_TAIL, "STARS.SUBSET", PS_DATA_ARRAY, "stars from analysis", subset);
 	    }
 	}
@@ -215,44 +61,60 @@
     return true;
 }
- 
-bool psastroProjectRawstars (psArray *stars, pmReadout *readout) {
-
-    pmCell *cell = readout->cell;
-    pmChip *chip = cell->chip;
-    pmFPA  *fpa  = chip->cell;
-
-    for (int i = 0; i < readout->stars->n; i++) {
-	pmAstromObj *star = readout->stars->data[i];
-	psCoordReadoutToCell (&star->cell, &star->pix, readout);
-	psCoordCellToChip (&star->chip, &star->cell, cell);
-	psCoordChipToFP (&star->FP, &star->chip, chip);
-	psCoordFPtoTP (&star->TP, &star->FP, fpa);
-	psCoordTPtoSky (&star->sky, &star->TP, fpa);
-    }
-    return true;
-}
- 
-bool psastroProjectRefstars (psArray *stars, pmReadout *readout) {
-
-    pmCell *cell = readout->cell;
-    pmChip *chip = cell->chip;
-    pmFPA  *fpa  = chip->cell;
-
-    for (int i = 0; i < stars->n; i++) {
-	pmAstromObj *star = stars->data[i];
-
-	psCoordSkytoTP (&star->TP, &star->sky, fpa);
-	psCoordTPtoFP (&star->FP, &star->TP, fpa);
-	psCoordFPtoChip (&star->chip, &star->FP, chip);
-	psCoordChipToCell (&star->cell, &star->chip, cell);
-	psCoordReadoutToCell (&star->pix, &star->cell, readout);
-    }
-    return true;
-}
- 
+
+pmFPA *pmFPACopyAstrom (pmFPA *inFPA) {
+
+    pmFPA *fpa = pmFPAAlloc (NULL, NULL);
+
+    // copy FPA astrometry data
+    fpa->toSky   = psProjectionCopy (inFPA->toSky);
+    fpa->toTPA   = psPlaneDistortCopy (inFPA->toTPA);
+    fpa->fromTPA = psPlaneDistortCopy (inFPA->fromTPA);
+
+    psArrayRealloc (fpa->chips, inFPA->chips->n);
+    for (int i = 0; i < inFPA->chips->n; i++) {
+	pmChip *inChip = inFPA->chips->data[i];
+	pmChip *chip = pmChipAlloc (fpa);
+
+	// copy Chip astrometry data
+	chip->toFPA = psPlaneTransformCopy (inChip->toFPA);
+	chip->fromFPA = psPlaneTransformCopy (inChip->fromFPA);
+
+	psArrayRealloc (chip->cells, inChip->cells->n);
+	for (int j = 0; j < inChip->cells->n; j++) {
+	    pmCell *inCell = inChip->cells->data[j];
+	    pmCell *cell = pmCellAlloc (chip);
+
+	    // cell.header is a view on inCell.header
+	    cell->header = psMemCopy (inCell->header);
+
+	    // copy Cell astrometry data 
+	    cell->toChip = psPlaneTransformCopy (inCell->toChip);
+
+	    psArrayRealloc (cell->readouts, inCell->readouts->n);
+	    for (int k = 0; k < inCell->readouts->n; k++) {
+		pmReadout *inReadout = inCell->readouts->data[k];
+		pmReadout *readout = pmReadoutAlloc (cell);
+
+		// copy Readout data
+		*readout = *inReadout;
+
+		cell->readouts->data[k] = readout;
+	    }
+	    chip->cells->data[j] = cell;
+	}
+	fpa->chips->data[i] = chip;
+    }
+    return (fpa);
+}
+
 // measure per-chip astrometry terms 
 bool psastroChipAstrom (pmFPA *fpa, psMetadata *config) {
 
-    pmFPA *raw = pmFPACopy (fpa);
+    bool status;
+    psArray *match;
+    pmAstromStats stats;
+
+    // save the raw astrometry for later reference
+    pmFPA *raw = pmFPACopyAstrom (fpa);
 
     // first pass: measure the per-chip solutions, modify the chip.toFPA terms
@@ -282,22 +144,40 @@
 		pmReadout *readout = cell->readouts->data[k];
 
-		// use the header & config info to project refstars on the focal plane
+		// pull out the SUBSET rawstars (a view)
+		psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "STARS.SUBSET");
+
+		// project the rawstars to the current best guess astrometry
+		psastroProjectRawstars (rawstars, readout);
+
+		// use the header & config info to project refstars onto the focal plane
 		psastroProjectRefstars (refstars, readout);
 
+		testWriteRaw ("ref.inp", refstars);
+		testWriteRaw ("raw.inp", rawstars);
+
+		// fprintf (stderr, "rawstars:\n");
+		// psastroDumpStars (rawstars);
+		// fprintf (stderr, "refstars:\n");
+		// psastroDumpStars (refstars);
+
 		// find initial offset / rotation
-		stat = pmAstromGridMatch (readout->stars, refstars, config);
+		stats = pmAstromGridMatch (rawstars, refstars, config);
 
 		// adjust the chip.toFPA terms only
-		pmAstromGridApply (chip->toFPA, stat);
+		pmAstromGridApply (chip->toFPA, stats);
+		chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA);
 
 		// use fit result to re-project rawstars
-		psastroProjectRawstars (readout->stars, readout);
-		psastroProjectRefstars (refstars,       readout);
+		psastroProjectRawstars (rawstars, readout);
+		psastroProjectRefstars (refstars, readout);
+
+		testWriteRaw ("ref.dat", refstars);
+		testWriteRaw ("raw.dat", rawstars);
     
 		// use small radius to match stars
-		match = pmAstromRadiusMatch (rawstars, refstars, options);
+		match = pmAstromRadiusMatch (rawstars, refstars, config);
 
 		// fit astrometric terms
-		pmAstromMatchedListFit (chip->toFPA, readout->stars, refstars, match, options);
+		pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, config);
 	    }
 	}
@@ -306,40 +186,161 @@
     // second stage: re-normalize the chip terms, passing the 
     // average rotation and offset values to the fpa.toSky
-	    
-# if (0)
-    // code disabled for now 
-    // this is not as trivial as it seems at first:
-    // applying a new rotation to the FPA implies changes to the 
-    // (x,y) reference coordinates for the chips as well as a rotation
-    // I need to calculate both the x,y offset effect and the rotation
-    // in one step.  what about the effect of higher order terms in either
-    // the chip.toFPA or the fpa.toTPA? 
-    // this calculation also needs to account for scale effects between the 
-    // fpa and the chips.
-
-    // calculate the average rotation and boresite offset relative to raw
-    for (int i = 0; i < fpa->chips->n; i++) {
-	pmChip *iChip = raw->chips->data[i];
-	pmChip *oChip = fpa->chips->data[i];
-
-	T1 = psPlaneTransformGetRotation (iChip->toFPA);
-	T2 = psPlaneTransformGetRotation (oChip->toFPA);
-	dT += T1 - T2;
-	dN ++;
-
+    if (RENORM) {
+
+	// this code is needed for the mosastro stage, with multiple chip solutions
+
+	double dX, dY, dT, dN;
+	dX = dY = dT = dN = 0;
+
+	psPlane origin, P1, P2;
+	origin.x = 0;
+	origin.y = 0;
+
+	// calculate the average rotation and boresite offset relative to raw
+	for (int i = 0; i < fpa->chips->n; i++) {
+	    pmChip *iChip = raw->chips->data[i];
+	    pmChip *oChip = fpa->chips->data[i];
+
+	    // offset of chip
+	    psCoordChipToFPA (&P1, &origin, iChip);
+	    psCoordChipToFPA (&P2, &origin, oChip);
+	    dX += (P2.x - P1.x);
+	    dY += (P2.y - P1.y);
+
+	    // get parity-independent rotations for old and new solutions
+	    double T1 = psPlaneTransformGetRotation (iChip->toFPA);
+	    double T2 = psPlaneTransformGetRotation (oChip->toFPA);
+	    dT += T2 - T1;
+	    dN ++;
+	}
+
+	dT /= dN;
+	dX /= dN;
+	dY /= dN;
+
+	// R(T)
+	double PC1_1 = fpa->toTPA->x->coeff[1][0][0][0];
+	double PC1_2 = fpa->toTPA->x->coeff[0][1][0][0];
+	double PC2_1 = fpa->toTPA->y->coeff[1][0][0][0];
+	double PC2_2 = fpa->toTPA->y->coeff[0][1][0][0];
+
+	// R(dT)
+	double dPC1_1 = +cos (dT);
+	double dPC1_2 = +sin (dT);
+	double dPC2_1 = -sin (dT);
+	double dPC2_2 = +cos (dT);
+
+	// R'(T) = R(T) * R(dT)
+	double pc1_1 = PC1_1*dPC1_1 + PC1_2*dPC2_1;
+	double pc1_2 = PC1_1*dPC1_2 + PC1_2*dPC2_2;
+	double pc2_1 = PC2_1*dPC1_1 + PC2_2*dPC2_1;
+	double pc2_2 = PC2_1*dPC1_2 + PC2_2*dPC2_2;
+
+	double det = 1.0 / (pc1_1*pc2_2 - pc1_2*pc2_1);
+
+	// R'(-T)  (matrix inverse, not just rotation inverse -- keeps parity)
+	double pi1_1 = +pc2_2 * det;
+	double pi1_2 = -pc1_2 * det;
+	double pi2_1 = -pc2_1 * det;
+	double pi2_2 = +pc1_1 * det;
+    
+	// apply the new modifcations in rotation and boresite
+	for (int i = 0; i < fpa->chips->n; i++) {
+	    pmChip *oChip = fpa->chips->data[i];
+
+	    // r(T)
+	    double pr1_1 = oChip->toFPA->x->coeff[1][0];
+	    double pr1_2 = oChip->toFPA->x->coeff[0][1];
+	    double pr2_1 = oChip->toFPA->y->coeff[1][0];
+	    double pr2_2 = oChip->toFPA->y->coeff[0][1];
+
+	    // ri'(T) = R(T) r(t)
+	    double ri1_1 = PC1_1*pr1_1 + PC1_2*pr2_1;
+	    double ri1_2 = PC1_1*pr1_2 + PC1_2*pr2_2;
+	    double ri2_1 = PC2_1*pr1_1 + PC2_2*pr2_1;
+	    double ri2_2 = PC2_1*pr1_2 + PC2_2*pr2_2;
+
+	    // r'(T) = R'(-T) ri'(T)
+	    oChip->toFPA->x->coeff[1][0] = pi1_1*ri1_1 + pi1_2*ri2_1;
+	    oChip->toFPA->x->coeff[0][1] = pi1_1*ri1_2 + pi1_2*ri2_2;
+	    oChip->toFPA->y->coeff[1][0] = pi2_1*ri1_1 + pi2_2*ri2_1;
+	    oChip->toFPA->y->coeff[0][1] = pi2_1*ri1_2 + pi2_2*ri2_2;
+
+	    double dx = PC1_1*oChip->toFPA->x->coeff[0][0] + PC1_2*oChip->toFPA->y->coeff[0][0] + dX;
+	    double dy = PC2_1*oChip->toFPA->x->coeff[0][0] + PC2_2*oChip->toFPA->y->coeff[0][0] + dY;
+
+	    oChip->toFPA->x->coeff[0][0] = pi1_1*dx + pi1_2*dy;
+	    oChip->toFPA->y->coeff[0][0] = pi2_1*dx + pi2_2*dy;
+	}
+
+	fpa->toTPA->x->coeff[0][0][0][0] -= dX;
+	fpa->toTPA->y->coeff[0][0][0][0] -= dY;
+
+	fpa->toTPA->x->coeff[1][0][0][0] = pc1_1;
+	fpa->toTPA->x->coeff[0][1][0][0] = pc1_2;
+	fpa->toTPA->y->coeff[1][0][0][0] = pc2_1;
+	fpa->toTPA->y->coeff[0][1][0][0] = pc2_2;
+    }
+    return true; 
+}
+
+// returns the rotation term, forcing positive parity
+double psPlaneTransformGetRotation (psPlaneTransform *map) {
+
+    if (map->x->nX < 1) return 0;
+    if (map->x->nY < 1) return 0;
+
+    if (map->y->nX < 1) return 0;
+    if (map->y->nY < 1) return 0;
+    
+    double pc1_1 = map->x->coeff[1][0];
+    double pc1_2 = map->x->coeff[0][1];
+    double pc2_1 = map->y->coeff[1][0];
+    double pc2_2 = map->y->coeff[0][1];
+
+    double px = SIGN (pc1_1);
+    double py = SIGN (pc2_2);
+
+    // both x and y terms imply an angle. take the average
+    double t1 = -atan2 (px*pc1_2, px*pc1_1);
+    double t2 = +atan2 (py*pc2_1, py*pc2_2);
+    
+    // careful near -pi,+pi boundary...
+    if (t1 - t2 > M_PI/2) t2 += 2*M_PI;
+    if (t2 - t1 > M_PI/2) t1 += 2*M_PI;
+
+    double theta = 0.5*(t1 + t2);
+    while (theta < M_PI) theta += 2*M_PI;
+    while (theta > M_PI) theta -= 2*M_PI;
+    
+    return (theta);
+}
+
+// elixir-style pseudo FITS table (header + ascii list)
+bool testWriteRaw (char *filename, psArray *sources) {
+
+    int i;
+
+    // re-open, add data to end of file
+    FILE *f = fopen (filename, "w");
+    if (f == NULL) {
+	psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename);
+	return false;
+    }
+
+    for (i = 0; i < sources->n; i++) {
 	
-    }
-
-    dT /= dN;
-
-    psPlaneTransformSetRotation (fpa->toTPA, dT);
-	
-    // apply the new modifcations in rotation and boresite
-    for (int i = 0; i < fpa->chips->n; i++) {
-	pmChip *iChip = raw->chips->data[i];
-	pmChip *oChip = fpa->chips->data[i];
-
-	psPlaneTransformSetRotation (oChip->toFPA, -dT);
-    }
-# endif 
-}
+	pmAstromObj *star = sources->data[i];
+
+	fprintf (f, "%8.2f %8.2f   %8.2f %8.2f   %8.2f %8.2f   %8.2f %8.2f   %8.2f %8.2f   %10.6f %10.6f   %8.2f %8.2f\n", 
+		 star->pix.x, star->pix.y, 
+		 star->cell.x, star->cell.y, 
+		 star->chip.x, star->chip.y, 
+		 star->FP.x, star->FP.y, 
+		 star->TP.x, star->TP.y, 
+		 star->sky.r*DEG_RAD, star->sky.d*DEG_RAD, 
+		 star->Mag, star->dMag);
+    }
+    fclose (f);
+    return true;
+}