Index: branches/cnb_branch_20080830/psastro/src/.cvsignore
===================================================================
--- branches/cnb_branch_20080830/psastro/src/.cvsignore	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/.cvsignore	(revision 20033)
@@ -15,2 +15,3 @@
 psastroModel
 gpcModel
+psastroModelFit
Index: branches/cnb_branch_20080830/psastro/src/Makefile.am
===================================================================
--- branches/cnb_branch_20080830/psastro/src/Makefile.am	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/Makefile.am	(revision 20033)
@@ -1,18 +1,23 @@
 
 lib_LTLIBRARIES = libpsastro.la
-libpsastro_la_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSASTRO_CFLAGS)
+libpsastro_la_CFLAGS = $(PSASTRO_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+libpsastro_la_LDFLAGS = $(PSASTRO_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 
-bin_PROGRAMS = psastro psastroModel gpcModel
+bin_PROGRAMS = psastro psastroModel psastroModelFit gpcModel
 
-psastro_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PPSTATS_CFLAGS) $(PSASTRO_CFLAGS)
-psastro_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PPSTATS_LIBS) $(PSASTRO_LIBS)
+psastro_CFLAGS = $(PSASTRO_CFLAGS) $(PPSTATS_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+psastro_LDFLAGS = $(PSASTRO_LIBS) $(PPSTATS_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 psastro_LDADD = libpsastro.la
 
-psastroModel_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSASTRO_CFLAGS)
-psastroModel_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSASTRO_LIBS)
+psastroModel_CFLAGS = $(PSASTRO_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+psastroModel_LDFLAGS = $(PSASTRO_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 psastroModel_LDADD = libpsastro.la
 
-gpcModel_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSASTRO_CFLAGS)
-gpcModel_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSASTRO_LIBS)
+psastroModelFit_CFLAGS = $(PSASTRO_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+psastroModelFit_LDFLAGS = $(PSASTRO_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
+psastroModelFit_LDADD = libpsastro.la
+
+gpcModel_CFLAGS = $(PSASTRO_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+gpcModel_LDFLAGS = $(PSASTRO_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 gpcModel_LDADD = libpsastro.la
 
@@ -38,4 +43,9 @@
 	psastroCleanup.c
 
+psastroModelFit_SOURCES = \
+	psastroModelFit.c \
+	psastroModelBoresite.c      \
+	psastroModelFitBoresite.c
+
 gpcModel_SOURCES = gpcModel.c
 
@@ -50,5 +60,4 @@
 	psastroConvert.c	    \
 	psastroChipAstrom.c         \
-	psastroOneChip.c	    \
 	psastroOneChipGrid.c	    \
 	psastroOneChipFit.c	    \
@@ -57,4 +66,5 @@
 	psastroTestFuncs.c          \
 	psastroLuminosityFunction.c \
+	psastroLuminosityFunctionPlot.c \
 	psastroRefstarSubset.c      \
 	psastroFixChips.c           \
@@ -62,5 +72,8 @@
 	psastroUseModel.c           \
 	psastroMosaicAstrom.c       \
+	psastroMosaicDistortion.c   \
+	psastroMosaicFPtoTP.c       \
 	psastroMosaicGradients.c    \
+	psastroMosaicCorrectDistortion.c   \
 	psastroMosaicChipAstrom.c   \
 	psastroMosaicOneChip.c      \
Index: branches/cnb_branch_20080830/psastro/src/psastro.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastro.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastro.c	(revision 20033)
@@ -9,5 +9,4 @@
 
     pmConfig *config = NULL;
-
     psTimerStart ("complete");
 
@@ -24,6 +23,6 @@
     // load identify the data sources
     if (!psastroParseCamera (config)) {
-	psErrorStackPrint(stderr, "error setting up the camera\n");
-	exit (1);
+        psErrorStackPrint(stderr, "error setting up the camera\n");
+        exit (1);
     }
 
@@ -31,18 +30,18 @@
     // select subset of stars for astrometry
     if (!psastroDataLoad (config)) {
-	psErrorStackPrint(stderr, "error loading input data\n");
-	exit (1);
+        psErrorStackPrint(stderr, "error loading input data\n");
+        exit (1);
     }
 
     // run the full astrometry analysis (chip and/or mosaic)
     if (!psastroAnalysis (config)) {
-	psErrorStackPrint(stderr, "failure in psastro analysis\n");
-	exit (1);
+        psErrorStackPrint(stderr, "failure in psastro analysis\n");
+        exit (1);
     }
-    
+
     // write out the results
     if (!psastroDataSave (config)) {
-	psErrorStackPrint(stderr, "failed to write out data\n");
-	exit (1);
+        psErrorStackPrint(stderr, "failed to write out data\n");
+        exit (1);
     }
 
Index: branches/cnb_branch_20080830/psastro/src/psastro.h
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastro.h	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastro.h	(revision 20033)
@@ -17,11 +17,11 @@
 // logN = offset + slope * logS
 typedef struct {
-    double mMin;			// minimum magnitude bin with data
-    double mMax;			// maximum magnitude bin with data
-    double offset;			// fitted line offset
-    double slope;			// fitted line slope
-    double mPeak;			// mag of peak bin 
-    int nPeak;				// # of stars in peak bin
-    int sPeak;				// sum of stars to peak bin
+    double mMin;                        // minimum magnitude bin with data
+    double mMax;                        // maximum magnitude bin with data
+    double offset;                      // fitted line offset
+    double slope;                       // fitted line slope
+    double mPeak;                       // mag of peak bin
+    int nPeak;                          // # of stars in peak bin
+    int sPeak;                          // sum of stars to peak bin
 } pmLumFunc;
 
@@ -35,4 +35,5 @@
 psArray          *pmSourceToAstromObj (psArray *sources);
 bool              psastroAstromGuess (int *nStars, pmConfig *config);
+bool              psastroAstromGuessCheck (pmConfig *config);
 
 psPlaneDistort   *psPlaneDistortIdentity ();
@@ -46,5 +47,6 @@
 bool              psastroChooseRefstars (pmConfig *config, psArray *refs);
 bool              psastroRefstarSubset (pmReadout *readout);
-pmLumFunc        *psastroLuminosityFunction (psArray *stars);
+pmLumFunc        *psastroLuminosityFunction (psArray *stars, pmLumFunc *rawFunc);
+bool              psastroLuminosityFunctionPlot(psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc);
 psArray          *psastroRemoveClumps (psArray *input, int scale);
 
@@ -55,37 +57,41 @@
 
 // mosaic fitting functions
-bool 		  psastroMosaicGradients (pmFPA *fpa, psMetadata *recipe);
-bool 		  psastroMosaicCommonScale (pmFPA *fpa, psMetadata *recipe);
-bool 		  psastroMosaicAstrom (pmConfig *config);
-bool 		  psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe, int iteration);
-bool 		  psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe, int iteration);
-bool 		  psastroMosaicSetAstrom (pmFPA *fpa);
-bool 		  psastroMosaicHeaders (pmConfig *config);
-bool 		  psastroMosaicRescaleChips (pmFPA *fpa);
-bool 		  psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, int iteration);
+bool              psastroMosaicDistortion (pmFPA *fpa, psMetadata *recipe, int pass);
+psPlaneTransform *psastroMosaicFitRotAndScale (pmFPA *fpa);
+bool              psastroMosaicApplyRotAndScale (pmFPA *fpa, psPlaneTransform *TPtoFP);
+bool              psastroMosaicDistortionFromGradients (pmFPA *fpa, psMetadata *recipe);
+bool              psastroMosaicCorrectDistortion (pmFPA *fpa, psPlaneTransform *TPtoFP);
+bool              psastroMosaicCommonScale (pmFPA *fpa, psMetadata *recipe);
+bool              psastroMosaicAstrom (pmConfig *config);
+bool              psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe, int iteration);
+bool              psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe, int iteration);
+bool              psastroMosaicSetAstrom (pmFPA *fpa);
+bool              psastroMosaicHeaders (pmConfig *config);
+bool              psastroMosaicRescaleChips (pmFPA *fpa);
+bool              psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, int iteration);
 
 // Return version strings.
-psString 	  psastroVersion(void);
-psString 	  psastroVersionLong(void);
+psString          psastroVersion(void);
+psString          psastroVersionLong(void);
 
 // demo plots
-bool 		  psastroPlotRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe);
-bool 		  psastroPlotRefstars (psArray *refstars, psMetadata *recipe);
-bool 		  psastroPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe);
+bool              psastroPlotRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe);
+bool              psastroPlotRefstars (psArray *refstars, psMetadata *recipe);
+bool              psastroPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe);
 
-bool 		  psastroDumpRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip);
-bool 		  psastroDumpRefstars (psArray *refstars, char *filename);
-bool 		  psastroDumpMatches (pmFPA *fpa, char *filename);
-bool 		  psastroDumpStars (psArray *stars, char *filename);
+bool              psastroDumpRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip);
+bool              psastroDumpRefstars (psArray *refstars, char *filename);
+bool              psastroDumpMatches (pmFPA *fpa, char *filename);
+bool              psastroDumpStars (psArray *stars, char *filename);
 bool              psastroDumpMatchedStars (char *filename, psArray *rawstars, psArray *refstars, psArray *match);
 bool              psastroDumpGradients (psArray *gradients, char *filename);
 
-bool		  psastroMosaicSetAstrom_tmp (pmFPA *fpa);
-int 		  psastroSortByMag (const void *a, const void *b);
+bool              psastroMosaicSetAstrom_tmp (pmFPA *fpa);
+int               psastroSortByMag (const void *a, const void *b);
 
 bool              psastroFixChips (pmConfig *config, psMetadata *recipe);
 bool              psastroFixChipsTest (pmConfig *config, psMetadata *recipe);
 bool              psastroUseModel (pmConfig *config, psMetadata *recipe);
-bool              psastroDumpCorners (char *filename, pmFPA *fpa);
+bool              psastroDumpCorners (char *filenameU, char *filenameD, pmFPA *fpa);
 
 
Index: branches/cnb_branch_20080830/psastro/src/psastroAnalysis.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroAnalysis.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroAnalysis.c	(revision 20033)
@@ -63,5 +63,4 @@
         chipastro = true;
     }
-
     if (chipastro) {
         if (!psastroChipAstrom (config)) {
@@ -81,4 +80,6 @@
     // psastroZeroPoint (config);
 
+    psastroAstromGuessCheck (config);
+
     // XXX how do we specify stack astrometry?
     // psastroStackAstrom (config, refs);
Index: branches/cnb_branch_20080830/psastro/src/psastroArguments.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroArguments.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroArguments.c	(revision 20033)
@@ -81,4 +81,11 @@
     }
 
+    // dump the configuration to a file?
+    if ((N = psArgumentGet (argc, argv, "-dumpconfig"))) {
+        psArgumentRemove (N, &argc, argv);
+        psMetadataAddStr (config->arguments, PS_LIST_TAIL, "DUMP_CONFIG", PS_META_REPLACE, "", argv[N]);
+        psArgumentRemove (N, &argc, argv);
+    }
+
     status = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list");
     if (!status) {
Index: branches/cnb_branch_20080830/psastro/src/psastroAstromGuess.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroAstromGuess.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroAstromGuess.c	(revision 20033)
@@ -1,3 +1,4 @@
 # include "psastroInternal.h"
+# define DEBUG 0
 
 // this function loads the header WCS astrometry terms into the fpa terms and applies the
@@ -28,6 +29,6 @@
     psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE);
     if (!recipe) {
-	psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!");
-	return false;
+        psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!");
+        return false;
     }
 
@@ -35,5 +36,5 @@
     bool useModel = psMetadataLookupBool (&status, config->arguments, "PSASTRO.USE.MODEL");
     if (!status) {
-	useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL");
+        useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL");
     }
 
@@ -41,6 +42,6 @@
     pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
     if (!input) {
-	psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
-	return false;
+        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
+        return false;
     }
 
@@ -48,14 +49,23 @@
     double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
     if (!status) {
-	psError(PS_ERR_IO, true, "Failed to lookup pixel scale"); 
-	return false; 
-    } 
+        psError(PS_ERR_IO, true, "Failed to lookup pixel scale");
+        return false;
+    }
+
+    psVector *cornerL = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerM = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerP = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerQ = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerR = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerD = psVectorAllocEmpty (100, PS_TYPE_F32);
 
     pmFPA *fpa = input->fpa;
+
+    if (DEBUG) psastroDumpCorners ("corners.up.guess1.dat", "corners.dn.guess1.dat", fpa);
 
     // load mosaic-level astrometry?
     bool bilevelAstrometry = false;
     if (!useModel) {
-	psastroAstromGuessSetFPA (fpa, &bilevelAstrometry);
+        psastroAstromGuessSetFPA (fpa, &bilevelAstrometry);
     }
 
@@ -65,14 +75,34 @@
         if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; }
 
-	if (!useModel) {
-	    if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue;
-	}
+        if (!useModel) {
+            if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue;
+        }
 
         if (newFPA) {
             newFPA = false;
-	    while (fpa->toSky->R <        0) fpa->toSky->R += 2.0*M_PI;
-	    while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI;
+            while (fpa->toSky->R <        0) fpa->toSky->R += 2.0*M_PI;
+            while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI;
             RAminSky = fpa->toSky->R - M_PI;
             RAmaxSky = fpa->toSky->R + M_PI;
+        }
+
+        // report and save the current best guess for the chip 0,0 pixel coordinates
+        {
+            psPlane ptCH, ptFP, ptTP;
+            psSphere ptSky;
+
+            ptCH.x = 0;
+            ptCH.y = 0;
+            psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
+            psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
+            psDeproject (&ptSky, &ptTP, fpa->toSky);
+            psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
+
+            psVectorAppend (cornerL, ptFP.x);
+            psVectorAppend (cornerM, ptFP.y);
+            psVectorAppend (cornerP, ptTP.x);
+            psVectorAppend (cornerQ, ptTP.y);
+            psVectorAppend (cornerR, ptSky.r);
+            psVectorAppend (cornerD, ptSky.d);
         }
 
@@ -86,21 +116,8 @@
                 if (! readout->data_exists) { continue; }
 
-		// report the current best guess for the cell 0,0 pixel coordinate
-		{ 
-		  psPlane ptCH, ptFP, ptTP;
-		  psSphere ptSky;
-
-		  ptCH.x = 0;
-		  ptCH.y = 0;
-		  psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
-		  psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
-		  psDeproject (&ptSky, &ptTP, fpa->toSky);
-		  psLogMsg ("psastro", 2, "0,0 pix for chip,cell %d,%d = %f,%f\n", view->chip, view->cell, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
-		}
-
                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
                 if (rawstars == NULL) { continue; }
 
-		*nStars += rawstars->n;
+                *nStars += rawstars->n;
                 for (int i = 0; i < rawstars->n; i++) {
                     pmAstromObj *raw = rawstars->data[i];
@@ -121,27 +138,29 @@
                 }
 
-		// dump or plot the resulting projected positions
-		if (psTraceGetLevel("psastro.dump") > 0) {
-		    psastroDumpRawstars (rawstars, fpa, chip);
-		}
-
-		if (psTraceGetLevel("psastro.plot") > 0) {
-		    psastroPlotRawstars (rawstars, fpa, chip, recipe);
-		}
+                // dump or plot the resulting projected positions
+                if (psTraceGetLevel("psastro.dump") > 0) {
+                    psastroDumpRawstars (rawstars, fpa, chip);
+                }
+
+                if (psTraceGetLevel("psastro.plot") > 0) {
+                    psastroPlotRawstars (rawstars, fpa, chip, recipe);
+                }
             }
         }
     }
+
+    if (DEBUG) psastroDumpCorners ("corners.up.guess2.dat", "corners.dn.guess2.dat", fpa);
 
     // how many total sources are available to us?
     psMetadataAddS32 (recipe, PS_LIST_TAIL, "NTOTSTAR",  PS_META_REPLACE, "", *nStars);
     if (*nStars == 0) {
-	psLogMsg ("psastro", 2, "no sources available for astrometry\n");
-	psFree (view);
-	return true;
-    }
-
-    psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", 
-	      DEG_RAD*RAmin, DEG_RAD*DECmin, 
-	      DEG_RAD*RAmax, DEG_RAD*DECmax);
+        psLogMsg ("psastro", 2, "no sources available for astrometry\n");
+        psFree (view);
+        return true;
+    }
+
+    psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n",
+              DEG_RAD*RAmin, DEG_RAD*DECmin,
+              DEG_RAD*RAmax, DEG_RAD*DECmax);
 
     psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MIN",  PS_META_REPLACE, "", RAmin);
@@ -149,4 +168,18 @@
     psMetadataAddF32 (recipe, PS_LIST_TAIL, "DEC_MIN", PS_META_REPLACE, "", DECmin);
     psMetadataAddF32 (recipe, PS_LIST_TAIL, "DEC_MAX", PS_META_REPLACE, "", DECmax);
+
+    psMetadataAddVector (input->fpa->analysis, PS_LIST_TAIL, "CORNER.L", PS_META_REPLACE, "corner pixel", cornerL);
+    psMetadataAddVector (input->fpa->analysis, PS_LIST_TAIL, "CORNER.M", PS_META_REPLACE, "corner pixel", cornerM);
+    psMetadataAddVector (input->fpa->analysis, PS_LIST_TAIL, "CORNER.P", PS_META_REPLACE, "corner pixel", cornerP);
+    psMetadataAddVector (input->fpa->analysis, PS_LIST_TAIL, "CORNER.Q", PS_META_REPLACE, "corner pixel", cornerQ);
+    psMetadataAddVector (input->fpa->analysis, PS_LIST_TAIL, "CORNER.R", PS_META_REPLACE, "corner pixel", cornerR);
+    psMetadataAddVector (input->fpa->analysis, PS_LIST_TAIL, "CORNER.D", PS_META_REPLACE, "corner pixel", cornerD);
+
+    psFree (cornerL);
+    psFree (cornerM);
+    psFree (cornerP);
+    psFree (cornerQ);
+    psFree (cornerR);
+    psFree (cornerD);
 
     psFree (view);
@@ -168,13 +201,13 @@
     pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
     if (bilevelAstrometry) {
-	if (!pmAstromReadBilevelChip (chip, hdu->header)) {
-	    psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); 
-	    return false;
-	} 
+        if (!pmAstromReadBilevelChip (chip, hdu->header)) {
+            psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
+            return false;
+        }
     } else {
-	if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) {
-	    psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); 
-	    return false;
-	} 
+        if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) {
+            psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
+            return false;
+        }
     }
     return true;
@@ -190,13 +223,164 @@
     // load mosaic-level astrometry?
     if (phu) {
-	char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
-	if (ctype) {
-	    *bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
-	}
+        char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
+        if (ctype) {
+            *bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
+        }
     }
     if (*bilevelAstrometry) {
-	pmAstromReadBilevelMosaic (fpa, phu->header);
-    } 
+        pmAstromReadBilevelMosaic (fpa, phu->header);
+    }
     psFree (view);
     return true;
 }
+
+// we made a guess at the beginning; how does the guess compare with the result?
+bool psastroAstromGuessCheck (pmConfig *config) {
+
+    bool status;
+
+    // select the input data sources
+    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
+    if (!input) {
+        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
+        return false;
+    }
+
+    pmFPA *fpa = input->fpa;
+
+    psVector *cornerLo = psMetadataLookupPtr (&status, input->fpa->analysis, "CORNER.L");
+    psVector *cornerMo = psMetadataLookupPtr (&status, input->fpa->analysis, "CORNER.M");
+    psVector *cornerPo = psMetadataLookupPtr (&status, input->fpa->analysis, "CORNER.P");
+    psVector *cornerQo = psMetadataLookupPtr (&status, input->fpa->analysis, "CORNER.Q");
+    psVector *cornerRo = psMetadataLookupPtr (&status, input->fpa->analysis, "CORNER.R");
+    psVector *cornerDo = psMetadataLookupPtr (&status, input->fpa->analysis, "CORNER.D");
+
+    if (cornerLo->n < 3) return true;
+
+    psVector *cornerLn = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerMn = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerPn = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerQn = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerRn = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerDn = psVectorAllocEmpty (100, PS_TYPE_F32);
+
+    if (DEBUG) psastroDumpCorners ("corners.up.guess3.dat", "corners.dn.guess3.dat", fpa);
+
+    pmChip *chip = NULL;
+    pmFPAview *view = pmFPAviewAlloc (0);
+
+    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
+        if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; }
+
+        psPlane ptCH, ptFP, ptTP;
+        psSphere ptSky;
+
+        ptCH.x = 0;
+        ptCH.y = 0;
+        psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
+        psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
+        psDeproject (&ptSky, &ptTP, fpa->toSky);
+        psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
+
+        // new corner locations based on the calibrated astrometry
+        psVectorAppend (cornerLn, ptFP.x);
+        psVectorAppend (cornerMn, ptFP.y);
+        psVectorAppend (cornerPn, ptTP.x);
+        psVectorAppend (cornerQn, ptTP.y);
+        psVectorAppend (cornerRn, ptSky.r);
+        psVectorAppend (cornerDn, ptSky.d);
+    }
+
+    // compare the old R,D values projected to the same tangent plane as the new R,D values:
+
+    psVector *cornerPs = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerQs = psVectorAllocEmpty (100, PS_TYPE_F32);
+
+    for (int i = 0; i < cornerRo->n; i++) {
+
+        psPlane ptTP;
+        psSphere ptSky;
+
+        ptSky.r = cornerRo->data.F32[i];
+        ptSky.d = cornerDo->data.F32[i];
+
+        psProject (&ptTP, &ptSky, fpa->toSky);
+        psVectorAppend (cornerPs, ptTP.x);
+        psVectorAppend (cornerQs, ptTP.y);
+    }
+
+    psPlaneTransform *map = psPlaneTransformAlloc (1, 1);
+    map->x->coeffMask[1][1] = PS_POLY_MASK_SET;
+    map->y->coeffMask[1][1] = PS_POLY_MASK_SET;
+
+    psVectorFitPolynomial2D (map->x, NULL, 0, cornerPn, NULL, cornerPs, cornerQs);
+    psVectorFitPolynomial2D (map->y, NULL, 0, cornerQn, NULL, cornerPs, cornerQs);
+
+    // apply the linear fit...
+    psVector *cornerPf = psPolynomial2DEvalVector (map->x, cornerPs, cornerQs);
+    psVector *cornerQf = psPolynomial2DEvalVector (map->y, cornerPs, cornerQs);
+
+    // ...and calculate the residual between Pn,Qn and Pf,Qf
+    psVector *cornerPd = (psVector *) psBinaryOp (NULL, cornerPn, "-", cornerPf);
+    psVector *cornerQd = (psVector *) psBinaryOp (NULL, cornerQn, "-", cornerQf);
+
+    psStats *statsP = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+    psStats *statsQ = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+
+    psVectorStats (statsP, cornerPd, NULL, NULL, 0);
+    psVectorStats (statsQ, cornerQd, NULL, NULL, 0);
+
+    float angle = atan2 (map->y->coeff[1][0], map->x->coeff[1][0]);
+    float scale = hypot (map->y->coeff[1][0], map->x->coeff[1][0]);
+
+    psLogMsg ("psastro", 3, "boresite offset  : %f,%f\n", map->x->coeff[0][0], map->y->coeff[0][0]);
+    psLogMsg ("psastro", 3, "boresite angle   : %f, scale: %f", angle*PS_DEG_RAD, scale);
+    psLogMsg ("psastro", 3, "boresite scatter : %f,%f\n", statsP->sampleStdev, statsQ->sampleStdev);
+
+    // write the elapsed time here; this will be updated in psastroMosaicAstrometry, if called
+    psMetadata *header = psMetadataLookupMetadata (&status, input->fpa->analysis, "PSASTRO.HEADER");
+    if (!header) {
+        header = psMetadataAlloc();
+        psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", header);
+        psFree (header);  // drop this reference
+    }
+
+    psMetadataAddF32 (header, PS_LIST_TAIL, "AST_R0", PS_META_REPLACE, "boresite offset in RA (TP units)", map->x->coeff[0][0]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "AST_D0", PS_META_REPLACE, "boresite offset in DEC (TP units)", map->y->coeff[0][0]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "AST_T0", PS_META_REPLACE, "boresite angle (degrees)", angle*PS_DEG_RAD);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "AST_S0", PS_META_REPLACE, "boresite scale correction", scale);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "AST_RS", PS_META_REPLACE, "boresite scatter in RA (TP units)", statsP->sampleStdev);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "AST_DS", PS_META_REPLACE, "boresite scatter in DEC (TP units)", statsQ->sampleStdev);
+
+    if (0) {
+        FILE *f = fopen ("corners.dat", "w");
+        for (int i = 0; i < cornerRo->n; i++) {
+            fprintf (f, "%10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f  |  %10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f\n",
+                     cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i],
+                     cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]);
+        }
+        fclose (f);
+    }
+
+    psFree (cornerPf);
+    psFree (cornerQf);
+    psFree (cornerPd);
+    psFree (cornerQd);
+
+    psFree (statsP);
+    psFree (statsQ);
+
+    psFree (cornerLn);
+    psFree (cornerMn);
+    psFree (cornerPn);
+    psFree (cornerQn);
+    psFree (cornerRn);
+    psFree (cornerDn);
+    psFree (cornerPs);
+    psFree (cornerQs);
+    psFree (map);
+    psFree (view);
+
+
+    return true;
+}
Index: branches/cnb_branch_20080830/psastro/src/psastroChipAstrom.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroChipAstrom.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroChipAstrom.c	(revision 20033)
@@ -115,19 +115,8 @@
     }
 
-    // psastroDumpCorners ("corners.chipAstrom.dat", input->fpa);
-
-# if (0)
-    if (!psastroFixChipsTest (config, recipe)) {
-        psError(PSASTRO_ERR_UNKNOWN, false, "failed to align problematic chips");
-        return false;
-    }
-# endif
-
     if (!psastroFixChips (config, recipe)) {
         psError(PSASTRO_ERR_UNKNOWN, false, "failed to align problematic chips");
         return false;
     }
-
-    // psastroDumpCorners ("corners.fixChips.dat", input->fpa);
 
     psFree (view);
Index: branches/cnb_branch_20080830/psastro/src/psastroChooseRefstars.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroChooseRefstars.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroChooseRefstars.c	(revision 20033)
@@ -51,5 +51,5 @@
         psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
         if (!chip->process || !chip->file_exists) { continue; }
-	if (!chip->fromFPA) { continue; }
+        if (!chip->fromFPA) { continue; }
 
         while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
@@ -62,12 +62,12 @@
                 if (! readout->data_exists) { continue; }
 
-		psRegion *extent = pmReadoutExtent (readout);
-		if (!extent) {
-		    psError(PSASTRO_ERR_CONFIG, true, "Can't find readout size!\n");
-		    return NULL;
-		}
+                psRegion *extent = pmReadoutExtent (readout);
+                if (!extent) {
+                    psError(PSASTRO_ERR_CONFIG, true, "Can't find readout size!\n");
+                    return NULL;
+                }
 
-		int Nx = abs(extent->x1 - extent->x0);
-		int Ny = abs(extent->y1 - extent->y0);
+                int Nx = abs(extent->x1 - extent->x0);
+                int Ny = abs(extent->y1 - extent->y0);
 
                 float minX = -fieldPadding*Nx;
@@ -98,24 +98,24 @@
                     psFree (ref);
 
-		    if (nMax && (refstars->n >= nMax)) break;
+                    if (nMax && (refstars->n >= nMax)) break;
                 }
                 psTrace ("psastro", 4, "Added %ld refstars\n", refstars->n);
 
-		psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars);
-		psFree (refstars);
-		psFree (extent);
+                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars);
+                psFree (refstars);
+                psFree (extent);
 
-		if (matchLumFunc) {
-		    // in this case, no PSASTRO.REFSTARS is added to readout->analysis
-		    if (!psastroRefstarSubset (readout)) {
-			psError(PSASTRO_ERR_DATA, false, "Can't determine an appropriate refstar subset\n");
-			psFree (index);
-			psFree (view);
-			return false;
-		    }
-		}
+                if (matchLumFunc) {
+                    // in this case, no PSASTRO.REFSTARS is added to readout->analysis
+                    if (!psastroRefstarSubset (readout)) {
+                        psError(PSASTRO_ERR_DATA, false, "Can't determine an appropriate refstar subset\n");
+                        psFree (index);
+                        psFree (view);
+                        return false;
+                    }
+                }
             }
         }
-	if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
+        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
     }
     if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
Index: branches/cnb_branch_20080830/psastro/src/psastroDemoDump.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroDemoDump.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroDemoDump.c	(revision 20033)
@@ -128,4 +128,6 @@
         if (!chip->process || !chip->file_exists) continue;
 	
+	char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
+
 	while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
             psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
@@ -153,6 +155,6 @@
 
 		    pmAstromObj *raw = rawstars->data[match->raw];
-		    fprintf (f, "%f %f  %f %f  %f %f  %f %f  %f   |   ",  
-			     DEG_RAD*raw->sky->r, DEG_RAD*raw->sky->d, 
+		    fprintf (f, "%s  %f %f  %f %f  %f %f  %f %f  %f   |   ",  
+			     chipName, DEG_RAD*raw->sky->r, DEG_RAD*raw->sky->d, 
 			     raw->TP->x, raw->TP->y, 
 			     raw->FP->x, raw->FP->y, 
@@ -193,10 +195,16 @@
 }
 
-bool psastroDumpCorners (char *filename, pmFPA *fpa) {
+bool psastroDumpCorners (char *filenameU, char *filenameD, pmFPA *fpa) {
 
   // XXX test output of chip corners based on model
-  FILE *f = fopen (filename, "w");
+  FILE *fu = fopen (filenameU, "w");
+  FILE *fd = fopen (filenameD, "w");
 
   pmFPAview *view = pmFPAviewAlloc (0);
+
+  float fpaAngle = PM_DEG_RAD * atan2 (fpa->toTPA->y->coeff[1][0], fpa->toTPA->x->coeff[1][0]);
+
+  fprintf (fu, "# boresite: %f, %f @ %f\n", fpa->toSky->R*PS_DEG_RAD, fpa->toSky->D*PS_DEG_RAD, fpaAngle);
+  fprintf (fd, "# boresite: %f, %f @ %f\n", fpa->toSky->R*PS_DEG_RAD, fpa->toSky->D*PS_DEG_RAD, fpaAngle);
 
   pmChip *chip = NULL;
@@ -209,38 +217,80 @@
 	psSphere ptSky;
 
+	// UP 0,0
 	ptCP.x = region->x0; ptCP.y = region->y0;
 	psPlaneTransformApply (&ptFP, chip->toFPA, &ptCP);
 	psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
 	psDeproject (&ptSky, &ptTP, fpa->toSky);
-	fprintf (f, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
-
+	fprintf (fu, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+
+	// DOWN 0,0
+	psProject (&ptTP, &ptSky, fpa->toSky);
+	psPlaneTransformApply (&ptFP, fpa->fromTPA, &ptTP);
+	psPlaneTransformApply (&ptCP, chip->fromFPA, &ptFP);
+	fprintf (fd, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+
+	// UP 1,0
 	ptCP.x = region->x1; ptCP.y = region->y0;
 	psPlaneTransformApply (&ptFP, chip->toFPA, &ptCP);
 	psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
 	psDeproject (&ptSky, &ptTP, fpa->toSky);
-	fprintf (f, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
-
+	fprintf (fu, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+	fprintf (fu, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+
+	// DOWN 1,0
+	psProject (&ptTP, &ptSky, fpa->toSky);
+	psPlaneTransformApply (&ptFP, fpa->fromTPA, &ptTP);
+	psPlaneTransformApply (&ptCP, chip->fromFPA, &ptFP);
+	fprintf (fd, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+	fprintf (fd, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+
+	// UP 1,1
 	ptCP.x = region->x1; ptCP.y = region->y1;
 	psPlaneTransformApply (&ptFP, chip->toFPA, &ptCP);
 	psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
 	psDeproject (&ptSky, &ptTP, fpa->toSky);
-	fprintf (f, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
-
+	fprintf (fu, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+	fprintf (fu, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+
+	// DOWN 1,1
+	psProject (&ptTP, &ptSky, fpa->toSky);
+	psPlaneTransformApply (&ptFP, fpa->fromTPA, &ptTP);
+	psPlaneTransformApply (&ptCP, chip->fromFPA, &ptFP);
+	fprintf (fd, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+	fprintf (fd, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+
+	// UP 0,1
 	ptCP.x = region->x0; ptCP.y = region->y1;
 	psPlaneTransformApply (&ptFP, chip->toFPA, &ptCP);
 	psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
 	psDeproject (&ptSky, &ptTP, fpa->toSky);
-	fprintf (f, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
-
+	fprintf (fu, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+	fprintf (fu, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+
+	// DOWN 0,1
+	psProject (&ptTP, &ptSky, fpa->toSky);
+	psPlaneTransformApply (&ptFP, fpa->fromTPA, &ptTP);
+	psPlaneTransformApply (&ptCP, chip->fromFPA, &ptFP);
+	fprintf (fd, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+	fprintf (fd, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+
+	// UP 0,0
 	ptCP.x = region->x0; ptCP.y = region->y0;
 	psPlaneTransformApply (&ptFP, chip->toFPA, &ptCP);
 	psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
 	psDeproject (&ptSky, &ptTP, fpa->toSky);
-	fprintf (f, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+	fprintf (fu, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
+
+	// DOWN 0,0
+	psProject (&ptTP, &ptSky, fpa->toSky);
+	psPlaneTransformApply (&ptFP, fpa->fromTPA, &ptTP);
+	psPlaneTransformApply (&ptCP, chip->fromFPA, &ptFP);
+	fprintf (fd, "%10.6f %10.6f  %8.1f %8.1f  %8.1f %8.1f  %8.1f %8.1f\n", ptSky.r, ptSky.d, ptTP.x, ptTP.y, ptFP.x, ptFP.y, ptCP.x, ptCP.y);
 
 	psFree (region);
   }
 
-  fclose (f);
+  fclose (fu);
+  fclose (fd);
   psFree (view);
   return true;
Index: branches/cnb_branch_20080830/psastro/src/psastroDemoPlot.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroDemoPlot.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroDemoPlot.c	(revision 20033)
@@ -101,4 +101,5 @@
     }
     xVec->n = yVec->n = zVec->n = n;
+
     pmKapaPlotVectorTriple_AutoLimits_OpenGraph (kapa, &graphdata, xVec, yVec, zVec, false);
 
@@ -124,4 +125,8 @@
     xVec->n = yVec->n = zVec->n = n;
     pmKapaPlotVectorTriple_AutoLimits_OpenGraph (kapa, &graphdata, xVec, yVec, zVec, false);
+
+    // flip x (East increase to left)
+    SWAP (graphdata.xmin, graphdata.xmax);
+    KapaSetLimits (kapa, &graphdata);
 
     // pause and wait for user input:
@@ -189,4 +194,8 @@
     pmKapaPlotVectorTriple_AutoLimits_OpenGraph (kapa, &graphdata, xVec, yVec, zVec, false);
 
+    // flip x (East increase to left)
+    SWAP (graphdata.xmin, graphdata.xmax);
+    KapaSetLimits (kapa, &graphdata);
+
     // pause and wait for user input:
     // continue, save (provide name), ??
Index: branches/cnb_branch_20080830/psastro/src/psastroEnforceChips.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroEnforceChips.c	(revision 20026)
+++ 	(revision )
@@ -1,152 +1,0 @@
-# include "psastroInternal.h"
-# define NONLIN_TOL 0.001 /* tolerance in pixels */
-
-XXX add a function to supply for all entries, as opposed to fixing the poor astrometry ones.
-bool psastroEnforceChips (pmConfig *config, psMetadata *recipe) {
-
-  bool status;
-
-  bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS");
-  if (!status) {
-      fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS");
-  }
-  if (!fixChips) return true;
-
-  // identify reference astrometry table
-  // if not defined, correction was not requested; skip step
-  pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.MODEL");
-  if (!astrom) return true;
-
-  // select the input data sources
-  pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
-  if (!input) {
-      psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
-      return false;
-  }
-
-  // make sure the astrometry model is loaded
-  // de-activate all files except PSASTRO.MODEL
-  // XXX do I need to supply the input->fpa->concepts to the astrom->fpa->concepts?
-  psFree (astrom->fpa->concepts);
-  astrom->fpa->concepts = psMemIncrRefCounter (input->fpa->concepts);
-  pmFPAfileActivate (config->files, false, NULL);
-  pmFPAfileActivate (config->files, true, "PSASTRO.MODEL");
-
-  pmFPAview *view = pmFPAviewAlloc (0);
-
-  // files associated with the science image
-  if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
-      psError (PS_ERR_IO, false, "Can't load the astrometry model file");
-      return false;
-  }
-
-  // loop over all chips, replace input astrometry elements with those from astrom 
-  pmChip *obsChip = NULL;
-  while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
-    psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process);
-    if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
-
-    // set the chip astrometry using the astrom file
-    pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
-
-    psFree (obsChip->toFPA);
-    psFree (obsChip->fromFPA);
-
-    // XXX can we do a straight copy?  I think so, but we may need to apply a scaling factor based
-    // on the toSky or toTPA values
-    obsChip->toFPA   = psMemIncrRefCounter (refChip->toFPA);
-    obsChip->fromFPA = psMemIncrRefCounter (refChip->fromFPA);
-
-    // XXX test to write out adjusted header
-    pmAstromWriteBilevelChip (obsChip->hdu->header, obsChip, NONLIN_TOL);
-  }
-
-  psFree (input->fpa->toSky);
-  psFree (input->fpa->toTPA);
-  psFree (input->fpa->fromTPA);
-  input->fpa->toSky   = psMemIncrRefCounter (astrom->fpa->toSky);
-  input->fpa->toTPA   = psMemIncrRefCounter (astrom->fpa->toTPA);
-  input->fpa->fromTPA = psMemIncrRefCounter (astrom->fpa->fromTPA);
-
-  psMetadata *updates = psMetadataAlloc();
-  pmAstromWriteBilevelMosaic (updates, input->fpa, NONLIN_TOL);
-  psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
-  psFree (updates);
-
-  psFree (view);
-  return true;
-}
-
-# if (0)
-
-  // physical pixel scale in microns per pixel
-  double pixelTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.TOLERANCE");
-  if (!status) {
-    psError(PS_ERR_IO, true, "Failed to lookup pixel tolerance"); 
-    psFree (view);
-    return false; 
-  } 
-  double angleTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.ANGLE.TOLERANCE");
-  if (!status) {
-    psError(PS_ERR_IO, true, "Failed to lookup angle tolerance"); 
-    psFree (view);
-    return false; 
-  } 
-
-  pmChip *obsChip = NULL;
-
-  // loop over all chips
-  while ((obsChip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
-    psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process);
-    if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
-
-    // set the chip astrometry using the astrom file
-    pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
-
-    // bad Astrometry test:  ref pixel or angle outside nominal
-
-    psPlane refPixel = {0.0, 0.0, 0.0, 0.0};
-    psPlane obsCoord, refCoord;
-
-    // find location of 0,0 pixel in focal plane coords for this chip
-    psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel);
-
-    // find location of 0,0 pixel in focal plane coords for ref chip
-    psPlaneTransformApply (&refCoord, refChip->toFPA, &refPixel);
-    
-    psPlane offPixel = {0.0, 0.0, 0.0, 0.0};
-    psPlane obsOffPt, refOffPt;
-
-    // find location of 0,0 pixel in focal plane coords for this chip
-    psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel);
-
-    // find location of 0,0 pixel in focal plane coords for ref chip
-    psPlaneTransformApply (&refOffPt, refChip->toFPA, &offPixel);
-    
-    double obsAngle = atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x);
-    double refAngle = atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x);
-
-    bool badAstrom = false;
-    badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol;
-    badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol;
-    badAstrom |= fabs(obsAngle   - refAngle)   > angleTol;
-
-    if (!badAstrom) continue;
-
-    psFree (obsChip->toFPA);
-    psFree (obsChip->fromFPA);
-
-    // XXX can we do a straight copy?  I think so, but we may need to apply a scaling factor based
-    // on the toSky or toTPA values
-    obsChip->toFPA   = psMemIncrRefCounter (refChip->toFPA);
-    obsChip->fromFPA = psMemIncrRefCounter (refChip->fromFPA);
-  }
-  psFree (view);
-  return true;
-}
-
-// XXX for this function to work, I need to:
-// 1 *) define badAstrom (ref pixel too far from nominal?) (ref angle too far off?)
-// 2) load the astrom fpa in the pmFPAfileIO stage 
-// 3) allow for this operation to be optional
-# endif
Index: branches/cnb_branch_20080830/psastro/src/psastroModelAdjust.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroModelAdjust.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroModelAdjust.c	(revision 20033)
@@ -1,4 +1,7 @@
 # include "psastroStandAlone.h"
 # define NONLIN_TOL 0.001 /* tolerance in pixels */
+# define DEBUG 0
+
+bool psastroModelAdjustBoresite (pmFPAfile *output, pmChip *refChip);
 
 bool psastroModelAdjust (pmConfig *config) {
@@ -13,4 +16,18 @@
     }
 
+    // if we have not measured the boresite position, no adjustment is needed
+    bool fitBoresite = psMetadataLookupBool (&status, recipe, "PSASTRO.MODEL.FIT.BORESITE");
+    if (!status) psAbort ("Can't find recipe option PSASTRO.MODEL.FIT.BORESITE");
+
+    // as an alternative to fit the boresite from a rotation sequence, we can set the boresite
+    // relative to the reference chip coordinates
+    bool setBoresite = psMetadataLookupBool (&status, recipe, "PSASTRO.MODEL.SET.BORESITE");
+    if (!status) psAbort ("Can't find recipe option PSASTRO.MODEL.SET.BORESITE");
+
+    if (fitBoresite && setBoresite) {
+	psError(PS_ERR_IO, true, "invalid to choose both FIT.BORESITE and SET.BORESITE"); 
+	return false; 
+    } 
+	
     pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSASTRO.OUT.MODEL");
     if (!status) psAbort ("Can't find output pmFPAfile PSASTRO.OUT.MODEL");
@@ -29,7 +46,91 @@
     if (!refChip->toFPA) psAbort ("invalid astrometry for reference chip");
 
-    psPlane *boreCH = psPlaneAlloc();
-    psPlane *boreFP = psPlaneAlloc();    
-    psPlane *boreTP = psPlaneAlloc();    
+    // save the TPA region for tranformation inversions below
+    // psRegion *tpaRegion = pmAstromFPInTP (output->fpa);
+    psRegion *fpaRegion = pmAstromFPAExtent (output->fpa);
+
+    if (DEBUG) psastroDumpCorners ("corners.up.raw.dat", "corners.dn.raw.dat", output->fpa);
+
+    if (setBoresite) {
+	float boreXchip = psMetadataLookupF32 (&status, recipe, "PSASTRO.MODEL.BORESITE.X");
+	float boreYchip = psMetadataLookupF32 (&status, recipe, "PSASTRO.MODEL.BORESITE.Y");
+	psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.X0", PS_META_REPLACE, "boresite parameter", boreXchip); 
+	psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.Y0", PS_META_REPLACE, "boresite parameter", boreYchip); 
+    }
+
+    if (fitBoresite || setBoresite) {
+	psastroModelAdjustBoresite (output, refChip);
+    } else {
+	// FPA.BORE.X0,Y0 should be 0,0 in the focal plane, not the chip.  Ask for the
+	// coordinates which make refChip->toFPA(x,y) = (0,0)
+	psPlane *PT = psPlaneTransformGetCenter (refChip->toFPA, NONLIN_TOL);
+	psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.X0", PS_META_REPLACE, "boresite parameter", PT->x); 
+	psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.Y0", PS_META_REPLACE, "boresite parameter", PT->y); 
+	psFree (PT);
+    }
+
+    // rotate the chip-to-FPA transforms to have 0.0 posangle for refChip; 
+    // compensate by rotating fpa to TPA transform
+
+    // get the current posangle of the ref chip
+    float chipAngle = atan2 (refChip->toFPA->y->coeff[1][0], refChip->toFPA->x->coeff[1][0]);
+    fprintf (stderr, "chipAngle: %f\n", chipAngle*PS_DEG_RAD);
+    // psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.POSANGLE", PS_META_REPLACE, "boresite parameter", posangle);
+
+    // rotate the chip transforms
+    for (int i = 0; i < output->fpa->chips->n; i++) {
+	pmChip *chip = output->fpa->chips->data[i];
+	if (!chip->toFPA) continue;
+	// skip chips without astrometry
+
+	// save the region of this chip for the inversion below
+	psRegion *region = pmChipPixels (chip);
+
+	psPlaneTransform *toFPA = psPlaneTransformRotate (NULL, chip->toFPA, chipAngle);
+	psFree (chip->toFPA);
+	chip->toFPA = toFPA;
+
+	// invert the new fromFPA transform to get the new toFPA transform
+	psPlaneTransform *fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50);
+	psFree (chip->fromFPA);
+	chip->fromFPA = fromFPA;
+
+	psFree (region);
+
+	// save the transformation in the header
+	pmAstromWriteBilevelChip (chip->hdu->header, chip, NONLIN_TOL);
+    }
+
+    // get the current posangle of the fpa
+    float fpaAngle = atan2 (output->fpa->toTPA->y->coeff[1][0], output->fpa->toTPA->x->coeff[1][0]);
+    fprintf (stderr, "fpaAngle: %f\n", fpaAngle*PS_DEG_RAD);
+    // psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.POSANGLE", PS_META_REPLACE, "boresite parameter", posangle);
+
+    // remove the fpa rotation to generate a rotation-free model
+    psPlaneTransform *toTPA = psPlaneTransformRotate (NULL, output->fpa->toTPA, fpaAngle);
+    psFree (output->fpa->toTPA);
+    output->fpa->toTPA = toTPA;
+
+    psFree (output->fpa->fromTPA);
+    output->fpa->fromTPA = psPlaneTransformInvert(NULL, output->fpa->toTPA, *fpaRegion, 50);
+
+    // the model now describes the unrotated focal-plane
+    if (DEBUG) psastroDumpCorners ("corners.up.rot.dat", "corners.dn.rot.dat", output->fpa);
+
+    psMetadata *header = output->fpa->hdu->header;
+    pmAstromWriteBilevelMosaic (header, output->fpa, NONLIN_TOL);
+
+    psFree (fpaRegion);
+
+    return true;
+}
+
+bool psastroModelAdjustBoresite (pmFPAfile *output, pmChip *refChip) {
+
+    bool status;
+
+    psPlane  *boreCH  = psPlaneAlloc();
+    psPlane  *boreFP  = psPlaneAlloc();    
+    psPlane  *boreTP  = psPlaneAlloc();    
     psSphere *boreSky = psSphereAlloc();    
 
@@ -46,91 +147,99 @@
 	// skip the chips without astrometry
 
-	psPlaneTransform *fromFPA = psPlaneTransformSetCenter (NULL, chip->fromFPA, boreFP->x, boreFP->y);
+	// save the FPA region of this chip for the inversion below
+	psRegion *region = pmChipPixels (chip);
+
+	// the current toFPA returns boreFP->x,y for the boresite; subtract this from the transformations
+	chip->toFPA->x->coeff[0][0] -= boreFP->x;
+	chip->toFPA->y->coeff[0][0] -= boreFP->y;
+
+	// psPlaneTransform *toFPA = psPlaneTransformSetCenter (NULL, chip->toFPA, -boreFP->x, -boreFP->y);
+	// psFree (chip->toFPA);
+	// chip->toFPA = toFPA;
+
+	// invert the new fromFPA transform to get the new toFPA transform
+	// the region used here is the region covered by the chip in the FPA
+	psPlaneTransform *fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50); 
 	psFree (chip->fromFPA);
 	chip->fromFPA = fromFPA;
 
-	// invert the new fromFPA transform to get the new toFPA transform
-	psRegion *region = pmChipPixels (chip);
-	psFree (chip->toFPA);
-	chip->toFPA = psPlaneTransformInvert(NULL, chip->fromFPA, *region, 50);
 	psFree (region);
     }
 
-    // XXX we have now shifted the (0,0) pixel of the focal plane to the true boresite from the
-    // reported boresite.  At this point, we need to shift the tangent plane to use the new
-    // center as well.  if the toTPA transform were not linear, we would need to modify fromFPA
-    // to yield 0,0 at the new boresite location (ie, find Po,Qo = toTPA(Lo,Mo), probably could modify the
-    // toTPA/fromTPA transforms to use the new ref pixel, but this would only give us a tranf
-
-    // save the old (L,M) = (0,0) TP coordinate
-    float Po = output->fpa->toTPA->x->coeff[0][0];
-    float Qo = output->fpa->toTPA->y->coeff[0][0];
-
-    // the new toTPA yields the same TP coordinates for FP coordinates offset by Lo,Mo
-    psPlaneTransform *toTPA = psPlaneTransformSetCenter (NULL, output->fpa->toTPA, boreFP->x, boreFP->y);
-    psFree (output->fpa->toTPA);
-    output->fpa->toTPA = toTPA;
-    
-    // we are going to shift P,Q to have toTPA(0,0) = Po, Qo.  
-    // find the sky coordinates of the 0,0 pixel for the new frame
-    boreTP->x = output->fpa->toTPA->x->coeff[0][0] - Po;
-    boreTP->y = output->fpa->toTPA->y->coeff[0][0] - Qo;
+    if (DEBUG) psastroDumpCorners ("corners.up.shf.dat", "corners.dn.shf.dat", output->fpa);
+
+    // we have now adjusted the chips to use the correct boresite position as the center of the focal-plane system.
+    // we now need to reconstruct the TP to FP transformation, starting from stars projected about this new boresite position.
+
+    // find the R,D of the new boresite (boreFP -> 0,0; 0,0 -> -boreFP)
+    boreFP->x = -boreFP->x;
+    boreFP->y = -boreFP->y;
+    psPlaneTransformApply (boreTP, output->fpa->toTPA, boreFP);
     psDeproject (boreSky, boreTP, output->fpa->toSky); // find the RA,DEC coord of the focal-plane coordinate
 
-    // adjust the new TP frame to yield the same old (L,M) = 0,0 TP coordinate:
-    output->fpa->toTPA->x->coeff[0][0] = Po;
-    output->fpa->toTPA->y->coeff[0][0] = Qo;
-
-    // set the projection to use the new (P,Q) = (0,0) position
-    output->fpa->toSky->R = boreSky->r;
-    output->fpa->toSky->D = boreSky->d;
-
+    psProjection *newSky = psProjectionAlloc (boreSky->r, boreSky->d, output->fpa->toSky->Xs, output->fpa->toSky->Ys, output->fpa->toSky->type);
+
+    // generate a collection of points on the sky using the old toTPA transformation and toSky projection, projected with the newSky projection
+    // this is the FPA coordinate range covered by the FP: 
     psRegion *fpaRegion = pmAstromFPAExtent (output->fpa);
-
-    psFree (output->fpa->fromTPA);
-    output->fpa->fromTPA = psPlaneTransformInvert(NULL, output->fpa->toTPA, *fpaRegion, 50);
-
-    // rotate the chip to FPA transforms to have 0.0 posangle; 
-    // compensate by rotating fpa to TPA transforms
-
-    // get the current posangle of the ref chip
-    float posangle = atan2 (refChip->toFPA->y->coeff[1][0], refChip->toFPA->x->coeff[1][0]);
-    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.POSANGLE", PS_META_REPLACE, "boresite parameter", posangle);
-
-    // rotate the chip transforms
-    for (int i = 0; i < output->fpa->chips->n; i++) {
-	pmChip *chip = output->fpa->chips->data[i];
-	if (!chip->toFPA) continue;
-	// skip chips without astrometry
-
-	psPlaneTransform *toFPA = psPlaneTransformRotate (NULL, chip->toFPA, posangle);
-	psFree (chip->toFPA);
-	chip->toFPA = toFPA;
-
-	// invert the new fromFPA transform to get the new toFPA transform
-	psRegion *region = pmChipPixels (chip);
-	psFree (chip->fromFPA);
-	chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50);
-	psFree (region);
-
-	// XXX temporary output dump
-	psMetadata *header = chip->hdu->header;
-	// XXX to make the output single-level, this needs to be in a loop *after* the fromTPA rotation below
-	// pmAstromWriteWCS (header, output->fpa, chip, NONLIN_TOL);
-	pmAstromWriteBilevelChip (header, chip, NONLIN_TOL);
-    }
-
-    psPlaneTransform *fromTPA = psPlaneTransformRotate (NULL, output->fpa->fromTPA, posangle);
-    psFree (output->fpa->fromTPA);
-    output->fpa->fromTPA = fromTPA;
-
-    psFree (output->fpa->toTPA);
-    output->fpa->toTPA = psPlaneTransformInvert(NULL, output->fpa->fromTPA, *fpaRegion, 50);
-
-    psMetadata *header = output->fpa->hdu->header;
-    pmAstromWriteBilevelMosaic (header, output->fpa, NONLIN_TOL);
-
+    float dx = (fpaRegion->x1 - fpaRegion->x0) / 50.0;
+    float dy = (fpaRegion->y1 - fpaRegion->y0) / 50.0;
+
+    psPlane fp, tp;
+    psSphere sky;
+
+    psVector *FPx = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *FPy = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *TPx = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *TPy = psVectorAllocEmpty (100, PS_TYPE_F32);
+
+    // XXX a test: boreFP->x,y, should transform to tp.x,y = 0,0
+    fp.x = boreFP->x;
+    fp.y = boreFP->y;
+    psPlaneTransformApply (&tp, output->fpa->toTPA, &fp);
+    psDeproject (&sky, &tp, output->fpa->toSky); // find the RA,DEC coord of the focal-plane coordinate
+    psProject (&tp, &sky, newSky); // find the RA,DEC coord of the focal-plane coordinate
+
+    int Npts = 0;
+    for (fp.x = fpaRegion->x0; fp.x <= fpaRegion->x1; fp.x += dx) {
+	for (fp.y = fpaRegion->y0; fp.y <= fpaRegion->y1; fp.y += dy) {
+	    psPlaneTransformApply (&tp, output->fpa->toTPA, &fp);
+	    psDeproject (&sky, &tp, output->fpa->toSky); // find the RA,DEC coord of the focal-plane coordinate
+	    psProject (&tp, &sky, newSky); // find the RA,DEC coord of the focal-plane coordinate
+
+	    // we are fitting points in the NEW FP system to points in the NEW TP system
+	    FPx->data.F32[Npts] = fp.x - boreFP->x;
+	    FPy->data.F32[Npts] = fp.y - boreFP->y;
+	    TPx->data.F32[Npts] = tp.x;
+	    TPy->data.F32[Npts] = tp.y;
+	    psVectorExtend (FPx, 100, 1);
+	    psVectorExtend (FPy, 100, 1);
+	    psVectorExtend (TPx, 100, 1);
+	    psVectorExtend (TPy, 100, 1);
+	    Npts ++;
+	}
+    }
     psFree (fpaRegion);
 
+    // fit both up and down transformations to the same points
+    psVectorFitPolynomial2D (output->fpa->toTPA->x, NULL, 0, TPx, NULL, FPx, FPy);
+    psVectorFitPolynomial2D (output->fpa->toTPA->y, NULL, 0, TPy, NULL, FPx, FPy);
+    psVectorFitPolynomial2D (output->fpa->fromTPA->x, NULL, 0, FPx, NULL, TPx, TPy);
+    psVectorFitPolynomial2D (output->fpa->fromTPA->y, NULL, 0, FPy, NULL, TPx, TPy);
+
+    psFree (output->fpa->toSky);
+    output->fpa->toSky = newSky;
+
+    if (DEBUG) psastroDumpCorners ("corners.up.bore.dat", "corners.dn.bore.dat", output->fpa);
+
+    psFree (FPx);
+    psFree (FPy);
+    psFree (TPx);
+    psFree (TPy);
+
+    psFree (boreCH);
+    psFree (boreFP);
+    psFree (boreTP);
+    psFree (boreSky);
     return true;
 }
Index: branches/cnb_branch_20080830/psastro/src/psastroModelAnalysis.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroModelAnalysis.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroModelAnalysis.c	(revision 20033)
@@ -1,3 +1,4 @@
 # include "psastroStandAlone.h"
+# define NONLIN_TOL 0.001
 
 bool psastroModelAnalysis (pmConfig *config) {
@@ -12,8 +13,5 @@
     }
 
-    pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSASTRO.OUT.MODEL");
-    if (!status) psAbort ("Can't find output pmFPAfile PSASTRO.OUT.MODEL");
-
-    // physical pixel scale in microns per pixel
+    // reference chip for boresite parameters
     char *refChip = psMetadataLookupStr (&status, recipe, "PSASTRO.MODEL.REF.CHIP");
     if (!refChip) {
@@ -21,4 +19,23 @@
 	return false; 
     } 
+
+    pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSASTRO.OUT.MODEL");
+    if (!status) psAbort ("Can't find output pmFPAfile PSASTRO.OUT.MODEL");
+
+    // measure the boresite position from a rotation sequence?
+    bool fitBoresite = psMetadataLookupBool (&status, recipe, "PSASTRO.MODEL.FIT.BORESITE");
+    if (!fitBoresite) {
+	psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.X0", PS_META_REPLACE, "boresite parameter", 0.0); 
+	psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.Y0", PS_META_REPLACE, "boresite parameter", 0.0); 
+	psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.RX", PS_META_REPLACE, "boresite parameter", 0.0); 
+	psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.RY", PS_META_REPLACE, "boresite parameter", 0.0); 
+	psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.T0", PS_META_REPLACE, "boresite parameter", 0.0); 
+	psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.P0", PS_META_REPLACE, "boresite parameter", 0.0); 
+	psMetadataAddStr (output->fpa->concepts, PS_LIST_TAIL, "FPA.REF.CHIP", PS_META_REPLACE, "boresite parameter", refChip);
+	return true; 
+    } 
+
+    char *outroot = psMetadataLookupStr (&status, config->arguments, "OUTPUT");
+    if (!status || !outroot) psAbort ("Can't find outroot on config->arguments");
 
     /* model analysis:
@@ -35,4 +52,5 @@
      * T_0 : reference angle for rotator
      * P_0 : orientation of boresite ellipse
+     *
      */
 
@@ -58,4 +76,14 @@
     output->fpa = NULL;
 
+    char filename[256];
+    snprintf (filename, 256, "%s.bore", outroot);
+    FILE *outfile = fopen (filename, "w");
+    if (!outfile) {
+        psError(PSASTRO_ERR_ARGUMENTS, true, "cannot open %s for output", filename);
+	return false;
+    }
+
+    float posBoundary = 0.0;
+
     // extract the relevant measured and reported values from the reference chip
     for (int i = 0; i < files->n; i++) {
@@ -63,5 +91,5 @@
 	pmFPAfile *input = file->data.V;
 
-	// reported rotator position angle (this should perhaps be ROT, not POS)
+	// reported rotator position angle
 	double POSANGLE = psMetadataLookupF64 (&status, input->fpa->concepts, "FPA.POSANGLE"); 
 	if (!status) psAbort ("missing FPA.POSANGLE");
@@ -71,5 +99,4 @@
 	if (!chip) psAbort ("invalid chip name for reference");
 
-	fprintf (stderr, "input %d : %zx : %zx : %zx\n", i, (size_t) input, (size_t) chip, (size_t) chip->toFPA);
 	if (!chip->toFPA) continue;
 
@@ -80,13 +107,29 @@
 
 	// we have two measurements of the posangle (may be parity flipped to different quadrants)
-	float posangle = PM_DEG_RAD * atan2 (chip->toFPA->y->coeff[1][0], chip->toFPA->x->coeff[1][0]);
-	posZero->data.F32[n] = POSANGLE - posangle;
-	while (posZero->data.F32[n] > 360.0) posZero->data.F32[n] -= 360.0;
-	while (posZero->data.F32[n] <   0.0) posZero->data.F32[n] += 360.0;
+	// atan2 returns values in the range 0-2pi
+	// all posZero values should be clustered in some region, but we need to flip over the 0,360 boundary correctly.
+	// push all to one side or the other
+	float chipAngle = PM_DEG_RAD * atan2 (chip->toFPA->y->coeff[1][0], chip->toFPA->x->coeff[1][0]);
+	float fpaAngle = PM_DEG_RAD * atan2 (input->fpa->toTPA->y->coeff[1][0], input->fpa->toTPA->x->coeff[1][0]);
+
+	posZero->data.F32[n] = POSANGLE - chipAngle - fpaAngle;
+	if (n == 0) {
+	    posBoundary = posZero->data.F32[n] + 180.0;
+	} else {
+	    while (posZero->data.F32[n] > posBoundary) posZero->data.F32[n] -= 360.0;
+	    while (posZero->data.F32[n] < posBoundary - 360.0) posZero->data.F32[n] += 360.0;
+	}
 
 	Po->data.F32[n] = POSANGLE * PM_RAD_DEG; // reported position angle
-	Xo->data.F32[n] = chip->fromFPA->x->coeff[0][0]; // reported boresite x position in ref chip coordinates
-	Yo->data.F32[n] = chip->fromFPA->y->coeff[0][0]; // reported boresite y position in ref chip coordinates
-	fprintf (stderr, "%d : %f %f : %f = %f - %f\n", i, Xo->data.F32[n], Yo->data.F32[n], posZero->data.F32[n], POSANGLE, posangle);
+	float xc = chip->fromFPA->x->coeff[0][0]; // reported boresite x position in ref chip coordinates
+	float yc = chip->fromFPA->y->coeff[0][0]; // reported boresite y position in ref chip coordinates
+	// XXX this can also be derived from toFPA via GetCenter....
+	
+	psPlane *PT = psPlaneTransformGetCenter (chip->toFPA, NONLIN_TOL);
+	Xo->data.F32[n] = PT->x; // reported boresite x position in ref chip coordinates
+	Yo->data.F32[n] = PT->y; // reported boresite y position in ref chip coordinates
+	psFree (PT);
+
+	fprintf (outfile, "%d : %f %f : %f = %f - %f - %f | %f %f\n", i, Xo->data.F32[n], Yo->data.F32[n], posZero->data.F32[n], POSANGLE, chipAngle, fpaAngle, xc, yc);
 	n ++;
     }
@@ -100,9 +143,11 @@
     psVectorStats (stats, posZero, NULL, NULL, 0);
 
-    fprintf (stderr, "pos zero %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+    fprintf (outfile, "# pos zero %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
     psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.POS_ZERO", PS_META_REPLACE, "offset between obs and meas posangle", stats->sampleMedian); 
+    fclose (outfile);
 
-    psVector *params = psastroModelFitBoresite (Xo, Yo, Po);
+    psVector *params = psastroModelFitBoresite (Xo, Yo, Po, outroot);
     if (params->n != 6) psAbort ("error");
+
 
     psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.X0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_X0]); 
Index: branches/cnb_branch_20080830/psastro/src/psastroModelBoresite.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroModelBoresite.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroModelBoresite.c	(revision 20033)
@@ -45,5 +45,5 @@
 	    dPAR[PAR_RY]  = +sntht*csphi;
 	    dPAR[PAR_P0]  = -PAR[PAR_RY]*sntht*snphi - PAR[PAR_RX]*cstht*csphi;
-	    dPAR[PAR_T0]  = -PAR[PAR_RY]*sntht*csphi - PAR[PAR_RX]*cstht*snphi;
+	    dPAR[PAR_T0]  = -PAR[PAR_RY]*cstht*csphi - PAR[PAR_RX]*sntht*snphi;
 	}
 	return (value);
Index: branches/cnb_branch_20080830/psastro/src/psastroModelDataSave.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroModelDataSave.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroModelDataSave.c	(revision 20033)
@@ -1,3 +1,4 @@
 # include "psastroStandAlone.h"
+# define NONLIN_TOL 0.001 /* tolerance in pixels */
 
 # define ESCAPE { \
@@ -22,6 +23,6 @@
 
     pmFPAview *view = pmFPAviewAlloc (0);
+    pmChip *chip = NULL;
 
-    pmChip *chip = NULL;
     while ((chip = pmFPAviewNextChip (view, output->fpa, 1)) != NULL) {
         psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
@@ -35,2 +36,5 @@
     return true;
 }
+
+
+
Index: branches/cnb_branch_20080830/psastro/src/psastroModelFit.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroModelFit.c	(revision 20033)
+++ branches/cnb_branch_20080830/psastro/src/psastroModelFit.c	(revision 20033)
@@ -0,0 +1,37 @@
+# include "psastroStandAlone.h"
+
+int main (int argc, char **argv) {
+
+    psTimerStart ("complete");
+
+    if (argc != 3) {
+	fprintf (stderr, "USAGE: psastroModelFit (input) (output)\n");
+	exit (1);
+    }
+
+    FILE *f = fopen (argv[1], "r");
+    if (f == NULL) {
+	fprintf (stderr, "problem opening data file %s\n", argv[1]);
+	exit (2);
+    }
+
+    char name[1024];
+    psVector *Xo = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *Yo = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *Po = psVectorAllocEmpty (100, PS_TYPE_F32);
+
+    float x, y, p;
+
+    while (fscanf (f, "%s %f %f %f", name, &x, &y, &p) != EOF) {
+	psVectorAppend (Xo, x);
+	psVectorAppend (Yo, y);
+	psVectorAppend (Po, p*PS_RAD_DEG);
+    }	
+
+    // psTraceSetLevel("psLib.math", 5);
+    psastroModelFitBoresite (Xo, Yo, Po, argv[2]);
+
+    psLogMsg ("psastro", 3, "complete psastroModelFit run: %f sec\n", psTimerMark ("complete"));
+
+    exit (EXIT_SUCCESS);
+}
Index: branches/cnb_branch_20080830/psastro/src/psastroModelFitBoresite.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroModelFitBoresite.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroModelFitBoresite.c	(revision 20033)
@@ -2,5 +2,5 @@
 
 // we now have a set of observed L,M values.  fit these to the boresite model
-psVector *psastroModelFitBoresite (psVector *Xo, psVector *Yo, psVector *Po) {
+psVector *psastroModelFitBoresite (psVector *Xo, psVector *Yo, psVector *Po, char *outroot) {
 
     assert (Xo->n > 2);
@@ -60,5 +60,5 @@
     params->data.F32[PAR_T0] = 0.0;
 
-    psMinimization *myMin = psMinimizationAlloc (15, 0.001);
+    psMinimization *myMin = psMinimizationAlloc (25, 0.001);
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
     
@@ -74,11 +74,19 @@
     psMinimizeLMChi2(myMin, covar, params, constraint, x, y, NULL, psastroModelBoresite);
 
-    fprintf (stderr, "fitted values:\n");
-    fprintf (stderr, "Xo:  %f\n", params->data.F32[PAR_X0]);
-    fprintf (stderr, "Yo:  %f\n", params->data.F32[PAR_Y0]);
-    fprintf (stderr, "RX:  %f\n", params->data.F32[PAR_RX]);
-    fprintf (stderr, "RY:  %f\n", params->data.F32[PAR_RY]);
-    fprintf (stderr, "P0:  %f\n", params->data.F32[PAR_P0]);
-    fprintf (stderr, "T0:  %f\n", params->data.F32[PAR_T0]);
+    char filename[256];
+    snprintf (filename, 256, "%s.pars", outroot);
+    FILE *outfile = fopen (filename, "w");
+    if (!outfile) {
+        psAbort("cannot open %s for output", filename);
+    }
+
+    fprintf (outfile, "# fitted values:\n");
+    fprintf (outfile, "Xo:  %f\n", params->data.F32[PAR_X0]);
+    fprintf (outfile, "Yo:  %f\n", params->data.F32[PAR_Y0]);
+    fprintf (outfile, "RX:  %f\n", params->data.F32[PAR_RX]);
+    fprintf (outfile, "RY:  %f\n", params->data.F32[PAR_RY]);
+    fprintf (outfile, "P0:  %f\n", params->data.F32[PAR_P0]);
+    fprintf (outfile, "T0:  %f\n", params->data.F32[PAR_T0]);
+    fclose (outfile);
 
     return params;
Index: branches/cnb_branch_20080830/psastro/src/psastroMosaicAstrom.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroMosaicAstrom.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroMosaicAstrom.c	(revision 20033)
@@ -1,10 +1,15 @@
 # include "psastroInternal.h"
 # define NONLIN_TOL 0.001 /* tolerance in pixels */
+
+bool psastroMosaicFit (pmFPA *fpa, psMetadata *recipe, const char *rootname, int pass);
 
 // XXX require this fpa to have multiple chip extensions and a PHU?
 bool psastroMosaicAstrom (pmConfig *config) {
 
+    bool status;
+    char filename[256];
+
     // select the current recipe
-    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE);
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE);
     if (!recipe) {
 	psError(PSASTRO_ERR_CONFIG, false, "Can't find PSASTRO recipe!\n");
@@ -34,92 +39,26 @@
     # endif
 
-    // given the existing per-chip astrometry, determine matches between raw and ref stars
-    // is this needed? yes, if we didn't do SingleChip astrometry first
-    if (!psastroMosaicSetMatch (fpa, recipe, 0)) {
-	psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic");
-	return false;
-    }
-    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.0.dat"); }
+    char *outroot = psMetadataLookupStr (&status, config->arguments, "OUTPUT");
+    if (!status || !outroot) psAbort ("Can't find outroot on config->arguments");
 
-    // fitted chips will follow the local plate-scale, hiding the distortion
-    // modify the chip->toFPA scaling to match knowledge about pixel scale,
-    // then recalculate raw and ref positions
-    if (!psastroMosaicCommonScale (fpa, recipe)) {
-	psError(PSASTRO_ERR_UNKNOWN, false, "failed to set a common scale for the chips");
-	return false;
-    }
-    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.1.dat"); }
-
-    // fit the distortion by fitting its gradient
-    // apply the new distortion terms up and down
-    // refit the per-chip terms with linear fits only
-    if (!psastroMosaicGradients (fpa, recipe)) {
-	psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure mosaic gradients");
-	return false;
-    }
-    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.2.dat"); }
-
-    // measure the astrometry for the chips under the distortion term
-    if (!psastroMosaicChipAstrom (fpa, recipe, 0)) {
-	psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure chip astrometry in mosaic mode");
-	return false;
-    }
-    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.3.dat"); }
-
-    // do a second pass on the distortion with improved chip positions and rotations
-    // first, re-perform the match with a slightly tighter circle
-    if (!psastroMosaicSetMatch (fpa, recipe, 1)) {
-	psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic (2nd pass)");
-	return false;
-    }
-    if (!psastroMosaicCommonScale (fpa, recipe)) {
-	psError(PSASTRO_ERR_UNKNOWN, false, "failed to set a common scale for the chips (2nd pass)");
-	return false;
-    }
-    if (!psastroMosaicGradients (fpa, recipe)) {
-	psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure mosaic gradients (2nd pass)");
-	return false;
-    }
-    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.4.dat"); }
-
-    if (!psastroMosaicChipAstrom (fpa, recipe, 1)) {
-	psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure chip astrometry in mosaic mode (2nd pass)");
-	return false;
-    }
-    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.5.dat"); }
-
-    // do a third pass on the distortion with improved chip positions and rotations
-    // first, re-perform the match with a slightly tighter circle
-    if (!psastroMosaicSetMatch (fpa, recipe, 2)) {
-	psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic (3rd pass)");
-	return false;
-    }
-    if (!psastroMosaicCommonScale (fpa, recipe)) {
-	psError(PSASTRO_ERR_UNKNOWN, false, "failed to set a common scale for the chips (3rd pass)");
-	return false;
-    }
-    if (!psastroMosaicGradients (fpa, recipe)) {
-	psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure mosaic gradients (3rd pass)");
-	return false;
-    }
-    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.6.dat"); }
-
-    if (!psastroMosaicChipAstrom (fpa, recipe, 2)) {
-	psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure chip astrometry in mosaic mode (3rd pass)");
-	return false;
-    }
-    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.7.dat"); }
+    if (!psastroMosaicFit (fpa, recipe, outroot, 0)) return false;
+    if (!psastroMosaicFit (fpa, recipe, outroot, 1)) return false;
+    if (!psastroMosaicFit (fpa, recipe, outroot, 2)) return false;
+    if (!psastroMosaicFit (fpa, recipe, outroot, 3)) return false;
 
     // now fit the chips under the common distortion with higher-order terms
     // first, re-perform the match with a slightly tighter circle
-    if (!psastroMosaicSetMatch (fpa, recipe, 3)) {
-	psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic (3rd pass)");
+    if (!psastroMosaicSetMatch (fpa, recipe, 4)) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic (4th pass)");
 	return false;
     }
-    if (!psastroMosaicChipAstrom (fpa, recipe, 3)) {
-	psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure chip astrometry in mosaic mode (3rd pass)");
+    if (!psastroMosaicChipAstrom (fpa, recipe, 4)) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure chip astrometry in mosaic mode (4th pass)");
 	return false;
     }
-    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.8.dat"); }
+    if (psTraceGetLevel("psastro.dump") > 0) { 
+	snprintf (filename, 256, "%s.10.dat", outroot);
+	psastroDumpMatches (fpa, filename); 
+    }
 
     // save WCS and analysis metadata in update header.
@@ -150,2 +89,63 @@
  * sky (ra, dec)
  */
+
+// 1: match 4,5
+// 2: match 6,7
+// 3: match 8,9
+bool psastroMosaicFit (pmFPA *fpa, psMetadata *recipe, const char *rootname, int pass) {
+
+    char filename[256];
+
+    // given the existing per-chip astrometry, determine matches between raw and ref stars
+    // is this needed? yes, if we didn't do SingleChip astrometry first
+    if (!psastroMosaicSetMatch (fpa, recipe, pass)) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic (pass %d)", pass);
+	return false;
+    }
+
+    if ((pass == 0) && (psTraceGetLevel("psastro.dump.psastroMosaicAstrom") > 1)) { 
+	snprintf (filename, 256, "%s.0.dat", rootname);
+	psastroDumpMatches (fpa, filename); 
+    }
+
+    // fitted chips will follow the local plate-scale, hiding the distortion
+    // modify the chip->toFPA scaling to match knowledge about pixel scale,
+    // then recalculate raw and ref positions
+    if (!psastroMosaicCommonScale (fpa, recipe)) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to set a common scale for the chips (pass %d)", pass);
+	return false;
+    }
+
+    if ((pass == 0) && (psTraceGetLevel("psastro.dump.psastroMosaicAstrom") > 1)) { 
+	snprintf (filename, 256, "%s.1.dat", rootname);
+	psastroDumpMatches (fpa, filename); 
+    }
+
+    // fit the distortion by fitting its gradient
+    // apply the new distortion terms up and down
+    // refit the per-chip terms with linear fits only
+    if (!psastroMosaicDistortion (fpa, recipe, pass)) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure mosaic gradients (pass %d)", pass);
+	return false;
+    }
+
+    snprintf (filename, 256, "%s.%d.dat", rootname, 2*pass + 2);
+    if (psTraceGetLevel("psastro.dump.psastroMosaicAstrom") > 1) { psastroDumpMatches (fpa, filename); }
+    
+    // measure the astrometry for the chips under the distortion term
+    if (!psastroMosaicChipAstrom (fpa, recipe, pass)) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure chip astrometry in mosaic mode (pass %d)", pass);
+	return false;
+    }
+
+    int traceLevel = psTraceGetLevel("psastro.dump.psastroMosaicAstrom");
+    if (traceLevel > 1) {
+	snprintf (filename, 256, "%s.%d.dat", rootname, 2*pass + 3);
+	psastroDumpMatches (fpa, filename); 
+    } else {
+	snprintf (filename, 256, "%s.dat", rootname);
+	psastroDumpMatches (fpa, filename); 
+    }
+
+    return true;
+}
Index: branches/cnb_branch_20080830/psastro/src/psastroMosaicCorrectDistortion.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroMosaicCorrectDistortion.c	(revision 20033)
+++ branches/cnb_branch_20080830/psastro/src/psastroMosaicCorrectDistortion.c	(revision 20033)
@@ -0,0 +1,81 @@
+# include "psastroInternal.h"
+
+bool psastroMosaicCorrectDistortion (pmFPA *fpa, psPlaneTransform *TPtoFP) {
+
+    // invert the linear TPtoFP transform
+    psPlaneTransform *FPtoTP = p_psPlaneTransformLinearInvert (TPtoFP);
+    if (FPtoTP == NULL) {
+        psError (PS_ERR_UNKNOWN, false, "failed to invert TPtoFP\n");
+        return false;
+    }
+
+    // store the new coeffs in a new structure
+    psPlaneTransform *toTPAnew = psPlaneTransformAlloc(fpa->toTPA->x->nX, fpa->toTPA->x->nY);
+
+    // set the new coeffs, or set the mask
+    for (int i = 0; i <= toTPAnew->y->nX; i++) {
+        for (int j = 0; j <= toTPAnew->y->nY; j++) {
+
+	    // init the coeffs
+	    toTPAnew->x->coeff[i][j] = 0.0;
+	    toTPAnew->y->coeff[i][j] = 0.0;
+
+	    // if both are masked, mask the outpu
+            if ((fpa->toTPA->x->coeffMask[i][j] & PS_POLY_MASK_SET) && (fpa->toTPA->x->coeffMask[i][j] & PS_POLY_MASK_SET)) {
+		toTPAnew->x->coeffMask[i][j] = PS_POLY_MASK_SET;
+		toTPAnew->y->coeffMask[i][j] = PS_POLY_MASK_SET;
+		continue;
+	    } 
+
+	    // set the X terms
+	    toTPAnew->x->coeff[i][j] += (fpa->toTPA->x->coeffMask[i][j] & PS_POLY_MASK_SET) ? 0.0 : FPtoTP->x->coeff[1][0] * fpa->toTPA->x->coeff[i][j];
+	    toTPAnew->x->coeff[i][j] += (fpa->toTPA->y->coeffMask[i][j] & PS_POLY_MASK_SET) ? 0.0 : FPtoTP->x->coeff[0][1] * fpa->toTPA->y->coeff[i][j];
+
+	    // set the Y terms
+	    toTPAnew->y->coeff[i][j] += (fpa->toTPA->x->coeffMask[i][j] & PS_POLY_MASK_SET) ? 0.0 : FPtoTP->y->coeff[1][0] * fpa->toTPA->x->coeff[i][j];
+	    toTPAnew->y->coeff[i][j] += (fpa->toTPA->y->coeffMask[i][j] & PS_POLY_MASK_SET) ? 0.0 : FPtoTP->y->coeff[0][1] * fpa->toTPA->y->coeff[i][j];
+        }
+    }
+
+    // adjust the 0,0 terms:
+    toTPAnew->x->coeff[0][0] += FPtoTP->x->coeff[0][0];
+    toTPAnew->y->coeff[0][0] += FPtoTP->y->coeff[0][0];
+
+    psFree (fpa->toTPA);
+    fpa->toTPA = toTPAnew;
+
+    // invert toTPA to determine fromTPA. choose an appropriate region based on the dimensions
+    // of the complete FPA
+    psRegion *region = pmAstromFPAExtent (fpa);
+
+    psFree (fpa->fromTPA);
+    fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, *region, 50);
+    psFree (region);
+
+    if (fpa->fromTPA == NULL) {
+        psError (PS_ERR_UNKNOWN, false, "failed to invert fpa->toTPA\n");
+	psFree (FPtoTP);
+        return false;
+    }
+
+    psFree (FPtoTP);
+    return true;
+}
+
+// we have three coordinate systems and two polynomial transformations:
+// TP is the raw tangent plane coordinate system (aligned with RA,DEC)
+// FP is the raw focal plane coordinate system (aligned with camera X,Y)
+// dFP is the distorted focal plane coordinate system
+
+// fpa->toTPA is a polynomial transformation from FP to dFP
+//   L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j and 
+//   M(x,y) = \sum_i \sum_j B_{i,j} x^i y^j 
+//   where (x,y) are the FP coords and L,M  are the dFP coords
+
+// FPtoTP is a linear transformation from dFP to TP (we recover this from TPtoFP by inversion)
+// P(L,M) = r_xo + r_xx L + r_xy M
+// Q(L,M) = r_yo + r_yx L + r_yy M
+
+// we need to merge them into a single transformation:
+// P(x,y) = r_xo + r_xx * \sum_i \sum_j A_{i,j} x^i y^j + r_xy * \sum_i \sum_j B_{i,j} x^i y^j 
+// Q(x,y) = r_yo + r_yx * \sum_i \sum_j A_{i,j} x^i y^j + r_yy * \sum_i \sum_j B_{i,j} x^i y^j 
Index: branches/cnb_branch_20080830/psastro/src/psastroMosaicDistortion.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroMosaicDistortion.c	(revision 20033)
+++ branches/cnb_branch_20080830/psastro/src/psastroMosaicDistortion.c	(revision 20033)
@@ -0,0 +1,55 @@
+# include "psastroInternal.h"
+# define DEBUG 0
+
+bool psastroMosaicDistortion (pmFPA *fpa, psMetadata *recipe, int pass) {
+
+    // fit a linear transformation from reference TP to observed FP coords
+    psPlaneTransform *TPtoFP = psastroMosaicFitRotAndScale (fpa);
+    if (!TPtoFP) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to fit a linear TP correction\n");
+        return false;
+    }
+
+    if (DEBUG && (pass == 0)) psastroDumpCorners ("corners.up.v1.dat", "corners.dn.v1.dat", fpa);
+
+    // Correct the current reference star TP coordinates to the nearly-FP
+    // system.  We note two points here: 1) the corrected TP coordinates are
+    // NOT consistent with the sky coordinates, 2) the remaining differnce
+    // between the new TP reference coords and the observerd FP coordinates is
+    // only distortion
+    if (!psastroMosaicApplyRotAndScale (fpa, TPtoFP)) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to apply the linear TP correction\n");
+	psFree (TPtoFP);
+        return false;
+    }
+
+    if (DEBUG && (pass == 0)) psastroDumpCorners ("corners.up.v2.dat", "corners.dn.v2.dat", fpa);
+
+    if (!psastroMosaicDistortionFromGradients (fpa, recipe)) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to fit the distortion field\n");
+	psFree (TPtoFP);
+        return false;
+    }
+
+    if (DEBUG && (pass == 0)) psastroDumpCorners ("corners.up.v3.dat", "corners.dn.v3.dat", fpa);
+
+    if (!psastroMosaicCorrectDistortion (fpa, TPtoFP)) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to correct the distortion for the linear fit\n");
+	psFree (TPtoFP);
+        return false;
+    }
+	
+    if (DEBUG && (pass == 0)) psastroDumpCorners ("corners.up.v4.dat", "corners.dn.v4.dat", fpa);
+
+    if (!psastroMosaicSetAstrom (fpa)) {
+	psError(PSASTRO_ERR_UNKNOWN, false, "failed to apply mosaic distortion terms\n");
+	psFree (TPtoFP);
+        return false;
+    }
+
+    if (DEBUG && (pass == 0)) psastroDumpCorners ("corners.up.v5.dat", "corners.dn.v5.dat", fpa);
+
+    psFree (TPtoFP);
+    return true;
+}
+
Index: branches/cnb_branch_20080830/psastro/src/psastroMosaicFPtoTP.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroMosaicFPtoTP.c	(revision 20033)
+++ branches/cnb_branch_20080830/psastro/src/psastroMosaicFPtoTP.c	(revision 20033)
@@ -0,0 +1,164 @@
+# include "psastroInternal.h"
+
+/************************************************/
+psPlaneTransform *psastroMosaicFitRotAndScale (pmFPA *fpa) {
+
+    pmChip *chip = NULL;
+    pmCell *cell = NULL;
+    pmReadout *readout = NULL;
+
+    // XXX first pass: fit and remove any linear fp->tp transformation
+    // accumulate FP(x,y) & TP(x,y) in a single set of vectors
+
+    psVector *X  = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *Y  = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *x  = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *y  = psVectorAllocEmpty(100, PS_TYPE_F32);
+
+    pmFPAview *view = pmFPAviewAlloc (0);
+
+    // int nPts = 0;
+    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
+        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+        if (!chip->process || !chip->file_exists) { continue; }
+	if (!chip->toFPA) { continue; }
+	
+	while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
+            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+            if (!cell->process || !cell->file_exists) { continue; }
+
+	    // process each of the readouts
+	    // XXX there can only be one readout per chip, right?
+	    while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
+		if (! readout->data_exists) { continue; }
+
+		// select the raw objects for this readout
+		psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
+		if (rawstars == NULL) { continue; }
+
+		// select the raw objects for this readout
+		psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
+		if (refstars == NULL) { continue; }
+
+		psArray *match = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH");
+		if (match == NULL) { continue; }
+
+		// take the matched stars, first fit
+		for (int i = 0; i < match->n; i++) {
+		    pmAstromMatch *pair = match->data[i];
+		    pmAstromObj *rawStar = rawstars->data[pair->raw];
+		    pmAstromObj *refStar = refstars->data[pair->ref];
+
+		    // independent variables
+		    X->data.F32[X->n] = refStar->TP->x;
+		    Y->data.F32[Y->n] = refStar->TP->y;
+
+		    // fitted values
+		    x->data.F32[x->n] = rawStar->FP->x;
+		    y->data.F32[y->n] = rawStar->FP->y;
+
+		    psVectorExtend (X, 100, 1);
+		    psVectorExtend (Y, 100, 1);
+		    psVectorExtend (x, 100, 1);
+		    psVectorExtend (y, 100, 1);
+		}
+	    }
+	}
+    }
+    // x->n = y->n = X->n = Y->n = nPts;
+
+    // linear fit without xy cross term
+    psPlaneTransform *map = psPlaneTransformAlloc (1, 1);
+    map->x->coeffMask[1][1] = PS_POLY_MASK_SET;
+    map->y->coeffMask[1][1] = PS_POLY_MASK_SET;
+
+    // constant errors
+    psVector *mask = psVectorAlloc (X->n, PS_TYPE_U8);
+    psVectorInit (mask, 0);
+
+    // the stats options supplied are used to perform the clip fitting
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    stats->clipIter = 1;
+
+    // fit TP-to-FP transformation
+
+    // we run 3 cycles clipping in each of x and y, with only one iteration each.
+    // XXX use the stats lookups functions to get the width and center
+    for (int i = 0; i < 3; i++) {
+	if (!psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, x, NULL, X, Y)) {
+            psError(PS_ERR_UNKNOWN, false, "failure in clip-fitting for x\n");
+	    psFree (map);
+	    map = NULL;
+	    goto escape;
+	}
+        psTrace ("psastro", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, x->n);
+
+        if (!psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, y, NULL, X, Y)) {
+            psError(PS_ERR_UNKNOWN, false, "failure in clip-fitting for y\n");
+	    psFree (map);
+	    map = NULL;
+	    goto escape;
+	}
+        psTrace ("psastro", 3, "y resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, y->n);
+    }
+
+escape:
+    psFree (x);
+    psFree (y);
+    psFree (X);
+    psFree (Y);
+    psFree (mask);
+    psFree (view);
+    psFree (stats);
+
+    return (map);
+}
+
+/*****************************/
+bool psastroMosaicApplyRotAndScale (pmFPA *fpa, psPlaneTransform *TPtoFP) {
+
+    pmChip *chip = NULL;
+    pmCell *cell = NULL;
+    pmReadout *readout = NULL;
+
+    pmFPAview *view = pmFPAviewAlloc (0);
+    psPlane newTP;
+
+    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
+        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+        if (!chip->process || !chip->file_exists) { continue; }
+	if (!chip->toFPA) { continue; }
+	
+	while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
+            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+            if (!cell->process || !cell->file_exists) { continue; }
+
+	    // process each of the readouts
+	    // XXX there can only be one readout per chip, right?
+	    while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
+		if (! readout->data_exists) { continue; }
+
+		// select the raw objects for this readout
+		psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
+		if (refstars == NULL) { continue; }
+
+		// take the matched stars, first fit
+		for (int i = 0; i < refstars->n; i++) {
+		    pmAstromObj *refStar = refstars->data[i];
+
+		    // Correct the current reference star TP coordinates to the nearly-FP
+		    // system.  We note two points here: 1) the corrected TP coordinates are
+		    // NOT consistent with the sky coordinates, 2) the remaining differnce
+		    // between the new TP reference coords and the observerd FP coordinates is
+		    // only distortion
+
+		    psPlaneTransformApply (&newTP, TPtoFP, refStar->TP);
+		    refStar->TP->x = newTP.x;
+		    refStar->TP->y = newTP.y;
+		}
+	    }
+	}
+    }
+    psFree (view);
+    return true;
+}
Index: branches/cnb_branch_20080830/psastro/src/psastroMosaicGradients.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroMosaicGradients.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroMosaicGradients.c	(revision 20033)
@@ -2,5 +2,5 @@
 static int nPass = 0;
 
-bool psastroMosaicGradients (pmFPA *fpa, psMetadata *recipe) {
+bool psastroMosaicDistortionFromGradients (pmFPA *fpa, psMetadata *recipe) {
 
     bool status;
@@ -9,4 +9,7 @@
     pmReadout *readout = NULL;
     psArray *gradients = NULL;
+
+    // Measure the gradient as a function of position.  This must be performed between the
+    // corrected ref->TP and the observed raw->FP, for which the distortion is a perturbation.
 
     pmFPAview *view = pmFPAviewAlloc (0);
@@ -59,4 +62,6 @@
     }
 
+    // Fit the gradient field and convert to the distortion terms.
+
     // allocate mosaic-level polynomial transformation and set masks needed by DVO
     int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.ORDER");
@@ -95,13 +100,7 @@
     }
 	
-    if (!psastroMosaicSetAstrom (fpa)) {
-	psError(PSASTRO_ERR_UNKNOWN, false, "failed to apply mosaic distortion terms\n");
-	psFree (gradients);
-	psFree (view);
-        return false;
-    }
-
     psFree (gradients);
     psFree (view);
     return true;
 }
+
Index: branches/cnb_branch_20080830/psastro/src/psastroMosaicOneChip.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroMosaicOneChip.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroMosaicOneChip.c	(revision 20033)
@@ -12,6 +12,4 @@
     char errorWord[64];
     char orderWord[64];
-
-    assert (iteration < 4);
 
     PS_ASSERT_PTR_NON_NULL(chip,    false);
@@ -86,7 +84,15 @@
 
     // XXX allow statitic to be set by the user
-    psStats *fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-    fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA");
-    fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER");
+    // only clip if we are fitting the chip parameters.
+    psStats *fitStats = NULL;
+//    if (order == 0) {
+//	fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+//	fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA");
+//	fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER");
+//    } else {
+	fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+	fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA");
+	fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER");
+//    }
 
     // need to pass in an update header, sent in from above
@@ -117,5 +123,5 @@
     // XXX should these result in errors or be handled another way?
     psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d", astError, astNstar);
-    if (astError > maxError) {
+    if ((maxError > 0) && (astError > maxError)) {
 	psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
 	validSolution = false;
Index: branches/cnb_branch_20080830/psastro/src/psastroOneChip.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroOneChip.c	(revision 20026)
+++ 	(revision )
@@ -1,223 +1,0 @@
-# include "psastroInternal.h"
-
-# define REQUIRED_RECIPE_VALUE(VALUE, NAME, TYPE)\
-  VALUE = psMetadataLookup##TYPE (&status, recipe, NAME); \
-  if (!status) { \
-   psAbort ("Failed to find %s in recipe", NAME); }
-
-bool psastroOneChip (pmFPA *fpa, pmChip *chip, psArray *refstars, psArray *rawstars, psMetadata *recipe, psMetadata *updates) {
-
-    bool status;
-    pmAstromStats *stats = NULL;
-
-    // default value for match/fit : radius is in pixels
-    REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32); 
-
-    // run the match/fit sequence NITER times
-    REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32); 
-
-    // correct radius to FP units (physical pixel scale in microns per pixel)
-    REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32); 
-    RADIUS *= pixelScale;
-
-    // select the desired chip order
-    REQUIRED_RECIPE_VALUE (int order, "PSASTRO.CHIP.ORDER", S32);
-
-    // allowed limits for valid solutions
-    REQUIRED_RECIPE_VALUE (float maxError, "PSASTRO.MAX.ERROR", F32);
-    REQUIRED_RECIPE_VALUE (int minNstar, "PSASTRO.MIN.NSTAR", S32);
-
-    // do we need to get a rough initial match?
-    REQUIRED_RECIPE_VALUE (bool gridSearch, "PSASTRO.GRID.SEARCH", Bool);
-
-    // do we need to get a rough initial match?
-    REQUIRED_RECIPE_VALUE (int maxNstar, "PSASTRO.GRID.NSTAR.MAX", S32);
-
-    if (gridSearch) {
-	// generate the bright subset of maxNstar entries (note: rawstars is sorted by S/N)
-	psArray *subset = psArrayAlloc (PS_MIN (maxNstar, rawstars->n));
-	for (int i = 0; (i < maxNstar) && (i < rawstars->n); i++) {
-	    subset->data[i] = psMemIncrRefCounter (rawstars->data[i]);
-	}
-
-	// XXX set clump scale from recipe
-	psArray *gridStars = psastroRemoveClumps (subset, 150);
-	psFree (subset);
-
-	psArray *refSubset = psastroRemoveClumps (refstars, 150);
-
-	psLogMsg ("psastro", 3, "grid search using %ld raw vs %ld ref stars\n", gridStars->n, refSubset->n);
-
-	// find initial offset / rotation / scale
-	pmAstromStats *gridStats = pmAstromGridMatch (gridStars, refSubset, recipe);
-	if (gridStats == NULL) {
-	    psLogMsg ("psastro", 3, "failed to find a grid match solution\n");
-	    psFree (gridStars);
-	    psFree (refSubset);
-	    return false;
-	}
-	psLogMsg ("psastro", 3, "basic grid search result - offset: %f,%f pixels, rotation: %f deg\n", gridStats->offset.x, gridStats->offset.y, DEG_RAD*gridStats->angle);
-
-	// tweak the position by finding peak of matches stars
-	stats = pmAstromGridTweak (gridStars, refSubset, recipe, gridStats);
-	if (stats == NULL) {
-	    psLogMsg ("psastro", 3, "failed to measure tweaked grid solution\n");
-	    psFree (gridStats);
-	    psFree (gridStars);
-	    psFree (refSubset);
-	    return false;
-	}
-	psLogMsg ("psastro", 3, "tweak grid search result - offset: %f,%f pixels, rotation: %f deg\n", stats->offset.x, stats->offset.y, DEG_RAD*stats->angle);
-
-	// adjust the chip.toFPA terms only
-	pmAstromGridApply (chip->toFPA, stats);
-	psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
-	psFree (gridStats);
-	psFree (gridStars);
-	psFree (refSubset);
-    }
-
-    psArray *match = NULL;
-    psStats *fitStats = NULL;
-    pmAstromFitResults *results = NULL;
-
-    for (int iter = 0; iter < nIter; iter++) {
-	
-	char name[128];
-
-	sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter);
-	float radius = psMetadataLookupF32 (&status, recipe, name);
-	radius *= pixelScale;
-	if (!status || (radius == 0.0)) {
-	    radius = RADIUS;
-	}
-
-
-	// use small radius to match stars
-	match = pmAstromRadiusMatchFP (rawstars, refstars, radius);
-	if (match == NULL) {
-	    psLogMsg ("psastro", 3, "failed to find radius-matched sources\n");
-	    psFree (stats);
-	    return false;
-	}
-
-	// modify the order to correspond to the actual number of matched stars:
-	if ((match->n < 11) && (order >= 3)) order = 2;
-	if ((match->n <  7) && (order >= 2)) order = 1;
-	if ((match->n <  4) && (order >= 1)) order = 0;
-	if (order < 1) {
-	    psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n); 
-	    psFree (match);
-	    psFree (stats);
-	    return false; 
-	} 
-
-	// create output toFPA; set masks appropriate to the Elixir DVO astrometry format
-	psFree (chip->toFPA);
-	chip->toFPA = psPlaneTransformAlloc (order, order);
-	for (int i = 0; i <= chip->toFPA->x->nX; i++) {
-	    for (int j = 0; j <= chip->toFPA->x->nY; j++) {
-		if (i + j > order) {
-		    chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET;
-		    chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET;
-		}
-	    }
-	}
-
-	// XXX allow statitic to be set by the user
-	// fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-	fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-	fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA");
-	fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER");
-
-	// improved fit for astrometric terms
-	results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
-	if (!results) {
-	    psLogMsg ("psastro", 3, "failed to perform the matched fit\n");
-	    psFree (match);
-	    psFree (stats);
-	    psFree (fitStats);
-	    return false;
-	}
-    
-	// determine fromFPA transformation and apply new transformation to raw & ref stars
-	psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
-    
-	if (iter < nIter - 1) {
-	    psFree (fitStats);
-	    psFree (results);
-	    psFree (match);
-	}
-
-	// toSky converts from FPA & TPA units (microns) to sky units (radians)
-	float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
-	// float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
-	float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
-	int astNstar = results->yStats->clippedNvalues;
-	psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
-    }
-
-
-
-    // toSky converts from FPA & TPA units (microns) to sky units (radians)
-    float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
-
-    // pixError is the average 1D scatter in pixels ('results' are in FPA units = microns)
-    // float pixError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale;
-    float pixError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) / pixelScale;
-
-    // astError is the average 1D scatter in arcsec ('results' are in FPA units = microns)
-    // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
-    float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
-    int astNstar = results->yStats->clippedNvalues;
-
-    bool validSolution = true;
-
-    // XXX should these result in errors or be handled another way?
-    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d", astError, astNstar);
-    if (astError > maxError) {
-        psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
-	validSolution = false;
-    }
-    if (astNstar < minNstar) {
-        psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar);
-	validSolution = false;
-    }
-
-    // DVO expects NASTRO = 0 if we fail to find a solution
-    psMetadataAddF32 (updates, PS_LIST_TAIL, "PERROR",   PS_META_REPLACE, "astrometry error (pixels)", pixError);
-    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
-    if (validSolution) {
-	psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
-	psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
-    } else {
-	psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
-	psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
-    }
-    psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
-
-    // XXX drop from here : determine fromFPA transformation and apply new transformation to raw & ref stars
-    // psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
-    
-    // XXX check if we correctly applied the new transformation:
-    if (psTraceGetLevel("psastro.dump") > 0) {
-	psastroDumpRawstars (rawstars, fpa, chip);
-	psastroDumpMatchedStars ("match.dat", rawstars, refstars, match);
-	psastroDumpStars (refstars, "refstars.cal.dat");
-    }
-
-    if (psTraceGetLevel("psastro.plot") > 0) {
-	psastroPlotOneChipFit (rawstars, refstars, match, recipe);
-    }
-
-    psFree (match);
-    psFree (stats);
-    psFree (results);
-    psFree (fitStats);
-    return validSolution;
-}
-
-// psastroWriteStars ("raw.1.dat", rawstars);
-// psastroWriteStars ("ref.1.dat", refstars);
-// psastroWriteTransform (chip->toFPA);
-
Index: branches/cnb_branch_20080830/psastro/src/psastroOneChipGrid.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroOneChipGrid.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroOneChipGrid.c	(revision 20033)
@@ -21,5 +21,5 @@
     psArray *subset = psArrayAlloc (PS_MIN (maxNstar, rawstars->n));
     for (int i = 0; (i < maxNstar) && (i < rawstars->n); i++) {
-	subset->data[i] = psMemIncrRefCounter (rawstars->data[i]);
+        subset->data[i] = psMemIncrRefCounter (rawstars->data[i]);
     }
 
@@ -35,8 +35,8 @@
     pmAstromStats *gridStats = pmAstromGridMatch (gridStars, refSubset, recipe);
     if (gridStats == NULL) {
-	psLogMsg ("psastro", 3, "failed to find a grid match solution\n");
-	psFree (gridStars);
-	psFree (refSubset);
-	return false;
+        psLogMsg ("psastro", 3, "failed to find a grid match solution\n");
+        psFree (gridStars);
+        psFree (refSubset);
+        return false;
     }
     psLogMsg ("psastro", 3, "basic grid search result - offset: %f,%f pixels, rotation: %f deg\n", gridStats->offset.x, gridStats->offset.y, DEG_RAD*gridStats->angle);
@@ -45,9 +45,9 @@
     stats = pmAstromGridTweak (gridStars, refSubset, recipe, gridStats);
     if (stats == NULL) {
-	psLogMsg ("psastro", 3, "failed to measure tweaked grid solution\n");
-	psFree (gridStats);
-	psFree (gridStars);
-	psFree (refSubset);
-	return false;
+        psLogMsg ("psastro", 3, "failed to measure tweaked grid solution\n");
+        psFree (gridStats);
+        psFree (gridStars);
+        psFree (refSubset);
+        return false;
     }
     psLogMsg ("psastro", 3, "tweak grid search result - offset: %f,%f pixels, rotation: %f deg\n", stats->offset.x, stats->offset.y, DEG_RAD*stats->angle);
Index: branches/cnb_branch_20080830/psastro/src/psastroParseCamera.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroParseCamera.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroParseCamera.c	(revision 20033)
@@ -33,4 +33,13 @@
     psFree (chips);
 
+    psString dump_file = psMetadataLookupStr(&status, config->arguments, "DUMP_CONFIG");
+    if (dump_file) {
+        pmConfigCamerasCull(config, NULL);
+        pmConfigRecipesCull(config, "PPIMAGE,PPSTATS,PSPHOT,MASKS,PSASTRO");
+
+        pmConfigDump(config, input->fpa, dump_file);
+    }
+
+
     psTrace("psastro", 1, "Done with psastroParseCamera...\n");
     return true;
Index: branches/cnb_branch_20080830/psastro/src/psastroRefstarSubset.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroRefstarSubset.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroRefstarSubset.c	(revision 20033)
@@ -21,5 +21,5 @@
   // is needed...
   psLogMsg ("psastro", 4, "measuring luminosity function for rawstars\n");
-  pmLumFunc *rawfunc = psastroLuminosityFunction (rawstars);
+  pmLumFunc *rawfunc = psastroLuminosityFunction (rawstars, NULL);
   if (rawfunc == NULL) {
     psLogMsg ("psastro", 4, "giving up on rawstars for this readout\n");
@@ -27,5 +27,5 @@
   }
   psLogMsg ("psastro", 4, "measuring luminosity function for refstars\n");
-  pmLumFunc *reffunc = psastroLuminosityFunction (refstars);
+  pmLumFunc *reffunc = psastroLuminosityFunction (refstars, rawfunc);
   if (reffunc == NULL) {
     psLogMsg ("psastro", 4, "giving up on refstars for this readout\n");
@@ -84,7 +84,7 @@
 /* this test is a bit sensitive to the total number of refstars or rawstars available
    watch out if:
-   - the fitted slopes are extremely different 
+   - the fitted slopes are extremely different
    - the average number of stars per bin is ~1
-   
+
    skip the cut in these cases?
 */
Index: branches/cnb_branch_20080830/psastro/src/psastroStandAlone.h
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroStandAlone.h	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroStandAlone.h	(revision 20033)
@@ -25,5 +25,5 @@
 bool 		  psastroModelDataSave (pmConfig *config);
 
-psVector         *psastroModelFitBoresite (psVector *Xo, psVector *Yo, psVector *Po);
+psVector         *psastroModelFitBoresite (psVector *Xo, psVector *Yo, psVector *Po, char *outroot);
 psF32 		  psastroModelBoresite (psVector *deriv, const psVector *params, const psVector *coord);
 
Index: branches/cnb_branch_20080830/psastro/src/psastroUseModel.c
===================================================================
--- branches/cnb_branch_20080830/psastro/src/psastroUseModel.c	(revision 20026)
+++ branches/cnb_branch_20080830/psastro/src/psastroUseModel.c	(revision 20033)
@@ -1,4 +1,5 @@
 # include "psastroInternal.h"
 # define NONLIN_TOL 0.001 /* tolerance in pixels */
+# define DEBUG 0
 
 // apply the generic astrometry model to this image.  this assumes the WCS terms either do not
@@ -35,6 +36,61 @@
   }
 
+  // XXX TEST: apply the input image RA & DEC to the astrometry model, save corners
+  if (0) {
+      // these externally supplied values are used to set the final transformation terms
+      double RA  = psMetadataLookupF64 (&status, input->fpa->concepts, "FPA.RA");
+      double DEC = psMetadataLookupF64 (&status, input->fpa->concepts, "FPA.DEC");
+      // double POS = PM_RAD_DEG * psMetadataLookupF64 (&status, input->fpa->concepts, "FPA.POSANGLE");
+
+      // get projection scale; center is supplied
+      // XXX should this be astrom or input??
+      float Xs = psMetadataLookupF32(&status, astrom->fpa->concepts, "XSCALE") * PM_RAD_DEG;
+      float Ys = psMetadataLookupF32(&status, astrom->fpa->concepts, "YSCALE") * PM_RAD_DEG;
+
+      // allocate a new toSky projection using the reported position
+      psFree (astrom->fpa->toSky);
+      astrom->fpa->toSky = psProjectionAlloc (RA, DEC, Xs, Ys, PS_PROJ_DIS);
+
+      // local view
+      pmFPAview *myView = pmFPAviewAlloc (0);
+
+      // loop over all chips, replace input astrometry elements with those from astrom 
+      pmChip *obsChip = NULL;
+      while ((obsChip = pmFPAviewNextChip (myView, input->fpa, 1)) != NULL) {
+	  psTrace ("psastro", 4, "Chip %d: %x %x\n", myView->chip, obsChip->file_exists, obsChip->process);
+	  if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
+
+	  // set the chip astrometry using the astrom file
+	  pmChip *refChip = pmFPAviewThisChip (myView, astrom->fpa);
+
+	  psFree (obsChip->toFPA);
+	  psFree (obsChip->fromFPA);
+
+	  // supply astrometry from model 
+	  obsChip->toFPA   = psMemIncrRefCounter (refChip->toFPA);
+	  obsChip->fromFPA = psMemIncrRefCounter (refChip->fromFPA);
+
+	  // XXX if we want to write out the result, update the header here.  this needs to be
+	  // updated with the correct HDU selection.  obsChip->hdu may not exist.
+	  // pmAstromWriteBilevelChip (obsChip->hdu->header, obsChip, NONLIN_TOL);
+      }
+
+      psFree (input->fpa->toSky);
+      psFree (input->fpa->toTPA);
+      psFree (input->fpa->fromTPA);
+      input->fpa->toSky   = psMemIncrRefCounter (astrom->fpa->toSky);
+      input->fpa->toTPA   = psMemIncrRefCounter (astrom->fpa->toTPA);
+      input->fpa->fromTPA = psMemIncrRefCounter (astrom->fpa->fromTPA);
+
+      // XXX this is temporarily hardwired because of model error
+      input->fpa->toSky->type = PS_PROJ_TAN;
+
+      if (DEBUG) psastroDumpCorners ("corners.up.ast1.dat", "corners.dn.ast1.dat", input->fpa);
+  }
+
   // set the model using the RA, DEC, POSANGLE of the input image.
   pmAstromModelSetTP (astrom, input->fpa->concepts);
+
+  if (DEBUG) psastroDumpCorners ("corners.up.ast2.dat", "corners.dn.ast2.dat", astrom->fpa);
 
   // loop over all chips, replace input astrometry elements with those from astrom 
@@ -57,7 +113,4 @@
     // updated with the correct HDU selection.  obsChip->hdu may not exist.
     // pmAstromWriteBilevelChip (obsChip->hdu->header, obsChip, NONLIN_TOL);
-
-    
-
   }
 
@@ -72,5 +125,5 @@
   input->fpa->toSky->type = PS_PROJ_TAN;
 
-  // psastroDumpCorners ("corners.useModel.dat", input->fpa);
+  if (DEBUG) psastroDumpCorners ("corners.up.inp.dat", "corners.dn.inp.dat", input->fpa);
 
   // updated with the correct HDU selection.  obsChip->hdu may not exist.