IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 9, 2009, 11:25:34 AM (17 years ago)
Author:
Paul Price
Message:

Merging cnb_branch_20090113. No major conflicts. Compiles, but not tested.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r21350 r21422  
    1 /** Diagnostic plots called from within pmAstrometry routines.
    2  * @author      Chris Beaumont
    3  * @date        October 2008
    4  */
    5 
    6 /* Include Files   */
    71#ifdef HAVE_CONFIG_H
    82#include <config.h>
    … …  
    1610#include <pslib.h>
    1711
    18 #include "pmKapaPlots.h"
    19 
    2012#include "pmHDU.h"
    2113#include "pmFPA.h"
     14#include "pmFPAfile.h"
    2215#include "pmAstrometryObjects.h"
     16#include "pmAstrometryVisual.h"
     17#include "pmFPAExtent.h"
    2318
    2419# if (HAVE_KAPA)
    2520# include <kapa.h>
    26 
    27 # define KAPAX 700
    28 # define KAPAY 700
    29 
     21# include "pmKapaPlots.h"
     22#include "pmVisual.h"
    3023
    3124//variables to determine when things are plotted
    … …  
    3326static bool plotGridMatch        = true;
    3427static bool plotTweak            = true;
     28static bool plotRawStars         = true;
     29static bool plotRefStars         = false;
     30static bool plotLumFunc          = true;
     31static bool plotRemoveClumps     = true;
     32static bool plotOneChipFit       = true;
     33static bool plotFixChips         = true;
     34static bool plotAstromGuessCheck = true;
     35static bool plotMosaicMatches    = true;
     36static bool plotCommonScale      = true;
     37static bool plotMosaicOneChip    = true;
    3538
    3639// variables to store plotting window indices
    3740static int kapa = -1;
    38 
    39 
    40 /* Utility Routine Prototypes */
    41 bool pmAstromVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec, bool clip);
     41static int kapa2 = -1;
     42
     43//helper prototypes
     44bool residPlot (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe,
     45                char *title);
    4246
    4347
    4448/* Initialization Routines  */
    4549
    46 
    47 /** start or stop plotting */
    4850bool pmAstromSetVisual (bool mode) {
    4951    isVisual = mode;
    … …  
    5254
    5355
    54 /** open name, and resize a window if necessary */
    55 bool pmAstromVisualInitWindow (int *kapid, char *name) {
    56     if (*kapid == -1) {
    57         *kapid = KapaOpenNamedSocket("kapa", name);
    58         if (*kapid == -1) {
    59             fprintf (stderr, "failure to open kapa; visual mode disabled for pmAstrom\n");
    60             isVisual = false;
    61             return false;
    62         }
    63         KapaResize (*kapid, KAPAX, KAPAY);
    64     }
    65     return true;
    66 }
    67 
    68 
    69 /** ask the user how to proceed */
    70 bool pmAstromVisualAskUser(bool *plotflag)
    71 {
    72     char key[10];
    73     fprintf (stdout, "[c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots?");
    74     if (!fgets(key, 8, stdin)) {
    75         psWarning("Unable to read option");
    76     }
    77     if (key[0] == 's') {
    78         *plotflag = false;
    79     }
    80     if (key[0] == 'a') {
    81         pmAstromSetVisual(false);
    82     }
    83     return true;
    84 }
    85 
    86 
    87 /** destroy windows at the end of a run*/
    8856bool pmAstromVisualClose()
    8957{
    9058    if(kapa != -1)
    9159        KiiClose(kapa);
     60    if(kapa2 != -1)
     61        KiiClose(kapa2);
    9262    return true;
    9363}
    9464
    9565
    96 /* plotting routines */
    97 
    98 /**
    99  * Plot the offset between every pair of reference and raw source locations. The peak of this
    100  * distribution nominally gives the offset, scale difference, and rotation of the two catalogs.
    101  * Overplots the location of this peak as determined by pmAstromGridAngle, as well as some profiles
    102  * along horizontal and vertical cuts through this peak.
    103  */
    104 bool pmAstromVisualPlotGridMatch (const psArray *raw, ///< raw stars
    105                                   const psArray *ref, ///< reference stars
    106                                   psImage *gridNP,    ///< a 2D histogram of raw-ref star distances
    107                                   double offsetX,     ///< The X location (FP coordinates) of the peak of gridNP
    108                                   double offsetY,     ///< the Y location (FP coordinates) of the peak of gridNP
    109                                   double maxOffpix,   ///< The half-width of gridNP in FP coordinates
    110                                   double Scale,       ///< The pixel size of gridNP in histogram-bin-coordinates
    111                                   double Offset       ///< The (x,y) location (histogram-bin coordinates) of the FP point (0,0) in gridNP
     66/* Plotting Routines */
     67bool pmAstromVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe)
     68{
     69    // make sure we want to plot this
     70    if (!plotRawStars || !isVisual) return true;
     71
     72    //set up plot region
     73    if (!pmVisualInitWindow (&kapa, "psastro:plots")){
     74        isVisual = false;
     75        return false;
     76    }
     77    Graphdata graphdata;
     78    KapaSection section;
     79
     80    KapaInitGraph (&graphdata);
     81    KapaClearPlots (kapa);
     82
     83    graphdata.color = KapaColorByName ("black");
     84    graphdata.ptype = 7;
     85    graphdata.size = 0.5;
     86    graphdata.style = 2;
     87
     88    section.dx = 0.4;
     89    section.dy = 0.4;
     90
     91    //initialize and populate plotting vectors
     92    bool status = false;
     93    float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MIN");
     94    float iMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MAX");
     95
     96    psVector *xVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
     97    psVector *yVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
     98    psVector *zVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
     99
     100    section.x = 0.0;
     101    section.y = 0.5;
     102    section.name = NULL;
     103    psStringAppend (&section.name, "a0");
     104    KapaSetSection (kapa, &section);
     105    psFree (section.name);
     106
     107    //Chip coordinates
     108    int n = 0;
     109    for (int i = 0; i < rawstars->n; i++) {
     110        pmAstromObj *raw = rawstars->data[i];
     111        if (!isfinite(raw->Mag)) continue;
     112        if (raw->Mag < iMagMin) continue;
     113        if (raw->Mag > iMagMax) continue;
     114
     115        xVec->data.F32[n] = raw->chip->x;
     116        yVec->data.F32[n] = raw->chip->y;
     117        zVec->data.F32[n] = raw->Mag;
     118        n++;
     119    }
     120    xVec->n = yVec->n = zVec->n = n;
     121
     122    KapaSendLabel(kapa, "Chip", KAPA_LABEL_XP);
     123    KapaSendLabel(kapa, "X", KAPA_LABEL_XM);
     124    KapaSendLabel(kapa, "Y", KAPA_LABEL_YM);
     125    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
     126
     127    //Focal Plane Coordinates
     128    section.x = 0.5;
     129    section.y = 0.5;
     130    section.name = NULL;
     131    psStringAppend (&section.name, "a1");
     132    KapaSetSection (kapa, &section);
     133    psFree (section.name);
     134
     135    n = 0;
     136    for (int i = 0; i < rawstars->n; i++) {
     137        pmAstromObj *raw = rawstars->data[i];
     138        if (!isfinite(raw->Mag)) continue;
     139        if (raw->Mag < iMagMin) continue;
     140        if (raw->Mag > iMagMax) continue;
     141
     142        xVec->data.F32[n] = raw->FP->x;
     143        yVec->data.F32[n] = raw->FP->y;
     144        zVec->data.F32[n] = raw->Mag;
     145        n++;
     146    }
     147    xVec->n = yVec->n = zVec->n = n;
     148
     149    KapaSendLabel (kapa, "Focal Plane", KAPA_LABEL_XP);
     150    KapaSendLabel (kapa, "L", KAPA_LABEL_XM);
     151    KapaSendLabel (kapa, "M", KAPA_LABEL_YM);
     152    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
     153
     154    // Tangent Plane Coordinates
     155    section.x = 0.0;
     156    section.y = 0.0;
     157    section.name = NULL;
     158    psStringAppend (&section.name, "a2");
     159    KapaSetSection (kapa, &section);
     160    psFree (section.name);
     161
     162    n = 0;
     163    for (int i = 0; i < rawstars->n; i++) {
     164        pmAstromObj *raw = rawstars->data[i];
     165        if (!isfinite(raw->Mag)) continue;
     166        if (raw->Mag < iMagMin) continue;
     167        if (raw->Mag > iMagMax) continue;
     168
     169        xVec->data.F32[n] = raw->TP->x;
     170        yVec->data.F32[n] = raw->TP->y;
     171        zVec->data.F32[n] = raw->Mag;
     172        n++;
     173    }
     174    xVec->n = yVec->n = zVec->n = n;
     175
     176    KapaSendLabel (kapa, "Tangential Plane", KAPA_LABEL_XP);
     177    KapaSendLabel (kapa, "P", KAPA_LABEL_XM);
     178    KapaSendLabel (kapa, "Q", KAPA_LABEL_YM);
     179    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
     180
     181    //sky coordinates
     182    section.x = 0.5;
     183    section.y = 0.0;
     184    section.name = NULL;
     185    psStringAppend (&section.name, "a3");
     186    KapaSetSection (kapa, &section);
     187    psFree (section.name);
     188
     189    n = 0;
     190    for (int i = 0; i < rawstars->n; i++) {
     191        pmAstromObj *raw = rawstars->data[i];
     192        if (!isfinite(raw->Mag)) continue;
     193        if (raw->Mag < iMagMin) continue;
     194        if (raw->Mag > iMagMax) continue;
     195
     196        xVec->data.F32[n] = DEG_RAD*raw->sky->r;
     197        yVec->data.F32[n] = DEG_RAD*raw->sky->d;
     198        zVec->data.F32[n] = raw->Mag;
     199        n++;
     200    }
     201    xVec->n = yVec->n = zVec->n = n;
     202
     203    KapaSendLabel (kapa, "Sky", KAPA_LABEL_XP);
     204    KapaSendLabel(kapa, "RA", KAPA_LABEL_XM);
     205    KapaSendLabel(kapa, "Dec", KAPA_LABEL_YM);
     206    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
     207
     208    // flip x (East increase to left)
     209    SWAP (graphdata.xmin, graphdata.xmax);
     210    KapaSetLimits (kapa, &graphdata);
     211
     212    // plot label
     213    section.x = 0.0;
     214    section.y = 0.0;
     215    section.dx = .95;
     216    section.dy = .95;
     217    section.name = NULL;
     218    psStringAppend (&section.name, "a5");
     219    KapaSetSection (kapa, &section);
     220    KapaSendLabel (kapa, "Raw Star Selection - Initial Astrometry", KAPA_LABEL_XP);
     221    psFree (section.name);
     222
     223    // pause and wait for user input:
     224    pmVisualAskUser(&plotRawStars, &isVisual);
     225
     226    psFree (xVec);
     227    psFree (yVec);
     228    psFree (zVec);
     229    return true;
     230}
     231
     232
     233bool pmAstromVisualPlotRefStars (psArray *refstars, psMetadata *recipe)
     234{
     235    //make sure we want to plot this
     236    if (!isVisual || !plotRefStars) return true;
     237
     238    //set up plotting variables
     239    if (!pmVisualInitWindow (&kapa, "psastro:plots")) {
     240        isVisual = false;
     241        return false;
     242    }
     243
     244    Graphdata graphdata;
     245    KapaInitGraph (&graphdata);
     246    KapaClearSections (kapa);
     247
     248    graphdata.color = KapaColorByName ("black");
     249    graphdata.ptype = 7;
     250    graphdata.size = 0.5;
     251    graphdata.style = 2;
     252
     253    //initialize and populate plot vectors
     254    bool status = false;
     255    float rMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MIN");
     256    float rMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MAX");
     257
     258    psVector *xVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
     259    psVector *yVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
     260    psVector *zVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
     261
     262    int n = 0;
     263    for (int i = 0; i < refstars->n; i++) {
     264        pmAstromObj *ref = refstars->data[i];
     265        if (!isfinite(ref->Mag)) continue;
     266        if (ref->Mag > rMagMax) continue;
     267        if (ref->Mag < rMagMin) continue;
     268
     269        xVec->data.F32[n] = DEG_RAD*ref->sky->r;
     270        yVec->data.F32[n] = DEG_RAD*ref->sky->d;
     271        zVec->data.F32[n] = ref->Mag;
     272        n++;
     273    }
     274    xVec->n = yVec->n = zVec->n = n;
     275
     276    KapaSendLabel (kapa, "RA", KAPA_LABEL_XM);
     277    KapaSendLabel (kapa, "Dec", KAPA_LABEL_YM);
     278    KapaSendLabel (kapa, "Reference Stars", KAPA_LABEL_XP);
     279    pmVisualTriplePlot(kapa, &graphdata, xVec, yVec, zVec, false);
     280
     281    // flip x (East increase to left)
     282    SWAP (graphdata.xmin, graphdata.xmax);
     283    KapaSetLimits (kapa, &graphdata);
     284
     285    // pause and wait for user input:
     286    pmVisualAskUser(&plotRefStars, &isVisual);
     287
     288    psFree (xVec);
     289    psFree (yVec);
     290    psFree (zVec);
     291    return true;
     292}
     293
     294
     295bool pmAstromVisualPlotLuminosityFunction (psVector *lnMag,   // Log(n) for each magnitude bin
     296                                          psVector *Mag,     // magnitude bins
     297                                          pmLumFunc *lumFunc,// Fit to the reference star luminosity function
     298                                          pmLumFunc *rawFunc // Fit to the raw star luminoisty function
     299                                          )
     300{
     301
     302    // make sure we want to plot this
     303    if ( !isVisual || !plotLumFunc ) return true;
     304
     305    //set up plotting variables
     306    if ( !pmVisualInitWindow (&kapa, "psastro:plots")){
     307        isVisual = false;
     308        return false;
     309    }
     310
     311    Graphdata graphdata;
     312    KapaSection section1 = {"s1", .1, .1, .85, .35};
     313    KapaSection section2 =  {"s2", .1, .5, .85, .35};
     314    KapaSection section;
     315    if(rawFunc == NULL)
     316        section = section1;
     317    else
     318        section = section2;
     319    if (rawFunc == NULL)
     320        KapaClearPlots(kapa);
     321    KapaInitGraph(&graphdata);
     322
     323    //Determine Plot Limits
     324    pmVisualScaleGraphdata(&graphdata, Mag, lnMag, false);
     325
     326    //Make a line for the fit
     327    float x[2] = {graphdata.xmin, graphdata.xmax};
     328    float y[2] = {lumFunc->offset + x[0] * lumFunc->slope,
     329                 lumFunc->offset + x[1] * lumFunc->slope};
     330
     331    //Plot Data
     332    KapaSetSection(kapa, &section);
     333    KapaSetLimits(kapa, &graphdata);
     334
     335    KapaSetFont (kapa, "helvetica", 14);
     336    KapaBox (kapa, &graphdata);
     337    KapaSendLabel (kapa, "Magnitude", KAPA_LABEL_XM);
     338    KapaSendLabel (kapa, "Log(N)", KAPA_LABEL_YM);
     339    if (rawFunc == NULL)
     340        KapaSendLabel (kapa, "Raw Star Luminosity Function", KAPA_LABEL_XP);
     341    else
     342        KapaSendLabel (kapa,
     343                       "Reference Star Luminosity Function, Shifted Raw Fit, and Cutoff",
     344                       KAPA_LABEL_XP);
     345    graphdata.color=KapaColorByName("black");
     346    graphdata.style = 1;
     347    KapaPrepPlot (kapa, lnMag->n, &graphdata);
     348    KapaPlotVector(kapa, lnMag->n,   Mag->data.F32, "x");
     349    KapaPlotVector(kapa, lnMag->n, lnMag->data.F32, "y");
     350
     351    //Overplot fit
     352    graphdata.style=0;
     353    KapaPrepPlot(kapa,2,&graphdata);
     354    KapaPlotVector(kapa, 2, x, "x");
     355    KapaPlotVector(kapa, 2, y, "y");
     356
     357    //If rawFunc was supplied, overplot the raw star fit + cutoff
     358    if( rawFunc != NULL) {
     359        double mRef = 0.5*(lumFunc->mMin + lumFunc->mMax);
     360        double logRho = mRef * lumFunc->slope + lumFunc->offset;
     361        double mRaw = (logRho - rawFunc->offset) / rawFunc->slope;
     362        double deltaM = mRef - mRaw;
     363        double mRefMax = rawFunc->mMax + deltaM;
     364
     365        float xraw[2] = {rawFunc->mMin + deltaM, rawFunc->mMax + deltaM};
     366        float yraw[2] = {rawFunc->offset + (rawFunc->slope) * rawFunc->mMin,
     367                        rawFunc->offset + (rawFunc->slope) * rawFunc->mMax};
     368        float x[2] = {mRefMax, mRefMax};
     369        float y[2] = {graphdata.ymin, graphdata.ymax};
     370        graphdata.color= KapaColorByName("red");
     371        KapaPrepPlot(kapa, 2, &graphdata);
     372        KapaPlotVector(kapa, 2, x, "x");
     373        KapaPlotVector(kapa, 2, y, "y");
     374        KapaPrepPlot (kapa, 2, &graphdata);
     375        KapaPlotVector (kapa, 2, xraw, "x");
     376        KapaPlotVector (kapa, 2, yraw, "y");
     377
     378        // pause and wait for user input:
     379        pmVisualAskUser (&plotLumFunc, &isVisual);
     380    }
     381    return true;
     382} // end of pmAstromVisualPlotLuminosityFunction
     383
     384
     385bool pmAstromVisualPlotRemoveClumps (psArray *input, // Array containing the field stars
     386                                    psImage *count, // A 2D histogram of the field star distribution
     387                                    int scale,      // The pixel size of the histogram
     388                                    float limit     // The minimum numuber of stars in a bin flagged as a clump
     389                                    )
     390{
     391
     392    //make sure we want to plot this
     393    if (!isVisual || !plotRemoveClumps) return true;
     394
     395    //set up plot variables
     396    if ( !pmVisualInitWindow (&kapa, "psastro:plots")) {
     397        isVisual = false;
     398        return false;
     399    }
     400
     401    KapaSection section;
     402    Graphdata graphdata;
     403    section.x = 0;
     404    section.dx = .9;
     405    section.y = 0;
     406    section.dy = .9;
     407    section.name = NULL;
     408    psStringAppend( &section.name, "a0");
     409    KapaInitGraph(&graphdata);
     410    KapaSetSection(kapa, &section);
     411    psFree(section.name);
     412
     413    graphdata.ptype = 7;
     414    graphdata.size = 0.5;
     415    graphdata.style = 2;
     416    graphdata.color = KapaColorByName ("black");
     417    KapaClearPlots(kapa);
     418
     419    //set up plot vectors
     420    float Xmin = +FLT_MAX;
     421    float Xmax = -FLT_MAX;
     422    float Ymin = +FLT_MAX;
     423    float Ymax = -FLT_MAX;
     424    psVector *xVec = psVectorAlloc (input->n, PS_TYPE_F32);
     425    psVector *yVec = psVectorAlloc (input->n, PS_TYPE_F32);
     426
     427    //determine boundaries for histogram bin calculation
     428    int n = 0;
     429    for (int i=0; i< input->n; i++) {
     430        pmAstromObj *obj = (pmAstromObj *)input->data[i];
     431        if (!isfinite(obj->FP->x)) continue;
     432        if (!isfinite(obj->FP->y)) continue;
     433        xVec->data.F32[n] = obj->FP->x;
     434        yVec->data.F32[n] = obj->FP->y;
     435        Xmin = PS_MIN (Xmin, xVec->data.F32[n]);
     436        Xmax = PS_MAX (Xmax, xVec->data.F32[n]);
     437        Ymin = PS_MIN (Ymin, yVec->data.F32[n]);
     438        Ymax = PS_MAX (Ymax, yVec->data.F32[n]);
     439        n++;
     440    }
     441    xVec->n = yVec->n = n;
     442
     443    //plot stars
     444    graphdata.xmax = Xmax;
     445    graphdata.xmin = Xmin;
     446    graphdata.ymax = Ymax;
     447    graphdata.ymin = Ymin;
     448    KapaSetLimits (kapa, &graphdata);
     449    KapaSetFont (kapa, "helvetica", 14);
     450
     451    KapaBox (kapa, &graphdata);
     452
     453    KapaSendLabel (kapa, "L (pixels)", KAPA_LABEL_XM);
     454    KapaSendLabel (kapa, "M (pixels)", KAPA_LABEL_YM);
     455    KapaSendLabel (kapa, "Regions Flagged as Clumps (Red Boxes)",
     456                   KAPA_LABEL_XP);
     457
     458    KapaPrepPlot (kapa, xVec->n, &graphdata);
     459    KapaPlotVector (kapa, xVec->n, xVec->data.F32, "x");
     460    KapaPlotVector (kapa, yVec->n, yVec->data.F32, "y");
     461
     462    graphdata.color = KapaColorByName ("red");
     463    graphdata.style = 1;
     464
     465    //overplot clumpy regions excluded from analysis
     466    for(int i = 0; i < count->numCols; i++) {
     467        for (int j = 0; j < count->numRows; j++) {
     468            if(count->data.U32[j][i] <= limit) continue; //not a clump
     469            float Xbot = (i - 5) * scale + Xmin;
     470            float Ybot = (j - 5) * scale + Ymin;
     471            if(Xbot < graphdata.xmin || Xbot > graphdata.xmax ||
     472               Ybot < graphdata.ymin || Ybot > graphdata.ymax) continue;
     473            float x[5] = {Xbot, Xbot + scale, Xbot + scale, Xbot, Xbot};
     474            float y[5] = {Ybot, Ybot, Ybot + scale, Ybot + scale, Ybot};
     475            KapaPrepPlot (kapa, 5, &graphdata);
     476            KapaPlotVector (kapa, 5, x, "x");
     477            KapaPlotVector (kapa, 5, y, "y");
     478        }
     479    }
     480
     481    //ask for user input and finish
     482    pmVisualAskUser (&plotRemoveClumps, &isVisual);
     483    psFree (xVec);
     484    psFree (yVec);
     485
     486    return true;
     487}
     488
     489
     490bool pmAstromVisualPlotOneChipFit (psArray *rawstars, // stars detected in the image
     491                                  psArray *refstars, // reference stars over the same region
     492                                  psArray *match,    // contains which rawstars match to which refstars
     493                                  psMetadata *recipe // data reduction recipe
    112494                                  )
     495{
     496
     497    //make sure we want to plot this
     498    if (!isVisual || !plotOneChipFit)
     499        return true;
     500
     501    if(!pmVisualInitWindow(&kapa, "psastro:plots") || !pmVisualInitWindow(&kapa2, "psastro:plots")) {
     502        isVisual = false;
     503        return false;
     504    }
     505
     506    //plot the residuals
     507    if (!residPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals (Chip Coordinates)")) {
     508        plotOneChipFit = false;
     509        return false;
     510    }
     511
     512    //ask for user input and finish
     513    pmVisualAskUser(&plotOneChipFit, &isVisual);
     514    return true;
     515}
     516
     517
     518bool pmAstromVisualPlotFixChips (pmFPAfile *input, // focal plane array file
     519                                psVector *xOld, // old X location of chip cornerss
     520                                psVector *yOld // old Y location of chip corners
     521                                )
     522{
     523    //make sure we want to plot this
     524    if(!isVisual || !plotFixChips) return true;
     525
     526    if(!pmVisualInitWindow(&kapa, "psastro:plots")) {
     527        isVisual = false;
     528        return false;
     529    }
     530
     531    KapaSection section = {"s1", .05, .05, .9, .9};
     532    Graphdata graphdata;
     533    KapaInitGraph( &graphdata);
     534    KapaClearPlots (kapa);
     535    graphdata.ptype = 2;
     536    graphdata.style = 2;
     537
     538    psVector *xNew = psVectorAllocEmpty (xOld->n, PS_TYPE_F32);
     539    psVector *yNew = psVectorAllocEmpty (yOld->n, PS_TYPE_F32);
     540
     541    // copy of the code in psastroFixChips that generated xOld, yOld, but for xNew, yNew
     542    pmFPAview *view = pmFPAviewAlloc (0);
     543
     544    pmChip *obsChip = NULL;
     545    while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
     546        if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
     547
     548        psRegion *region = pmChipPixels(obsChip);
     549        psPlane ptCP, ptFP;
     550
     551        ptCP.x = region->x0; ptCP.y = region->y0;
     552        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     553        psVectorAppend (xNew, ptFP.x);
     554        psVectorAppend (yNew, ptFP.y);
     555
     556        ptCP.x = region->x0; ptCP.y = region->y1;
     557        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     558        psVectorAppend (xNew, ptFP.x);
     559        psVectorAppend (yNew, ptFP.y);
     560
     561        ptCP.x = region->x1; ptCP.y = region->y1;
     562        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     563        psVectorAppend (xNew, ptFP.x);
     564        psVectorAppend (yNew, ptFP.y);
     565
     566        ptCP.x = region->x1; ptCP.y = region->y0;
     567        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     568        psVectorAppend (xNew, ptFP.x);
     569        psVectorAppend (yNew, ptFP.y);
     570
     571        psFree (region);
     572    }
     573
     574    //set up graph
     575    pmVisualScaleGraphdata(&graphdata, xOld, yOld, true);
     576    pmVisualInitGraph(kapa, &section, &graphdata);
     577    KapaSendLabel (kapa, "L (FP)", KAPA_LABEL_XM);
     578    KapaSendLabel (kapa, "M (FP)", KAPA_LABEL_YM);
     579    KapaSendLabel (kapa, "Chip corners before (black) and after (red) FixChips", KAPA_LABEL_XP);
     580    KapaPrepPlot (kapa, xOld->n, &graphdata);
     581    KapaPlotVector (kapa, xOld->n, xOld->data.F32, "x");
     582    KapaPlotVector (kapa, xOld->n, yOld->data.F32, "y");
     583
     584    graphdata.ptype = 1;
     585    graphdata.color = KapaColorByName("red");
     586    KapaPrepPlot (kapa, xNew->n, &graphdata);
     587    KapaPlotVector (kapa, xNew->n, xNew->data.F32, "x");
     588    KapaPlotVector (kapa, xNew->n, yNew->data.F32, "y");
     589
     590    pmVisualAskUser(&plotFixChips, &isVisual);
     591    psFree(xNew);
     592    psFree(yNew);
     593    psFree(view);
     594    return true;
     595}
     596
     597
     598bool pmAstromVisualPlotAstromGuessCheck (psVector *cornerPo, // P coordinates of chip corners before fitting
     599                                        psVector *cornerQo, // Q coordinates of chip corners before fitting
     600                                        psVector *cornerPn, // P coordinates of chip corners after fitting
     601                                        psVector *cornerQn, // Q coordinates of chip corners after fitting
     602                                        psVector *cornerPd, // P coordinate residuals of fit from old to new coordinates
     603                                        psVector *cornerQd  // Q coordinate residuals of fit from old to new coordinates
     604                                        )
     605{
     606
     607    //make sure we want to plot this
     608    if (!isVisual || !plotAstromGuessCheck) return true;
     609
     610    //set up graph window
     611    if ( !pmVisualInitWindow (&kapa, "psastro:plots")) {
     612        isVisual = false;
     613        return false;
     614    }
     615
     616    Graphdata graphdata;
     617    KapaSection section;
     618    KapaInitGraph(&graphdata);
     619    KapaClearPlots (kapa);
     620
     621    graphdata.color = KapaColorByName ("black");
     622    graphdata.ptype = 2;
     623    graphdata.size = 0.5;
     624    graphdata.style = 2;
     625
     626    section.dx = 0.4;
     627    section.dy = 0.4;
     628
     629    //Old Corners
     630    section.x = 0.30;
     631    section.y = 0.50;
     632    section.name = NULL;
     633    psStringAppend (&section.name, "a0");
     634    KapaSetSection (kapa, &section);
     635    psFree(section.name);
     636
     637    pmVisualScaleGraphdata (&graphdata, cornerPo, cornerPo, true);
     638    KapaSetLimits (kapa, &graphdata);
     639    KapaBox (kapa, &graphdata);
     640    KapaSendLabel (kapa, "P (Pixels)", KAPA_LABEL_XM);
     641    KapaSendLabel (kapa, "Q (Pixels)", KAPA_LABEL_YM);
     642    KapaSendLabel (kapa,
     643                   "Fiducial Points in the Tangent Plane. Black: Initial Astrometry. Red: Final Astrometry",
     644                   KAPA_LABEL_XP);
     645    KapaPrepPlot (kapa, cornerPo->n, &graphdata);
     646    KapaPlotVector (kapa, cornerPo->n, cornerPo->data.F32, "x");
     647    KapaPlotVector (kapa, cornerQo->n, cornerQo->data.F32, "y");
     648
     649    // New Corners
     650    graphdata.color = KapaColorByName("red");
     651    graphdata.ptype = 7;
     652    graphdata.size = 1.5;
     653    KapaPrepPlot (kapa, cornerPn->n, &graphdata);
     654    KapaPlotVector (kapa, cornerPn->n, cornerPn->data.F32, "x");
     655    KapaPlotVector (kapa, cornerQn->n, cornerQn->data.F32, "y");
     656
     657    // Residuals
     658    psVector *xResid = psVectorAlloc(cornerPn->n, PS_DATA_F32);
     659    psVector *yResid = psVectorAlloc(cornerQn->n, PS_DATA_F32);
     660    for(int i=0; i < cornerPn->n; i++) {
     661        xResid->data.F32[i] = (cornerPd->data.F32[i]);
     662        yResid->data.F32[i] = (cornerQd->data.F32[i]);
     663    }
     664
     665    graphdata.color = KapaColorByName("black");
     666    graphdata.size=0.5;
     667    section.x = 0.3;
     668    section.y = 0.0;
     669    section.name = NULL;
     670    psStringAppend (&section.name, "a1");
     671    KapaSetSection (kapa, &section);
     672    psFree(section.name);
     673
     674    pmVisualScaleGraphdata (&graphdata, xResid, yResid, true);
     675    KapaSetLimits (kapa, &graphdata);
     676    KapaBox (kapa, &graphdata);
     677    KapaSendLabel (kapa, "dP", KAPA_LABEL_XM);
     678    KapaSendLabel (kapa, "dQ", KAPA_LABEL_YM);
     679    KapaSendLabel (kapa,
     680                   "Residual of the Fit from the Initial Astrometry to the Final Astrometry",
     681                   KAPA_LABEL_XP);
     682    KapaPrepPlot (kapa, cornerPd->n, &graphdata);
     683    KapaPlotVector (kapa, cornerPd->n, xResid->data.F32, "x");
     684    KapaPlotVector (kapa, cornerQd->n, yResid->data.F32, "y");
     685
     686    psFree(xResid);
     687    psFree(yResid);
     688    pmVisualAskUser (&plotAstromGuessCheck, &isVisual);
     689    return true;
     690}
     691
     692
     693bool pmAstromVisualPlotCommonScale (pmFPA *fpa,         // the fpa
     694                                   psVector *oldScale  // the old pixel scale of each chip in the fpa
     695                                   )
     696{
     697    //make sure we want to plot this
     698    if (!isVisual || !plotCommonScale) return true;
     699
     700    if (!pmVisualInitWindow(&kapa, "psastro:plots")){
     701        isVisual = false;
     702        return false;
     703    }
     704
     705    KapaSection section = {"s1", .05, .05, .9, .9};
     706    Graphdata graphdata;
     707    psPlane ptCH, ptFP;
     708    ptCH.x = 0;
     709    ptCH.y = 0;
     710    psVector *xVec = psVectorAlloc (oldScale->n, PS_TYPE_F32);
     711    psVector *yVec = psVectorAlloc (oldScale->n, PS_TYPE_F32);
     712
     713    int nobj = 0;
     714
     715    //project each chip corner to the Focal Plane
     716    for(int i = 0; i < fpa->chips->n; i++) {
     717        pmChip *chip = fpa->chips->data[i];
     718        if (!chip->process || !chip->file_exists) { continue; }
     719        if (!chip->toFPA) { continue; }
     720
     721        psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
     722        xVec->data.F32[nobj] = ptFP.x;
     723        yVec->data.F32[nobj] = ptFP.y;
     724        nobj++;
     725        if (nobj == oldScale->n) break;
     726    }
     727
     728    //set up plot window
     729    KapaInitGraph (&graphdata);
     730    KapaClearPlots (kapa);
     731    KapaSetSection (kapa, &section);
     732    KapaSetFont (kapa, "helvetica", 14);
     733    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, oldScale, false);
     734    KapaSendLabel (kapa, "L (FP)", KAPA_LABEL_XM);
     735    KapaSendLabel (kapa, "M (FP)", KAPA_LABEL_YM);
     736    KapaSendLabel (kapa, "Old Pixel Scale of FPA Chips (Not to Scale)", KAPA_LABEL_XP);
     737
     738    pmVisualAskUser (&plotCommonScale, &isVisual);
     739
     740    psFree(xVec);
     741    psFree(yVec);
     742
     743    return true;
     744}
     745
     746
     747bool pmAstromVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars,
     748                                     psArray *match, psMetadata *recipe)
     749{
     750
     751    //make sure we want to plot this
     752    if (!isVisual || !plotMosaicOneChip) return true;
     753
     754    if(!pmVisualInitWindow(&kapa, "psastro:plots") || !pmVisualInitWindow(&kapa2, "psastro:plots")) {
     755        isVisual = false;
     756        return false;
     757    }
     758
     759    //plot the residuals
     760    if (!residPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals - Mosaic Mode")) {
     761        isVisual = false;
     762        return false;
     763    }
     764
     765    //ask for user input and finish
     766    pmVisualAskUser(&plotMosaicOneChip, &isVisual);
     767
     768    return true;
     769}
     770
     771
     772bool pmAstromVisualPlotMosaicMatches( psArray *rawstars, psArray *refstars,
     773                                    psArray *match, int iteration,
     774                                    psMetadata *recipe)
     775{
     776    //make sure we want to plot this
     777    if (!isVisual || !plotMosaicMatches) return true;
     778
     779    char title[60];
     780    sprintf(title, "Matches found during psastroMosaicSetMatch iteration %d", iteration);
     781
     782
     783    if(!pmVisualInitWindow(&kapa, "psastro:plots") || !pmVisualInitWindow(&kapa2, "psastro:plots")) {
     784        isVisual = false;
     785        return false;
     786    }
     787
     788    if (!residPlot(rawstars, refstars, match, recipe, title)){
     789        isVisual = false;
     790        return false;
     791    }
     792
     793    //ask for user input
     794    pmVisualAskUser (&plotMosaicMatches, &isVisual);
     795    return true;
     796}
     797
     798
     799bool pmAstromVisualPlotGridMatch (const psArray *raw,
     800                                  const psArray *ref,
     801                                  psImage *gridNP,
     802                                  double offsetX,
     803                                  double offsetY,
     804                                  double maxOffpix,
     805                                  double Scale,
     806                                  double Offset)
    113807{
    114808    //make sure we want to plot this
    115809    if (!isVisual || !plotGridMatch) return true;
    116     if (!pmAstromVisualInitWindow(&kapa, "pmAstrom:plots"))
    117         return false;
     810    if (!pmVisualInitWindow(&kapa, "psastro:plots")){
     811        isVisual = false;
     812        return false;
     813    }
    118814
    119815    KapaSection section = {"s1", 0.05, 0.05, .75, .75};
    … …  
    123819    Graphdata graphdata;
    124820    int nplot = raw->n * ref->n;                      // number of points to plot
    125     float dXplot[nplot];                              // x data points
    126     float dYplot[nplot];                              // y data points
     821    psVector *dXplot = psVectorAlloc (nplot, PS_TYPE_F32); // x data points
     822    psVector *dYplot = psVectorAlloc (nplot, PS_TYPE_F32); // y data points
     823
    127824    pmAstromObj *ob1; pmAstromObj *ob2;               // shortcuts to the data in raw and ref
    128825    psU32 **NP = gridNP->data.U32;                    // shortcut to the gridNP data
    … …  
    168865            dX = ob1->FP->x - ob2->FP->x;
    169866            dY = ob1->FP->y - ob2->FP->y;
    170             dXplot[(i * ref->n) + j] = dX;
    171             dYplot[(i * ref->n) + j] = dY;
     867            dXplot->data.F32[(i * ref->n) + j] = dX;
     868            dYplot->data.F32[(i * ref->n) + j] = dY;
    172869        }
    173870    }
    … …  
    191888    //Plot the offsets
    192889    KapaPrepPlot(kapa, nplot, &graphdata);
    193     KapaPlotVector (kapa, nplot, dXplot, "x");
    194     KapaPlotVector (kapa, nplot, dYplot, "y");
     890    KapaPlotVector (kapa, nplot, dXplot->data.F32, "x");
     891    KapaPlotVector (kapa, nplot, dYplot->data.F32, "y");
    195892
    196893    //Overplot bounding box, peak of distribution
    … …  
    255952    KapaPlotVector (kapa, 2, yslice, "y");
    256953
    257     pmAstromVisualAskUser(&plotGridMatch);
     954    pmVisualAskUser(&plotGridMatch, &isVisual);
     955    psFree(dXplot);
     956    psFree(dYplot);
    258957    return true;
    259958} // end of pmAstromVisualPlotGridMatch
    260959
    261960
    262 /** Plot the refinements made within pmAstromGridTweak.
    263  * After pmAstromGridMatch finds the best rotaion/scale/offset between raw and reference stars
    264  * within a coarse grid of rotations/scales, pmAstromGridTweak computes a higher precision
    265  * estimate of the offset. It computes the 2 point correlation function between raw and ref
    266  * stars along horizontal and vertical cuts through the first-guess offset. It finds the peak
    267  * of these two profiles and adjusts the offset accordingly. This procedure plots the profiles.
    268  */
    269 bool pmAstromVisualPlotTweak (psVector *xHist, ///< Smoothed Horizontal cut through the histogram
    270                               psVector *yHist, ///< Smoothed Vertical cut throug the histogram
    271                               int xBin,        ///< X Bin index of the histogram peak
    272                               int yBin         ///< Y bin index of the histogram peak
     961bool pmAstromVisualPlotTweak (psVector *xHist, // Smoothed Horizontal cut through the histogram
     962                              psVector *yHist, // Smoothed Vertical cut throug the histogram
     963                              int xBin,        // X Bin index of the histogram peak
     964                              int yBin         // Y bin index of the histogram peak
    273965    )
    274966{
    275967    //make sure we want to plot this
    276968    if (!isVisual || !plotTweak) return true;
    277     if (!pmAstromVisualInitWindow(&kapa, "pmAstrom:plots")) {
     969    if (!pmVisualInitWindow(&kapa, "psastro:plots")) {
     970        isVisual = false;
    278971        return false;
    279972    }
    … …  
    299992
    300993    // plot the X histogram
    301     pmAstromVisualScaleGraphdata(&graphdata, xIndices, xHist, false);
     994    pmVisualScaleGraphdata(&graphdata, xIndices, xHist, false);
    302995    KapaSetSection(kapa, &section1);
    303996    KapaSetLimits (kapa, &graphdata);
    … …  
    3261019
    3271020    //plot the Y histogram
    328     pmAstromVisualScaleGraphdata(&graphdata, yIndices, yHist, false);
     1021    pmVisualScaleGraphdata(&graphdata, yIndices, yHist, false);
    3291022    KapaSetSection(kapa, &section2);
    3301023    KapaSetLimits (kapa, &graphdata);
    … …  
    3581051                   KAPA_LABEL_XP);
    3591052
    360     pmAstromVisualAskUser(&plotTweak);
     1053    pmVisualAskUser(&plotTweak, &isVisual);
    3611054
    3621055    psFree(xIndices);
    … …  
    3661059
    3671060
    368 /** pmAstromScaleGraphdata
    369  * Scale the graphdata structure based on x and y coordinates. Use sigma clipping to
    370  * prevent outliers from making te plot region too big.
    371  */
    372 bool pmAstromVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec, bool clip) {
    373 
    374     graphdata->xmin = +FLT_MAX;
    375     graphdata->xmax = -FLT_MAX;
    376     graphdata->ymin = +FLT_MAX;
    377     graphdata->ymax = -FLT_MAX;
    378 
    379     //determine standard deviation of xVec and yVec
    380     psStats *statsX = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    381     psStats *statsY = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    382     psVectorStats (statsX, xVec, NULL, NULL, 0);
    383     psVectorStats (statsY, yVec, NULL, NULL, 0);
    384 
    385     float xhi  = statsX->sampleMedian + 3 *statsX->sampleStdev;
    386     float xlo = statsX->sampleMedian - 3 *statsX->sampleStdev;
    387     float yhi = statsY->sampleMedian + 3 *statsY->sampleStdev;
    388     float ylo = statsY->sampleMedian - 3 *statsY->sampleStdev;
    389 
    390     // don't sigma clip
    391     if (!clip) {
    392         xhi = +FLT_MAX;
    393         xlo = -FLT_MAX;
    394         yhi = +FLT_MAX;
    395         ylo = -FLT_MAX;
    396     }
    397 
    398     // abort if there is no good data
    399     if (!isfinite(xhi) || !isfinite(xlo) || !isfinite(yhi) || !isfinite(ylo)) {
    400         graphdata->xmin = -1;
    401         graphdata->ymin  = -1;
    402         graphdata->xmax = 1;
    403         graphdata->ymax = 1;
    404         psFree(statsX);
    405         psFree(statsY);
    406         return false;
    407     }
    408 
    409     for(int i = 0; i < xVec->n; i++) {
    410         if (!isfinite(xVec->data.F32[i])) continue;
    411         if (xVec->data.F32[i] > xhi || xVec->data.F32[i] < xlo) continue;
    412         graphdata->xmin = PS_MIN (graphdata->xmin, xVec->data.F32[i]);
    413         graphdata->xmax = PS_MAX (graphdata->xmax, xVec->data.F32[i]);
    414     }
    415 
    416     for (int i = 0; i < yVec->n; i++) {
    417         if (!isfinite(xVec->data.F32[i])) continue;
    418         if (yVec->data.F32[i] > yhi || yVec->data.F32[i] < ylo) continue;
    419         graphdata->ymin = PS_MIN (graphdata->ymin, yVec->data.F32[i]);
    420         graphdata->ymax = PS_MAX (graphdata->ymax, yVec->data.F32[i]);
    421     }
    422 
    423     // add a whitespace border
    424     float range = graphdata->xmax - graphdata->xmin;
    425     if (range == 0) range = 1;
    426     graphdata->xmin -= .05 * range;
    427     graphdata->xmax += .05 * range;
    428 
    429     range = graphdata->ymax - graphdata->ymin;
    430     if (range == 0) range = 1;
    431     graphdata->ymin -= .05 * range;
    432     graphdata->ymax += .05 * range;
    433 
    434     psFree (statsX);
    435     psFree (statsY);
     1061bool residPlot (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe,
     1062                        char *title) {
     1063
     1064
     1065    //initialize graph information
     1066    Graphdata graphdata;
     1067    KapaSection section;
     1068
     1069    KapaInitGraph (&graphdata);
     1070    KapaClearPlots (kapa);
     1071
     1072    graphdata.color = KapaColorByName ("black");
     1073    graphdata.ptype = 7;
     1074    graphdata.size = 0.5;
     1075    graphdata.style = 2;
     1076
     1077    section.dx = 0.4;
     1078    section.dy = 0.4;
     1079
     1080    //initialize and populate the plotting vectors
     1081    bool status = false;
     1082    float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MIN");
     1083    float iMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MAX");
     1084    float rMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MIN");
     1085    float rMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MAX");
     1086
     1087    psVector *xVec = psVectorAlloc (match->n, PS_TYPE_F32);
     1088    psVector *yVec = psVectorAlloc (match->n, PS_TYPE_F32);
     1089    psVector *zVec = psVectorAlloc (match->n, PS_TYPE_F32);
     1090
     1091    // X vs dX
     1092    section.x = 0.0;
     1093    section.y = 0.5;
     1094    section.name = NULL;
     1095    psStringAppend (&section.name, "a0");
     1096    KapaSetSection (kapa, &section);
     1097    psFree (section.name);
     1098
     1099    int n = 0;
     1100    for (int i = 0; i < match->n; i++) {
     1101        pmAstromMatch *pair = match->data[i];
     1102        pmAstromObj *raw = rawstars->data[pair->raw];
     1103        pmAstromObj *ref = refstars->data[pair->ref];
     1104
     1105        if (!isfinite(raw->Mag)) continue;
     1106        if (raw->Mag < iMagMin) continue;
     1107        if (raw->Mag > iMagMax) continue;
     1108        if (ref->Mag < rMagMin) continue;
     1109        if (ref->Mag > rMagMax) continue;
     1110
     1111        xVec->data.F32[n] = raw->chip->x;
     1112        yVec->data.F32[n] = raw->chip->x - ref->chip->x;
     1113        zVec->data.F32[n] = raw->Mag;
     1114        n++;
     1115    }
     1116    xVec->n = yVec->n = zVec->n = n;
     1117
     1118    KapaSendLabel (kapa, "X", KAPA_LABEL_XM);
     1119    KapaSendLabel (kapa, "dX", KAPA_LABEL_YM);
     1120    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
     1121
     1122    // X vs dY
     1123    section.x = 0.5;
     1124    section.y = 0.5;
     1125    section.name = NULL;
     1126    psStringAppend (&section.name, "a1");
     1127    KapaSetSection (kapa, &section);
     1128    psFree (section.name);
     1129
     1130    n = 0;
     1131    for (int i = 0; i < match->n; i++) {
     1132        pmAstromMatch *pair = match->data[i];
     1133        pmAstromObj *raw = rawstars->data[pair->raw];
     1134        pmAstromObj *ref = refstars->data[pair->ref];
     1135
     1136        if (!isfinite(raw->Mag)) continue;
     1137        if (raw->Mag < iMagMin) continue;
     1138        if (raw->Mag > iMagMax) continue;
     1139        if (ref->Mag < rMagMin) continue;
     1140        if (ref->Mag > rMagMax) continue;
     1141
     1142        xVec->data.F32[n] = raw->chip->x;
     1143        yVec->data.F32[n] = raw->chip->y - ref->chip->y;
     1144        zVec->data.F32[n] = raw->Mag;
     1145        n++;
     1146    }
     1147    xVec->n = yVec->n = zVec->n = n;
     1148
     1149    KapaSendLabel (kapa, "X", KAPA_LABEL_XM);
     1150    KapaSendLabel (kapa, "dY", KAPA_LABEL_YM);
     1151    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
     1152
     1153    // Y vs dX
     1154    section.x = 0.0;
     1155    section.y = 0.0;
     1156    section.name = NULL;
     1157    psStringAppend (&section.name, "a2");
     1158    KapaSetSection (kapa, &section);
     1159    psFree (section.name);
     1160
     1161    n = 0;
     1162    for (int i = 0; i < match->n; i++) {
     1163        pmAstromMatch *pair = match->data[i];
     1164        pmAstromObj *raw = rawstars->data[pair->raw];
     1165        pmAstromObj *ref = refstars->data[pair->ref];
     1166
     1167        if (!isfinite(raw->Mag)) continue;
     1168        if (raw->Mag < iMagMin) continue;
     1169        if (raw->Mag > iMagMax) continue;
     1170        if (ref->Mag < rMagMin) continue;
     1171        if (ref->Mag > rMagMax) continue;
     1172
     1173        xVec->data.F32[n] = raw->chip->y;
     1174        yVec->data.F32[n] = raw->chip->x - ref->chip->x;
     1175        zVec->data.F32[n] = raw->Mag;
     1176        n++;
     1177    }
     1178    xVec->n = yVec->n = zVec->n = n;
     1179
     1180    KapaSendLabel (kapa, "Y", KAPA_LABEL_XM);
     1181    KapaSendLabel (kapa, "dX", KAPA_LABEL_YM);
     1182    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
     1183
     1184    // Y vs dY
     1185    section.x = 0.5;
     1186    section.y = 0.0;
     1187    section.name = NULL;
     1188    psStringAppend (&section.name, "a3");
     1189    KapaSetSection (kapa, &section);
     1190    psFree (section.name);
     1191
     1192    n = 0;
     1193    for (int i = 0; i < match->n; i++) {
     1194        pmAstromMatch *pair = match->data[i];
     1195        pmAstromObj *raw = rawstars->data[pair->raw];
     1196        pmAstromObj *ref = refstars->data[pair->ref];
     1197
     1198        if (!isfinite(raw->Mag)) continue;
     1199        if (raw->Mag < iMagMin) continue;
     1200        if (raw->Mag > iMagMax) continue;
     1201        if (ref->Mag < rMagMin) continue;
     1202        if (ref->Mag > rMagMax) continue;
     1203
     1204        xVec->data.F32[n] = raw->chip->y;
     1205        yVec->data.F32[n] = raw->chip->y - ref->chip->y;
     1206        zVec->data.F32[n] = raw->Mag;
     1207        n++;
     1208    }
     1209    xVec->n = yVec->n = zVec->n = n;
     1210
     1211    KapaSendLabel (kapa, "Y", KAPA_LABEL_XM);
     1212    KapaSendLabel (kapa, "dY", KAPA_LABEL_YM);
     1213    pmKapaPlotVectorTriple_AutoLimits_OpenGraph (kapa, &graphdata, xVec, yVec, zVec, false);
     1214
     1215    section.x = 0.0;
     1216    section.y = 0.0;
     1217    section.dx = 0.95;
     1218    section.dy = 0.95;
     1219    section.name = NULL;
     1220    psStringAppend (&section.name, "a5");
     1221    KapaSetSection (kapa, &section);
     1222    KapaSendLabel (kapa, title, KAPA_LABEL_XP);
     1223    psFree (section.name);
     1224
     1225    //second window
     1226
     1227    KapaInitGraph (&graphdata);
     1228    KapaClearPlots (kapa2);
     1229
     1230    graphdata.color = KapaColorByName ("black");
     1231    graphdata.ptype = 2;
     1232    graphdata.style = 2;
     1233
     1234    psFree (xVec);
     1235    psFree (yVec);
     1236    psFree (zVec);
     1237
     1238    xVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
     1239    yVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
     1240    zVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
     1241
     1242    // X vs Y by mag (raw)
     1243    n = 0;
     1244    for (int i = 0; i < rawstars->n; i++) {
     1245        pmAstromObj *raw = rawstars->data[i];
     1246        if (!isfinite(raw->Mag)) continue;
     1247        if (raw->Mag < iMagMin) continue;
     1248        if (raw->Mag > iMagMax) continue;
     1249
     1250        xVec->data.F32[n] = raw->chip->x;
     1251        yVec->data.F32[n] = raw->chip->y;
     1252        zVec->data.F32[n] = raw->Mag;
     1253        n++;
     1254    }
     1255    xVec->n = yVec->n = zVec->n = n;
     1256
     1257    KapaSendLabel (kapa2, "X", KAPA_LABEL_XM);
     1258    KapaSendLabel (kapa2, "Y", KAPA_LABEL_YM);
     1259    KapaSendLabel (kapa2,
     1260                   "Chip Coordinates. Black = Raw Stars. Red = Ref Stars. Blue = Matched Stars"
     1261                   , KAPA_LABEL_XP);
     1262    pmVisualTriplePlot (kapa2, &graphdata, xVec, yVec, zVec, false);
     1263
     1264    // X vs Y by mag (ref)
     1265    psFree (xVec);
     1266    psFree (yVec);
     1267    psFree (zVec);
     1268
     1269    xVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
     1270    yVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
     1271    zVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
     1272
     1273    graphdata.color = KapaColorByName ("red");
     1274    graphdata.ptype = 7;
     1275    graphdata.style = 2;
     1276
     1277    n = 0;
     1278    for (int i = 0; i < refstars->n; i++) {
     1279        pmAstromObj *ref = refstars->data[i];
     1280        if (!isfinite(ref->Mag)) continue;
     1281        if (ref->Mag < rMagMin) continue;
     1282        if (ref->Mag > rMagMax) continue;
     1283
     1284        xVec->data.F32[n] = ref->chip->x;
     1285        yVec->data.F32[n] = ref->chip->y;
     1286        zVec->data.F32[n] = ref->Mag;
     1287        n++;
     1288    }
     1289    xVec->n = yVec->n = zVec->n = n;
     1290    pmVisualTripleOverplot (kapa2, &graphdata, xVec, yVec, zVec, false);
     1291
     1292    //rescale the graph to include all points
     1293    float xmin = graphdata.xmin;
     1294    float ymin = graphdata.ymin;
     1295    float xmax = graphdata.xmax;
     1296    float ymax = graphdata.ymax;
     1297    pmVisualScaleGraphdata(&graphdata, xVec, yVec, true);
     1298    graphdata.xmin = PS_MIN(xmin, graphdata.xmin);
     1299    graphdata.ymin = PS_MIN(ymin, graphdata.ymin);
     1300    graphdata.xmax = PS_MAX(xmax, graphdata.xmax);
     1301    graphdata.ymax = PS_MAX(ymax, graphdata.ymax);
     1302    KapaSetLimits (kapa2, &graphdata);
     1303
     1304    //overplot matched stars in blue
     1305    psFree (xVec);
     1306    psFree (yVec);
     1307    psFree (zVec);
     1308
     1309    xVec = psVectorAlloc (match->n, PS_TYPE_F32);
     1310    yVec = psVectorAlloc (match->n, PS_TYPE_F32);
     1311    zVec = psVectorAlloc (match->n, PS_TYPE_F32);
     1312
     1313    graphdata.color = KapaColorByName ("blue");
     1314    n = 0;
     1315    for (int i = 0; i < match->n; i++) {
     1316        pmAstromMatch *pair = match->data[i];
     1317        pmAstromObj *raw = rawstars->data[pair->raw];
     1318        pmAstromObj *ref = refstars->data[pair->ref];
     1319        if (raw->Mag < iMagMin) continue;
     1320        if (raw->Mag > iMagMax) continue;
     1321        if (ref->Mag < rMagMin) continue;
     1322        if (ref->Mag > rMagMax) continue;
     1323
     1324        xVec->data.F32[n] = raw->chip->x;
     1325        yVec->data.F32[n] = raw->chip->y;
     1326        zVec->data.F32[n] = iMagMin;
     1327        n++;
     1328    }
     1329    xVec->n = yVec->n = zVec->n = n;
     1330    pmVisualTripleOverplot (kapa2, &graphdata, xVec, yVec, zVec, false);
     1331
     1332    psFree (xVec);
     1333    psFree (yVec);
     1334    psFree (zVec);
    4361335    return true;
    4371336}
    4381337
    4391338
     1339
     1340
    4401341# else
    4411342
    4421343bool pmAstromSetVisual(bool mode) { return true; }
    443 bool pmAstromVisualInitWindow (int *kapid, char *name) { return true; }
    444 bool pmAstromVisualAskUser (bool *plotflag) { return true; }
    4451344bool pmAstromVisualClose() { return true; }
    4461345bool pmAstromVisualPlotGridMatch (const psArray *raw, const psArray *ref, psImage *gridNP, double offsetX, double offsetY, double maxOffpix, double Scale, double Offset) { return true; }
    4471346bool pmAstromVisualPlotTweak (psVector *xHist, psVector *yHist, int xBin, int yBin) {return true;}
     1347bool pmAstromVisualPlotLuminosityFunction (psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc) {return true;}
     1348bool pmAstromVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe) {return true;}
     1349bool pmAstromVisualPlotRefStars (psArray *refstars, psMetadata *recipe) {return true;}
     1350bool pmAstromVisualPlotRemoveClumps (psArray *input, psImage *count, int scale, float limit) {return true;}
     1351bool pmAstromVisualPlotFixChips (pmFPAfile *input, psVector *xOld, psVector *yOld) {return true;}
     1352bool pmAstromVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) {return true;}
     1353bool pmAstromVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, psVector *cornerPn, psVector *cornerQn, psVector *cornerPd, psVector *cornerQd) {return true;}
     1354bool pmAstromVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) {return true;}
     1355bool pmAstromVisualPlotCommonScale (pmFPA *fpa, psVector *oldScale) {return true;}
     1356bool pmAstromVisualPlotMosaicMatches (psArray *rawstars, psArray *refstars, psArray *match, int iteration, psMetadata *recipe) {return true;}
    4481357
    4491358# endif
Note: See TracChangeset for help on using the changeset viewer.