Index: /trunk/psModules/src/astrom/pmAstrometryModel.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryModel.c (revision 26259) +++ /trunk/psModules/src/astrom/pmAstrometryModel.c (revision 26260) @@ -725,5 +725,5 @@ double X = Xo + RX*cos(POS - To)*cos(Po) + RY*sin(POS - To)*sin(Po); double Y = Yo + RY*sin(POS - To)*cos(Po) - RX*cos(POS - To)*sin(Po); - psLogMsg ("psModules.astrom", 4, "Boresite coords on reference chip: %f, %f\n", X, Y); + psLogMsg ("psModules.astrom", 4, "Boresite coords on reference chip: %f, %f pix = %f, %f sky\n", X, Y, PM_DEG_RAD*RA, PM_DEG_RAD*DEC); // get reference chip from name Index: /trunk/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 26259) +++ /trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 26260) @@ -37,4 +37,7 @@ #include "pmKapaPlots.h" #include "pmAstrometryVisual.h" + +// XXX this is defined in pmPSFtry.h, which makes no sense +float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction); #define PM_ASTROMETRYOBJECTS_DEBUG 1 @@ -283,5 +286,6 @@ return results; } - psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", results->xStats->clippedMean, results->xStats->clippedStdev, results->xStats->clippedNvalues, x->n); + // psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", results->xStats->clippedMean, results->xStats->clippedStdev, results->xStats->clippedNvalues, x->n); + psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", results->xStats->robustMedian, results->xStats->robustStdev, results->xStats->clippedNvalues, x->n); if (!psVectorClipFitPolynomial2D (map->y, results->yStats, mask, 0xff, y, wt, X, Y)) { @@ -296,8 +300,73 @@ return results; } - psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", results->yStats->clippedMean, results->yStats->clippedStdev, results->yStats->clippedNvalues, y->n); + // psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", results->yStats->clippedMean, results->yStats->clippedStdev, results->yStats->clippedNvalues, y->n); + psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", results->yStats->robustMedian, results->yStats->robustStdev, results->yStats->clippedNvalues, y->n); } results->xStats->clipIter = stats->clipIter; results->yStats->clipIter = stats->clipIter; + + // *** calculate the 90%-ile and the systematic scatter for each direction. + + // generate the X residual vector + psVector *xFit = psPolynomial2DEvalVector (map->x, X, Y); + if (!xFit) abort(); + psVector *xRes = (psVector *) psBinaryOp (NULL, x, "-", xFit); + if (!xRes) abort(); + psFree (xFit); + + psVector *yFit = psPolynomial2DEvalVector (map->y, X, Y); + if (!yFit) abort(); + psVector *yRes = (psVector *) psBinaryOp (NULL, y, "-", yFit); + if (!yRes) abort(); + psFree (yFit); + + // extract a high-quality subset (unmasked, S/N > XXX) and position errors + // XXX for now, generate a position error based on the magnitude error + psVector *xErr = psVectorAllocEmpty (match->n, PS_TYPE_F32); + psVector *yErr = psVectorAllocEmpty (match->n, PS_TYPE_F32); + psVector *xResGood = psVectorAllocEmpty (match->n, PS_TYPE_F32); + psVector *yResGood = psVectorAllocEmpty (match->n, PS_TYPE_F32); + + for (int i = 0; i < match->n; i++) { + if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; + pmAstromMatch *pair = match->data[i]; + pmAstromObj *rawStar = raw->data[pair->raw]; + if (!isfinite(rawStar->dMag)) continue; + if (rawStar->dMag > 0.02) continue; + + // two likely failure values: NAN or 0.0 --> use dMag in this case + float xErrValue, yErrValue; + if (isfinite(rawStar->chip->xErr) && (rawStar->chip->xErr > 0.0)) { + xErrValue = rawStar->chip->xErr; + } else { + xErrValue = PS_MAX(0.005, rawStar->dMag); + } + if (isfinite(rawStar->chip->yErr) && (rawStar->chip->yErr > 0.0)) { + yErrValue = rawStar->chip->yErr; + } else { + yErrValue = PS_MAX(0.005, rawStar->dMag); + } + + psVectorAppend (xErr, xErrValue); + psVectorAppend (yErr, yErrValue); + psVectorAppend (xResGood, xRes->data.F32[i]); + psVectorAppend (yResGood, yRes->data.F32[i]); + } + + results->dXsys = psVectorSystematicError (xResGood, xErr, 0.05); + results->dYsys = psVectorSystematicError (yResGood, yErr, 0.05); + + results->dXrange = pmAstromVectorRange (xResGood, 0.1, 0.9, results->xStats->clippedStdev); + results->dYrange = pmAstromVectorRange (yResGood, 0.1, 0.9, results->yStats->clippedStdev); + + psTrace ("psModules.astrom", 3, "dXsys: %f, dXrange: %f\n", results->dXsys, results->dXrange); + psTrace ("psModules.astrom", 3, "dYsys: %f, dYrange: %f\n", results->dYsys, results->dYrange); + + psFree (xErr); + psFree (yErr); + psFree (xRes); + psFree (yRes); + psFree (xResGood); + psFree (yResGood); psFree (x); @@ -311,4 +380,138 @@ } +// set the bin closest to the corresponding value. if USE_END is +/- 1, +// out-of-range saturates on lower/upper bin REGARDLESS of actual value +#define PS_BIN_FOR_VALUE(RESULT, VECTOR, VALUE, USE_END) { \ + psVectorBinaryDisectResult result; \ + psScalar tmpScalar; \ + tmpScalar.type.type = PS_TYPE_F32; \ + tmpScalar.data.F32 = (VALUE); \ + RESULT = psVectorBinaryDisect (&result, VECTOR, &tmpScalar); \ + switch (result) { \ + case PS_BINARY_DISECT_PASS: \ + break; \ + case PS_BINARY_DISECT_OUTSIDE_RANGE: \ + psTrace("psModules.astrom", 6, "selected bin outside range"); \ + if (USE_END == -1) { RESULT = 0; } \ + if (USE_END == +1) { RESULT = VECTOR->n - 1; } \ + break; \ + case PS_BINARY_DISECT_INVALID_INPUT: \ + case PS_BINARY_DISECT_INVALID_TYPE: \ + psAbort ("programming error"); \ + break; \ + } } + +# define PS_BIN_INTERPOLATE(RESULT, VECTOR, BOUNDS, BIN, VALUE) { \ + float dX, dY, Xo, Yo, Xt; \ + if (BIN == BOUNDS->n - 1) { \ + dX = 0.5*(BOUNDS->data.F32[BIN+1] - BOUNDS->data.F32[BIN-1]); \ + dY = VECTOR->data.F32[BIN] - VECTOR->data.F32[BIN-1]; \ + Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \ + Yo = VECTOR->data.F32[BIN]; \ + } else { \ + dX = 0.5*(BOUNDS->data.F32[BIN+2] - BOUNDS->data.F32[BIN]); \ + dY = VECTOR->data.F32[BIN+1] - VECTOR->data.F32[BIN]; \ + Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \ + Yo = VECTOR->data.F32[BIN]; \ + } \ + if (dY != 0) { \ + Xt = (VALUE - Yo)*dX/dY + Xo; \ + } else { \ + Xt = Xo; \ + } \ + Xt = PS_MIN (BOUNDS->data.F32[BIN+1], PS_MAX(BOUNDS->data.F32[BIN], Xt)); \ + psTrace("psModules.astrom", 6, "(Xo, Yo, dX, dY, Xt, Yt) is (%.2f %.2f %.2f %.2f %.2f %.2f)\n", \ + Xo, Yo, dX, dY, Xt, VALUE); \ + RESULT = Xt; } + +float pmAstromVectorRange (psVector *myVector, float minFrac, float maxFrac, float stdevGuess) { + + psStats *stats = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max + psHistogram *histogram = NULL; // Histogram of the data + psHistogram *cumulative = NULL; // Cumulative histogram of the data + float min = NAN, max = NAN; // Mimimum and maximum values + + // Get the minimum and maximum values + if (!psVectorStats(stats, myVector, NULL, NULL, 0)) { + psFree(stats); + return NAN; + } + min = stats->min; + max = stats->max; + if (isnan(min) || isnan(max)) { + psFree(stats); + return NAN; + } + + psTrace("psModules.astrom", 5, "Data min/max is (%.2f, %.2f)\n", min, max); + + // If all data points have the same value, then we set the appropriate members of stats and return. + if (fabs(max - min) <= FLT_EPSILON) { + psFree (stats); + return 0.0; + } + + // Define the histogram bin size. + float binSize = 0.001; + long numBins = PS_MIN(100000, (max - min) / binSize); // Number of bins + psTrace("psModules.astrom", 5, "Numbins is %ld\n", numBins); + psTrace("psModules.astrom", 5, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max); + + // allocate the histogram containers + histogram = psHistogramAlloc(min, max, numBins); + cumulative = psHistogramAlloc(min, max, numBins); + + if (!psVectorHistogram(histogram, myVector, NULL, NULL, 0)) { + // if psVectorHistogram returns false, we have a programming error + psAbort ("Unable to generate histogram"); + } + if (psTraceGetLevel("psModules.astrom") >= 8) { + PS_VECTOR_PRINT_F32(histogram->bounds); + PS_VECTOR_PRINT_F32(histogram->nums); + } + + // Convert the specific histogram to a cumulative histogram + // The cumulative histogram data points correspond to the UPPER bound value (N < Bin[i+1]) + cumulative->nums->data.F32[0] = histogram->nums->data.F32[0]; + for (long i = 1; i < histogram->nums->n; i++) { + cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i]; + cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i]; + } + if (psTraceGetLevel("psModules.astrom") >= 8) { + PS_VECTOR_PRINT_F32(cumulative->bounds); + PS_VECTOR_PRINT_F32(cumulative->nums); + } + + // Find the bin which contains the first data point above the limit + long totalDataPoints = cumulative->nums->data.F32[numBins - 1]; + psTrace("psModules.astrom", 6, "Total data points is %ld\n", totalDataPoints); + + // find bin which is the lower bound of the limit value (value[bin] < f < value[bin+1] + long binMin; + PS_BIN_FOR_VALUE(binMin, cumulative->nums, minFrac * totalDataPoints, 0); + psTrace("psModules.astrom", 6, "The bin is %ld (%.4f to %.4f)\n", binMin, cumulative->bounds->data.F32[binMin], cumulative->bounds->data.F32[binMin+1]); + + // Linear interpolation to the limit value in bin units + float valueMin; + PS_BIN_INTERPOLATE (valueMin, cumulative->nums, cumulative->bounds, binMin, totalDataPoints * minFrac); + psTrace("psModules.astrom", 6, "limit value is %f\n", valueMin); + + // find bin which is the lower bound of the limit value (value[bin] < f < value[bin+1] + long binMax; + PS_BIN_FOR_VALUE(binMax, cumulative->nums, maxFrac * totalDataPoints, 0); + psTrace("psModules.astrom", 6, "The bin is %ld (%.4f to %.4f)\n", binMax, cumulative->bounds->data.F32[binMax], cumulative->bounds->data.F32[binMax+1]); + + // Linear interpolation to the limit value in bin units + float valueMax; + PS_BIN_INTERPOLATE (valueMax, cumulative->nums, cumulative->bounds, binMax, totalDataPoints * maxFrac); + psTrace("psModules.astrom", 6, "limit value is %f\n", valueMax); + + // Clean up + psFree(histogram); + psFree(cumulative); + psFree(stats); + + return (valueMax - valueMin); +} /****************************************************************************** @@ -634,17 +837,17 @@ } -# if 0 - char line[16]; - psFits *fits = psFitsOpen ("grid.image.fits", "w"); - psFitsWriteImage (fits, NULL, gridNP, 0, NULL); - psFitsClose (fits); - fprintf (stderr, "wrote grid image, press return to continue\n"); - fgets (line, 15, stdin); -# endif + if (psTraceGetLevel("psModules.astrom") >= 5) { + char line[16]; + psFits *fits = psFitsOpen ("grid.image.fits", "w"); + psFitsWriteImage (fits, NULL, gridNP, 0, NULL); + psFitsClose (fits); + fprintf (stderr, "wrote grid image, press return to continue\n"); + fgets (line, 15, stdin); + } // only check bins with at least 1/2 of max bin // XXX requiring at least 3 matches in bin int minNpts = PS_MAX (0.5*imStats->max, 5); - psTrace("psModule.astrom", 5, "minNpts: %d, min: %d, max: %d, median: %f, stdev: %f", minNpts, (int)(imStats->min), (int)(imStats->max), imStats->sampleMedian, imStats->sampleStdev); + psTrace("psModule.astrom", 4, "minNpts: %d, min: %d, max: %d, median: %f, stdev: %f", minNpts, (int)(imStats->min), (int)(imStats->max), imStats->sampleMedian, imStats->sampleStdev); // find the 'best' bin @@ -687,5 +890,6 @@ // XXX this function is crashing - // pmAstromVisualPlotGridMatch(raw, ref, gridNP, stats->offset.x, stats->offset.y, maxOffpix, Scale, Offset); + pmAstromVisualPlotGridMatch(raw, ref, gridNP, stats->offset.x, stats->offset.y, maxOffpix, Scale, Offset); + pmAstromVisualPlotGridMatchOverlay(raw, ref); psFree (imStats); @@ -962,2 +1166,91 @@ */ +/*****************************************************************************/ +static void pmAstromMatchInfoFree (pmAstromMatchInfo *info) +{ + if (info == NULL) return; + return; +} + + +/*****************************************************************************/ +pmAstromMatchInfo *pmAstromMatchInfoAlloc() +{ + pmAstromMatchInfo *info = psAlloc (sizeof(pmAstromMatchInfo)); + psMemSetDeallocator(info, (psFreeFunc) pmAstromMatchInfoFree); + + info->match = NULL; + info->radius = NAN; + + return (info); +} + +// generate a unique set of matches (choose closest match) +psArray *pmAstromRadiusMatchUniq (psArray *rawstars, psArray *refstars, psArray *matches) { + + // I have the matches between the refstars and the rawstars. + // For each refstar, find the single match which has the smallest radius and reject + // all others. + + // create an array of refstars->n arrays, each containing all of the matches for the + // given refstar. + + psArray *refstarMatches = psArrayAlloc (refstars->n); + + for (int i = 0; i < matches->n; i++) { + + pmAstromMatch *match = matches->data[i]; + + // accumulate this refstar match on the array for this refstar (create if needed) + psArray *refSet = refstarMatches->data[match->ref]; + if (!refSet) { + refstarMatches->data[match->ref] = psArrayAllocEmpty (8); + refSet = refstarMatches->data[match->ref]; + } + + pmAstromMatchInfo *matchInfo = pmAstromMatchInfoAlloc(); + + pmAstromObj *refStar = refstars->data[match->ref]; + pmAstromObj *rawStar = rawstars->data[match->raw]; + + matchInfo->match = match; // reference to the match of interest + matchInfo->radius = hypot (refStar->FP->x - rawStar->FP->x, refStar->FP->y - rawStar->FP->y); + + psArrayAdd (refSet, 8, matchInfo); // matchInfo->match is just a reference + psFree (matchInfo); + } + + // we now have a set of matches for each refstar and their distances; create a new set + // keeping only the closest entry for each match + + psArray *unique = psArrayAllocEmpty (PS_MAX(16, matches->n / 2)); + for (int i = 0; i < refstars->n; i++) { + + psArray *refSet = refstarMatches->data[i]; + if (!refSet) continue; + if (refSet->n == 0) continue; // not certain how this can happen... + + if (refSet->n == 1) { + pmAstromMatchInfo *matchInfo = refSet->data[0]; + psArrayAdd (unique, 32, matchInfo->match); + continue; + } + + pmAstromMatchInfo *matchInfo = refSet->data[0]; + float minRadius = matchInfo->radius; + pmAstromMatch *minMatch = matchInfo->match; + for (int j = 1; j < refSet->n; j++) { + pmAstromMatchInfo *matchInfo = refSet->data[j]; + if (minRadius < matchInfo->radius) continue; + minMatch = matchInfo->match; + minRadius = matchInfo->radius; + } + + psArrayAdd (unique, 32, minMatch); // minMatch is just a reference to a match on matches, + } + + psLogMsg ("psModules.astrom", 3, "generate unique matches to reference stars: %ld matches -> %ld matches\n", matches->n, unique->n); + psFree (refstarMatches); + + return unique; +} Index: /trunk/psModules/src/astrom/pmAstrometryObjects.h =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 26259) +++ /trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 26260) @@ -54,8 +54,19 @@ typedef struct { - int raw; ///< What is this? - int ref; ///< What is this? + int raw; ///< reference to the rawstar entry + int ref; ///< reference to the refstar entry } pmAstromMatch; + + +/* + * The pmAstromMatchInfo structure is used to generate a unique set of matches + */ +typedef struct +{ + pmAstromMatch *match; ///< reference to the match + float radius; ///< distance between the object +} +pmAstromMatchInfo; @@ -85,4 +96,8 @@ int nMatch; ///< double nSigma; ///< + double dXsys; ///< systematic error in X + double dYsys; ///< systematic error in Y + double dXrange; ///< 10% - 90% range X residuals (unmasked, high S/N) + double dYrange; ///< 10% - 90% range Y residuals (unmasked, high S/N) } pmAstromFitResults; @@ -121,4 +136,5 @@ ); +psArray *pmAstromRadiusMatchUniq (psArray *refstars, psArray *rawstars, psArray *matches); pmAstromStats *pmAstromStatsAlloc(void); @@ -343,4 +359,6 @@ ); +float pmAstromVectorRange (psVector *myVector, float minFrac, float maxFrac, float stdevGuess); + /// @} #endif // PM_ASTROMETRY_OBJECTS_H Index: /trunk/psModules/src/astrom/pmAstrometryVisual.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryVisual.c (revision 26259) +++ /trunk/psModules/src/astrom/pmAstrometryVisual.c (revision 26260) @@ -450,5 +450,5 @@ graphdata.color = KapaColorByName ("red"); - graphdata.style = 1; + graphdata.style = 0; //overplot clumpy regions excluded from analysis @@ -905,6 +905,8 @@ KapaPlotVector (kapa, gridNP->numCols, horizontalIndices, "x"); KapaPlotVector (kapa, gridNP->numCols, horizHistSlice, "y"); + float xslice[2] = {offsetX - Scale / 2., offsetX - Scale / 2.}; float yslice[2] = {-5, 100}; + graphdata.style = 0; graphdata.color = KapaColorByName("red"); KapaPrepPlot(kapa, 2, &graphdata); @@ -927,6 +929,8 @@ KapaPlotVector (kapa, gridNP->numRows, vertHistSlice, "x"); KapaPlotVector (kapa, gridNP->numRows, verticalIndices, "y"); + yslice[0] = yslice[1] = offsetY - Scale / 2.; xslice[0] = -5; xslice[1] = 100; + graphdata.style = 0; graphdata.color = KapaColorByName("red"); KapaPrepPlot(kapa, 2, &graphdata); @@ -940,4 +944,105 @@ } // end of pmAstromVisualPlotGridMatch + +bool pmAstromVisualPlotGridMatchOverlay (const psArray *raw, + const psArray *ref) +{ + //make sure we want to plot this + if (!pmVisualIsVisual() || !plotGridMatch) return true; + if (!pmVisualInitWindow(&kapa2, "psastro:plots")){ + return false; + } + + Graphdata graphdata; + psVector *xPlot = psVectorAlloc (PS_MAX(raw->n, ref->n), PS_TYPE_F32); // x data points + psVector *yPlot = psVectorAlloc (PS_MAX(raw->n, ref->n), PS_TYPE_F32); // y data points + psVector *zPlot = psVectorAlloc (PS_MAX(raw->n, ref->n), PS_TYPE_F32); // y data points + + // set up plot information + KapaClearPlots(kapa2); + KapaInitGraph(&graphdata); + + KapaSetFont(kapa2, "helvetica", 14); + KapaBox(kapa2, &graphdata); + KapaSendLabel (kapa2, "X (FP)", KAPA_LABEL_XM); + KapaSendLabel (kapa2, "Y (FP)", KAPA_LABEL_YM); + KapaSendLabel (kapa2, "pmAstromGridAngle residuals. Box: Correlation Peak.", KAPA_LABEL_XP); + + // plot the REF data. (also calculate the plot ranges, accumulate the plot vectors) + graphdata.xmin = +INT_MAX; + graphdata.xmax = -INT_MAX; + graphdata.ymin = +INT_MAX; + graphdata.ymax = -INT_MAX; + for (int i = 0; i < ref->n; i++) { + pmAstromObj *obj = ref->data[i]; + graphdata.xmin = PS_MIN(graphdata.xmin, obj->FP->x); + graphdata.xmax = PS_MAX(graphdata.xmax, obj->FP->x); + graphdata.ymin = PS_MIN(graphdata.ymin, obj->FP->y); + graphdata.ymax = PS_MAX(graphdata.ymax, obj->FP->y); + xPlot->data.F32[i] = obj->FP->x; + yPlot->data.F32[i] = obj->FP->y; + zPlot->data.F32[i] = obj->Mag; + } + xPlot->n = yPlot->n = zPlot->n = ref->n; + KapaSetLimits(kapa2, &graphdata); + + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + psVectorStats (stats, zPlot, NULL, NULL, 0); + float zero = stats->sampleMedian + 3.0; + float range = 6.0; + + for (int i = 0; i < zPlot->n; i++) { + float value = (zero - zPlot->data.F32[i]) / range; + zPlot->data.F32[i] = PS_MAX(0.0, PS_MIN(1.0, value)); + } + + // the point size will be scaled from the z vector + graphdata.style = 2; + graphdata.ptype = 7; + graphdata.size = -1; + graphdata.color = KapaColorByName ("black"); + + KapaPrepPlot (kapa2, xPlot->n, &graphdata); + KapaPlotVector (kapa2, xPlot->n, xPlot->data.F32, "x"); + KapaPlotVector (kapa2, yPlot->n, yPlot->data.F32, "y"); + KapaPlotVector (kapa2, zPlot->n, zPlot->data.F32, "z"); + + // plot the RAW data (keep previous limits) + for (int i = 0; i < raw->n; i++) { + pmAstromObj *obj = raw->data[i]; + xPlot->data.F32[i] = obj->FP->x; + yPlot->data.F32[i] = obj->FP->y; + zPlot->data.F32[i] = obj->Mag; + } + xPlot->n = yPlot->n = zPlot->n = raw->n; + + psStatsInit(stats); + psVectorStats (stats, zPlot, NULL, NULL, 0); + zero = stats->sampleMedian + 3.0; + range = 6.0; + + for (int i = 0; i < zPlot->n; i++) { + float value = (zero - zPlot->data.F32[i]) / range; + zPlot->data.F32[i] = PS_MAX(0.0, PS_MIN(1.0, value)); + } + + // the point size will be scaled from the z vector + graphdata.style = 2; + graphdata.ptype = 7; + graphdata.size = -1; + graphdata.color = KapaColorByName ("red"); + + KapaPrepPlot (kapa2, xPlot->n, &graphdata); + KapaPlotVector (kapa2, xPlot->n, xPlot->data.F32, "x"); + KapaPlotVector (kapa2, yPlot->n, yPlot->data.F32, "y"); + KapaPlotVector (kapa2, zPlot->n, zPlot->data.F32, "z"); + + pmVisualAskUser(&plotGridMatch); + psFree(xPlot); + psFree(yPlot); + psFree(zPlot); + psFree(stats); + return true; +} bool pmAstromVisualPlotTweak (psVector *xHist, // Smoothed Horizontal cut through the histogram Index: /trunk/psModules/src/astrom/pmAstrometryVisual.h =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryVisual.h (revision 26259) +++ /trunk/psModules/src/astrom/pmAstrometryVisual.h (revision 26260) @@ -45,4 +45,7 @@ ); + +bool pmAstromVisualPlotGridMatchOverlay (const psArray *raw, + const psArray *ref); /** Index: /trunk/psModules/src/objects/pmPSFtryMetric.c =================================================================== --- /trunk/psModules/src/objects/pmPSFtryMetric.c (revision 26259) +++ /trunk/psModules/src/objects/pmPSFtryMetric.c (revision 26260) @@ -76,4 +76,5 @@ pmSourceVisualPlotPSFMetric (psfTry); + pmSourceVisualPlotPSFMetricSubpix (psfTry); psFree (stats); Index: /trunk/psModules/src/objects/pmPSFtryModel.c =================================================================== --- /trunk/psModules/src/objects/pmPSFtryModel.c (revision 26259) +++ /trunk/psModules/src/objects/pmPSFtryModel.c (revision 26260) @@ -72,9 +72,5 @@ } - // XXX set the min number of needed source more carefully (depends on psfTrendMode?) - int orderMax = PS_MAX (options->psfTrendNx, options->psfTrendNy); - if ((sources->n < 15) && (orderMax >= 3)) orderMax = 2; - if ((sources->n < 11) && (orderMax >= 2)) orderMax = 1; - if ((sources->n < 8) && (orderMax >= 1)) orderMax = 0; + // hard limit on minimum number of stars if ((sources->n < 3)) { psError (PS_ERR_UNKNOWN, true, "failed to determine PSF parameters"); @@ -82,11 +78,23 @@ } - int orderMin; + // this is a bit tricky, because we have two cases (MAP vs POLY), and they have a different + // definition for 'order' (order_MAP = order_POLY + 1). in addition, we have a + // user-specified MAX order, which we should respect, regardless of the mode + + // set the max order (0 = constant) which the number of psf stars can support: + // rule of thumb: require 3 stars per 'cell' (order+1)^2 + int MaxOrderForStars = 0; + if (sources->n >= 12) MaxOrderForStars = 1; // 4 cells + if (sources->n >= 27) MaxOrderForStars = 2; // 9 cells + if (sources->n >= 48) MaxOrderForStars = 3; // 16 cells + if (sources->n > 75) MaxOrderForStars = 4; // 25 cells + + int orderMax = PS_MAX (options->psfTrendNx, options->psfTrendNy); + int orderMin = 0; if (options->psfTrendMode == PM_TREND_MAP) { - orderMin = 1; - orderMax = PS_MAX(orderMax, 1); - } else { - orderMin = 0; - } + MaxOrderForStars ++; + orderMin ++; + } + orderMax = PS_MIN (orderMax, MaxOrderForStars); // save the raw source mask (generated by pmPSFtryFitEXT) Index: /trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO.c (revision 26259) +++ /trunk/psModules/src/objects/pmSourceIO.c (revision 26260) @@ -980,4 +980,11 @@ if (!tableHeader) psAbort("cannot read table header"); + char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION"); + if (!xtension) psAbort("cannot read table type"); + if (strcmp (xtension, "BINTABLE")) { + psWarning ("no binary table in extension %s, skipping\n", dataname); + return false; + } + char *exttype = psMetadataLookupStr (NULL, tableHeader, "EXTTYPE"); if (!exttype) psAbort("cannot read table type"); Index: /trunk/psModules/src/objects/pmSourceIO_MatchedRefs.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_MatchedRefs.c (revision 26259) +++ /trunk/psModules/src/objects/pmSourceIO_MatchedRefs.c (revision 26260) @@ -120,4 +120,6 @@ psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP", PS_DATA_F32, "x coord on chip", raw->chip->x); psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP", PS_DATA_F32, "y coord on chip", raw->chip->y); + psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP_FIT",PS_DATA_F32, "x fitted coord on chip", ref->chip->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP_FIT",PS_DATA_F32, "y fitted coord on chip", ref->chip->y); psMetadataAdd (row, PS_LIST_TAIL, "X_FPA", PS_DATA_F32, "x coord on focal plane", raw->FP->x); psMetadataAdd (row, PS_LIST_TAIL, "Y_FPA", PS_DATA_F32, "y coord on focal plane", raw->FP->y); Index: /trunk/psModules/src/objects/pmSourceVisual.c =================================================================== --- /trunk/psModules/src/objects/pmSourceVisual.c (revision 26259) +++ /trunk/psModules/src/objects/pmSourceVisual.c (revision 26260) @@ -73,4 +73,167 @@ float range; + range = graphdata.xmax - graphdata.xmin; + graphdata.xmax += 0.05*range; + graphdata.xmin -= 0.05*range; + range = graphdata.ymax - graphdata.ymin; + graphdata.ymax += 0.05*range; + graphdata.ymin -= 0.05*range; + + // better choice for range? + // graphdata.xmin = -17.0; + // graphdata.xmax = -9.0; + graphdata.ymin = -0.51; + graphdata.ymax = +0.51; + + KapaSetLimits (kapa1, &graphdata); + + KapaSetFont (kapa1, "helvetica", 14); + KapaBox (kapa1, &graphdata); + KapaSendLabel (kapa1, "PSF Mag", KAPA_LABEL_XM); + KapaSendLabel (kapa1, "Ap Mag - PSF Mag", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + graphdata.etype |= 0x01; + + KapaPrepPlot (kapa1, n, &graphdata); + KapaPlotVector (kapa1, n, x->data.F32, "x"); + KapaPlotVector (kapa1, n, y->data.F32, "y"); + KapaPlotVector (kapa1, n, dy->data.F32, "dym"); + KapaPlotVector (kapa1, n, dy->data.F32, "dyp"); + + psFree (x); + psFree (y); + psFree (dy); + + pmVisualAskUser(NULL); + return true; +} + +bool pmSourceVisualPlotPSFMetricSubpix (pmPSFtry *psfTry) { + + KapaSection section; // put the positive profile in one and the residuals in another? + Graphdata graphdata; + + if (!pmVisualIsVisual()) return true; + + if (kapa1 == -1) { + kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots"); + if (kapa1 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + return false; + } + } + + KapaClearSections (kapa1); + KapaInitGraph (&graphdata); + + int n; + float range; + psVector *x = psVectorAllocEmpty (psfTry->sources->n, PS_TYPE_F32); + psVector *y = psVectorAllocEmpty (psfTry->sources->n, PS_TYPE_F32); + psVector *dy = psVectorAllocEmpty(psfTry->sources->n, PS_TYPE_F32); + + // section a: fractional-x pixel + section.dx = 1.0; + section.dy = 0.5; + section.x = 0.0; + section.y = 0.0; + section.name = NULL; + psStringAppend (§ion.name, "a1"); + KapaSetSection (kapa1, §ion); + psFree (section.name); + + graphdata.xmin = +32.0; + graphdata.xmax = -32.0; + graphdata.ymin = +32.0; + graphdata.ymax = -32.0; + + // construct the plot vectors + n = 0; + for (int i = 0; i < psfTry->sources->n; i++) { + if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue; + + pmSource *source = psfTry->sources->data[i]; + x->data.F32[n] = source->modelEXT->params->data.F32[PM_PAR_XPOS] - (int)source->modelEXT->params->data.F32[PM_PAR_XPOS]; + + y->data.F32[n] = psfTry->metric->data.F32[i]; + dy->data.F32[n] = psfTry->metricErr->data.F32[i]; + graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); + graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]); + graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]); + graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]); + n++; + } + x->n = y->n = dy->n = n; + + range = graphdata.xmax - graphdata.xmin; + graphdata.xmax += 0.05*range; + graphdata.xmin -= 0.05*range; + range = graphdata.ymax - graphdata.ymin; + graphdata.ymax += 0.05*range; + graphdata.ymin -= 0.05*range; + + // better choice for range? + // graphdata.xmin = -17.0; + // graphdata.xmax = -9.0; + graphdata.ymin = -0.51; + graphdata.ymax = +0.51; + + KapaSetLimits (kapa1, &graphdata); + + KapaSetFont (kapa1, "helvetica", 14); + KapaBox (kapa1, &graphdata); + KapaSendLabel (kapa1, "PSF Mag", KAPA_LABEL_XM); + KapaSendLabel (kapa1, "Ap Mag - PSF Mag", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + graphdata.etype |= 0x01; + + KapaPrepPlot (kapa1, n, &graphdata); + KapaPlotVector (kapa1, n, x->data.F32, "x"); + KapaPlotVector (kapa1, n, y->data.F32, "y"); + KapaPlotVector (kapa1, n, dy->data.F32, "dym"); + KapaPlotVector (kapa1, n, dy->data.F32, "dyp"); + + // *** section b: fractional-x pixel + section.dx = 1.0; + section.dy = 0.5; + section.x = 0.0; + section.y = 0.5; + section.name = NULL; + psStringAppend (§ion.name, "a2"); + KapaSetSection (kapa1, §ion); + psFree (section.name); + + graphdata.xmin = +32.0; + graphdata.xmax = -32.0; + graphdata.ymin = +32.0; + graphdata.ymax = -32.0; + + // construct the plot vectors + n = 0; + for (int i = 0; i < psfTry->sources->n; i++) { + if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue; + + pmSource *source = psfTry->sources->data[i]; + x->data.F32[n] = source->modelEXT->params->data.F32[PM_PAR_YPOS] - (int)source->modelEXT->params->data.F32[PM_PAR_YPOS]; + + y->data.F32[n] = psfTry->metric->data.F32[i]; + dy->data.F32[n] = psfTry->metricErr->data.F32[i]; + graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); + graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]); + graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]); + graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]); + n++; + } + x->n = y->n = dy->n = n; + range = graphdata.xmax - graphdata.xmin; graphdata.xmax += 0.05*range; Index: /trunk/psModules/src/objects/pmSourceVisual.h =================================================================== --- /trunk/psModules/src/objects/pmSourceVisual.h (revision 26259) +++ /trunk/psModules/src/objects/pmSourceVisual.h (revision 26260) @@ -19,4 +19,6 @@ bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask); bool pmSourceVisualPlotPSFMetric (pmPSFtry *try); +bool pmSourceVisualPlotPSFMetricSubpix (pmPSFtry *try); + bool pmSourceVisualShowModelFit (pmSource *source); bool pmSourceVisualShowModelFits (pmPSF *psf, psArray *sources, psImageMaskType maskVal);