Index: trunk/psastro/src/psastroUtils.c
===================================================================
--- trunk/psastro/src/psastroUtils.c	(revision 5510)
+++ 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;
+}