IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 9, 2008, 3:07:38 PM (18 years ago)
Author:
Paul Price
Message:

Merging from Chris Beaumont branch (cns_branch_20081009) to get luminosity function plots.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/astrom/pmAstrometryObjects.c

    r17437 r20036  
    88*  @author EAM, IfA
    99*
    10 *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2008-04-11 07:42:05 $
     10*  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2008-10-10 01:07:38 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    1717#include <config.h>
    1818#endif
     19
    1920
    2021/******************************************************************************/
    … …  
    2829#include <unistd.h>   // for unlink
    2930#include <pslib.h>
     31#include <kapa.h> // cnb: do I need an IFDEF?
    3032
    3133#include "pmHDU.h"
    3234#include "pmFPA.h"
    3335#include "pmAstrometryObjects.h"
     36#include "pmKapaPlots.h"
    3437
    3538#define PM_ASTROMETRYOBJECTS_DEBUG 1
     39
    3640/******************************************************************************
    3741pmAstromObjSortByMag(**a, **b): sort by mag (descending)
    … …  
    243247    // need to use the stats lookups functions to get the width and center
    244248    for (int i = 0; i < nIter; i++) {
    245         if (!psVectorClipFitPolynomial2D (map->x, results->xStats, mask, 0xff, x, wt, X, Y)) {
     249        if (!psVectorClipFitPolynomial2D (map->x, results->xStats, mask, 0xff, x, wt, X, Y)) {
    246250            psError(PS_ERR_UNKNOWN, false, "failure in clip-fitting for x\n");
    247             psFree (x);
    248             psFree (y);
    249             psFree (X);
    250             psFree (Y);
    251             psFree (wt);
    252             psFree (mask);
    253 
    254             return results;
    255         }
     251            psFree (x);
     252            psFree (y);
     253            psFree (X);
     254            psFree (Y);
     255            psFree (wt);
     256            psFree (mask);
     257
     258            return results;
     259        }
    256260        psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", results->xStats->clippedMean, results->xStats->clippedStdev, results->xStats->clippedNvalues, x->n);
    257261
    258262        if (!psVectorClipFitPolynomial2D (map->y, results->yStats, mask, 0xff, y, wt, X, Y)) {
    259263            psError(PS_ERR_UNKNOWN, false, "failure in clip-fitting for y\n");
    260             psFree (x);
    261             psFree (y);
    262             psFree (X);
    263             psFree (Y);
    264             psFree (wt);
    265             psFree (mask);
    266 
    267             return results;
    268         }
     264            psFree (x);
     265            psFree (y);
     266            psFree (X);
     267            psFree (Y);
     268            psFree (wt);
     269            psFree (mask);
     270
     271            return results;
     272        }
    269273        psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", results->yStats->clippedMean, results->yStats->clippedStdev, results->yStats->clippedNvalues, y->n);
    270274    }
    … …  
    514518    const psMetadata *config)
    515519{
     520
    516521    PS_ASSERT_PTR_NON_NULL(raw, NULL);
    517522    PS_ASSERT_PTR_NON_NULL(ref, NULL);
    … …  
    559564    psF32 **D2 = gridD2->data.F32;
    560565
     566    // vectors to hold dX and dY
     567    int nplot = raw->n * ref->n;
     568    float dXplot[nplot];
     569    float dYplot[nplot];
     570
    561571    // accumulate grids for focal plane (L,M) matches
    562572    for (int i = 0; i < raw->n; i++) {
    … …  
    567577            dY = ob1->FP->y - ob2->FP->y;
    568578
     579            dXplot[(i * ref->n) + j] = dX;
     580            dYplot[(i * ref->n) + j] = dY;
    569581            // 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);
    570582            // find bin coordinates for this delta-delta
    … …  
    673685    // fprintf (stderr, "sigma: nMatch: %d, nTest: %d, nTen: %d\n", stats->nMatch, stats->nTest, sort->data.U32[sort->n - 10]);
    674686
    675     psFree (sort);
     687    /*** plot DXs and DYs ***/
     688    KapaSection section = {"s1", 0.00, 0.00, 1.0, 1.0};
     689    //KapaSection left = {"s2", 0.05, 0.05, 0.35, 0.8};
     690    //KapaSection right = {"s3", 0.55, 0.05, 0.35, 0.8};
     691    //KapaSection yprof = {"s4", 0.45, 0.05, 0.10, 0.9};
     692    //KapaSection xprof = {"s5", .05, .90, .35, .1};
     693
     694    Graphdata graphdata;
     695    int kapa = pmKapaOpen(true);
     696    KapaClearPlots(kapa);
     697    KapaInitGraph(&graphdata);
     698    KapaSetSection(kapa, &section);
     699    graphdata.xmin = stats->offset.x - 1.5 * maxOffpix;
     700    graphdata.xmax = stats->offset.x + 1.5 * maxOffpix;
     701    graphdata.ymin = stats->offset.y - 1.5 * maxOffpix;
     702    graphdata.ymax = stats->offset.y + 1.5 * maxOffpix;
     703
     704    KapaSetLimits(kapa, &graphdata);
     705    KapaSetFont(kapa, "helvetica", 14);
     706    KapaBox(kapa, &graphdata);
     707    KapaSendLabel (kapa, "X offset", KAPA_LABEL_XM);
     708    KapaSendLabel (kapa, "Y offset", KAPA_LABEL_YM);
     709    KapaSendLabel (kapa, "pmAstromGridAngle residuals. Big Box: Serach Region. Small Box: Correlation Peak.",
     710                   KAPA_LABEL_XP);
     711    graphdata.style = 2;
     712    graphdata.ptype = 0;
     713    graphdata.size = 0.4;
     714    graphdata.color = KapaColorByName ("black");
     715    KapaPrepPlot(kapa, nplot, &graphdata);
     716    KapaPlotVector (kapa, nplot, dXplot, "x");
     717    KapaPlotVector (kapa, nplot, dYplot, "y");
     718
     719    //Overplot bounding box, peak of distribution
     720    float xbound[5] = { -maxOffpix, maxOffpix, maxOffpix, -maxOffpix, -maxOffpix};
     721    float ybound[5] = { -maxOffpix, -maxOffpix, maxOffpix, maxOffpix, -maxOffpix};
     722    float xbin[5] = {stats->offset.x - 0.5 * Scale, stats->offset.x + 0.5 * Scale, stats->offset.x + 0.5 * Scale,
     723                     stats->offset.x - 0.5 * Scale, stats->offset.x - 0.5 * Scale};
     724    float ybin[5] = {stats->offset.y - 0.5 * Scale, stats->offset.y - 0.5 * Scale, stats->offset.y + 0.5 * Scale,
     725                     stats->offset.y + 0.5 * Scale, stats->offset.y - 0.5 * Scale};
     726    graphdata.color = KapaColorByName("red");
     727    graphdata.style = 0;
     728    graphdata.size = 1.0;
     729    KapaPrepPlot(kapa, 5, &graphdata);
     730    KapaPlotVector (kapa, 5, xbound, "x");
     731    KapaPlotVector (kapa, 5, ybound, "y");
     732    KapaPrepPlot(kapa, 5, &graphdata);
     733    KapaPlotVector (kapa, 5, xbin, "x");
     734    KapaPlotVector (kapa, 5, ybin, "y");
     735
     736    KapaPNG(kapa, "dXdY.png");
     737    char c;
     738    fscanf(stdin, "%c", &c);
     739    /*** end ***/
     740
     741  psFree (sort);
    676742    psFree (listNP);
    677743    psFree (gridNP);
    … …  
    730796    minStat->minMetric = 1e10;
    731797    for (double scale = minScale; scale <= maxScale; scale += delScale) {
    732         for (double angle = minAngle; angle <= maxAngle; angle += delAngle) {
    733             rot = pmAstromRotateObj (raw, center, angle, scale);
     798        for (double angle = minAngle; angle <= maxAngle; angle += delAngle) {
     799            rot = pmAstromRotateObj (raw, center, angle, scale);
    734800
    735801# if 0
    736             FILE *f1 = fopen ("raw.dat", "w");
    737             for (int i = 0; i < rot->n; i++) {
    738                 pmAstromObj *obj = rot->data[i];
    739                 fprintf (f1, "%8.2f %8.2f   %6.2f\n", obj->FP->x, obj->FP->y, obj->Mag);
    740             }
    741             fclose (f1);
    742             FILE *f2 = fopen ("ref.dat", "w");
    743             for (int i = 0; i < ref->n; i++) {
    744                 pmAstromObj *obj = ref->data[i];
    745                 fprintf (f2, "%8.2f %8.2f   %6.2f\n", obj->FP->x, obj->FP->y, obj->Mag);
    746             }
    747             fclose (f2);
    748             fprintf (stderr, "type return");
    749             char c;
    750             fscanf (stdin, "%c", &c);
     802            FILE *f1 = fopen ("raw.dat", "w");
     803            for (int i = 0; i < rot->n; i++) {
     804                pmAstromObj *obj = rot->data[i];
     805                fprintf (f1, "%8.2f %8.2f   %6.2f\n", obj->FP->x, obj->FP->y, obj->Mag);
     806            }
     807            fclose (f1);
     808            FILE *f2 = fopen ("ref.dat", "w");
     809            for (int i = 0; i < ref->n; i++) {
     810                pmAstromObj *obj = ref->data[i];
     811                fprintf (f2, "%8.2f %8.2f   %6.2f\n", obj->FP->x, obj->FP->y, obj->Mag);
     812            }
     813            fclose (f2);
     814            fprintf (stderr, "type return");
     815            char c;
     816            fscanf (stdin, "%c", &c);
    751817# endif
    752 
    753             newStat = pmAstromGridAngle (rot, ref, config);
    754             newStat->angle  = angle;
    755             newStat->scale  = scale;
    756             newStat->center = center;
    757 
    758             if (isfinite(newStat->minMetric) && (newStat->minMetric > 0.0) && (newStat->minMetric < minStat->minMetric)) {
    759                 *minStat = *newStat;
    760                 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) *",
    761                           minStat->offset.x, minStat->offset.y, PS_DEG_RAD*minStat->angle, minStat->scale, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric));
    762             } else {
    763                 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)",
    764                           newStat->offset.x, newStat->offset.y, PS_DEG_RAD*newStat->angle, newStat->scale, newStat->nMatch, newStat->nSigma, newStat->minVar, log10(newStat->minMetric));
    765 
    766             }
    767             psFree (newStat);
    768             psFree (rot);
    769         }
     818            newStat = pmAstromGridAngle (rot, ref, config);
     819            newStat->angle  = angle;
     820            newStat->scale  = scale;
     821            newStat->center = center;
     822
     823            if (isfinite(newStat->minMetric) && (newStat->minMetric > 0.0) && (newStat->minMetric < minStat->minMetric)) {
     824                *minStat = *newStat;
     825                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) *",
     826                          minStat->offset.x, minStat->offset.y, PS_DEG_RAD*minStat->angle, minStat->scale, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric));
     827            } else {
     828                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)",
     829                          newStat->offset.x, newStat->offset.y, PS_DEG_RAD*newStat->angle, newStat->scale, newStat->nMatch, newStat->nSigma, newStat->minVar, log10(newStat->minMetric));
     830
     831            }
     832            psFree (newStat);
     833            psFree (rot);
     834        }
    770835    }
    771836    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)",
Note: See TracChangeset for help on using the changeset viewer.