Changeset 20033

Timestamp:

Oct 9, 2008, 2:39:01 PM (18 years ago)

Author:

beaumont

Message:

merged mainline into branch. resolved conflicts. added plots.

Location:

branches/cnb_branch_20080830/psastro/src

Files:

• . (modified) (1 prop)
• .cvsignore (modified) (1 diff)
• Makefile.am (modified) (5 diffs)
• psastro.c (modified) (3 diffs)
• psastro.h (modified) (4 diffs)
• psastroAnalysis.c (modified) (2 diffs)
• psastroArguments.c (modified) (1 diff)
• psastroAstromGuess.c (modified) (11 diffs)
• psastroChipAstrom.c (modified) (1 diff)
• psastroChooseRefstars.c (modified) (3 diffs)
• psastroDemoDump.c (modified) (4 diffs)
• psastroDemoPlot.c (modified) (3 diffs)
• psastroEnforceChips.c (deleted)
• psastroModelAdjust.c (modified) (4 diffs)
• psastroModelAnalysis.c (modified) (9 diffs)
• psastroModelBoresite.c (modified) (1 diff)
• psastroModelDataSave.c (modified) (3 diffs)
• psastroModelFit.c (added)
• psastroModelFitBoresite.c (modified) (3 diffs)
• psastroMosaicAstrom.c (modified) (3 diffs)
• psastroMosaicCorrectDistortion.c (added)
• psastroMosaicDistortion.c (added)
• psastroMosaicFPtoTP.c (added)
• psastroMosaicGradients.c (modified) (4 diffs)
• psastroMosaicOneChip.c (modified) (3 diffs)
• psastroOneChip.c (deleted)
• psastroOneChipGrid.c (modified) (3 diffs)
• psastroParseCamera.c (modified) (1 diff)
• psastroRefstarSubset.c (modified) (3 diffs)
• psastroStandAlone.h (modified) (1 diff)
• psastroUseModel.c (modified) (4 diffs)

branches/cnb_branch_20080830/psastro/src

@@ -15 +15 @@
 psastroModel
 gpcModel
+psastroModelFit

@@ -15 +15 @@
 psastroModel
 gpcModel
+psastroModelFit

branches/cnb_branch_20080830/psastro/src/.cvsignore

@@ -15 +15 @@
 psastroModel
 gpcModel
+psastroModelFit

branches/cnb_branch_20080830/psastro/src/Makefile.am

@@ -1 +1 @@
 
 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 +43 @@
         psastroCleanup.c
 
+psastroModelFit_SOURCES = \
+        psastroModelFit.c \
+        psastroModelBoresite.c      \
+        psastroModelFitBoresite.c
+
 gpcModel_SOURCES = gpcModel.c
 
@@ -50 +60 @@
         psastroConvert.c            \
         psastroChipAstrom.c         \
-        psastroOneChip.c            \
         psastroOneChipGrid.c        \
         psastroOneChipFit.c         \
@@ -57 +66 @@
         psastroTestFuncs.c          \
         psastroLuminosityFunction.c \
+        psastroLuminosityFunctionPlot.c \
         psastroRefstarSubset.c      \
         psastroFixChips.c           \
@@ -62 +72 @@
         psastroUseModel.c           \
         psastroMosaicAstrom.c       \
+        psastroMosaicDistortion.c   \
+        psastroMosaicFPtoTP.c       \
         psastroMosaicGradients.c    \
+        psastroMosaicCorrectDistortion.c   \
         psastroMosaicChipAstrom.c   \
         psastroMosaicOneChip.c      \

branches/cnb_branch_20080830/psastro/src/psastro.c

@@ -9 +9 @@
 
     pmConfig *config = NULL;
-
     psTimerStart ("complete");
 
@@ -24 +23 @@
     // 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 +30 @@
     // 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);
     }
 

branches/cnb_branch_20080830/psastro/src/psastro.h

@@ -17 +17 @@
 // 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 +35 @@
 psArray          *pmSourceToAstromObj (psArray *sources);
 bool              psastroAstromGuess (int *nStars, pmConfig *config);
+bool              psastroAstromGuessCheck (pmConfig *config);
 
 psPlaneDistort   *psPlaneDistortIdentity ();
@@ -46 +47 @@
 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 +57 @@
 
 // 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);
 
 

branches/cnb_branch_20080830/psastro/src/psastroAnalysis.c

@@ -63 +63 @@
         chipastro = true;
     }
-
     if (chipastro) {
         if (!psastroChipAstrom (config)) {
@@ -81 +80 @@
     // psastroZeroPoint (config);
 
+    psastroAstromGuessCheck (config);
+
     // XXX how do we specify stack astrometry?
     // psastroStackAstrom (config, refs);

branches/cnb_branch_20080830/psastro/src/psastroArguments.c

@@ -81 +81 @@
     }
 
+    // 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) {

branches/cnb_branch_20080830/psastro/src/psastroAstromGuess.c

@@ -1 +1 @@
 # include "psastroInternal.h"
+# define DEBUG 0
 
 // this function loads the header WCS astrometry terms into the fpa terms and applies the
@@ -28 +29 @@
     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 +36 @@
     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 +42 @@
     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 +49 @@
     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 +75 @@
         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 +116 @@
                 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 +138 @@
                 }
 
-                // 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 +168 @@
     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 +201 @@
     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 +223 @@
     // 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;
+}

branches/cnb_branch_20080830/psastro/src/psastroChipAstrom.c

@@ -115 +115 @@
     }
 
-    // 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);

branches/cnb_branch_20080830/psastro/src/psastroChooseRefstars.c

@@ -51 +51 @@
         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 +62 @@
                 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 +98 @@
                     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;

branches/cnb_branch_20080830/psastro/src/psastroDemoDump.c

@@ -128 +128 @@
         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 +155 @@
 
                     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 +195 @@
 }
 
-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 +217 @@
         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;

branches/cnb_branch_20080830/psastro/src/psastroDemoPlot.c

@@ -101 +101 @@
     }
     xVec->n = yVec->n = zVec->n = n;
+
     pmKapaPlotVectorTriple_AutoLimits_OpenGraph (kapa, &graphdata, xVec, yVec, zVec, false);
 
@@ -124 +125 @@
     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 +194 @@
     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), ??

branches/cnb_branch_20080830/psastro/src/psastroModelAdjust.c

@@ -1 +1 @@
 # 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 +16 @@
     }
 
+    // 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 +46 @@
     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 +147 @@
         // 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;
 }

branches/cnb_branch_20080830/psastro/src/psastroModelAnalysis.c

@@ -1 +1 @@
 # include "psastroStandAlone.h"
+# define NONLIN_TOL 0.001
 
 bool psastroModelAnalysis (pmConfig *config) {
@@ -12 +13 @@
     }
 
-    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 +19 @@
         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 +52 @@
      * T_0 : reference angle for rotator
      * P_0 : orientation of boresite ellipse
+     *
      */
 
@@ -58 +76 @@
     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 +91 @@
         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 +99 @@
         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 +107 @@
 
         // 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 +143 @@
     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]); 

branches/cnb_branch_20080830/psastro/src/psastroModelBoresite.c

@@ -45 +45 @@
             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);

branches/cnb_branch_20080830/psastro/src/psastroModelDataSave.c

@@ -1 +1 @@
 # include "psastroStandAlone.h"
+# define NONLIN_TOL 0.001 /* tolerance in pixels */
 
 # define ESCAPE { \
@@ -22 +23 @@
 
     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 +36 @@
     return true;
 }
+
+
+

branches/cnb_branch_20080830/psastro/src/psastroModelFitBoresite.c

@@ -2 +2 @@
 
 // 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 +60 @@
     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 +74 @@
     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;

branches/cnb_branch_20080830/psastro/src/psastroMosaicAstrom.c

@@ -1 +1 @@
 # 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 +39 @@
     # 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 +89 @@
  * 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;
+}

branches/cnb_branch_20080830/psastro/src/psastroMosaicGradients.c

@@ -2 +2 @@
 static int nPass = 0;
 
-bool psastroMosaicGradients (pmFPA *fpa, psMetadata *recipe) {
+bool psastroMosaicDistortionFromGradients (pmFPA *fpa, psMetadata *recipe) {
 
     bool status;
@@ -9 +9 @@
     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 +62 @@
     }
 
+    // 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 +100 @@
     }
         
-    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;
 }
+

branches/cnb_branch_20080830/psastro/src/psastroMosaicOneChip.c

@@ -12 +12 @@
     char errorWord[64];
     char orderWord[64];
-
-    assert (iteration < 4);
 
     PS_ASSERT_PTR_NON_NULL(chip,    false);
@@ -86 +84 @@
 
     // 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 +123 @@
     // 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;

branches/cnb_branch_20080830/psastro/src/psastroOneChipGrid.c

@@ -21 +21 @@
     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 +35 @@
     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 +45 @@
     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);

branches/cnb_branch_20080830/psastro/src/psastroParseCamera.c

@@ -33 +33 @@
     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;

branches/cnb_branch_20080830/psastro/src/psastroRefstarSubset.c

@@ -21 +21 @@
   // 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 +27 @@
   }
   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 +84 @@
 /* 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?
 */

branches/cnb_branch_20080830/psastro/src/psastroStandAlone.h

@@ -25 +25 @@
 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);
 

branches/cnb_branch_20080830/psastro/src/psastroUseModel.c

@@ -1 +1 @@
 # 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 +36 @@
   }
 
+  // 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 +113 @@
     // updated with the correct HDU selection.  obsChip->hdu may not exist.
     // pmAstromWriteBilevelChip (obsChip->hdu->header, obsChip, NONLIN_TOL);
-
-    
-
   }
 
@@ -72 +125 @@
   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.

@@ -15 +15 @@
 psastroModel
 gpcModel
+psastroModelFit