Index: /trunk/psastro/src/psastroAnalysis.c
===================================================================
--- /trunk/psastro/src/psastroAnalysis.c	(revision 19512)
+++ /trunk/psastro/src/psastroAnalysis.c	(revision 19513)
@@ -81,4 +81,6 @@
     // psastroZeroPoint (config);
 
+    psastroAstromGuessCheck (config);
+
     // XXX how do we specify stack astrometry?
     // psastroStackAstrom (config, refs);
Index: /trunk/psastro/src/psastroAstromGuess.c
===================================================================
--- /trunk/psastro/src/psastroAstromGuess.c	(revision 19512)
+++ /trunk/psastro/src/psastroAstromGuess.c	(revision 19513)
@@ -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
@@ -52,5 +53,14 @@
     } 
 
+    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?
@@ -77,4 +87,24 @@
         }
 
+	// 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);
+	}
+
         // apply the new WCS guess data to all of the data in the readouts
         while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
@@ -85,17 +115,4 @@
             while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
                 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");
@@ -132,4 +149,6 @@
         }
     }
+
+    if (DEBUG) psastroDumpCorners ("corners.up.guess2.dat", "corners.dn.guess2.dat", fpa);
 
     // how many total sources are available to us?
@@ -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);
@@ -201,2 +234,151 @@
     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;
+    }
+
+    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);
+
+    pmFPA *fpa = input->fpa;
+
+    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);
+    }
+
+    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");
+
+    // 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;
+}