Index: trunk/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 20010) +++ trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 20036) @@ -8,6 +8,6 @@ * @author EAM, IfA * -* @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ -* @date $Date: 2008-04-11 07:42:05 $ +* @version $Revision: 1.40 $ $Name: not supported by cvs2svn $ +* @date $Date: 2008-10-10 01:07:38 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -17,4 +17,5 @@ #include #endif + /******************************************************************************/ @@ -28,10 +29,13 @@ #include // for unlink #include +#include // cnb: do I need an IFDEF? #include "pmHDU.h" #include "pmFPA.h" #include "pmAstrometryObjects.h" +#include "pmKapaPlots.h" #define PM_ASTROMETRYOBJECTS_DEBUG 1 + /****************************************************************************** pmAstromObjSortByMag(**a, **b): sort by mag (descending) @@ -243,28 +247,28 @@ // need to use the stats lookups functions to get the width and center for (int i = 0; i < nIter; i++) { - if (!psVectorClipFitPolynomial2D (map->x, results->xStats, mask, 0xff, x, wt, X, Y)) { + if (!psVectorClipFitPolynomial2D (map->x, results->xStats, mask, 0xff, x, wt, X, Y)) { psError(PS_ERR_UNKNOWN, false, "failure in clip-fitting for x\n"); - psFree (x); - psFree (y); - psFree (X); - psFree (Y); - psFree (wt); - psFree (mask); - - return results; - } + psFree (x); + psFree (y); + psFree (X); + psFree (Y); + psFree (wt); + psFree (mask); + + 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); if (!psVectorClipFitPolynomial2D (map->y, results->yStats, mask, 0xff, y, wt, X, Y)) { psError(PS_ERR_UNKNOWN, false, "failure in clip-fitting for y\n"); - psFree (x); - psFree (y); - psFree (X); - psFree (Y); - psFree (wt); - psFree (mask); - - return results; - } + psFree (x); + psFree (y); + psFree (X); + psFree (Y); + psFree (wt); + psFree (mask); + + 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); } @@ -514,4 +518,5 @@ const psMetadata *config) { + PS_ASSERT_PTR_NON_NULL(raw, NULL); PS_ASSERT_PTR_NON_NULL(ref, NULL); @@ -559,4 +564,9 @@ psF32 **D2 = gridD2->data.F32; + // vectors to hold dX and dY + int nplot = raw->n * ref->n; + float dXplot[nplot]; + float dYplot[nplot]; + // accumulate grids for focal plane (L,M) matches for (int i = 0; i < raw->n; i++) { @@ -567,4 +577,6 @@ dY = ob1->FP->y - ob2->FP->y; + dXplot[(i * ref->n) + j] = dX; + dYplot[(i * ref->n) + j] = dY; // fprintf (f, "dX,dY: %8.2f %8.2f : %8.2f %8.2f : %8.2f %8.2f\n", dX, dY, ob1->FP->x, ob2->FP->x, ob1->FP->y, ob2->FP->y); // find bin coordinates for this delta-delta @@ -673,5 +685,59 @@ // fprintf (stderr, "sigma: nMatch: %d, nTest: %d, nTen: %d\n", stats->nMatch, stats->nTest, sort->data.U32[sort->n - 10]); - psFree (sort); + /*** plot DXs and DYs ***/ + KapaSection section = {"s1", 0.00, 0.00, 1.0, 1.0}; + //KapaSection left = {"s2", 0.05, 0.05, 0.35, 0.8}; + //KapaSection right = {"s3", 0.55, 0.05, 0.35, 0.8}; + //KapaSection yprof = {"s4", 0.45, 0.05, 0.10, 0.9}; + //KapaSection xprof = {"s5", .05, .90, .35, .1}; + + Graphdata graphdata; + int kapa = pmKapaOpen(true); + KapaClearPlots(kapa); + KapaInitGraph(&graphdata); + KapaSetSection(kapa, §ion); + graphdata.xmin = stats->offset.x - 1.5 * maxOffpix; + graphdata.xmax = stats->offset.x + 1.5 * maxOffpix; + graphdata.ymin = stats->offset.y - 1.5 * maxOffpix; + graphdata.ymax = stats->offset.y + 1.5 * maxOffpix; + + KapaSetLimits(kapa, &graphdata); + KapaSetFont(kapa, "helvetica", 14); + KapaBox(kapa, &graphdata); + KapaSendLabel (kapa, "X offset", KAPA_LABEL_XM); + KapaSendLabel (kapa, "Y offset", KAPA_LABEL_YM); + KapaSendLabel (kapa, "pmAstromGridAngle residuals. Big Box: Serach Region. Small Box: Correlation Peak.", + KAPA_LABEL_XP); + graphdata.style = 2; + graphdata.ptype = 0; + graphdata.size = 0.4; + graphdata.color = KapaColorByName ("black"); + KapaPrepPlot(kapa, nplot, &graphdata); + KapaPlotVector (kapa, nplot, dXplot, "x"); + KapaPlotVector (kapa, nplot, dYplot, "y"); + + //Overplot bounding box, peak of distribution + float xbound[5] = { -maxOffpix, maxOffpix, maxOffpix, -maxOffpix, -maxOffpix}; + float ybound[5] = { -maxOffpix, -maxOffpix, maxOffpix, maxOffpix, -maxOffpix}; + float xbin[5] = {stats->offset.x - 0.5 * Scale, stats->offset.x + 0.5 * Scale, stats->offset.x + 0.5 * Scale, + stats->offset.x - 0.5 * Scale, stats->offset.x - 0.5 * Scale}; + float ybin[5] = {stats->offset.y - 0.5 * Scale, stats->offset.y - 0.5 * Scale, stats->offset.y + 0.5 * Scale, + stats->offset.y + 0.5 * Scale, stats->offset.y - 0.5 * Scale}; + graphdata.color = KapaColorByName("red"); + graphdata.style = 0; + graphdata.size = 1.0; + KapaPrepPlot(kapa, 5, &graphdata); + KapaPlotVector (kapa, 5, xbound, "x"); + KapaPlotVector (kapa, 5, ybound, "y"); + KapaPrepPlot(kapa, 5, &graphdata); + KapaPlotVector (kapa, 5, xbin, "x"); + KapaPlotVector (kapa, 5, ybin, "y"); + + KapaPNG(kapa, "dXdY.png"); + char c; + fscanf(stdin, "%c", &c); + /*** end ***/ + + psFree (sort); psFree (listNP); psFree (gridNP); @@ -730,42 +796,41 @@ minStat->minMetric = 1e10; for (double scale = minScale; scale <= maxScale; scale += delScale) { - for (double angle = minAngle; angle <= maxAngle; angle += delAngle) { - rot = pmAstromRotateObj (raw, center, angle, scale); + for (double angle = minAngle; angle <= maxAngle; angle += delAngle) { + rot = pmAstromRotateObj (raw, center, angle, scale); # if 0 - FILE *f1 = fopen ("raw.dat", "w"); - for (int i = 0; i < rot->n; i++) { - pmAstromObj *obj = rot->data[i]; - fprintf (f1, "%8.2f %8.2f %6.2f\n", obj->FP->x, obj->FP->y, obj->Mag); - } - fclose (f1); - FILE *f2 = fopen ("ref.dat", "w"); - for (int i = 0; i < ref->n; i++) { - pmAstromObj *obj = ref->data[i]; - fprintf (f2, "%8.2f %8.2f %6.2f\n", obj->FP->x, obj->FP->y, obj->Mag); - } - fclose (f2); - fprintf (stderr, "type return"); - char c; - fscanf (stdin, "%c", &c); + FILE *f1 = fopen ("raw.dat", "w"); + for (int i = 0; i < rot->n; i++) { + pmAstromObj *obj = rot->data[i]; + fprintf (f1, "%8.2f %8.2f %6.2f\n", obj->FP->x, obj->FP->y, obj->Mag); + } + fclose (f1); + FILE *f2 = fopen ("ref.dat", "w"); + for (int i = 0; i < ref->n; i++) { + pmAstromObj *obj = ref->data[i]; + fprintf (f2, "%8.2f %8.2f %6.2f\n", obj->FP->x, obj->FP->y, obj->Mag); + } + fclose (f2); + fprintf (stderr, "type return"); + char c; + fscanf (stdin, "%c", &c); # endif - - newStat = pmAstromGridAngle (rot, ref, config); - newStat->angle = angle; - newStat->scale = scale; - newStat->center = center; - - if (isfinite(newStat->minMetric) && (newStat->minMetric > 0.0) && (newStat->minMetric < minStat->minMetric)) { - *minStat = *newStat; - psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric) *", - minStat->offset.x, minStat->offset.y, PS_DEG_RAD*minStat->angle, minStat->scale, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric)); - } else { - psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)", - newStat->offset.x, newStat->offset.y, PS_DEG_RAD*newStat->angle, newStat->scale, newStat->nMatch, newStat->nSigma, newStat->minVar, log10(newStat->minMetric)); - - } - psFree (newStat); - psFree (rot); - } + newStat = pmAstromGridAngle (rot, ref, config); + newStat->angle = angle; + newStat->scale = scale; + newStat->center = center; + + if (isfinite(newStat->minMetric) && (newStat->minMetric > 0.0) && (newStat->minMetric < minStat->minMetric)) { + *minStat = *newStat; + psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric) *", + minStat->offset.x, minStat->offset.y, PS_DEG_RAD*minStat->angle, minStat->scale, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric)); + } else { + psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)", + newStat->offset.x, newStat->offset.y, PS_DEG_RAD*newStat->angle, newStat->scale, newStat->nMatch, newStat->nSigma, newStat->minVar, log10(newStat->minMetric)); + + } + psFree (newStat); + psFree (rot); + } } psLogMsg ("psModule.astrom.grid.match", 4, "grid best - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)",