Changeset 25027 for branches/pap/psastro

Timestamp:

Aug 7, 2009, 4:08:25 PM (17 years ago)

Author:

Paul Price

Message:

Merging trunk (r25026) to get up-to-date on old branch.

Location:

branches/pap

Files:

• . (modified) (1 prop)
• psastro (modified) (1 prop)
• psastro/src/Makefile.am (modified) (1 diff)
• psastro/src/psastro.h (modified) (3 diffs)
• psastro/src/psastroAnalysis.c (modified) (1 diff)
• psastro/src/psastroAstromGuess.c (modified) (1 diff)
• psastro/src/psastroChipAstrom.c (modified) (1 diff)
• psastro/src/psastroChooseGlintStars.c (copied) (copied from trunk/psastro/src/psastroChooseGlintStars.c )
• psastro/src/psastroChooseRefstars.c (modified) (3 diffs)
• psastro/src/psastroExtractGhosts.c (modified) (1 diff)
• psastro/src/psastroExtractStars.c (modified) (1 diff)
• psastro/src/psastroFindChip.c (modified) (1 diff, 1 prop)
• psastro/src/psastroFixChipsTest.c (modified) (3 diffs)
• psastro/src/psastroLoadGhosts.c (modified) (4 diffs)
• psastro/src/psastroLoadGlints.c (copied) (copied from trunk/psastro/src/psastroLoadGlints.c )
• psastro/src/psastroLoadRefstars.c (modified) (3 diffs)
• psastro/src/psastroMaskUpdates.c (modified) (5 diffs)
• psastro/src/psastroMaskUtils.c (modified) (1 diff)
• psastro/src/psastroModelAnalysis.c (modified) (1 diff)
• psastro/src/psastroModelFitBoresite.c (modified) (1 diff)
• psastro/src/psastroMosaicOneChip.c (modified) (1 diff)
• psastro/src/psastroOneChipFit.c (modified) (1 diff)
• psastro/src/psastroZeroPoint.c (modified) (5 diffs)

branches/pap/psastro/src/Makefile.am

@@ -89 +89 @@
         psastroMaskUpdates.c        \
         psastroMaskUtils.c          \
+        psastroChooseGlintStars.c   \
+        psastroLoadGlints.c   \
         psastroLoadGhosts.c         \
         psastroGhostUtils.c         \

branches/pap/psastro/src/psastro.h

@@ -73 +73 @@
 bool              psastroMaskUpdates (pmConfig *config);
 
+bool              psastroLoadGlints (pmConfig *config);
+bool              psastroChooseGlintStars (pmConfig *config, psArray *refs, const char *source);
+
 psArray          *psastroLoadRefstars (pmConfig *config, const char *source);
 bool              psastroChipAstrom (pmConfig *config);
@@ -144 +147 @@
 bool              psastroZeroPoint (pmConfig *config);
 bool              psastroZeroPointReadout(pmReadout *readout, float zeropt, float exptime);
-bool              psastroZeroPointFromRecipe (float *zeropt, float *exptime, pmFPA *fpa, psMetadata *recipe);
+bool              psastroZeroPointFromRecipe (float *zeropt, float *exptime, float *ghostMaxMag, pmFPA *fpa, psMetadata *recipe);
 
 // masking functions
@@ -162 +165 @@
 pmChip           *psastroFindChip (double *xChip, double *yChip, pmFPA *fpa, double xFPA, double yFPA);
 bool              psastroChipBounds (pmFPA *fpa);
+bool              psastroFindChipInXrange (pmFPA *fpa, int nChip, double xFPA, double yFPA);
+bool              psastroFindChipInYrange (pmFPA *fpa, int nChip, double xFPA, double yFPA);
+bool              psastroFindChipYedges (double *yFPAs, double *yFPAe, pmFPA *fpa, int nChip);
+bool              psastroFindChipXedges (double *yFPAs, double *yFPAe, pmFPA *fpa, int nChip);
+bool              psastroFPAtoChip (double *xChip, double *yChip, pmFPA *fpa, int nChip, double xFPA, double yFPA);
+
 
 // psastroExtract functions

branches/pap/psastro/src/psastroAnalysis.c

@@ -86 +86 @@
         return false;
     }
+
+    if (!psastroChooseGlintStars (config, refs, "PSASTRO.INPUT")) {
+        psError (PSASTRO_ERR_UNKNOWN, false, "failed to select glint stars from reference star list\n");
+        psFree(refs);
+        return false;
+    }
     psFree (refs);  // refs of interest are saved on readout->analysis
 

branches/pap/psastro/src/psastroAstromGuess.c

@@ -382 +382 @@
     psStats *statsQ = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
 
-    psVectorStats (statsP, cornerPd, NULL, cornerMK, 1);
-    psVectorStats (statsQ, cornerQd, NULL, cornerMK, 1);
+    if (!psVectorStats (statsP, cornerPd, NULL, cornerMK, 1)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+        return false;
+    }
+    if (!psVectorStats (statsQ, cornerQd, NULL, cornerMK, 1)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+        return false;
+    }
 
     float angle = atan2 (map->y->coeff[1][0], map->x->coeff[1][0]);

branches/pap/psastro/src/psastroChipAstrom.c

@@ -130 +130 @@
     }
 
+    if (numGoodChips == 0) {
+        psError(PSASTRO_ERR_UNKNOWN, false, "Failed to fit any chips");
+        return false;
+    }
+
     if (!psastroFixChips (config, recipe)) {
         psError(PSASTRO_ERR_UNKNOWN, false, "failed to align problematic chips");

branches/pap/psastro/src/psastroChooseRefstars.c

@@ -6 +6 @@
  *
  *  @author IfA
- *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
+ *  @version $Revision: 1.21 $
  *  @date $Date: 2009-02-07 02:03:34 $
  *  Copyright 2009 Institute for Astronomy, University of Hawaii
@@ -49 +49 @@
     pmFPAview *view = pmFPAviewAlloc (0);
     pmFPA *fpa = input->fpa;
+
+    // XXX kind of a hack -- think this through a bit more clearly:
+    if (psMetadataLookupPtr (&status, fpa->analysis, "MATCHED_REFS")) {
+        // we loaded a set of matched references from an earlier astrometry
+        // analysis.  however, we are re-doing the astrometry here, so remove
+        // that prior set of matched references
+        psMetadataRemoveKey (fpa->analysis, "MATCHED_REFS");
+    }
 
     // sort by mag
@@ -114 +122 @@
                 psTrace ("psastro", 4, "Added %ld refstars\n", refstars->n);
 
+# if (0)
+                // XXX test
+                FILE *outfile = fopen ("refstars.dat", "w");
+                assert (outfile);
+                for (int nn = 0; nn < refstars->n; nn++) {
+                          pmAstromObj *ref = refstars->data[nn];
+                          fprintf (outfile, "%lf %lf\n", ref->sky->r*PS_DEG_RAD, ref->sky->d*PS_DEG_RAD);
+                }
+                fclose (outfile);
+# endif
+
                 psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars);
                 psFree (refstars);

branches/pap/psastro/src/psastroExtractGhosts.c

@@ -53 +53 @@
 
     // really error-out here?  or just skip?
-    if (!psastroZeroPointFromRecipe (&zeropt, &exptime, fpa, recipe)) {
+    if (!psastroZeroPointFromRecipe (&zeropt, &exptime, NULL, fpa, recipe)) {
         psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe");
         return false;

branches/pap/psastro/src/psastroExtractStars.c

@@ -65 +65 @@
 
     // really error-out here?  or just skip?
-    if (!psastroZeroPointFromRecipe (&zeropt, &exptime, fpa, recipe)) {
+    if (!psastroZeroPointFromRecipe (&zeropt, &exptime, NULL, fpa, recipe)) {
         psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe");
         return false;

branches/pap/psastro/src/psastroFindChip.c

@@ -110 +110 @@
 }
 
+// identify chips which land on this column (FP coords)
+bool psastroFindChipInXrange (pmFPA *fpa, int nChip, double xFPA, double yFPA) {
+
+    if (!chipXmin || !chipXmax) {
+        psAbort ("chip bounds not set");
+    }
+
+    if (xFPA <  chipXmin->data.F32[nChip]) return false;
+    if (xFPA >= chipXmax->data.F32[nChip]) return false;
+    return true;
+}
+
+// identify chips which land on this row (FP coords)
+bool psastroFindChipInYrange (pmFPA *fpa, int nChip, double xFPA, double yFPA) {
+
+    if (!chipYmin || !chipYmax) {
+        psAbort ("chip bounds not set");
+    }
+
+    if (yFPA <  chipYmin->data.F32[nChip]) return false;
+    if (yFPA >= chipYmax->data.F32[nChip]) return false;
+    return true;
+}
+
+// return the FPA coordinates of the Y edges of the chip
+bool psastroFindChipYedges (double *yFPAs, double *yFPAe, pmFPA *fpa, int nChip) {
+
+    *yFPAs = chipYmin->data.F32[nChip];
+    *yFPAe = chipYmax->data.F32[nChip];
+    return true;
+}
+
+// return the FPA coordinates of the X edges of the chip
+bool psastroFindChipXedges (double *yFPAs, double *yFPAe, pmFPA *fpa, int nChip) {
+
+    *yFPAs = chipXmin->data.F32[nChip];
+    *yFPAe = chipXmax->data.F32[nChip];
+    return true;
+}
+
+// convert FPA to Chip coordinates
+bool psastroFPAtoChip (double *xChip, double *yChip, pmFPA *fpa, int nChip, double xFPA, double yFPA) {
+
+    pmChip *chip = fpa->chips->data[nChip];
+
+    psPlane ptCH, ptFP;
+    ptFP.x = xFPA;
+    ptFP.y = yFPA;
+    psPlaneTransformApply (&ptCH, chip->fromFPA, &ptFP);
+
+    *xChip = ptCH.x;
+    *yChip = ptCH.y;
+    return true;
+}
+
 bool psastroExtractFreeChipBounds () {
   

branches/pap/psastro/src/psastroFixChipsTest.c

@@ -140 +140 @@
     psStats *stats;
     stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    psVectorStats (stats, dX, NULL, NULL, 0);
+    if (!psVectorStats (stats, dX, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+        return false;
+    }
     Xo = stats->robustMedian;
     fprintf (stderr, "offset x: %f +/- %f\n", stats->robustMedian, stats->robustStdev);
@@ -146 +149 @@
 
     stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    psVectorStats (stats, dY, NULL, NULL, 0);
+    if (!psVectorStats (stats, dY, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+        return false;
+    }
     Yo = stats->robustMedian;
     fprintf (stderr, "offset y: %f +/- %f\n", stats->robustMedian, stats->robustStdev);
@@ -168 +174 @@
     }
     stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    psVectorStats (stats, dT, NULL, NULL, 0);
+    if (!psVectorStats (stats, dT, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+        return false;
+    }
     float To = stats->robustMedian;
     fprintf (stderr, "offset t: %f +/- %f\n", stats->robustMedian, stats->robustStdev);

branches/pap/psastro/src/psastroLoadGhosts.c

@@ -46 +46 @@
     pmCell *cell = NULL;
     pmReadout *readout = NULL;
-    float zeropt, exptime;
+    float zeropt, exptime, MAX_MAG;
     psMetadata *md = NULL;
     psPolynomial2D *centerX = NULL;
@@ -71 +71 @@
     if (!strcasecmp(ghostFile, "NONE")) return true;
 
-    double MAX_MAG = psMetadataLookupF32 (&status, recipe, "GHOST_MAX_MAG");
-
     if (!pmConfigFileRead (&ghostModel, ghostFile, "GHOST MODEL")) {
         psError(PSASTRO_ERR_CONFIG, true, "Trouble loading ghost model");
@@ -97 +95 @@
 
     // really error-out here?  or just skip?
-    if (!psastroZeroPointFromRecipe (&zeropt, &exptime, fpa, recipe)) {
+    if (!psastroZeroPointFromRecipe (&zeropt, &exptime, &MAX_MAG, fpa, recipe)) {
         psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe");
         goto escape;
@@ -162 +160 @@
                     // which encompass the full size of the ghost
                     pmChip *ghostChip = psastroFindChip (&ghost->chip->x, &ghost->chip->y, fpa, -ghost->srcFP->x, -ghost->srcFP->y);
-                    fprintf (stderr, "raw chip position: %f, %f ", ghost->chip->x, ghost->chip->y);
+                    // fprintf (stderr, "raw chip position: %f, %f ", ghost->chip->x, ghost->chip->y);
 
                     ghostChip = psastroFindChip (&ghost->chip->x, &ghost->chip->y, fpa, ghost->FP->x, ghost->FP->y);
-                    fprintf (stderr, "-> model chip position: %f, %f\n", ghost->chip->x, ghost->chip->y);
+                    // fprintf (stderr, "-> model chip position: %f, %f\n", ghost->chip->x, ghost->chip->y);
 
                     if (!ghostChip) goto skip;

branches/pap/psastro/src/psastroLoadRefstars.c

@@ -128 +128 @@
     if (!strcmp (outformat, "PS1_DEV_0")) {
       refstars = psastroReadGetstar_PS1_DEV_0 (fits);
+
+      // XXX test
+      // FILE *outfile = fopen ("refstars.dat", "w");
+      // assert (outfile);
+      // for (int nn = 0; nn < refstars->n; nn++) {
+      //          pmAstromObj *ref = refstars->data[nn];
+      //          fprintf (outfile, "%lf %lf\n", ref->sky->r*PS_DEG_RAD, ref->sky->d*PS_DEG_RAD);
+      // }
+      // fclose (outfile);
     }
     if (refstars == NULL) {
@@ -181 +190 @@
             ref->sky->d   = RAD_DEG*psMetadataLookupF32 (&status, row, "DEC");
             ref->Mag      = 0.001*psMetadataLookupS32 (&status, row, "MAG");  // ELIXIR uses millimags
+            ref->Color    = 0.0;
         } else {
             ref->sky->r   = RAD_DEG*psMetadataLookupF64 (&status, row, "RA");
             ref->sky->d   = RAD_DEG*psMetadataLookupF64 (&status, row, "DEC");
             ref->Mag      = psMetadataLookupF32 (&status, row, "MAG"); // PANSTARRS uses mags
+            ref->Color    = 0.0;
         }
 
@@ -216 +227 @@
         psMetadata *row = psFitsReadTableRow(fits, i); // Table row
 
-        ref->sky->r   = RAD_DEG*psMetadataLookupF32 (&status, row, "RA");
-        ref->sky->d   = RAD_DEG*psMetadataLookupF32 (&status, row, "DEC");
+        ref->sky->r   = RAD_DEG*psMetadataLookupF64 (&status, row, "RA");
+        ref->sky->d   = RAD_DEG*psMetadataLookupF64 (&status, row, "DEC");
         ref->Mag      = psMetadataLookupF32 (&status, row, "MAG");
+        float MagC1   = psMetadataLookupF32 (&status, row, "MAG_C1");
+        float MagC2   = psMetadataLookupF32 (&status, row, "MAG_C2");
+        if (!isnan(MagC1) && !isnan(MagC2)) {
+            ref->Color = MagC1 - MagC2;
+        } else {
+            // XXX save the color and the slope in the table header?
+            ref->Color = 0.0;
+        }
 
         // XXX VERY temporary hack to avoid M31 bulge

branches/pap/psastro/src/psastroMaskUpdates.c

@@ -28 +28 @@
     pmCell *cell = NULL;
     pmReadout *readout = NULL;
-    float zeropt, exptime;
+    float zeropt, exptime, GHOST_MAX_MAG;
 
     psImageMaskType ghostMaskValue = pmConfigMaskGet("GHOST", config); // Mask value for ghost pixels
+    psImageMaskType glintMaskValue = pmConfigMaskGet("GHOST", config); // Mask value for glint pixels (overload ghost)
     psImageMaskType spikeMaskValue = pmConfigMaskGet("SPIKE", config); // Mask value for ghost pixels
     psImageMaskType starMaskValue  = pmConfigMaskGet("STARCORE", config); // Mask value for ghost pixels
@@ -53 +54 @@
     }
     bool COUNT_GHOSTS = psMetadataLookupF32 (&status, recipe, "REFSTAR_COUNT_GHOSTS");
-    double GHOST_MAX_MAG = psMetadataLookupF32 (&status, recipe, "GHOST_MAX_MAG");
     int nGhosts = 0;
+
+    // convert star positions to glint positions and add to the fpa->analysis data
+    if (!psastroLoadGlints (config)) {
+        psError(PSASTRO_ERR_CONFIG, false, "Error loading glints");
+        return false;
+    }
 
     psLogMsg ("psastro", PS_LOG_INFO, "generating a bright-star mask");
@@ -79 +85 @@
 
     // really error-out here?  or just skip?
-    if (!psastroZeroPointFromRecipe (&zeropt, &exptime, fpa, recipe)) {
+    if (!psastroZeroPointFromRecipe (&zeropt, &exptime, &GHOST_MAX_MAG, fpa, recipe)) {
         psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe");
         return false;
@@ -221 +227 @@
                     }
 
+                    // This masking option was needed for persistent charge trails in GPC1; it
+                    // has since been replaced with 'burntool', which is applied upon readout
+                    // by the camera software, and therefore is aware of the image sequence.
                     if (REFSTAR_MASK_BLEED) {
-
                         // convert x,y chip coordinates to cells in maskChip
                         pmCell *refCell = pmCellInChip (refChip, ref->chip->x, ref->chip->y);
@@ -245 +253 @@
                 }
 
-                // select the raw objects for this readout (loaded in psastroExtract.c)
+                // select the ghost object for this readout (loaded in psastroExtractGhosts.c).
+                // These differ from the reference stars since the star position is not
+                // contained by the readout; instead, the ghost position is predicted based on
+                // the ghost model, and the ghost positions associated with a given readout are
+                // supplied here.
                 psArray *ghosts = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GHOSTS");
-                if (ghosts == NULL) { continue; }
-
-                // identify the bright stars of interest
-                for (int i = 0; i < ghosts->n; i++) {
-                    psastroGhost *ghost = ghosts->data[i];
-                    // XXX bright vs faint ghost bits? (OR with SUSPECT) 
-                    psastroMaskEllipticalAnnulus (readoutMask->mask, ghostMaskValue, ghost->chip->x, ghost->chip->y, ghost->inner, ghost->outer);
-                }
-
-                // select the raw objects for this readout, flag is they fall in a mask
-                psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS");
-                if (rawstars == NULL) return false;
+                if (ghosts) { 
+                    // mask the ghosts on this readout
+                    for (int i = 0; i < ghosts->n; i++) {
+                        psastroGhost *ghost = ghosts->data[i];
+                        // XXX bright vs faint ghost bits? (OR with SUSPECT) 
+                        psastroMaskEllipticalAnnulus (readoutMask->mask, ghostMaskValue, ghost->chip->x, ghost->chip->y, ghost->inner, ghost->outer);
+                    }
+                }
+
+                // Select the glint mask regions for this readout (loaded in
+                // psastroExtractGlints.c).  These glint regions are defined as rectangular
+                // boxes and are generated for each chip based on the position of the bright
+                // stars beyond the edge of the focal plane.  This masking is currently very
+                // GPC1-specific
+                psArray *glints = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GLINTS");
+                if (glints) {
+                    // mask the glints on this readout
+                    for (int i = 0; i < glints->n; i++) {
+                        psRegion *glint = glints->data[i];
+                        psastroMaskRectangle (readoutMask->mask, glintMaskValue, glint->x0, glint->y0, glint->x1, glint->y1);
+                    }
+                }
+
+                // this probably should move into a function of its own:
+                {
+                    // select the raw objects for this readout, flag is they fall in a mask
+                    psArray *inSources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
+                    if (inSources == NULL) continue;
                 
-                // XXX finish this:
-                for (int i = 0; i < rawstars->n; i++) {
-                    pmAstromObj *raw = rawstars->data[i];
-                    psImageMaskType value = readoutMask->mask->data.PS_TYPE_IMAGE_MASK_DATA[(int)(raw->chip->x)][(int)(raw->chip->y)];
-                    if (value) continue;
+                    // create a replacement output array:
+                    // psArray *outSources = psAllocArrayEmpty(100);
+
+                    // XXX finish this: raise a bit for stars that land on certain types of masks;
+                    // others (eg, bright star core) should be ignored.
+                    for (int i = 0; i < inSources->n; i++) {
+                        pmSource *source = inSources->data[i];
+
+                        int xChip = source->peak->x;
+                        int yChip = source->peak->y;
+
+                        bool onChip = true;
+                        onChip &= (xChip >= 0);
+                        onChip &= (xChip < readoutMask->mask->numCols);
+                        onChip &= (yChip >= 0);
+                        onChip &= (yChip < readoutMask->mask->numRows);
+                        if (!onChip) {
+                            // if the source is off the edge of the chip, raise a different bit?
+                            source->mode |= PM_SOURCE_MODE_OFF_CHIP;
+                            continue;
+                        }
+
+                        psImageMaskType value = readoutMask->mask->data.PS_TYPE_IMAGE_MASK_DATA[yChip][xChip];
+                        if (value & ghostMaskValue) {
+                            source->mode |= PM_SOURCE_MODE_ON_GHOST;
+                        }
+                        // XXX note that for now, glint and ghost are identical
+                        pmSourceMode PM_SOURCE_MODE_ON_GLINT = PM_SOURCE_MODE_ON_GHOST;
+                        if (value & glintMaskValue) {
+                            source->mode |= PM_SOURCE_MODE_ON_GLINT;
+                        }
+                        if (value & spikeMaskValue) {
+                            source->mode |= PM_SOURCE_MODE_ON_SPIKE;
+                        }
+                    }
                 }
             }

branches/pap/psastro/src/psastroMaskUtils.c

@@ -181 +181 @@
 bool psastroMaskEllipticalAnnulus (psImage *mask, psImageMaskType value, float x0, float y0, psEllipseAxes eInner, psEllipseAxes eOuter) {
 
+    // skip the masking if the outer ellipse is nonsensical
+    psEllipseShape sOuter = psEllipseAxesToShape (eOuter);
+    if (isnan(sOuter.sx) || isnan(sOuter.sy) || isnan(sOuter.sxy)) return false;
+
     psEllipseShape sInner = psEllipseAxesToShape (eInner);
-    psEllipseShape sOuter = psEllipseAxesToShape (eOuter);
+    if (isnan(sInner.sx) || isnan(sInner.sy) || isnan(sInner.sxy)) {
+        // use a solid ellipse if the inner ellipse is nonsensical
+        sInner.sx = 0.1;
+        sInner.sy = 0.1;
+        sInner.sxy = 0.0;
+    }
 
     // phi is the coordinate along the elliptical path

branches/pap/psastro/src/psastroModelAnalysis.c

@@ -153 +153 @@
 
     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
-    psVectorStats (stats, posZero, NULL, NULL, 0);
+    if (!psVectorStats (stats, posZero, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+        return false;
+    }
 
     fprintf (outfile, "# pos zero %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);

branches/pap/psastro/src/psastroModelFitBoresite.c

@@ -62 +62 @@
 
     // center (Xo) = mean(Xo), RX = range / 2
-    psVectorStats (stats, Xo, NULL, NULL, 0);
+    if (!psVectorStats (stats, Xo, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+        return NULL;
+    }
     params->data.F32[PAR_X0] = stats->sampleMean;
     params->data.F32[PAR_RX] = (stats->max - stats->min) / 2.0;
 
     // center (Yo) = mean(Yo), RY = range / 2
-    psVectorStats (stats, Yo, NULL, NULL, 0);
+    if (!psVectorStats (stats, Yo, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+        return NULL;
+    }
     params->data.F32[PAR_Y0] = stats->sampleMean;
     params->data.F32[PAR_RY] = (stats->max - stats->min) / 2.0;

branches/pap/psastro/src/psastroMosaicOneChip.c

@@ -58 +58 @@
 
     // modify the order to correspond to the actual number of matched stars:
-    if ((match->n < 17) && (order >= 3)) order = 2;
-    if ((match->n < 13) && (order >= 2)) order = 1;
-    if ((match->n <  9) && (order >= 1)) order = 0;
-    if ((match->n <  3) || (order < 0) || (order > 3)) {
+    int Ndof_min = 3;
+    int order_max = 0.5*(sqrt(4*match->n - 4*Ndof_min + 1) - 3);
+    order = PS_MIN (order, order_max);
+
+    // if ((match->n < 17) && (order >= 3)) order = 2;
+    // if ((match->n < 13) && (order >= 2)) order = 1;
+    // if ((match->n <  9) && (order >= 1)) order = 0;
+
+    if (order < 0) {
         psLogMsg ("psastro", 3, "insufficient stars (%ld) or invalid order (%d)", match->n, order);
         return false;

branches/pap/psastro/src/psastroOneChipFit.c

@@ -64 +64 @@
         // modify the order to correspond to the actual number of matched stars:
         int Ndof_min = 3;
-        int order_max = 0.5*(3 + sqrt(4*match->n - 4*Ndof_min + 1));
+        int order_max = 0.5*(sqrt(4*match->n - 4*Ndof_min + 1) - 3);
         order = PS_MIN (order, order_max);
+
+        // order 0 : Ro -> nterms = 1 * 2;
+        // order 1 : Ro, Rx, Ry -> nterms = 3 * 2;
+        // order 2 : Ro, Rx, Ry, Rxx, Rxy, Ryy -> nterms = 6 * 2;
+        // order 3 : Ro, Rx, Ry, Rxx, Rxy, Ryy, Rxxx, Rxxy, Rxyy, Ryyy -> nterms = 10 * 2
+        // 2*(N+1)*(N+2)/2 = (N+1)*(N+2) = nterms;
+        // (order+1)(order+2) + ndof = nvalues
+        // order^2 + 3*order + 2 + ndof = nvalue;
+        // order^2 + 3*order + 2 + ndof - nvalue = 0;
+        // 2*order = -3 +/- sqrt (9 - 4*(2 - nvalue + ndof));
+        // 2*order = -3 +/- sqrt (9 - 8 + 4*nvalue - 4*ndof);
+        // 2*order = (sqrt (1 + 4*nvalue - 4*ndof) - 3);
 
         // if ((match->n < 11) && (order >= 3)) order = 2;

branches/pap/psastro/src/psastroZeroPoint.c

@@ -50 +50 @@
 
     // really error-out here?  or just skip?
-    if (!psastroZeroPointFromRecipe (&zeropt, &exptime, fpa, recipe)) {
+    if (!psastroZeroPointFromRecipe (&zeropt, &exptime, NULL, fpa, recipe)) {
         psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe");
         return false;
@@ -126 +126 @@
     // this analysis has too few data points to use the robust median method
     psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-    psVectorStats (stats, dMag, NULL, NULL, 0);
+    if (!psVectorStats (stats, dMag, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+        return false;
+    }
     fprintf (stderr, "zero point %f +/- %f using %d stars; transparency %f\n", stats->clippedMean, stats->clippedStdev, Npts, zeropt - stats->clippedMean);
 
@@ -145 +148 @@
   return false; }
 
-bool psastroZeroPointFromRecipe (float *zeropt, float *exptime, pmFPA *fpa, psMetadata *recipe) {
+bool psastroZeroPointFromRecipe (float *zeropt, float *exptime, float *ghostMaxMag, pmFPA *fpa, psMetadata *recipe) {
 
     bool status;
@@ -170 +173 @@
         char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER");
         if (!status) {
-            psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
+            // psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
             continue;
         }
@@ -181 +184 @@
         *zeropt = psMetadataLookupF32 (&status, refItem->data.md, "ZEROPT");
         if (!status) {
-            psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
+            psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing ZEROPT");
             continue;
+        }
+        if (ghostMaxMag) {
+            *ghostMaxMag = psMetadataLookupF32 (&status, refItem->data.md, "GHOST_MAX_MAG");
+            if (!status) {
+                psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing GHOST_MAX_MAG");
+                continue;
+            }
         }
         psFree (iter);