IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 3, 2010, 8:50:52 AM (16 years ago)
Author:
eugene
Message:

updates from trunk

Location:
branches/simtest_nebulous_branches
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • branches/simtest_nebulous_branches

  • branches/simtest_nebulous_branches/psphot

  • branches/simtest_nebulous_branches/psphot/src

    • Property svn:ignore

      •  

        old new  
        1818psphotVersionDefinitions.h
        1919psphotMomentsStudy
         20psphotPetrosianStudy
         21psphotForced
         22psphotMakePSF
         23psphotStack

  • branches/simtest_nebulous_branches/psphot/src/psphotVisual.c

    r24636 r27840  
    1515# include <kapa.h>
    1616
     17bool pmVisualLimitsFromVectors (Graphdata *graphdata, psVector *xVec, psVector *yVec);
     18
    1719// functions used to visualize the analysis as it goes
    1820// these are invoked by the -visual options
    1921
    20 static int kapa = -1;
     22static int kapa1 = -1;
    2123static int kapa2 = -1;
    2224static int kapa3 = -1;
    2325
     26int psphotKapaChannel (int channel) {
     27
     28    switch (channel) {
     29      case 1:
     30        if (kapa1 == -1) {
     31            kapa1 = KapaOpenNamedSocket ("kapa", "psphot:images");
     32            if (kapa1 == -1) {
     33                fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     34                pmVisualSetVisual(false);
     35            }
     36        }
     37        return kapa1;
     38      case 2:
     39        if (kapa2 == -1) {
     40            kapa2 = KapaOpenNamedSocket ("kapa", "psphot:plots");
     41            if (kapa2 == -1) {
     42                fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     43                pmVisualSetVisual(false);
     44            }
     45        }
     46        return kapa2;
     47      case 3:
     48        if (kapa3 == -1) {
     49            kapa3 = KapaOpenNamedSocket ("kapa", "psphot:stamps");
     50            if (kapa3 == -1) {
     51                fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     52                pmVisualSetVisual(false);
     53            }
     54        }
     55        return kapa3;
     56      default:
     57        psAbort ("unknown kapa channel");
     58    }
     59    psAbort ("unknown kapa channel");
     60}
     61
     62bool psphotVisualEraseOverlays (int channel, char *overlay) {
     63
     64    int myKapa = psphotKapaChannel (channel);
     65    if (!(strcasecmp (overlay, "all"))) {
     66      KiiEraseOverlay (myKapa, "red");
     67      KiiEraseOverlay (myKapa, "green");
     68      KiiEraseOverlay (myKapa, "blue");
     69      KiiEraseOverlay (myKapa, "yellow");
     70      return true;
     71    }
     72    KiiEraseOverlay (myKapa, overlay);
     73    return true;
     74}
    2475
    2576bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel) {
    … …  
    69120}
    70121
    71 bool psphotVisualScaleImage (int kapaFD, psImage *inImage, const char *name, int channel) {
     122bool psphotVisualScaleImage (int kapaFD, psImage *inImage, psImage *inMask, const char *name, int channel) {
    72123
    73124    KiiImage image;
    … …  
    79130    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    80131    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
    81     if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) {
     132    if (!psImageBackground(stats, NULL, inImage, inMask, 0xffff, rng)) {
    82133        fprintf (stderr, "failed to get background values\n");
    83134        return false;
    … …  
    131182    if (!pmVisualIsVisual()) return true;
    132183
    133     if (kapa == -1) {
    134         kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
    135         if (kapa == -1) {
    136             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    137             pmVisualSetVisual(false);
    138             return false;
    139         }
    140     }
    141 
    142     // psphotVisualShowMask (kapa, readout->mask, "mask", 2);
    143     psphotVisualScaleImage (kapa, readout->variance, "variance", 1);
    144     psphotVisualScaleImage (kapa, readout->image, "image", 0);
    145 
    146     // pause and wait for user input:
    147     // continue, save (provide name), ??
    148     char key[10];
    149     fprintf (stdout, "[c]ontinue? ");
    150     if (!fgets(key, 8, stdin)) {
    151         psWarning("Unable to read option");
    152     }
     184    int kapa = psphotKapaChannel (1);
     185    if (kapa == -1) return false;
     186
     187    psphotVisualShowMask (kapa, readout->mask, "mask", 2);
     188    psphotVisualScaleImage (kapa, readout->variance, readout->mask, "variance", 1);
     189    psphotVisualScaleImage (kapa, readout->image, readout->mask, "image", 0);
     190
     191    pmVisualAskUser(NULL);
    153192    return true;
    154193}
    … …  
    160199    if (!pmVisualIsVisual()) return true;
    161200
    162     if (kapa == -1) {
    163         kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
    164         if (kapa == -1) {
    165             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    166             pmVisualSetVisual(false);
    167             return false;
    168         }
    169     }
     201    int kapa = psphotKapaChannel (1);
     202    if (kapa == -1) return false;
    170203
    171204    bool status = false;
    … …  
    178211    }
    179212
    180     psphotVisualScaleImage (kapa, backgnd->image, "backgnd", 2);
    181     psphotVisualScaleImage (kapa, readout->image, "backsub", 0);
    182 
    183     // pause and wait for user input:
    184     // continue, save (provide name), ??
    185     char key[10];
    186     fprintf (stdout, "[c]ontinue? ");
    187     if (!fgets(key, 8, stdin)) {
    188         psWarning("Unable to read option");
    189     }
    190     return true;
    191 }
    192 
    193 bool psphotVisualShowSignificance (psImage *image) {
     213    psphotVisualScaleImage (kapa, backgnd->image, readout->mask, "backgnd", 2);
     214    psphotVisualScaleImage (kapa, readout->image, readout->mask, "backsub", 0);
     215
     216    pmVisualAskUser(NULL);
     217    return true;
     218}
     219
     220bool psphotVisualShowSignificance (psImage *image, float min, float max) {
    194221
    195222    if (!pmVisualIsVisual()) return true;
    196223
    197     if (kapa == -1) {
    198         kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
    199         if (kapa == -1) {
    200             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    201             pmVisualSetVisual(false);
    202             return false;
    203         }
    204     }
    205 
    206     // XXX test: image->data.F32[10][10] = 10000;
    207     psphotVisualRangeImage (kapa, image, "signif", 2, -1.0, 25.0*25.0);
    208 
    209     // pause and wait for user input:
    210     // continue, save (provide name), ??
    211     char key[10];
    212     fprintf (stdout, "[c]ontinue? ");
    213     if (!fgets(key, 8, stdin)) {
    214         psWarning("Unable to read option");
    215     }
     224    int kapa = psphotKapaChannel (1);
     225    if (kapa == -1) return false;
     226
     227    psphotVisualRangeImage (kapa, image, "signif", 2, min, max);
     228
     229    pmVisualAskUser(NULL);
    216230    return true;
    217231}
    … …  
    224238    if (!pmVisualIsVisual()) return true;
    225239
    226     if (kapa == -1) {
    227         fprintf (stderr, "kapa not opened, skipping\n");
    228         return false;
    229     }
     240    int kapa = psphotKapaChannel (1);
     241    if (kapa == -1) return false;
    230242
    231243    psArray *peaks = detections->peaks;
    … …  
    233245    // note: this uses the Ohana allocation tools:
    234246    // ALLOCATE (overlay, KiiOverlay, 3*peaks->n + 1);
    235     ALLOCATE (overlay, KiiOverlay, peaks->n);
     247    ALLOCATE (overlay, KiiOverlay, peaks->n + 2);
    236248
    237249    Noverlay = 0;
    … …  
    249261        overlay[Noverlay].text = NULL;
    250262        Noverlay ++;
    251 
    252 # if (0)
    253         overlay[Noverlay].type = KII_OVERLAY_BOX;
    254         overlay[Noverlay].x = peak->x;
    255         overlay[Noverlay].y = peak->y;
    256         overlay[Noverlay].dx = 1.0;
    257         overlay[Noverlay].dy = 1.0;
    258         overlay[Noverlay].angle = 0.0;
    259         overlay[Noverlay].text = NULL;
    260         Noverlay ++;
    261 
    262         overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    263         overlay[Noverlay].x = peak->xf;
    264         overlay[Noverlay].y = peak->yf;
    265         overlay[Noverlay].dx = 2.0;
    266         overlay[Noverlay].dy = 2.0;
    267         overlay[Noverlay].angle = 0.0;
    268         overlay[Noverlay].text = NULL;
    269         Noverlay ++;
    270 # endif
    271     }
    272 
    273 # if (0)
    274     overlay[Noverlay].type = KII_OVERLAY_BOX;
    275     overlay[Noverlay].x = 10.0;
    276     overlay[Noverlay].y = 10.0;
    277     overlay[Noverlay].dx = 0.5;
    278     overlay[Noverlay].dy = 0.5;
    279     overlay[Noverlay].angle = 0.0;
    280     overlay[Noverlay].text = NULL;
    281     Noverlay ++;
    282 # endif
     263    }
    283264
    284265    KiiLoadOverlay (kapa, overlay, Noverlay, "red");
    285266    FREE (overlay);
    286267
    287     // pause and wait for user input:
    288     // continue, save (provide name), ??
    289     char key[10];
    290     fprintf (stdout, "[c]ontinue? ");
    291     if (!fgets(key, 8, stdin)) {
    292         psWarning("Unable to read option");
    293     }
     268    pmVisualAskUser(NULL);
    294269    return true;
    295270}
    … …  
    302277    if (!pmVisualIsVisual()) return true;
    303278
    304     if (kapa == -1) {
    305         fprintf (stderr, "kapa not opened, skipping\n");
    306         return false;
    307     }
     279    int kapa = psphotKapaChannel (1);
     280    if (kapa == -1) return false;
    308281
    309282    psArray *footprints = detections->footprints;
    … …  
    325298
    326299        // draw the top
     300        // XXX need to allow top (and bottom) to have more than one span
    327301        span = footprint->spans->data[0];
    328302        overlay[Noverlay].type = KII_OVERLAY_LINE;
    … …  
    399373    FREE (overlay);
    400374
    401     // pause and wait for user input:
    402     // continue, save (provide name), ??
    403     char key[10];
    404     fprintf (stdout, "[c]ontinue? ");
    405     if (!fgets(key, 8, stdin)) {
    406         psWarning("Unable to read option");
    407     }
     375    pmVisualAskUser(NULL);
    408376    return true;
    409377}
    … …  
    419387    if (!pmVisualIsVisual()) return true;
    420388
    421     if (kapa == -1) {
    422         fprintf (stderr, "kapa not opened, skipping\n");
    423         return false;
    424     }
     389    int kapa = psphotKapaChannel (1);
     390    if (kapa == -1) return false;
     391
     392    // XXX mark the different source classes with different color/shape dots
     393    // XXX are moments S/N and peak S/N consistent?
    425394
    426395    // note: this uses the Ohana allocation tools:
    … …  
    448417        overlay[Noverlay].dx = 2.0*axes.major;
    449418        overlay[Noverlay].dy = 2.0*axes.minor;
    450         overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD;  // XXXXXXXX the axes angle is negative to display of object on kapa
     419
     420        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
     421
    451422        overlay[Noverlay].text = NULL;
    452423        Noverlay ++;
    … …  
    456427    FREE (overlay);
    457428
    458     // pause and wait for user input:
    459     // continue, save (provide name), ??
    460     char key[10];
    461     fprintf (stdout, "[c]ontinue? ");
    462     if (!fgets(key, 8, stdin)) {
    463         psWarning("Unable to read option");
    464     }
    465 
    466     return true;
    467 }
    468 
    469 bool psphotVisualPlotMoments (psMetadata *recipe, psArray *sources) {
     429    pmVisualAskUser(NULL);
     430    return true;
     431}
     432
     433bool psphotVisualPlotMoments (psMetadata *recipe, psMetadata *analysis, psArray *sources) {
    470434
    471435    bool status;
    … …  
    475439    if (!pmVisualIsVisual()) return true;
    476440
    477     if (kapa3 == -1) {
    478         kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
    479         if (kapa3 == -1) {
    480             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    481             pmVisualSetVisual(false);
    482             return false;
    483         }
    484     }
    485 
    486     KapaClearPlots (kapa3);
     441    int myKapa = psphotKapaChannel (2);
     442    if (myKapa == -1) return false;
     443
     444    KapaClearPlots (myKapa);
    487445    KapaInitGraph (&graphdata);
    488     KapaSetFont (kapa3, "courier", 14);
     446    KapaSetFont (myKapa, "courier", 14);
     447
     448    section.bg = KapaColorByName ("none"); // XXX probably should be 'none'
    489449
    490450    float SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM");
    491451
    492452    // select the max psfX,Y values for the plot limits
    493     float Xmin = 0.0, Xmax = 0.0;
    494     float Ymin = 0.0, Ymax = 0.0;
     453    float Xmin = 1000.0, Xmax = 0.0;
     454    float Ymin = 1000.0, Ymax = 0.0;
    495455    {
    496         int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
     456        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
    497457        for (int n = 0; n < nRegions; n++) {
    498458
    499459            char regionName[64];
    500460            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
    501             psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
    502 
    503             float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
     461            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
     462
     463            float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
    504464            float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
    505465            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
    506466            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
    507467
    508             float X0 = psfX - 4.0*psfdX;
    509             float X1 = psfX + 4.0*psfdX;
    510             float Y0 = psfY - 4.0*psfdY;
    511             float Y1 = psfY + 4.0*psfdY;
    512 
    513             if (isfinite(X0)) { Xmin = PS_MAX(Xmin, X0); }
    514             if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
    515             if (isfinite(Y0)) { Ymin = PS_MAX(Ymin, Y0); }
    516             if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
    517         }
    518     }
     468            float X0 = psfX - 4.0*psfdX;
     469            float X1 = psfX + 4.0*psfdX;
     470            float Y0 = psfY - 4.0*psfdY;
     471            float Y1 = psfY + 4.0*psfdY;
     472
     473            if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); }
     474            if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
     475            if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); }
     476            if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
     477        }
     478    }
     479    Xmin = PS_MAX(Xmin, -0.1);
     480    Ymin = PS_MAX(Ymin, -0.1);
     481
     482    // XXX test: hardwire plot limits
     483    // Xmin = -0.1; Ymin = -0.1;
     484    // Xmax = 20.1; Ymax = 20.1;
    519485
    520486    // storage vectors for data to be plotted
    … …  
    564530    section.y  = 0.00;
    565531    section.name = psStringCopy ("MxxMyy");
    566     KapaSetSection (kapa3, &section);
     532    KapaSetSection (myKapa, &section);
    567533    psFree (section.name);
    568534
    … …  
    572538    graphdata.xmax = Xmax;
    573539    graphdata.ymax = Ymax;
    574     KapaSetLimits (kapa3, &graphdata);
    575 
    576     KapaBox (kapa3, &graphdata);
    577     KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM);
    578     KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM);
     540    KapaSetLimits (myKapa, &graphdata);
     541
     542    graphdata.padXm = NAN;
     543    graphdata.padYm = NAN;
     544    graphdata.padXp = 0.5;
     545    graphdata.padYp = 0.5;
     546    KapaBox (myKapa, &graphdata);
     547
     548    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
     549    KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM);
    579550
    580551    graphdata.color = KapaColorByName ("black");
    … …  
    582553    graphdata.size = 0.3;
    583554    graphdata.style = 2;
    584     KapaPrepPlot (kapa3, nF, &graphdata);
    585     KapaPlotVector (kapa3, nF, xFaint->data.F32, "x");
    586     KapaPlotVector (kapa3, nF, yFaint->data.F32, "y");
     555    KapaPrepPlot (myKapa, nF, &graphdata);
     556
     557    KapaPlotVector (myKapa, nF, xFaint->data.F32, "x");
     558    KapaPlotVector (myKapa, nF, yFaint->data.F32, "y");
    587559
    588560    graphdata.color = KapaColorByName ("red");
    … …  
    590562    graphdata.size = 0.5;
    591563    graphdata.style = 2;
    592     KapaPrepPlot (kapa3, nB, &graphdata);
    593     KapaPlotVector (kapa3, nB, xBright->data.F32, "x");
    594     KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
     564    KapaPrepPlot (myKapa, nB, &graphdata);
     565    KapaPlotVector (myKapa, nB, xBright->data.F32, "x");
     566    KapaPlotVector (myKapa, nB, yBright->data.F32, "y");
    595567
    596568    // second section: MagMyy
    … …  
    598570    section.dy = 0.25;
    599571    section.x  = 0.00;
    600     section.y  = 0.80;
     572    section.y  = 0.75;
    601573    section.name = psStringCopy ("MagMyy");
    602     KapaSetSection (kapa3, &section);
     574    KapaSetSection (myKapa, &section);
    603575    psFree (section.name);
    604576
    … …  
    608580    graphdata.ymin = Ymin;
    609581    graphdata.ymax = Ymax;
    610     KapaSetLimits (kapa3, &graphdata);
    611 
     582    KapaSetLimits (myKapa, &graphdata);
     583
     584    graphdata.padXm = 0.5;
     585    graphdata.padYm = NAN;
     586    graphdata.padXp = NAN;
     587    graphdata.padYp = 0.5;
    612588    strcpy (graphdata.labels, "0210");
    613     KapaBox (kapa3, &graphdata);
    614     KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_XP);
    615     KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM);
     589    KapaBox (myKapa, &graphdata);
     590    KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_XP);
     591    KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM);
    616592
    617593    graphdata.color = KapaColorByName ("black");
    … …  
    619595    graphdata.size = 0.3;
    620596    graphdata.style = 2;
    621     KapaPrepPlot (kapa3, nF, &graphdata);
    622     KapaPlotVector (kapa3, nF, mFaint->data.F32, "x");
    623     KapaPlotVector (kapa3, nF, yFaint->data.F32, "y");
     597    KapaPrepPlot (myKapa, nF, &graphdata);
     598    KapaPlotVector (myKapa, nF, mFaint->data.F32, "x");
     599    KapaPlotVector (myKapa, nF, yFaint->data.F32, "y");
    624600
    625601    graphdata.color = KapaColorByName ("red");
    … …  
    627603    graphdata.size = 0.5;
    628604    graphdata.style = 2;
    629     KapaPrepPlot (kapa3, nB, &graphdata);
    630     KapaPlotVector (kapa3, nB, mBright->data.F32, "x");
    631     KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
     605    KapaPrepPlot (myKapa, nB, &graphdata);
     606    KapaPlotVector (myKapa, nB, mBright->data.F32, "x");
     607    KapaPlotVector (myKapa, nB, yBright->data.F32, "y");
    632608
    633609    // third section: MagMxx
    634610    section.dx = 0.25;
    635611    section.dy = 0.75;
    636     section.x  = 0.80;
     612    section.x  = 0.75;
    637613    section.y  = 0.00;
    638614    section.name = psStringCopy ("MagMxx");
    639     KapaSetSection (kapa3, &section);
     615    KapaSetSection (myKapa, &section);
    640616    psFree (section.name);
    641617
    … …  
    645621    graphdata.ymin =  -7.9;
    646622    graphdata.ymax = -17.1;
    647     KapaSetLimits (kapa3, &graphdata);
    648 
     623    KapaSetLimits (myKapa, &graphdata);
     624
     625    graphdata.padXm = NAN;
     626    graphdata.padYm = 0.5;
     627    graphdata.padXp = 0.5;
     628    graphdata.padYp = NAN;
    649629    strcpy (graphdata.labels, "2001");
    650     KapaBox (kapa3, &graphdata);
    651     KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM);
    652     KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_YP);
     630    KapaBox (myKapa, &graphdata);
     631    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
     632    KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_YP);
    653633
    654634    graphdata.color = KapaColorByName ("black");
    … …  
    656636    graphdata.size = 0.3;
    657637    graphdata.style = 2;
    658     KapaPrepPlot (kapa3, nF, &graphdata);
    659     KapaPlotVector (kapa3, nF, xFaint->data.F32, "x");
    660     KapaPlotVector (kapa3, nF, mFaint->data.F32, "y");
     638    KapaPrepPlot (myKapa, nF, &graphdata);
     639    KapaPlotVector (myKapa, nF, xFaint->data.F32, "x");
     640    KapaPlotVector (myKapa, nF, mFaint->data.F32, "y");
    661641
    662642    graphdata.color = KapaColorByName ("red");
    … …  
    664644    graphdata.size = 0.5;
    665645    graphdata.style = 2;
    666     KapaPrepPlot (kapa3, nB, &graphdata);
    667     KapaPlotVector (kapa3, nB, xBright->data.F32, "x");
    668     KapaPlotVector (kapa3, nB, mBright->data.F32, "y");
     646    KapaPrepPlot (myKapa, nB, &graphdata);
     647    KapaPlotVector (myKapa, nB, xBright->data.F32, "x");
     648    KapaPlotVector (myKapa, nB, mBright->data.F32, "y");
    669649
    670650    // draw N circles to outline the clumps
    671651    {
    672         KapaSelectSection (kapa3, "MxxMyy");
     652        KapaSelectSection (myKapa, "MxxMyy");
    673653
    674654        // draw a circle centered on psfX,Y with size of the psf limit
    … …  
    676656        psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
    677657
    678         int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
     658        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
    679659        float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
    680660
    … …  
    682662        graphdata.style = 0;
    683663
    684         graphdata.xmin = Xmin;
    685         graphdata.ymin = Ymin;
    686         graphdata.xmax = Xmax;
    687         graphdata.ymax = Ymax;
    688         KapaSetLimits (kapa3, &graphdata);
     664        graphdata.xmin = Xmin;
     665        graphdata.ymin = Ymin;
     666        graphdata.xmax = Xmax;
     667        graphdata.ymax = Ymax;
     668        KapaSetLimits (myKapa, &graphdata);
    689669
    690670        for (int n = 0; n < nRegions; n++) {
    … …  
    692672            char regionName[64];
    693673            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
    694             psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
     674            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
    695675
    696676            float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
    … …  
    705685                yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
    706686            }
    707             KapaPrepPlot (kapa3, xLimit->n, &graphdata);
    708             KapaPlotVector (kapa3, xLimit->n, xLimit->data.F32, "x");
    709             KapaPlotVector (kapa3, yLimit->n, yLimit->data.F32, "y");
     687            KapaPrepPlot (myKapa, xLimit->n, &graphdata);
     688            KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x");
     689            KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y");
    710690        }
    711691        psFree (xLimit);
    712692        psFree (yLimit);
    713693    }
    714 
    715 # if (0)
    716     // *** make a histogram of the source counts in the x and y directions
    717     psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0);
    718     psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0);
    719     psVectorHistogram (nX, xFaint, NULL, NULL, 0);
    720     psVectorHistogram (nY, yFaint, NULL, NULL, 0);
    721     psVector *dX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
    722     psVector *vX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
    723     psVector *dY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
    724     psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
    725     for (int i = 0; i < nX->nums->n; i++) {
    726         dX->data.F32[i] = nX->nums->data.S32[i];
    727         vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]);
    728     }
    729     for (int i = 0; i < nY->nums->n; i++) {
    730         dY->data.F32[i] = nY->nums->data.S32[i];
    731         vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]);
    732     }
    733 
    734     graphdata.color = KapaColorByName ("black");
    735     graphdata.ptype = 0;
    736     graphdata.size = 0.0;
    737     graphdata.style = 0;
    738     KapaPrepPlot (kapa3, dX->n, &graphdata);
    739     KapaPlotVector (kapa3, dX->n, dX->data.F32, "x");
    740     KapaPlotVector (kapa3, vX->n, vX->data.F32, "y");
    741 
    742     psFree (nX);
    743     psFree (dX);
    744     psFree (vX);
    745 
    746     psFree (nY);
    747     psFree (dY);
    748     psFree (vY);
    749 # endif
    750694
    751695    psFree (xBright);
    … …  
    756700    psFree (mFaint);
    757701
    758     // pause and wait for user input:
    759     // continue, save (provide name), ??
    760     char key[10];
    761     fprintf (stdout, "[c]ontinue? ");
    762     if (!fgets(key, 8, stdin)) {
    763         psWarning("Unable to read option");
    764     }
     702    pmVisualAskUser(NULL);
    765703    return true;
    766704}
    767705
    768706// assumes 'kapa' value is checked and set
    769 bool psphotVisualShowRoughClass_Single (psArray *sources, pmSourceType type, pmSourceMode mode, char *color) {
     707bool psphotVisualShowRoughClass_Single (int myKapa, psArray *sources, pmSourceType type, pmSourceMode mode, char *color) {
    770708
    771709    int Noverlay;
    … …  
    802740        overlay[Noverlay].dx = 2.0*axes.major;
    803741        overlay[Noverlay].dy = 2.0*axes.minor;
    804         overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD;
     742        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
    805743        overlay[Noverlay].text = NULL;
    806744        Noverlay ++;
    807745    }
    808746
    809     KiiLoadOverlay (kapa, overlay, Noverlay, color);
     747    KiiLoadOverlay (myKapa, overlay, Noverlay, color);
    810748    FREE (overlay);
    811749
    … …  
    817755    if (!pmVisualIsVisual()) return true;
    818756
    819     if (kapa == -1) {
    820         fprintf (stderr, "kapa not opened, skipping\n");
    821         return false;
    822     }
    823 
    824     KiiEraseOverlay (kapa, "yellow"); // moments
    825 
    826     psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_STAR, 0, "red");
    827     psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_EXTENDED, 0, "blue");
    828     psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_DEFECT, 0, "blue");
    829     psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_SATURATED, 0, "red");
    830     psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_PSFSTAR, "yellow");
    831     psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_SATSTAR, "green");
    832 
    833     // pause and wait for user input:
    834     // continue, save (provide name), ??
    835     char key[10];
     757    int myKapa = psphotKapaChannel (1);
     758    if (myKapa == -1) return false;
     759
     760    KiiEraseOverlay (myKapa, "yellow"); // moments
     761
     762    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_STAR, 0, "red");
     763    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_EXTENDED, 0, "blue");
     764    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_DEFECT, 0, "blue");
     765    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_SATURATED, 0, "red");
     766    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_PSFSTAR, "yellow");
     767    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_SATSTAR, "green");
     768
    836769    fprintf (stdout, "red: STAR or SAT AREA; blue: EXTENDED or DEFECT; green: SATSTAR; yellow: PSFSTAR\n");
    837     fprintf (stdout, "[c]ontinue? ");
    838     if (!fgets(key, 8, stdin)) {
    839         psWarning("Unable to read option");
    840     }
     770    pmVisualAskUser(NULL);
    841771    return true;
    842772}
    … …  
    846776    if (!pmVisualIsVisual()) return true;
    847777
    848     if (kapa2 == -1) {
    849         kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars");
    850         if (kapa2 == -1) {
    851             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    852             pmVisualSetVisual(false);
    853             return false;
    854         }
    855     }
     778    int myKapa = psphotKapaChannel (3);
     779    if (myKapa == -1) return false;
    856780
    857781    int DX = 64;
    … …  
    898822
    899823    psImage *psfLogFlux = (psImage *) psUnaryOp (NULL, psfMosaic, "log");
    900     psphotVisualRangeImage (kapa2, psfLogFlux, "psf_mosaic",    0, -2.0, 3.0);
    901     psphotVisualRangeImage (kapa2, funMosaic, "psf_analytical", 1, -10.0, 100.0);
    902     psphotVisualRangeImage (kapa2, resMosaic, "psf_residual",   2, -10.0, 100.0);
     824    psphotVisualRangeImage (myKapa, psfLogFlux, "psf_mosaic",    0, -2.0, 3.0);
     825    psphotVisualRangeImage (myKapa, funMosaic, "psf_analytical", 1, -10.0, 100.0);
     826    psphotVisualRangeImage (myKapa, resMosaic, "psf_residual",   2, -10.0, 100.0);
    903827
    904828    psFree (psfMosaic);
    … …  
    908832    psFree (modelRef);
    909833
    910     // pause and wait for user input:
    911     // continue, save (provide name), ??
    912     char key[10];
    913     fprintf (stdout, "[c]ontinue? ");
    914     if (!fgets(key, 8, stdin)) {
    915         psWarning("Unable to read option");
    916     }
     834    pmVisualAskUser(NULL);
    917835    return true;
    918836}
    … …  
    924842    if (!pmVisualIsVisual()) return true;
    925843
    926     if (kapa2 == -1) {
    927         kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars");
    928         if (kapa2 == -1) {
    929             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    930             pmVisualSetVisual(false);
    931             return false;
    932         }
    933     }
     844    int myKapa = psphotKapaChannel (3);
     845    if (myKapa == -1) return false;
    934846
    935847    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    1017929            if (Xo == 0) {
    1018930                // place source alone on this row
    1019                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     931                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     932                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    1020933                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
    1021934
    1022                 psphotSubWithTest (source, false, maskVal); // remove source (force)
     935                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    1023936                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
    1024                 psphotSetState (source, false, maskVal); // reset source Add/Sub state to recorded
     937
     938                if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    1025939
    1026940                Yo += DY;
    … …  
    1032946                Xo = 0;
    1033947
    1034                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     948                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     949                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    1035950                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
    1036951
    1037                 psphotSubWithTest (source, false, maskVal); // remove source (force)
     952                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    1038953                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
    1039                 psphotSetState (source, false, maskVal); // replace source (has been subtracted)
     954
     955                if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    1040956
    1041957                Xo = DX;
    … …  
    1044960        } else {
    1045961            // extend this row
    1046             psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     962            bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     963            if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    1047964            psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
    1048965
    1049             psphotSubWithTest (source, false, maskVal); // remove source (force)
     966            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    1050967            psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
    1051             psphotSetState (source, false, maskVal); // replace source (has been subtracted)
     968            if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    1052969
    1053970            Xo += DX;
    … …  
    1056973    }
    1057974
    1058     psphotVisualRangeImage (kapa2, outpos, "psfpos", 0, -0.05, 0.95);
    1059     psphotVisualRangeImage (kapa2, outsub, "psfsub", 1, -0.05, 0.95);
    1060 
    1061     // pause and wait for user input:
    1062     // continue, save (provide name), ??
    1063     char key[10];
    1064     fprintf (stdout, "[c]ontinue? ");
    1065     if (!fgets(key, 8, stdin)) {
    1066         psWarning("Unable to read option");
    1067     }
    1068 
     975    psphotVisualRangeImage (myKapa, outpos, "psfpos", 0, -0.05, 0.95);
     976    psphotVisualRangeImage (myKapa, outsub, "psfsub", 1, -0.05, 0.95);
     977
     978    pmVisualAskUser(NULL);
    1069979    psFree (outpos);
    1070980    psFree (outsub);
    … …  
    1084994    if (!pmVisualIsVisual()) return true;
    1085995
    1086     if (kapa2 == -1) {
    1087         kapa2 = KapaOpenNamedSocket ("kapa", "psphot:images");
    1088         if (kapa2 == -1) {
    1089             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1090             pmVisualSetVisual(false);
    1091             return false;
    1092         }
    1093     }
     996    int myKapa = psphotKapaChannel (3);
     997    if (myKapa == -1) return false;
    1094998
    1095999    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    11191023        pmSource *source = sources->data[i];
    11201024
    1121         bool keep = false;
    1122         keep |= (source->mode & PM_SOURCE_MODE_SATSTAR);
    1123         if (!keep) continue;
     1025        // only show "real" saturated stars (not defects)
     1026        if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;;
     1027        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;;
    11241028
    11251029        // how does this subimage get placed into the output image?
    … …  
    11641068        pmSource *source = sources->data[i];
    11651069
    1166         bool keep = false;
    1167         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    1168             nSAT ++;
    1169             keep = true;
    1170         }
    1171         if (!keep) continue;
     1070        // only show "real" saturated stars (not defects)
     1071        if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;;
     1072        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;;
     1073        nSAT ++;
    11721074
    11731075        if (Xo + DX > NX) {
    … …  
    11751077            if (Xo == 0) {
    11761078                // place source alone on this row
    1177                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     1079                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     1080                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    11781081                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
    1179                 psphotSetState (source, true, maskVal); // reset source Add/Sub state to recorded
     1082                if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    11801083
    11811084                Yo += DY;
    … …  
    11861089                Yo += dY;
    11871090                Xo = 0;
    1188                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     1091
     1092                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     1093                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    11891094                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
    1190                 psphotSetState (source, true, maskVal); // replace source (has been subtracted)
     1095                if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    11911096
    11921097                Xo = DX;
    … …  
    11951100        } else {
    11961101            // extend this row
    1197             psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     1102            bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     1103            if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    11981104            psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
    1199             psphotSetState (source, true, maskVal); // replace source (has been subtracted)
     1105            if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    12001106
    12011107            Xo += DX;
    … …  
    12041110    }
    12051111
    1206     psphotVisualScaleImage (kapa2, outsat, "satstar", 2);
    1207 
    1208     // pause and wait for user input:
    1209     // continue, save (provide name), ??
    1210     char key[10];
    1211     fprintf (stdout, "[c]ontinue? ");
    1212     if (!fgets(key, 8, stdin)) {
    1213         psWarning("Unable to read option");
    1214     }
    1215 
     1112    psphotVisualScaleImage (myKapa, outsat, NULL, "satstar", 2);
     1113
     1114    pmVisualAskUser(NULL);
    12161115    psFree (outsat);
    12171116    return true;
    12181117}
    12191118
     1119static void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1)
     1120{
     1121    float x[2], y[2];
     1122    x[0] = x0;
     1123    x[1] = x1;
     1124    y[0] = y0;
     1125    y[1] = y1;
     1126    KapaPrepPlot   (myKapa, 2, graphdata);
     1127    KapaPlotVector (myKapa, 2, x, "x");
     1128    KapaPlotVector (myKapa, 2, y, "y");
     1129}
     1130
    12201131bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal) {
    12211132
    12221133    Graphdata graphdata;
    12231134
    1224     bool state = !(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
    1225     psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     1135    bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     1136    if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    12261137
    12271138    int nPts = source->pixels->numRows * source->pixels->numCols;
    … …  
    12401151        for (int ix = 0; ix < source->pixels->numCols; ix++) {
    12411152            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
    1242                 // rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
    1243                 rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
     1153                rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
     1154                // rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
    12441155                Rb->data.F32[nb] = log10(rb->data.F32[nb]);
    12451156                fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]);
    12461157                nb++;
    12471158            } else {
    1248                 // rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
    1249                 rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ;
     1159                rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
     1160                // rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ;
    12501161                Rg->data.F32[ng] = log10(rg->data.F32[ng]);
    12511162                fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]);
    … …  
    12551166    }
    12561167
     1168    // generate model profiles (major and minor axis):
     1169    // create a model with theta = 0.0 so major and minor axes are equiv to x and y:
     1170    psEllipseShape rawShape, rotShape;
     1171
     1172    rawShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
     1173    rawShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
     1174    rawShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
     1175
     1176    psEllipseAxes axes = psEllipseShapeToAxes (rawShape, 20.0);
     1177
     1178    axes.theta = 0.0;
     1179
     1180    rotShape = psEllipseAxesToShape (axes);
     1181
     1182    psVector *params = psVectorAlloc(source->modelPSF->params->n, PS_TYPE_F32);
     1183    for (int i = 0; i < source->modelPSF->params->n; i++) {
     1184        params->data.F32[i] = source->modelPSF->params->data.F32[i];
     1185    }
     1186    params->data.F32[PM_PAR_SXX] = rotShape.sx * M_SQRT2;
     1187    params->data.F32[PM_PAR_SYY] = rotShape.sy * M_SQRT2;
     1188    params->data.F32[PM_PAR_SXY] = rotShape.sxy;
     1189    params->data.F32[PM_PAR_XPOS] = 0.0;
     1190    params->data.F32[PM_PAR_YPOS] = 0.0;
     1191
     1192    psVector *rmod = psVectorAlloc(300, PS_TYPE_F32);
     1193    psVector *fmaj = psVectorAlloc(300, PS_TYPE_F32);
     1194    psVector *fmin = psVectorAlloc(300, PS_TYPE_F32);
     1195
     1196    psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     1197
     1198    float r = 0.0;
     1199    for (int i = 0; i < rmod->n; i++) {
     1200        r = i*0.1;
     1201        rmod->data.F32[i] = r;
     1202
     1203        coord->data.F32[1] = r;
     1204        coord->data.F32[0] = 0.0;
     1205        fmaj->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
     1206
     1207        coord->data.F32[0] = r;
     1208        coord->data.F32[1] = 0.0;
     1209        fmin->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
     1210    }
     1211    psFree (coord);
     1212    psFree (params);
     1213
     1214    float FWHM_MAJOR = 2.0*source->modelPSF->modelRadius (source->modelPSF->params, 0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
     1215    float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
     1216    if (FWHM_MAJOR < FWHM_MINOR) PS_SWAP (FWHM_MAJOR, FWHM_MINOR);
     1217
     1218    psEllipseMoments emoments;
     1219    emoments.x2 = source->moments->Mxx;
     1220    emoments.xy = source->moments->Mxy;
     1221    emoments.y2 = source->moments->Myy;
     1222    axes = psEllipseMomentsToAxes (emoments, 20.0);
     1223    float MOMENTS_MAJOR = 2.355*axes.major;
     1224    float MOMENTS_MINOR = 2.355*axes.minor;
     1225
     1226    float logHM = log10(0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
     1227
    12571228    // reset source Add/Sub state to recorded
    1258     psphotSetState (source, state, maskVal);
     1229    if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    12591230
    12601231    KapaInitGraph (&graphdata);
    … …  
    12911262    KapaPlotVector (myKapa, nb, fb->data.F32, "y");
    12921263
     1264    graphdata.color = KapaColorByName ("blue");
     1265    graphdata.ptype = 0;
     1266    graphdata.size = 0.0;
     1267    graphdata.style = 0;
     1268    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1269    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1270    KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
     1271    plotline (myKapa, &graphdata, 0.0, logHM, 30.0, logHM);
     1272    plotline (myKapa, &graphdata, 0.5*FWHM_MINOR, 0.0, 0.5*FWHM_MINOR, 5.0);
     1273    graphdata.ltype = 1;
     1274    plotline (myKapa, &graphdata, 0.5*MOMENTS_MINOR, 0.0, 0.5*MOMENTS_MINOR, 5.0);
     1275    graphdata.ltype = 0;
     1276       
     1277    graphdata.color = KapaColorByName ("green");
     1278    graphdata.ptype = 0;
     1279    graphdata.size = 0.0;
     1280    graphdata.style = 0;
     1281    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1282    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1283    KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
     1284    plotline (myKapa, &graphdata, 0.5*FWHM_MAJOR, 0.0, 0.5*FWHM_MAJOR, 5.0);
     1285    graphdata.ltype = 1;
     1286    plotline (myKapa, &graphdata, 0.5*MOMENTS_MAJOR, 0.0, 0.5*MOMENTS_MAJOR, 5.0);
     1287    graphdata.ltype = 0;
     1288       
     1289    for (int i = 0; i < rmod->n; i++) {
     1290        rmod->data.F32[i] = log10(rmod->data.F32[i]);
     1291    }
     1292
    12931293    // ** loglog **
    12941294    KapaSelectSection (myKapa, "loglog");
    … …  
    12991299    graphdata.ymin = -0.05;
    13001300    graphdata.ymax = +5.05;
     1301    graphdata.color = KapaColorByName ("black");
    13011302    KapaSetLimits (myKapa, &graphdata);
    13021303
    … …  
    13211322    KapaPlotVector (myKapa, nb, Rb->data.F32, "x");
    13221323    KapaPlotVector (myKapa, nb, fb->data.F32, "y");
     1324
     1325    graphdata.color = KapaColorByName ("blue");
     1326    graphdata.ptype = 0;
     1327    graphdata.size = 0.0;
     1328    graphdata.style = 0;
     1329    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1330    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1331    KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
     1332
     1333    graphdata.color = KapaColorByName ("green");
     1334    graphdata.ptype = 0;
     1335    graphdata.size = 0.0;
     1336    graphdata.style = 0;
     1337    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1338    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1339    KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
     1340
     1341    psFree (rmod);
     1342    psFree (fmin);
     1343    psFree (fmaj);
    13231344
    13241345    psFree (rg);
    … …  
    13371358    if (!pmVisualIsVisual()) return true;
    13381359
    1339     if (kapa3 == -1) {
    1340         kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
    1341         if (kapa3 == -1) {
    1342             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1343             pmVisualSetVisual(false);
    1344             return false;
    1345         }
    1346     }
     1360    int myKapa = psphotKapaChannel (2);
     1361    if (myKapa == -1) return false;
    13471362
    13481363    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    13511366    assert (maskVal);
    13521367
    1353     KapaClearPlots (kapa3);
     1368    section.bg  = KapaColorByName ("none"); // XXX probably should be 'none'
     1369
     1370    KapaClearPlots (myKapa);
    13541371    // first section : mag vs CR nSigma
    13551372    section.dx = 1.0;
    … …  
    13591376    section.name = NULL;
    13601377    psStringAppend (&section.name, "linlog");
    1361     KapaSetSection (kapa3, &section);
     1378    KapaSetSection (myKapa, &section);
    13621379    psFree (section.name);
    13631380
    … …  
    13691386    section.name = NULL;
    13701387    psStringAppend (&section.name, "loglog");
    1371     KapaSetSection (kapa3, &section);
     1388    KapaSetSection (myKapa, &section);
    13721389    psFree (section.name);
    13731390
    … …  
    13781395        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
    13791396
    1380         psphotVisualPlotRadialProfile (kapa3, source, maskVal);
     1397        psphotVisualPlotRadialProfile (myKapa, source, maskVal);
    13811398
    13821399        // pause and wait for user input:
    … …  
    13881405        }
    13891406        if (key[0] == 'e') {
    1390             KapaClearPlots (kapa3);
     1407            KapaClearPlots (myKapa);
    13911408        }
    13921409        if (key[0] == 's') {
    … …  
    14071424    psEllipseAxes axes;
    14081425
     1426    // XXX skip this for now: it is not very clear
     1427    return true;
     1428
    14091429    if (!pmVisualIsVisual()) return true;
    14101430
    1411     if (kapa == -1) {
    1412         fprintf (stderr, "kapa not opened, skipping\n");
    1413         return false;
    1414     }
     1431    int myKapa = psphotKapaChannel (1);
     1432    if (myKapa == -1) return false;
    14151433
    14161434    // note: this uses the Ohana allocation tools:
    … …  
    14851503    }
    14861504
    1487     KiiLoadOverlay (kapa, overlayE, NoverlayE, "red");
    1488     KiiLoadOverlay (kapa, overlayO, NoverlayO, "yellow");
     1505    KiiLoadOverlay (myKapa, overlayE, NoverlayE, "red");
     1506    KiiLoadOverlay (myKapa, overlayO, NoverlayO, "yellow");
    14891507    FREE (overlayE);
    14901508    FREE (overlayO);
    14911509
    1492     // pause and wait for user input:
    1493     // continue, save (provide name), ??
    1494     char key[10];
    14951510    fprintf (stdout, "even bits (0x0001, 0x0004, ... : red\n");
    14961511    fprintf (stdout, "odd bits (0x0002, 0x0008, ... : yellow\n");
    1497     fprintf (stdout, "[c]ontinue? ");
    1498     if (!fgets(key, 8, stdin)) {
    1499         psWarning("Unable to read option");
    1500     }
    1501 
    1502     return true;
    1503 }
    1504 
    1505 bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) {
    1506 
    1507     int Noverlay, NOVERLAY;
     1512    pmVisualAskUser(NULL);
     1513
     1514    return true;
     1515}
     1516
     1517bool psphotVisualShowSourceSize_Single (int myKapa, psArray *sources, pmSourceMode mode, bool keep, float scale, char *color) {
     1518
     1519    int Noverlay;
    15081520    KiiOverlay *overlay;
    15091521
    1510     if (!pmVisualIsVisual()) return true;
    1511 
    1512     if (kapa == -1) {
    1513         fprintf (stderr, "kapa not opened, skipping\n");
    1514         return false;
    1515     }
     1522    psEllipseMoments emoments;
     1523    psEllipseAxes axes;
    15161524
    15171525    // note: this uses the Ohana allocation tools:
     1526    ALLOCATE (overlay, KiiOverlay, sources->n);
     1527
    15181528    Noverlay = 0;
    1519     NOVERLAY = 100;
    1520     ALLOCATE (overlay, KiiOverlay, sources->n);
    1521 
    1522     // mark CRs with red boxes
    15231529    for (int i = 0; i < sources->n; i++) {
    15241530
    … …  
    15261532        if (source == NULL) continue;
    15271533
    1528         if (!(source->mode & PM_SOURCE_MODE_CR_LIMIT)) continue;
    1529 
    1530         overlay[Noverlay].type = KII_OVERLAY_BOX;
    1531         overlay[Noverlay].x = source->peak->xf;
    1532         overlay[Noverlay].y = source->peak->yf;
    1533 
    1534         overlay[Noverlay].dx = 4;
    1535         overlay[Noverlay].dy = 4;
    1536         overlay[Noverlay].angle = 0;
     1534        if (mode) {
     1535            if (keep) {
     1536                if (!(source->mode & mode)) continue;
     1537            } else {
     1538                if (source->mode & mode) continue;
     1539            }
     1540        }
     1541
     1542        pmMoments *moments = source->moments;
     1543        if (moments == NULL) continue;
     1544
     1545        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     1546        overlay[Noverlay].x = moments->Mx;
     1547        overlay[Noverlay].y = moments->My;
     1548
     1549        emoments.x2 = moments->Mxx;
     1550        emoments.y2 = moments->Myy;
     1551        emoments.xy = moments->Mxy;
     1552
     1553        axes = psEllipseMomentsToAxes (emoments, 20.0);
     1554
     1555        overlay[Noverlay].dx = scale*2.0*axes.major;
     1556        overlay[Noverlay].dy = scale*2.0*axes.minor;
     1557        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
    15371558        overlay[Noverlay].text = NULL;
    15381559        Noverlay ++;
    1539         CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    1540     }
    1541     KiiLoadOverlay (kapa, overlay, Noverlay, "red");
    1542 
    1543 
    1544     Noverlay = 0;
    1545     for (int i = 0; i < sources->n; i++) {
    1546 
    1547         pmSource *source = sources->data[i];
    1548         if (source == NULL) continue;
    1549 
    1550         // mark EXTs with yellow circles
    1551         if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
    1552 
    1553         overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    1554         overlay[Noverlay].x = source->peak->xf;
    1555         overlay[Noverlay].y = source->peak->yf;
    1556 
    1557         overlay[Noverlay].dx = 10;
    1558         overlay[Noverlay].dy = 10;
    1559         overlay[Noverlay].angle = 0;
    1560         overlay[Noverlay].text = NULL;
    1561         Noverlay ++;
    1562         CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    1563     }
    1564 
    1565     KiiLoadOverlay (kapa, overlay, Noverlay, "blue");
     1560    }
     1561
     1562    KiiLoadOverlay (myKapa, overlay, Noverlay, color);
    15661563    FREE (overlay);
    15671564
    1568     psphotVisualShowMask (kapa, readout->mask, "mask", 2);
    1569 
    1570     // pause and wait for user input:
    1571     // continue, save (provide name), ??
    1572     char key[10];
    1573     fprintf (stdout, "CR: 4pix red BOX; EXT: 10pix blue circle\n");
    1574     fprintf (stdout, "[c]ontinue? ");
    1575     if (!fgets(key, 8, stdin)) {
    1576         psWarning("Unable to read option");
    1577     }
    1578 
    1579     return true;
    1580 }
    1581 
    1582 bool psphotVisualPlotSourceSize (psArray *sources) {
    1583 
     1565    return true;
     1566}
     1567
     1568bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) {
     1569
     1570    if (!pmVisualIsVisual()) return true;
     1571
     1572    int myKapa = psphotKapaChannel (1);
     1573    if (myKapa == -1) return false;
     1574
     1575    KiiEraseOverlay (myKapa, "red");
     1576    KiiEraseOverlay (myKapa, "blue");
     1577    KiiEraseOverlay (myKapa, "green");
     1578    KiiEraseOverlay (myKapa, "yellow");
     1579
     1580    psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_EXT_LIMIT | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_SATSTAR, 0, 1.0, "green");
     1581    psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_EXT_LIMIT, 1, 1.0, "blue");
     1582    psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_CR_LIMIT, 1, 1.0, "red");
     1583    psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_DEFECT, 1, 2.0, "red");
     1584    psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_SATSTAR, 1, 1.0, "yellow");
     1585
     1586    fprintf (stdout, "red: CR; blue: EXTENDED; green: PSF-like; yellow: SATSTAR\n");
     1587    pmVisualAskUser(NULL);
     1588    return true;
     1589}
     1590
     1591bool psphotVisualPlotSourceSize (psMetadata *recipe, psMetadata *analysis, psArray *sources) {
     1592
     1593    bool status;
    15841594    Graphdata graphdata;
    15851595    KapaSection section;
    … …  
    15871597    if (!pmVisualIsVisual()) return true;
    15881598
    1589     if (kapa3 == -1) {
    1590         kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
    1591         if (kapa3 == -1) {
    1592             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1593             pmVisualSetVisual(false);
    1594             return false;
    1595         }
    1596     }
    1597 
    1598     KapaClearPlots (kapa3);
     1599    int myKapa = psphotKapaChannel (2);
     1600    if (myKapa == -1) return false;
     1601
     1602    KapaClearPlots (myKapa);
    15991603    KapaInitGraph (&graphdata);
    1600 
    1601     // first section : mag vs CR nSigma
    1602     section.dx = 1.0;
    1603     section.dy = 0.5;
    1604     section.x = 0.0;
    1605     section.y = 0.0;
    1606     section.name = NULL;
    1607     psStringAppend (&section.name, "a1");
    1608     KapaSetSection (kapa3, &section);
     1604    KapaSetFont (myKapa, "courier", 14);
     1605
     1606    section.bg  = KapaColorByName ("none"); // XXX probably should be 'none'
     1607
     1608    // select the max psfX,Y values for the plot limits
     1609    float Xmin = 1000.0, Xmax = 0.0;
     1610    float Ymin = 1000.0, Ymax = 0.0;
     1611    {
     1612        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
     1613        for (int n = 0; n < nRegions; n++) {
     1614
     1615            char regionName[64];
     1616            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
     1617            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
     1618
     1619            float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
     1620            float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
     1621            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
     1622            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
     1623
     1624            float X0 = psfX - 10.0*psfdX;
     1625            float X1 = psfX + 10.0*psfdX;
     1626            float Y0 = psfY - 10.0*psfdY;
     1627            float Y1 = psfY + 10.0*psfdY;
     1628
     1629            if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); }
     1630            if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
     1631            if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); }
     1632            if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
     1633        }
     1634    }
     1635    Xmin = PS_MAX(Xmin, -0.1);
     1636    Ymin = PS_MAX(Ymin, -0.1);
     1637
     1638    // storage vectors for data to be plotted
     1639    psVector *xSAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1640    psVector *ySAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1641    psVector *mSAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1642    psVector *sSAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1643
     1644    psVector *xPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1645    psVector *yPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1646    psVector *mPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1647    psVector *sPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1648
     1649    psVector *xEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1650    psVector *yEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1651    psVector *mEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1652    psVector *sEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1653
     1654    psVector *xDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1655    psVector *yDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1656    psVector *mDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1657    psVector *sDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1658
     1659    psVector *xLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1660    psVector *yLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1661    psVector *mLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1662    psVector *sLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1663
     1664    psVector *xCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1665    psVector *yCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1666    psVector *mCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1667    psVector *sCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1668
     1669    // construct the vectors
     1670    int nSAT = 0;
     1671    int nEXT = 0;
     1672    int nPSF = 0;
     1673    int nDEF = 0;
     1674    int nLOW = 0;
     1675    int nCR  = 0;
     1676    for (int i = 0; i < sources->n; i++) {
     1677        pmSource *source = sources->data[i];
     1678        if (source->moments == NULL) continue;
     1679
     1680        // only plot the measured sources...
     1681        if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) continue;
     1682
     1683        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) {
     1684            xCR->data.F32[nCR] = source->moments->Mxx;
     1685            yCR->data.F32[nCR] = source->moments->Myy;
     1686            mCR->data.F32[nCR] = -2.5*log10(source->moments->Sum);
     1687            sCR->data.F32[nCR] = source->extNsigma;
     1688            nCR++;
     1689        }
     1690        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
     1691            xSAT->data.F32[nSAT] = source->moments->Mxx;
     1692            ySAT->data.F32[nSAT] = source->moments->Myy;
     1693            mSAT->data.F32[nSAT] = -2.5*log10(source->moments->Sum);
     1694            sSAT->data.F32[nSAT] = source->extNsigma;
     1695            nSAT++;
     1696        }
     1697        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     1698            xEXT->data.F32[nEXT] = source->moments->Mxx;
     1699            yEXT->data.F32[nEXT] = source->moments->Myy;
     1700            mEXT->data.F32[nEXT] = -2.5*log10(source->moments->Sum);
     1701            sEXT->data.F32[nEXT] = source->extNsigma;
     1702            nEXT++;
     1703            continue;
     1704        }
     1705        if (source->mode & PM_SOURCE_MODE_DEFECT) {
     1706            xDEF->data.F32[nDEF] = source->moments->Mxx;
     1707            yDEF->data.F32[nDEF] = source->moments->Myy;
     1708            mDEF->data.F32[nDEF] = -2.5*log10(source->moments->Sum);
     1709            sDEF->data.F32[nDEF] = source->extNsigma;
     1710            nDEF++;
     1711            continue;
     1712        }
     1713        if (source->errMag > 0.1) {
     1714            xLOW->data.F32[nLOW] = source->moments->Mxx;
     1715            yLOW->data.F32[nLOW] = source->moments->Myy;
     1716            mLOW->data.F32[nLOW] = -2.5*log10(source->moments->Sum);
     1717            sLOW->data.F32[nLOW] = source->extNsigma;
     1718            nLOW++;
     1719            continue;
     1720        }
     1721        xPSF->data.F32[nPSF] = source->moments->Mxx;
     1722        yPSF->data.F32[nPSF] = source->moments->Myy;
     1723        mPSF->data.F32[nPSF] = -2.5*log10(source->moments->Sum);
     1724        sPSF->data.F32[nPSF] = source->extNsigma;
     1725        nPSF++;
     1726    }
     1727
     1728    xSAT->n = nSAT;
     1729    ySAT->n = nSAT;
     1730    mSAT->n = nSAT;
     1731    sSAT->n = nSAT;
     1732
     1733    xPSF->n = nPSF;
     1734    yPSF->n = nPSF;
     1735    mPSF->n = nPSF;
     1736    sPSF->n = nPSF;
     1737
     1738    xEXT->n = nEXT;
     1739    yEXT->n = nEXT;
     1740    mEXT->n = nEXT;
     1741    sEXT->n = nEXT;
     1742
     1743    xCR->n = nCR;
     1744    yCR->n = nCR;
     1745    mCR->n = nCR;
     1746    sCR->n = nCR;
     1747
     1748    xDEF->n = nDEF;
     1749    yDEF->n = nDEF;
     1750    mDEF->n = nDEF;
     1751    sDEF->n = nDEF;
     1752
     1753    xLOW->n = nLOW;
     1754    yLOW->n = nLOW;
     1755    mLOW->n = nLOW;
     1756    sLOW->n = nLOW;
     1757
     1758    // four sections: MxxMyy, MagMxx, MagMyy, MagSigma
     1759
     1760    // first section: MxxMyy
     1761    section.dx = 0.75;
     1762    section.dy = 0.60;
     1763    section.x  = 0.00;
     1764    section.y  = 0.00;
     1765    section.name = psStringCopy ("MxxMyy");
     1766    KapaSetSection (myKapa, &section);
    16091767    psFree (section.name);
     1768
     1769    graphdata.color = KapaColorByName ("black");
     1770    graphdata.xmin = Xmin;
     1771    graphdata.ymin = Ymin;
     1772    graphdata.xmax = Xmax;
     1773    graphdata.ymax = Ymax;
     1774    KapaSetLimits (myKapa, &graphdata);
     1775
     1776    graphdata.padXm = NAN;
     1777    graphdata.padYm = NAN;
     1778    graphdata.padXp = 0.5;
     1779    graphdata.padYp = 0.5;
     1780    KapaBox (myKapa, &graphdata);
     1781    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
     1782    KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM);
     1783
     1784    graphdata.color = KapaColorByName ("black");
     1785    graphdata.ptype = 0;
     1786    graphdata.size = 0.5;
     1787    graphdata.style = 2;
     1788    KapaPrepPlot   (myKapa, nPSF, &graphdata);
     1789    KapaPlotVector (myKapa, nPSF, xPSF->data.F32, "x");
     1790    KapaPlotVector (myKapa, nPSF, yPSF->data.F32, "y");
     1791
     1792    graphdata.color = KapaColorByName ("blue");
     1793    graphdata.ptype = 0;
     1794    graphdata.size = 0.5;
     1795    graphdata.style = 2;
     1796    KapaPrepPlot   (myKapa, nEXT, &graphdata);
     1797    KapaPlotVector (myKapa, nEXT, xEXT->data.F32, "x");
     1798    KapaPlotVector (myKapa, nEXT, yEXT->data.F32, "y");
     1799
     1800    graphdata.color = KapaColorByName ("red");
     1801    graphdata.ptype = 0;
     1802    graphdata.size = 0.5;
     1803    graphdata.style = 2;
     1804    KapaPrepPlot   (myKapa, nDEF, &graphdata);
     1805    KapaPlotVector (myKapa, nDEF, xDEF->data.F32, "x");
     1806    KapaPlotVector (myKapa, nDEF, yDEF->data.F32, "y");
     1807
     1808    graphdata.color = KapaColorByName ("red");
     1809    graphdata.ptype = 7;
     1810    graphdata.size = 1.0;
     1811    graphdata.style = 2;
     1812    KapaPrepPlot   (myKapa, nCR, &graphdata);
     1813    KapaPlotVector (myKapa, nCR, xCR->data.F32, "x");
     1814    KapaPlotVector (myKapa, nCR, yCR->data.F32, "y");
     1815
     1816    graphdata.color = KapaColorByName ("blue");
     1817    graphdata.ptype = 7;
     1818    graphdata.size = 1.0;
     1819    graphdata.style = 2;
     1820    KapaPrepPlot   (myKapa, nSAT, &graphdata);
     1821    KapaPlotVector (myKapa, nSAT, xSAT->data.F32, "x");
     1822    KapaPlotVector (myKapa, nSAT, ySAT->data.F32, "y");
     1823
     1824    graphdata.color = KapaColorByName ("black");
     1825    graphdata.ptype = 7;
     1826    graphdata.size = 1.0;
     1827    graphdata.style = 2;
     1828    KapaPrepPlot   (myKapa, nLOW, &graphdata);
     1829    KapaPlotVector (myKapa, nLOW, xLOW->data.F32, "x");
     1830    KapaPlotVector (myKapa, nLOW, yLOW->data.F32, "y");
     1831
     1832    // second section: MagMyy
     1833    section.dx = 0.75;
     1834    section.dy = 0.20;
     1835    section.x  = 0.00;
     1836    section.y  = 0.80;
     1837    section.name = psStringCopy ("MagMyy");
     1838    KapaSetSection (myKapa, &section);
     1839    psFree (section.name);
     1840
     1841    graphdata.color = KapaColorByName ("black");
     1842    graphdata.xmin = -17.1;
     1843    graphdata.xmax =  -6.9;
     1844    graphdata.ymin = Ymin;
     1845    graphdata.ymax = Ymax;
     1846    KapaSetLimits (myKapa, &graphdata);
     1847
     1848    graphdata.padXm = 0.5;
     1849    graphdata.padYm = NAN;
     1850    graphdata.padXp = NAN;
     1851    graphdata.padYp = 0.5;
     1852    strcpy (graphdata.labels, "0210");
     1853    KapaBox (myKapa, &graphdata);
     1854    KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_XP);
     1855    KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM);
     1856
     1857    graphdata.color = KapaColorByName ("black");
     1858    graphdata.ptype = 0;
     1859    graphdata.size = 0.5;
     1860    graphdata.style = 2;
     1861    KapaPrepPlot   (myKapa, nPSF, &graphdata);
     1862    KapaPlotVector (myKapa, nPSF, mPSF->data.F32, "x");
     1863    KapaPlotVector (myKapa, nPSF, yPSF->data.F32, "y");
     1864
     1865    graphdata.color = KapaColorByName ("blue");
     1866    graphdata.ptype = 0;
     1867    graphdata.size = 0.5;
     1868    graphdata.style = 2;
     1869    KapaPrepPlot   (myKapa, nEXT, &graphdata);
     1870    KapaPlotVector (myKapa, nEXT, mEXT->data.F32, "x");
     1871    KapaPlotVector (myKapa, nEXT, yEXT->data.F32, "y");
     1872
     1873    graphdata.color = KapaColorByName ("red");
     1874    graphdata.ptype = 0;
     1875    graphdata.size = 0.5;
     1876    graphdata.style = 2;
     1877    KapaPrepPlot   (myKapa, nDEF, &graphdata);
     1878    KapaPlotVector (myKapa, nDEF, mDEF->data.F32, "x");
     1879    KapaPlotVector (myKapa, nDEF, yDEF->data.F32, "y");
     1880
     1881    graphdata.color = KapaColorByName ("red");
     1882    graphdata.ptype = 7;
     1883    graphdata.size = 1.0;
     1884    graphdata.style = 2;
     1885    KapaPrepPlot   (myKapa, nCR, &graphdata);
     1886    KapaPlotVector (myKapa, nCR, mCR->data.F32, "x");
     1887    KapaPlotVector (myKapa, nCR, yCR->data.F32, "y");
     1888
     1889    graphdata.color = KapaColorByName ("blue");
     1890    graphdata.ptype = 7;
     1891    graphdata.size = 1.0;
     1892    graphdata.style = 2;
     1893    KapaPrepPlot   (myKapa, nSAT, &graphdata);
     1894    KapaPlotVector (myKapa, nSAT, mSAT->data.F32, "x");
     1895    KapaPlotVector (myKapa, nSAT, ySAT->data.F32, "y");
     1896
     1897    // third section: MagMxx
     1898    section.dx = 0.25;
     1899    section.dy = 0.60;
     1900    section.x  = 0.75;
     1901    section.y  = 0.00;
     1902    section.name = psStringCopy ("MagMxx");
     1903    KapaSetSection (myKapa, &section);
     1904    psFree (section.name);
     1905
     1906    graphdata.color = KapaColorByName ("black");
     1907    graphdata.xmin = Xmin;
     1908    graphdata.xmax = Xmax;
     1909    graphdata.ymin =  -6.9;
     1910    graphdata.ymax = -17.1;
     1911    KapaSetLimits (myKapa, &graphdata);
     1912
     1913    graphdata.padXm = NAN;
     1914    graphdata.padYm = 0.5;
     1915    graphdata.padXp = 0.5;
     1916    graphdata.padYp = NAN;
     1917    strcpy (graphdata.labels, "2001");
     1918    KapaBox (myKapa, &graphdata);
     1919    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
     1920    KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_YP);
     1921
     1922    graphdata.color = KapaColorByName ("black");
     1923    graphdata.ptype = 0;
     1924    graphdata.size = 0.5;
     1925    graphdata.style = 2;
     1926    KapaPrepPlot   (myKapa, nPSF, &graphdata);
     1927    KapaPlotVector (myKapa, nPSF, xPSF->data.F32, "x");
     1928    KapaPlotVector (myKapa, nPSF, mPSF->data.F32, "y");
     1929
     1930    graphdata.color = KapaColorByName ("blue");
     1931    graphdata.ptype = 0;
     1932    graphdata.size = 0.5;
     1933    graphdata.style = 2;
     1934    KapaPrepPlot   (myKapa, nEXT, &graphdata);
     1935    KapaPlotVector (myKapa, nEXT, xEXT->data.F32, "x");
     1936    KapaPlotVector (myKapa, nEXT, mEXT->data.F32, "y");
     1937
     1938    graphdata.color = KapaColorByName ("red");
     1939    graphdata.ptype = 0;
     1940    graphdata.size = 0.5;
     1941    graphdata.style = 2;
     1942    KapaPrepPlot   (myKapa, nDEF, &graphdata);
     1943    KapaPlotVector (myKapa, nDEF, xDEF->data.F32, "x");
     1944    KapaPlotVector (myKapa, nDEF, mDEF->data.F32, "y");
     1945
     1946    graphdata.color = KapaColorByName ("red");
     1947    graphdata.ptype = 7;
     1948    graphdata.size = 1.0;
     1949    graphdata.style = 2;
     1950    KapaPrepPlot   (myKapa, nCR, &graphdata);
     1951    KapaPlotVector (myKapa, nCR, xCR->data.F32, "x");
     1952    KapaPlotVector (myKapa, nCR, mCR->data.F32, "y");
     1953
     1954    graphdata.color = KapaColorByName ("blue");
     1955    graphdata.ptype = 7;
     1956    graphdata.size = 1.0;
     1957    graphdata.style = 2;
     1958    KapaPrepPlot   (myKapa, nSAT, &graphdata);
     1959    KapaPlotVector (myKapa, nSAT, xSAT->data.F32, "x");
     1960    KapaPlotVector (myKapa, nSAT, mSAT->data.F32, "y");
     1961
     1962    // fourth section: MagSigma
     1963    section.dx = 0.75;
     1964    section.dy = 0.20;
     1965    section.x  = 0.00;
     1966    section.y  = 0.60;
     1967    section.name = psStringCopy ("MagSigma");
     1968    KapaSetSection (myKapa, &section);
     1969    psFree (section.name);
     1970
     1971    graphdata.color = KapaColorByName ("black");
     1972    graphdata.xmax =  -6.9;
     1973    graphdata.xmin = -17.1;
     1974    graphdata.ymin = -20.1;
     1975    graphdata.ymax = +20.1;
     1976    KapaSetLimits (myKapa, &graphdata);
     1977
     1978    graphdata.padXm = 0.5;
     1979    graphdata.padYm = NAN;
     1980    graphdata.padXp = 0.5;
     1981    graphdata.padYp = 0.5;
     1982    strcpy (graphdata.labels, "0100");
     1983    KapaBox (myKapa, &graphdata);
     1984    // KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_XM);
     1985    KapaSendLabel (myKapa, "EXT&ss&c", KAPA_LABEL_YM);
     1986
     1987    graphdata.color = KapaColorByName ("black");
     1988    graphdata.ptype = 0;
     1989    graphdata.size = 0.5;
     1990    graphdata.style = 2;
     1991    KapaPrepPlot   (myKapa, nPSF, &graphdata);
     1992    KapaPlotVector (myKapa, nPSF, mPSF->data.F32, "x");
     1993    KapaPlotVector (myKapa, nPSF, sPSF->data.F32, "y");
     1994
     1995    graphdata.color = KapaColorByName ("blue");
     1996    graphdata.ptype = 0;
     1997    graphdata.size = 0.5;
     1998    graphdata.style = 2;
     1999    KapaPrepPlot   (myKapa, nEXT, &graphdata);
     2000    KapaPlotVector (myKapa, nEXT, mEXT->data.F32, "x");
     2001    KapaPlotVector (myKapa, nEXT, sEXT->data.F32, "y");
     2002
     2003    graphdata.color = KapaColorByName ("red");
     2004    graphdata.ptype = 0;
     2005    graphdata.size = 0.5;
     2006    graphdata.style = 2;
     2007    KapaPrepPlot   (myKapa, nDEF, &graphdata);
     2008    KapaPlotVector (myKapa, nDEF, mDEF->data.F32, "x");
     2009    KapaPlotVector (myKapa, nDEF, sDEF->data.F32, "y");
     2010
     2011    graphdata.color = KapaColorByName ("red");
     2012    graphdata.ptype = 7;
     2013    graphdata.size = 1.0;
     2014    graphdata.style = 2;
     2015    KapaPrepPlot   (myKapa, nCR, &graphdata);
     2016    KapaPlotVector (myKapa, nCR, mCR->data.F32, "x");
     2017    KapaPlotVector (myKapa, nCR, sCR->data.F32, "y");
     2018
     2019    graphdata.color = KapaColorByName ("blue");
     2020    graphdata.ptype = 7;
     2021    graphdata.size = 1.0;
     2022    graphdata.style = 2;
     2023    KapaPrepPlot   (myKapa, nSAT, &graphdata);
     2024    KapaPlotVector (myKapa, nSAT, mSAT->data.F32, "x");
     2025    KapaPlotVector (myKapa, nSAT, sSAT->data.F32, "y");
     2026
     2027    // draw N circles to outline the clumps
     2028    {
     2029        KapaSelectSection (myKapa, "MxxMyy");
     2030
     2031        // draw a circle centered on psfX,Y with size of the psf limit
     2032        psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
     2033        psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
     2034
     2035        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
     2036        float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
     2037
     2038        graphdata.color = KapaColorByName ("blue");
     2039        graphdata.style = 0;
     2040
     2041        graphdata.xmin = Xmin;
     2042        graphdata.ymin = Ymin;
     2043        graphdata.xmax = Xmax;
     2044        graphdata.ymax = Ymax;
     2045        KapaSetLimits (myKapa, &graphdata);
     2046
     2047        for (int n = 0; n < nRegions; n++) {
     2048
     2049            char regionName[64];
     2050            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
     2051            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
     2052
     2053            float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
     2054            float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
     2055            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
     2056            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
     2057            float Rx = psfdX * PSF_CLUMP_NSIGMA;
     2058            float Ry = psfdY * PSF_CLUMP_NSIGMA;
     2059
     2060            for (int i = 0; i < xLimit->n; i++) {
     2061                xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX;
     2062                yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
     2063            }
     2064            KapaPrepPlot (myKapa, xLimit->n, &graphdata);
     2065            KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x");
     2066            KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y");
     2067        }
     2068        psFree (xLimit);
     2069        psFree (yLimit);
     2070    }
     2071
     2072    psFree (xSAT);
     2073    psFree (ySAT);
     2074    psFree (mSAT);
     2075    psFree (sSAT);
     2076
     2077    psFree (xEXT);
     2078    psFree (yEXT);
     2079    psFree (mEXT);
     2080    psFree (sEXT);
     2081
     2082    psFree (xPSF);
     2083    psFree (yPSF);
     2084    psFree (mPSF);
     2085    psFree (sPSF);
     2086
     2087    psFree (xDEF);
     2088    psFree (yDEF);
     2089    psFree (mDEF);
     2090    psFree (sDEF);
     2091
     2092    psFree (xLOW);
     2093    psFree (yLOW);
     2094    psFree (mLOW);
     2095    psFree (sLOW);
     2096
     2097    psFree (xCR);
     2098    psFree (yCR);
     2099    psFree (mCR);
     2100    psFree (sCR);
     2101
     2102    pmVisualAskUser(NULL);
     2103    return true;
     2104}
     2105
     2106bool psphotVisualShowResidualImage (pmReadout *readout) {
     2107
     2108    if (!pmVisualIsVisual()) return true;
     2109
     2110    int myKapa = psphotKapaChannel (1);
     2111    if (myKapa == -1) return false;
     2112
     2113    psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", 1);
     2114
     2115    pmVisualAskUser(NULL);
     2116    return true;
     2117}
     2118
     2119bool psphotVisualPlotApResid (psArray *sources, float mean, float error) {
     2120
     2121    Graphdata graphdata;
     2122    float lineX[2], lineY[2];
     2123
     2124    if (!pmVisualIsVisual()) return true;
     2125
     2126    int myKapa = psphotKapaChannel (2);
     2127    if (myKapa == -1) return false;
     2128
     2129    KapaClearPlots (myKapa);
     2130    KapaInitGraph (&graphdata);
    16102131
    16112132    psVector *x = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    16122133    psVector *y = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    1613 
    1614     graphdata.xmin = +32.0;
    1615     graphdata.xmax = -32.0;
    1616     graphdata.ymin = +32.0;
    1617     graphdata.ymax = -32.0;
    1618 
    1619     // construct the plot vectors
    1620     int n = 0;
    1621     for (int i = 0; i < sources->n; i++) {
    1622         pmSource *source = sources->data[i];
    1623         if (!source) continue;
    1624         if (source->type != PM_SOURCE_TYPE_STAR) continue;
    1625         if (!isfinite (source->crNsigma)) continue;
    1626 
    1627         x->data.F32[n] = -2.5*log10(source->peak->flux);
    1628         y->data.F32[n] = source->crNsigma;
    1629         graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
    1630         graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
    1631         graphdata.ymin = -0.5;
    1632         graphdata.ymax = 10.0;
    1633 
    1634         n++;
    1635     }
    1636     x->n = y->n = n;
    1637 
    1638     float range;
    1639     range = graphdata.xmax - graphdata.xmin;
    1640     graphdata.xmax += 0.05*range;
    1641     graphdata.xmin -= 0.05*range;
    1642 
    1643     // XXX set the plot range to match the image
    1644     KapaSetLimits (kapa3, &graphdata);
    1645 
    1646     KapaSetFont (kapa3, "helvetica", 14);
    1647     KapaBox (kapa3, &graphdata);
    1648     KapaSendLabel (kapa3, "Peak as Mag", KAPA_LABEL_XM);
    1649     KapaSendLabel (kapa3, "CR N Sigma", KAPA_LABEL_YM);
    1650 
    1651     graphdata.color = KapaColorByName ("black");
    1652     graphdata.ptype = 2;
    1653     graphdata.size = 0.5;
    1654     graphdata.style = 2;
    1655     KapaPrepPlot (kapa3, n, &graphdata);
    1656     KapaPlotVector (kapa3, n, x->data.F32, "x");
    1657     KapaPlotVector (kapa3, n, y->data.F32, "y");
    1658 
    1659     // second section : mag vs EXT nSigma
    1660     section.dx = 1.0;
    1661     section.dy = 0.5;
    1662     section.x = 0.0;
    1663     section.y = 0.5;
    1664     section.name = NULL;
    1665     psStringAppend (&section.name, "a2");
    1666     KapaSetSection (kapa3, &section);
    1667     psFree (section.name);
    1668 
    1669     graphdata.xmin = +32.0;
    1670     graphdata.xmax = -32.0;
    1671     graphdata.ymin = +32.0;
    1672     graphdata.ymax = -32.0;
    1673 
    1674     // construct the plot vectors
    1675     n = 0;
    1676     for (int i = 0; i < sources->n; i++) {
    1677         pmSource *source = sources->data[i];
    1678         if (!source) continue;
    1679         if (source->type != PM_SOURCE_TYPE_STAR) continue;
    1680         if (!isfinite (source->extNsigma)) continue;
    1681 
    1682         x->data.F32[n] = -2.5*log10(source->peak->flux);
    1683         y->data.F32[n] = source->extNsigma;
    1684         graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
    1685         graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
    1686         graphdata.ymin = -0.5;
    1687         graphdata.ymax = 10.0;
    1688 
    1689         n++;
    1690     }
    1691     x->n = y->n = n;
    1692 
    1693     range = graphdata.xmax - graphdata.xmin;
    1694     graphdata.xmax += 0.05*range;
    1695     graphdata.xmin -= 0.05*range;
    1696 
    1697     // XXX set the plot range to match the image
    1698     KapaSetLimits (kapa3, &graphdata);
    1699 
    1700     KapaSetFont (kapa3, "helvetica", 14);
    1701     KapaBox (kapa3, &graphdata);
    1702     KapaSendLabel (kapa3, "EXT N Sigma", KAPA_LABEL_YM);
    1703 
    1704     graphdata.color = KapaColorByName ("black");
    1705     graphdata.ptype = 2;
    1706     graphdata.size = 0.5;
    1707     graphdata.style = 2;
    1708     KapaPrepPlot (kapa3, n, &graphdata);
    1709     KapaPlotVector (kapa3, n, x->data.F32, "x");
    1710     KapaPlotVector (kapa3, n, y->data.F32, "y");
    1711 
    1712     psFree (x);
    1713     psFree (y);
    1714 
    1715     // pause and wait for user input:
    1716     // continue, save (provide name), ??
    1717     char key[10];
    1718     fprintf (stdout, "[c]ontinue? ");
    1719     if (!fgets(key, 8, stdin)) {
    1720         psWarning("Unable to read option");
    1721     }
    1722     return true;
    1723 }
    1724 
    1725 bool psphotVisualShowResidualImage (pmReadout *readout) {
    1726 
    1727     if (!pmVisualIsVisual()) return true;
    1728 
    1729     if (kapa == -1) {
    1730         kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
    1731         if (kapa == -1) {
    1732             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1733             pmVisualSetVisual(false);
    1734             return false;
    1735         }
    1736     }
    1737 
    1738     psphotVisualScaleImage (kapa, readout->image, "resid", 1);
    1739 
    1740     // pause and wait for user input:
    1741     // continue, save (provide name), ??
    1742     char key[10];
    1743     fprintf (stdout, "[c]ontinue? ");
    1744     if (!fgets(key, 8, stdin)) {
    1745         psWarning("Unable to read option");
    1746     }
    1747     return true;
    1748 }
    1749 
    1750 bool psphotVisualPlotApResid (psArray *sources) {
    1751 
    1752     Graphdata graphdata;
    1753 
    1754     if (!pmVisualIsVisual()) return true;
    1755 
    1756     if (kapa3 == -1) {
    1757         kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
    1758         if (kapa3 == -1) {
    1759             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1760             pmVisualSetVisual(false);
    1761             return false;
    1762         }
    1763     }
    1764 
    1765     KapaClearPlots (kapa3);
    1766     KapaInitGraph (&graphdata);
    1767 
    1768     psVector *x = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    1769     psVector *y = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2134    psVector *dy = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    17702135
    17712136    graphdata.xmin = +32.0;
    … …  
    17852150        x->data.F32[n] = source->psfMag;
    17862151        y->data.F32[n] = source->apMag - source->psfMag;
     2152        dy->data.F32[n] = source->errMag;
    17872153        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
    17882154        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
    … …  
    17922158        n++;
    17932159    }
    1794     x->n = y->n = n;
     2160    x->n = y->n = dy->n = n;
    17952161
    17962162    float range;
    … …  
    18022168    graphdata.ymin -= 0.05*range;
    18032169
    1804     // XXX set the plot range to match the image
    1805     KapaSetLimits (kapa3, &graphdata);
    1806 
    1807     KapaSetFont (kapa3, "helvetica", 14);
    1808     KapaBox (kapa3, &graphdata);
    1809     KapaSendLabel (kapa3, "PSF Mag", KAPA_LABEL_XM);
    1810     KapaSendLabel (kapa3, "Ap Mag - PSF Mag", KAPA_LABEL_YM);
     2170    // XXX test
     2171    graphdata.xmin = -17.0;
     2172    graphdata.xmax =  -9.0;
     2173    graphdata.ymin = -0.31;
     2174    graphdata.ymax = +0.31;
     2175
     2176    KapaSetLimits (myKapa, &graphdata);
     2177
     2178    KapaSetFont (myKapa, "helvetica", 14);
     2179    KapaBox (myKapa, &graphdata);
     2180    KapaSendLabel (myKapa, "PSF Mag", KAPA_LABEL_XM);
     2181    KapaSendLabel (myKapa, "Ap Mag - PSF Mag", KAPA_LABEL_YM);
    18112182
    18122183    graphdata.color = KapaColorByName ("black");
    … …  
    18142185    graphdata.size = 0.5;
    18152186    graphdata.style = 2;
    1816     KapaPrepPlot (kapa3, n, &graphdata);
    1817     KapaPlotVector (kapa3, n, x->data.F32, "x");
    1818     KapaPlotVector (kapa3, n, y->data.F32, "y");
     2187    graphdata.etype |= 0x01;
     2188    KapaPrepPlot (myKapa, n, &graphdata);
     2189    KapaPlotVector (myKapa, n, x->data.F32, "x");
     2190    KapaPlotVector (myKapa, n, y->data.F32, "y");
     2191    KapaPlotVector (myKapa, n, dy->data.F32, "dym");
     2192    KapaPlotVector (myKapa, n, dy->data.F32, "dyp");
     2193
     2194    graphdata.color = KapaColorByName ("blue");
     2195    graphdata.ptype = 0;
     2196    graphdata.size = 0.5;
     2197    graphdata.style = 0;
     2198    graphdata.etype = 0;
     2199    lineX[0] = graphdata.xmin;
     2200    lineX[1] = graphdata.xmax;
     2201    lineY[0] = lineY[1] = mean;
     2202    KapaPrepPlot (myKapa, 2, &graphdata);
     2203    KapaPlotVector (myKapa, 2, lineX, "x");
     2204    KapaPlotVector (myKapa, 2, lineY, "y");
     2205
     2206    graphdata.color = KapaColorByName ("red");
     2207    graphdata.ptype = 0;
     2208    graphdata.size = 0.5;
     2209    graphdata.style = 0;
     2210    graphdata.etype = 0;
     2211    lineX[0] = graphdata.xmin;
     2212    lineX[1] = graphdata.xmax;
     2213    lineY[0] = lineY[1] = mean + error;
     2214    KapaPrepPlot (myKapa, 2, &graphdata);
     2215    KapaPlotVector (myKapa, 2, lineX, "x");
     2216    KapaPlotVector (myKapa, 2, lineY, "y");
     2217
     2218    graphdata.color = KapaColorByName ("red");
     2219    graphdata.ptype = 0;
     2220    graphdata.size = 0.5;
     2221    graphdata.style = 0;
     2222    graphdata.etype = 0;
     2223    lineX[0] = graphdata.xmin;
     2224    lineX[1] = graphdata.xmax;
     2225    lineY[0] = lineY[1] = mean - error;
     2226    KapaPrepPlot (myKapa, 2, &graphdata);
     2227    KapaPlotVector (myKapa, 2, lineX, "x");
     2228    KapaPlotVector (myKapa, 2, lineY, "y");
    18192229
    18202230    psFree (x);
    18212231    psFree (y);
    1822 
    1823     // pause and wait for user input:
    1824     // continue, save (provide name), ??
    1825     char key[10];
    1826     fprintf (stdout, "[c]ontinue? ");
    1827     if (!fgets(key, 8, stdin)) {
    1828         psWarning("Unable to read option");
    1829     }
     2232    psFree (dy);
     2233
     2234    pmVisualAskUser(NULL);
     2235    return true;
     2236}
     2237
     2238bool psphotVisualPlotChisq (psArray *sources) {
     2239
     2240    Graphdata graphdata;
     2241
     2242    if (!pmVisualIsVisual()) return true;
     2243
     2244    int myKapa = psphotKapaChannel (2);
     2245    if (myKapa == -1) return false;
     2246
     2247    KapaClearPlots (myKapa);
     2248    KapaInitGraph (&graphdata);
     2249
     2250    psVector *x = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2251    psVector *y = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2252
     2253    graphdata.xmin = +32.0;
     2254    graphdata.xmax = -32.0;
     2255    graphdata.ymin = +32.0;
     2256    graphdata.ymax = -32.0;
     2257
     2258    FILE *f = fopen ("chisq.dat", "w");
     2259
     2260    // construct the plot vectors
     2261    int n = 0;
     2262    for (int i = 0; i < sources->n; i++) {
     2263        pmSource *source = sources->data[i];
     2264        if (!source) continue;
     2265        if (source->type != PM_SOURCE_TYPE_STAR) continue;
     2266        if (!source->moments) continue;
     2267        if (!isfinite(source->moments->Sum)) continue;
     2268        if (!source->modelPSF) continue;
     2269        if (!isfinite(source->modelPSF->chisq)) continue;
     2270
     2271        x->data.F32[n] = -2.5*log10(source->moments->Sum);
     2272        y->data.F32[n] = source->modelPSF->chisq / source->modelPSF->nDOF;
     2273        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
     2274        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
     2275        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]);
     2276        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]);
     2277
     2278        fprintf (f, "%d %d %f %f\n", i, n, x->data.F32[n], y->data.F32[n]);
     2279
     2280        n++;
     2281    }
     2282    x->n = y->n = n;
     2283    fclose (f);
     2284
     2285    float range;
     2286    range = graphdata.xmax - graphdata.xmin;
     2287    graphdata.xmax += 0.05*range;
     2288    graphdata.xmin -= 0.05*range;
     2289    range = graphdata.ymax - graphdata.ymin;
     2290    graphdata.ymax += 0.05*range;
     2291    graphdata.ymin -= 0.05*range;
     2292
     2293    // XXX test
     2294    graphdata.xmin = -17.0;
     2295    graphdata.xmax =  -3.0;
     2296    graphdata.ymin =  -0.1;
     2297    graphdata.ymax = +10.1;
     2298
     2299    KapaSetLimits (myKapa, &graphdata);
     2300
     2301    KapaSetFont (myKapa, "helvetica", 14);
     2302    KapaBox (myKapa, &graphdata);
     2303    KapaSendLabel (myKapa, "PSF Mag", KAPA_LABEL_XM);
     2304    KapaSendLabel (myKapa, "ChiSq", KAPA_LABEL_YM);
     2305
     2306    graphdata.color = KapaColorByName ("black");
     2307    graphdata.ptype = 2;
     2308    graphdata.size = 0.5;
     2309    graphdata.style = 2;
     2310    KapaPrepPlot (myKapa, n, &graphdata);
     2311    KapaPlotVector (myKapa, n, x->data.F32, "x");
     2312    KapaPlotVector (myKapa, n, y->data.F32, "y");
     2313
     2314    psFree (x);
     2315    psFree (y);
     2316
     2317    pmVisualAskUser(NULL);
     2318    return true;
     2319}
     2320
     2321bool psphotVisualShowPetrosians (psArray *sources) {
     2322
     2323    int Noverlay, NOVERLAY;
     2324    KiiOverlay *overlay;
     2325
     2326    if (!pmVisualIsVisual()) return true;
     2327
     2328    int kapa = psphotKapaChannel (1);
     2329    if (kapa == -1) return false;
     2330
     2331    Noverlay = 0;
     2332    NOVERLAY = 100;
     2333    ALLOCATE (overlay, KiiOverlay, NOVERLAY);
     2334
     2335    for (int i = 0; i < sources->n; i++) {
     2336        pmSource *source = sources->data[i];
     2337
     2338        if (!source) continue;
     2339        if (!source->extpars) continue;
     2340        if (!source->extpars->petProfile) continue;
     2341
     2342        float petrosianRadius = source->extpars->petrosianRadius;
     2343        psEllipseAxes *axes = &source->extpars->axes;
     2344
     2345        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     2346        overlay[Noverlay].x = source->peak->xf;
     2347        overlay[Noverlay].y = source->peak->yf;
     2348        overlay[Noverlay].dx = 1.0*petrosianRadius;
     2349        overlay[Noverlay].dy = 1.0*petrosianRadius*axes->minor/axes->major;
     2350        overlay[Noverlay].angle = axes->theta * PS_DEG_RAD;
     2351        overlay[Noverlay].text = NULL;
     2352        Noverlay ++;
     2353        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     2354
     2355        // overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     2356        // overlay[Noverlay].x = source->peak->xf;
     2357        // overlay[Noverlay].y = source->peak->yf;
     2358        // overlay[Noverlay].dx = 2.0*petrosianRadius;
     2359        // overlay[Noverlay].dy = 2.0*petrosianRadius*axes->minor/axes->major;
     2360        // overlay[Noverlay].angle = axes->theta * PS_DEG_RAD;
     2361        // overlay[Noverlay].text = NULL;
     2362        // Noverlay ++;
     2363        // CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     2364    }
     2365
     2366    KiiLoadOverlay (kapa, overlay, Noverlay, "red");
     2367    FREE (overlay);
     2368
     2369    pmVisualAskUser(NULL);
    18302370    return true;
    18312371}
    … …  
    18502390
    18512391# endif
     2392
     2393# if (0)
     2394// *** make a histogram of the source counts in the x and y directions
     2395psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0);
     2396psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0);
     2397psVectorHistogram (nX, xFaint, NULL, NULL, 0);
     2398psVectorHistogram (nY, yFaint, NULL, NULL, 0);
     2399psVector *dX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
     2400psVector *vX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
     2401psVector *dY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
     2402psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
     2403for (int i = 0; i < nX->nums->n; i++) {
     2404    dX->data.F32[i] = nX->nums->data.S32[i];
     2405    vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]);
     2406}
     2407for (int i = 0; i < nY->nums->n; i++) {
     2408    dY->data.F32[i] = nY->nums->data.S32[i];
     2409    vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]);
     2410}
     2411
     2412graphdata.color = KapaColorByName ("black");
     2413graphdata.ptype = 0;
     2414graphdata.size = 0.0;
     2415graphdata.style = 0;
     2416KapaPrepPlot (myKapa, dX->n, &graphdata);
     2417KapaPlotVector (myKapa, dX->n, dX->data.F32, "x");
     2418KapaPlotVector (myKapa, vX->n, vX->data.F32, "y");
     2419
     2420psFree (nX);
     2421psFree (dX);
     2422psFree (vX);
     2423
     2424psFree (nY);
     2425psFree (dY);
     2426psFree (vY);
     2427
     2428# endif
     2429
Note: See TracChangeset for help on using the changeset viewer.