IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 2, 2009, 3:12:47 PM (17 years ago)
Author:
eugene
Message:

check in changes from gene@development branch : extensive changes to moments calculation, psf model generation, aperture residuals; improvements to extended source analysis

Location:
trunk/psphot/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src

    • Property svn:ignore

      •  

        old new  
        1818psphotVersionDefinitions.h
        1919psphotMomentsStudy
         20psphotPetrosianStudy

  • trunk/psphot/src/psphotVisual.c

    r24636 r25755  
    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}
    2461
    2562bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel) {
    … …  
    131168    if (!pmVisualIsVisual()) return true;
    132169
    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     }
     170    int kapa = psphotKapaChannel (1);
     171    if (kapa == -1) return false;
    141172
    142173    // psphotVisualShowMask (kapa, readout->mask, "mask", 2);
    … …  
    144175    psphotVisualScaleImage (kapa, readout->image, "image", 0);
    145176
    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     }
     177    pmVisualAskUser(NULL);
    153178    return true;
    154179}
    … …  
    160185    if (!pmVisualIsVisual()) return true;
    161186
    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     }
     187    int kapa = psphotKapaChannel (1);
     188    if (kapa == -1) return false;
    170189
    171190    bool status = false;
    … …  
    181200    psphotVisualScaleImage (kapa, readout->image, "backsub", 0);
    182201
    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     }
     202    pmVisualAskUser(NULL);
    190203    return true;
    191204}
    … …  
    195208    if (!pmVisualIsVisual()) return true;
    196209
    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;
     210    int kapa = psphotKapaChannel (1);
     211    if (kapa == -1) return false;
     212
    207213    psphotVisualRangeImage (kapa, image, "signif", 2, -1.0, 25.0*25.0);
    208214
    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     }
     215    pmVisualAskUser(NULL);
    216216    return true;
    217217}
    … …  
    224224    if (!pmVisualIsVisual()) return true;
    225225
    226     if (kapa == -1) {
    227         fprintf (stderr, "kapa not opened, skipping\n");
    228         return false;
    229     }
     226    int kapa = psphotKapaChannel (1);
     227    if (kapa == -1) return false;
    230228
    231229    psArray *peaks = detections->peaks;
    … …  
    233231    // note: this uses the Ohana allocation tools:
    234232    // ALLOCATE (overlay, KiiOverlay, 3*peaks->n + 1);
    235     ALLOCATE (overlay, KiiOverlay, peaks->n);
     233    ALLOCATE (overlay, KiiOverlay, peaks->n + 2);
    236234
    237235    Noverlay = 0;
    … …  
    271269    }
    272270
    273 # if (0)
     271# if (1)
    274272    overlay[Noverlay].type = KII_OVERLAY_BOX;
    275273    overlay[Noverlay].x = 10.0;
    276274    overlay[Noverlay].y = 10.0;
    277     overlay[Noverlay].dx = 0.5;
    278     overlay[Noverlay].dy = 0.5;
     275    overlay[Noverlay].dx = 1.0;
     276    overlay[Noverlay].dy = 1.0;
    279277    overlay[Noverlay].angle = 0.0;
    280278    overlay[Noverlay].text = NULL;
    … …  
    285283    FREE (overlay);
    286284
    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     }
     285    pmVisualAskUser(NULL);
    294286    return true;
    295287}
    … …  
    302294    if (!pmVisualIsVisual()) return true;
    303295
    304     if (kapa == -1) {
    305         fprintf (stderr, "kapa not opened, skipping\n");
    306         return false;
    307     }
     296    int kapa = psphotKapaChannel (1);
     297    if (kapa == -1) return false;
    308298
    309299    psArray *footprints = detections->footprints;
    … …  
    325315
    326316        // draw the top
     317        // XXX need to allow top (and bottom) to have more than one span
    327318        span = footprint->spans->data[0];
    328319        overlay[Noverlay].type = KII_OVERLAY_LINE;
    … …  
    399390    FREE (overlay);
    400391
    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     }
     392    pmVisualAskUser(NULL);
    408393    return true;
    409394}
    … …  
    419404    if (!pmVisualIsVisual()) return true;
    420405
    421     if (kapa == -1) {
    422         fprintf (stderr, "kapa not opened, skipping\n");
    423         return false;
    424     }
     406    int kapa = psphotKapaChannel (1);
     407    if (kapa == -1) return false;
     408
     409    // XXX mark the different source classes with different color/shape dots
     410    // XXX are moments S/N and peak S/N consistent?
    425411
    426412    // note: this uses the Ohana allocation tools:
    … …  
    448434        overlay[Noverlay].dx = 2.0*axes.major;
    449435        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
     436
     437        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
     438
    451439        overlay[Noverlay].text = NULL;
    452440        Noverlay ++;
    … …  
    456444    FREE (overlay);
    457445
    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 
     446    pmVisualAskUser(NULL);
    466447    return true;
    467448}
    … …  
    475456    if (!pmVisualIsVisual()) return true;
    476457
    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);
     458    int myKapa = psphotKapaChannel (2);
     459    if (myKapa == -1) return false;
     460
     461    KapaClearPlots (myKapa);
    487462    KapaInitGraph (&graphdata);
    488     KapaSetFont (kapa3, "courier", 14);
     463    KapaSetFont (myKapa, "courier", 14);
    489464
    490465    float SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM");
    491466
    492467    // 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;
     468    float Xmin = 1000.0, Xmax = 0.0;
     469    float Ymin = 1000.0, Ymax = 0.0;
    495470    {
    496471        int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
    … …  
    511486            float Y1 = psfY + 4.0*psfdY;
    512487
    513             if (isfinite(X0)) { Xmin = PS_MAX(Xmin, X0); }
     488            if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); }
    514489            if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
    515             if (isfinite(Y0)) { Ymin = PS_MAX(Ymin, Y0); }
     490            if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); }
    516491            if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
    517492        }
    518493    }
     494    Xmin = PS_MAX(Xmin, -0.1);
     495    Ymin = PS_MAX(Ymin, -0.1);
    519496
    520497    // storage vectors for data to be plotted
    … …  
    564541    section.y  = 0.00;
    565542    section.name = psStringCopy ("MxxMyy");
    566     KapaSetSection (kapa3, &section);
     543    KapaSetSection (myKapa, &section);
    567544    psFree (section.name);
    568545
    … …  
    572549    graphdata.xmax = Xmax;
    573550    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);
     551    KapaSetLimits (myKapa, &graphdata);
     552
     553    KapaBox (myKapa, &graphdata);
     554    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
     555    KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM);
    579556
    580557    graphdata.color = KapaColorByName ("black");
    … …  
    582559    graphdata.size = 0.3;
    583560    graphdata.style = 2;
    584     KapaPrepPlot (kapa3, nF, &graphdata);
    585     KapaPlotVector (kapa3, nF, xFaint->data.F32, "x");
    586     KapaPlotVector (kapa3, nF, yFaint->data.F32, "y");
     561    KapaPrepPlot (myKapa, nF, &graphdata);
     562    KapaPlotVector (myKapa, nF, xFaint->data.F32, "x");
     563    KapaPlotVector (myKapa, nF, yFaint->data.F32, "y");
    587564
    588565    graphdata.color = KapaColorByName ("red");
    … …  
    590567    graphdata.size = 0.5;
    591568    graphdata.style = 2;
    592     KapaPrepPlot (kapa3, nB, &graphdata);
    593     KapaPlotVector (kapa3, nB, xBright->data.F32, "x");
    594     KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
     569    KapaPrepPlot (myKapa, nB, &graphdata);
     570    KapaPlotVector (myKapa, nB, xBright->data.F32, "x");
     571    KapaPlotVector (myKapa, nB, yBright->data.F32, "y");
    595572
    596573    // second section: MagMyy
    … …  
    600577    section.y  = 0.80;
    601578    section.name = psStringCopy ("MagMyy");
    602     KapaSetSection (kapa3, &section);
     579    KapaSetSection (myKapa, &section);
    603580    psFree (section.name);
    604581
    … …  
    608585    graphdata.ymin = Ymin;
    609586    graphdata.ymax = Ymax;
    610     KapaSetLimits (kapa3, &graphdata);
     587    KapaSetLimits (myKapa, &graphdata);
    611588
    612589    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);
     590    KapaBox (myKapa, &graphdata);
     591    KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_XP);
     592    KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM);
    616593
    617594    graphdata.color = KapaColorByName ("black");
    … …  
    619596    graphdata.size = 0.3;
    620597    graphdata.style = 2;
    621     KapaPrepPlot (kapa3, nF, &graphdata);
    622     KapaPlotVector (kapa3, nF, mFaint->data.F32, "x");
    623     KapaPlotVector (kapa3, nF, yFaint->data.F32, "y");
     598    KapaPrepPlot (myKapa, nF, &graphdata);
     599    KapaPlotVector (myKapa, nF, mFaint->data.F32, "x");
     600    KapaPlotVector (myKapa, nF, yFaint->data.F32, "y");
    624601
    625602    graphdata.color = KapaColorByName ("red");
    … …  
    627604    graphdata.size = 0.5;
    628605    graphdata.style = 2;
    629     KapaPrepPlot (kapa3, nB, &graphdata);
    630     KapaPlotVector (kapa3, nB, mBright->data.F32, "x");
    631     KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
     606    KapaPrepPlot (myKapa, nB, &graphdata);
     607    KapaPlotVector (myKapa, nB, mBright->data.F32, "x");
     608    KapaPlotVector (myKapa, nB, yBright->data.F32, "y");
    632609
    633610    // third section: MagMxx
    … …  
    637614    section.y  = 0.00;
    638615    section.name = psStringCopy ("MagMxx");
    639     KapaSetSection (kapa3, &section);
     616    KapaSetSection (myKapa, &section);
    640617    psFree (section.name);
    641618
    … …  
    645622    graphdata.ymin =  -7.9;
    646623    graphdata.ymax = -17.1;
    647     KapaSetLimits (kapa3, &graphdata);
     624    KapaSetLimits (myKapa, &graphdata);
    648625
    649626    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);
     627    KapaBox (myKapa, &graphdata);
     628    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
     629    KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_YP);
    653630
    654631    graphdata.color = KapaColorByName ("black");
    … …  
    656633    graphdata.size = 0.3;
    657634    graphdata.style = 2;
    658     KapaPrepPlot (kapa3, nF, &graphdata);
    659     KapaPlotVector (kapa3, nF, xFaint->data.F32, "x");
    660     KapaPlotVector (kapa3, nF, mFaint->data.F32, "y");
     635    KapaPrepPlot (myKapa, nF, &graphdata);
     636    KapaPlotVector (myKapa, nF, xFaint->data.F32, "x");
     637    KapaPlotVector (myKapa, nF, mFaint->data.F32, "y");
    661638
    662639    graphdata.color = KapaColorByName ("red");
    … …  
    664641    graphdata.size = 0.5;
    665642    graphdata.style = 2;
    666     KapaPrepPlot (kapa3, nB, &graphdata);
    667     KapaPlotVector (kapa3, nB, xBright->data.F32, "x");
    668     KapaPlotVector (kapa3, nB, mBright->data.F32, "y");
     643    KapaPrepPlot (myKapa, nB, &graphdata);
     644    KapaPlotVector (myKapa, nB, xBright->data.F32, "x");
     645    KapaPlotVector (myKapa, nB, mBright->data.F32, "y");
    669646
    670647    // draw N circles to outline the clumps
    671648    {
    672         KapaSelectSection (kapa3, "MxxMyy");
     649        KapaSelectSection (myKapa, "MxxMyy");
    673650
    674651        // draw a circle centered on psfX,Y with size of the psf limit
    … …  
    686663        graphdata.xmax = Xmax;
    687664        graphdata.ymax = Ymax;
    688         KapaSetLimits (kapa3, &graphdata);
     665        KapaSetLimits (myKapa, &graphdata);
    689666
    690667        for (int n = 0; n < nRegions; n++) {
    … …  
    705682                yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
    706683            }
    707             KapaPrepPlot (kapa3, xLimit->n, &graphdata);
    708             KapaPlotVector (kapa3, xLimit->n, xLimit->data.F32, "x");
    709             KapaPlotVector (kapa3, yLimit->n, yLimit->data.F32, "y");
     684            KapaPrepPlot (myKapa, xLimit->n, &graphdata);
     685            KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x");
     686            KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y");
    710687        }
    711688        psFree (xLimit);
    712689        psFree (yLimit);
    713690    }
    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
    750691
    751692    psFree (xBright);
    … …  
    756697    psFree (mFaint);
    757698
    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     }
     699    pmVisualAskUser(NULL);
    765700    return true;
    766701}
    767702
    768703// assumes 'kapa' value is checked and set
    769 bool psphotVisualShowRoughClass_Single (psArray *sources, pmSourceType type, pmSourceMode mode, char *color) {
     704bool psphotVisualShowRoughClass_Single (int myKapa, psArray *sources, pmSourceType type, pmSourceMode mode, char *color) {
    770705
    771706    int Noverlay;
    … …  
    802737        overlay[Noverlay].dx = 2.0*axes.major;
    803738        overlay[Noverlay].dy = 2.0*axes.minor;
    804         overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD;
     739        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
    805740        overlay[Noverlay].text = NULL;
    806741        Noverlay ++;
    807742    }
    808743
    809     KiiLoadOverlay (kapa, overlay, Noverlay, color);
     744    KiiLoadOverlay (myKapa, overlay, Noverlay, color);
    810745    FREE (overlay);
    811746
    … …  
    817752    if (!pmVisualIsVisual()) return true;
    818753
    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];
     754    int myKapa = psphotKapaChannel (1);
     755    if (myKapa == -1) return false;
     756
     757    KiiEraseOverlay (myKapa, "yellow"); // moments
     758
     759    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_STAR, 0, "red");
     760    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_EXTENDED, 0, "blue");
     761    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_DEFECT, 0, "blue");
     762    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_SATURATED, 0, "red");
     763    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_PSFSTAR, "yellow");
     764    psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_SATSTAR, "green");
     765
    836766    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     }
     767    pmVisualAskUser(NULL);
    841768    return true;
    842769}
    … …  
    846773    if (!pmVisualIsVisual()) return true;
    847774
    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     }
     775    int myKapa = psphotKapaChannel (3);
     776    if (myKapa == -1) return false;
    856777
    857778    int DX = 64;
    … …  
    898819
    899820    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);
     821    psphotVisualRangeImage (myKapa, psfLogFlux, "psf_mosaic",    0, -2.0, 3.0);
     822    psphotVisualRangeImage (myKapa, funMosaic, "psf_analytical", 1, -10.0, 100.0);
     823    psphotVisualRangeImage (myKapa, resMosaic, "psf_residual",   2, -10.0, 100.0);
    903824
    904825    psFree (psfMosaic);
    … …  
    908829    psFree (modelRef);
    909830
    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     }
     831    pmVisualAskUser(NULL);
    917832    return true;
    918833}
    … …  
    924839    if (!pmVisualIsVisual()) return true;
    925840
    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     }
     841    int myKapa = psphotKapaChannel (3);
     842    if (myKapa == -1) return false;
    934843
    935844    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    1017926            if (Xo == 0) {
    1018927                // place source alone on this row
    1019                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     928                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     929                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    1020930                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
    1021931
    1022                 psphotSubWithTest (source, false, maskVal); // remove source (force)
     932                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    1023933                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
    1024                 psphotSetState (source, false, maskVal); // reset source Add/Sub state to recorded
     934
     935                if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    1025936
    1026937                Yo += DY;
    … …  
    1032943                Xo = 0;
    1033944
    1034                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     945                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     946                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    1035947                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
    1036948
    1037                 psphotSubWithTest (source, false, maskVal); // remove source (force)
     949                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    1038950                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
    1039                 psphotSetState (source, false, maskVal); // replace source (has been subtracted)
     951
     952                if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    1040953
    1041954                Xo = DX;
    … …  
    1044957        } else {
    1045958            // extend this row
    1046             psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     959            bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     960            if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    1047961            psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
    1048962
    1049             psphotSubWithTest (source, false, maskVal); // remove source (force)
     963            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    1050964            psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
    1051             psphotSetState (source, false, maskVal); // replace source (has been subtracted)
     965            if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    1052966
    1053967            Xo += DX;
    … …  
    1056970    }
    1057971
    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 
     972    psphotVisualRangeImage (myKapa, outpos, "psfpos", 0, -0.05, 0.95);
     973    psphotVisualRangeImage (myKapa, outsub, "psfsub", 1, -0.05, 0.95);
     974
     975    pmVisualAskUser(NULL);
    1069976    psFree (outpos);
    1070977    psFree (outsub);
    … …  
    1084991    if (!pmVisualIsVisual()) return true;
    1085992
    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     }
     993    int myKapa = psphotKapaChannel (3);
     994    if (myKapa == -1) return false;
    1094995
    1095996    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    11191020        pmSource *source = sources->data[i];
    11201021
    1121         bool keep = false;
    1122         keep |= (source->mode & PM_SOURCE_MODE_SATSTAR);
    1123         if (!keep) continue;
     1022        // only show "real" saturated stars (not defects)
     1023        if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;;
     1024        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;;
    11241025
    11251026        // how does this subimage get placed into the output image?
    … …  
    11641065        pmSource *source = sources->data[i];
    11651066
    1166         bool keep = false;
    1167         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    1168             nSAT ++;
    1169             keep = true;
    1170         }
    1171         if (!keep) continue;
     1067        // only show "real" saturated stars (not defects)
     1068        if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;;
     1069        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;;
     1070        nSAT ++;
    11721071
    11731072        if (Xo + DX > NX) {
    … …  
    11751074            if (Xo == 0) {
    11761075                // place source alone on this row
    1177                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     1076                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     1077                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    11781078                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
    1179                 psphotSetState (source, true, maskVal); // reset source Add/Sub state to recorded
     1079                if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    11801080
    11811081                Yo += DY;
    … …  
    11861086                Yo += dY;
    11871087                Xo = 0;
    1188                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     1088
     1089                bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     1090                if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    11891091                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
    1190                 psphotSetState (source, true, maskVal); // replace source (has been subtracted)
     1092                if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    11911093
    11921094                Xo = DX;
    … …  
    11951097        } else {
    11961098            // extend this row
    1197             psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     1099            bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     1100            if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    11981101            psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
    1199             psphotSetState (source, true, maskVal); // replace source (has been subtracted)
     1102            if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    12001103
    12011104            Xo += DX;
    … …  
    12041107    }
    12051108
    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 
     1109    psphotVisualScaleImage (myKapa, outsat, "satstar", 2);
     1110
     1111    pmVisualAskUser(NULL);
    12161112    psFree (outsat);
    12171113    return true;
    … …  
    12221118    Graphdata graphdata;
    12231119
    1224     bool state = !(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
    1225     psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     1120    bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
     1121    if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    12261122
    12271123    int nPts = source->pixels->numRows * source->pixels->numCols;
    … …  
    12401136        for (int ix = 0; ix < source->pixels->numCols; ix++) {
    12411137            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) ;
     1138                rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
     1139                // rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
    12441140                Rb->data.F32[nb] = log10(rb->data.F32[nb]);
    12451141                fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]);
    12461142                nb++;
    12471143            } 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) ;
     1144                rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
     1145                // rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ;
    12501146                Rg->data.F32[ng] = log10(rg->data.F32[ng]);
    12511147                fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]);
    … …  
    12561152
    12571153    // reset source Add/Sub state to recorded
    1258     psphotSetState (source, state, maskVal);
     1154    if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    12591155
    12601156    KapaInitGraph (&graphdata);
    … …  
    13371233    if (!pmVisualIsVisual()) return true;
    13381234
    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     }
     1235    int myKapa = psphotKapaChannel (2);
     1236    if (myKapa == -1) return false;
    13471237
    13481238    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    13511241    assert (maskVal);
    13521242
    1353     KapaClearPlots (kapa3);
     1243    KapaClearPlots (myKapa);
    13541244    // first section : mag vs CR nSigma
    13551245    section.dx = 1.0;
    … …  
    13591249    section.name = NULL;
    13601250    psStringAppend (&section.name, "linlog");
    1361     KapaSetSection (kapa3, &section);
     1251    KapaSetSection (myKapa, &section);
    13621252    psFree (section.name);
    13631253
    … …  
    13691259    section.name = NULL;
    13701260    psStringAppend (&section.name, "loglog");
    1371     KapaSetSection (kapa3, &section);
     1261    KapaSetSection (myKapa, &section);
    13721262    psFree (section.name);
    13731263
    … …  
    13781268        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
    13791269
    1380         psphotVisualPlotRadialProfile (kapa3, source, maskVal);
     1270        psphotVisualPlotRadialProfile (myKapa, source, maskVal);
    13811271
    13821272        // pause and wait for user input:
    … …  
    13881278        }
    13891279        if (key[0] == 'e') {
    1390             KapaClearPlots (kapa3);
     1280            KapaClearPlots (myKapa);
    13911281        }
    13921282        if (key[0] == 's') {
    … …  
    14071297    psEllipseAxes axes;
    14081298
     1299    // XXX skip this for now: it is not very clear
     1300    return true;
     1301
    14091302    if (!pmVisualIsVisual()) return true;
    14101303
    1411     if (kapa == -1) {
    1412         fprintf (stderr, "kapa not opened, skipping\n");
    1413         return false;
    1414     }
     1304    int myKapa = psphotKapaChannel (1);
     1305    if (myKapa == -1) return false;
    14151306
    14161307    // note: this uses the Ohana allocation tools:
    … …  
    14851376    }
    14861377
    1487     KiiLoadOverlay (kapa, overlayE, NoverlayE, "red");
    1488     KiiLoadOverlay (kapa, overlayO, NoverlayO, "yellow");
     1378    KiiLoadOverlay (myKapa, overlayE, NoverlayE, "red");
     1379    KiiLoadOverlay (myKapa, overlayO, NoverlayO, "yellow");
    14891380    FREE (overlayE);
    14901381    FREE (overlayO);
    14911382
    1492     // pause and wait for user input:
    1493     // continue, save (provide name), ??
    1494     char key[10];
    14951383    fprintf (stdout, "even bits (0x0001, 0x0004, ... : red\n");
    14961384    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;
     1385    pmVisualAskUser(NULL);
     1386
     1387    return true;
     1388}
     1389
     1390bool psphotVisualShowSourceSize_Single (int myKapa, psArray *sources, pmSourceMode mode, bool keep, float scale, char *color) {
     1391
     1392    int Noverlay;
    15081393    KiiOverlay *overlay;
    15091394
    1510     if (!pmVisualIsVisual()) return true;
    1511 
    1512     if (kapa == -1) {
    1513         fprintf (stderr, "kapa not opened, skipping\n");
    1514         return false;
    1515     }
     1395    psEllipseMoments emoments;
     1396    psEllipseAxes axes;
    15161397
    15171398    // note: this uses the Ohana allocation tools:
     1399    ALLOCATE (overlay, KiiOverlay, sources->n);
     1400
    15181401    Noverlay = 0;
    1519     NOVERLAY = 100;
    1520     ALLOCATE (overlay, KiiOverlay, sources->n);
    1521 
    1522     // mark CRs with red boxes
    15231402    for (int i = 0; i < sources->n; i++) {
    15241403
    … …  
    15261405        if (source == NULL) continue;
    15271406
    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;
     1407        // if (source->type != type) continue;
     1408        if (mode) {
     1409            if (keep) {
     1410                if (!(source->mode & mode)) continue;
     1411            } else {
     1412                if (source->mode & mode) continue;
     1413            }
     1414        }
     1415
     1416        pmMoments *moments = source->moments;
     1417        if (moments == NULL) continue;
     1418
     1419        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     1420        overlay[Noverlay].x = moments->Mx;
     1421        overlay[Noverlay].y = moments->My;
     1422
     1423        emoments.x2 = moments->Mxx;
     1424        emoments.y2 = moments->Myy;
     1425        emoments.xy = moments->Mxy;
     1426
     1427        axes = psEllipseMomentsToAxes (emoments, 20.0);
     1428
     1429        overlay[Noverlay].dx = scale*2.0*axes.major;
     1430        overlay[Noverlay].dy = scale*2.0*axes.minor;
     1431        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
    15371432        overlay[Noverlay].text = NULL;
    15381433        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");
     1434    }
     1435
     1436    KiiLoadOverlay (myKapa, overlay, Noverlay, color);
    15661437    FREE (overlay);
    15671438
    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 
     1439    return true;
     1440}
     1441
     1442bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) {
     1443
     1444    if (!pmVisualIsVisual()) return true;
     1445
     1446    int myKapa = psphotKapaChannel (1);
     1447    if (myKapa == -1) return false;
     1448
     1449    KiiEraseOverlay (myKapa, "red");
     1450    KiiEraseOverlay (myKapa, "blue");
     1451    KiiEraseOverlay (myKapa, "green");
     1452    KiiEraseOverlay (myKapa, "yellow");
     1453
     1454    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");
     1455    psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_EXT_LIMIT, 1, 1.0, "blue");
     1456    psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_CR_LIMIT, 1, 1.0, "red");
     1457    psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_DEFECT, 1, 2.0, "red");
     1458    psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_SATSTAR, 1, 1.0, "yellow");
     1459
     1460    fprintf (stdout, "red: CR; blue: EXTENDED; green: PSF-like; yellow: SATSTAR\n");
     1461    pmVisualAskUser(NULL);
     1462    return true;
     1463}
     1464
     1465bool psphotVisualPlotSourceSize (psMetadata *recipe, psArray *sources) {
     1466
     1467    bool status;
    15841468    Graphdata graphdata;
    15851469    KapaSection section;
    … …  
    15871471    if (!pmVisualIsVisual()) return true;
    15881472
    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;
     1473    int myKapa = psphotKapaChannel (2);
     1474    if (myKapa == -1) return false;
     1475
     1476    KapaClearPlots (myKapa);
     1477    KapaInitGraph (&graphdata);
     1478    KapaSetFont (myKapa, "courier", 14);
     1479
     1480    // select the max psfX,Y values for the plot limits
     1481    float Xmin = 1000.0, Xmax = 0.0;
     1482    float Ymin = 1000.0, Ymax = 0.0;
     1483    {
     1484        int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
     1485        for (int n = 0; n < nRegions; n++) {
     1486
     1487            char regionName[64];
     1488            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
     1489            psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
     1490
     1491            float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
     1492            float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
     1493            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
     1494            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
     1495
     1496            float X0 = psfX - 10.0*psfdX;
     1497            float X1 = psfX + 10.0*psfdX;
     1498            float Y0 = psfY - 10.0*psfdY;
     1499            float Y1 = psfY + 10.0*psfdY;
     1500
     1501            if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); }
     1502            if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
     1503            if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); }
     1504            if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
    15951505        }
    15961506    }
    1597 
    1598     KapaClearPlots (kapa3);
     1507    Xmin = PS_MAX(Xmin, -0.1);
     1508    Ymin = PS_MAX(Ymin, -0.1);
     1509
     1510    // storage vectors for data to be plotted
     1511    psVector *xSAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1512    psVector *ySAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1513    psVector *mSAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1514    psVector *sSAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1515
     1516    psVector *xPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1517    psVector *yPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1518    psVector *mPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1519    psVector *sPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1520
     1521    psVector *xEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1522    psVector *yEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1523    psVector *mEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1524    psVector *sEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1525
     1526    psVector *xDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1527    psVector *yDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1528    psVector *mDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1529    psVector *sDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1530
     1531    psVector *xCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1532    psVector *yCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1533    psVector *mCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1534    psVector *sCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1535
     1536    // construct the vectors
     1537    int nSAT = 0;
     1538    int nEXT = 0;
     1539    int nPSF = 0;
     1540    int nDEF = 0;
     1541    int nCR  = 0;
     1542    for (int i = 0; i < sources->n; i++) {
     1543        pmSource *source = sources->data[i];
     1544        if (source->moments == NULL) continue;
     1545
     1546        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) {
     1547            xCR->data.F32[nCR] = source->moments->Mxx;
     1548            yCR->data.F32[nCR] = source->moments->Myy;
     1549            mCR->data.F32[nCR] = -2.5*log10(source->moments->Sum);
     1550            sCR->data.F32[nCR] = source->extNsigma;
     1551            nCR++;
     1552        }
     1553        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
     1554            xSAT->data.F32[nSAT] = source->moments->Mxx;
     1555            ySAT->data.F32[nSAT] = source->moments->Myy;
     1556            mSAT->data.F32[nSAT] = -2.5*log10(source->moments->Sum);
     1557            sSAT->data.F32[nSAT] = source->extNsigma;
     1558            nSAT++;
     1559        }
     1560        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     1561            xEXT->data.F32[nEXT] = source->moments->Mxx;
     1562            yEXT->data.F32[nEXT] = source->moments->Myy;
     1563            mEXT->data.F32[nEXT] = -2.5*log10(source->moments->Sum);
     1564            sEXT->data.F32[nEXT] = source->extNsigma;
     1565            nEXT++;
     1566            continue;
     1567        }
     1568        if (source->mode & PM_SOURCE_MODE_DEFECT) {
     1569            xDEF->data.F32[nDEF] = source->moments->Mxx;
     1570            yDEF->data.F32[nDEF] = source->moments->Myy;
     1571            mDEF->data.F32[nDEF] = -2.5*log10(source->moments->Sum);
     1572            sDEF->data.F32[nDEF] = source->extNsigma;
     1573            nDEF++;
     1574            continue;
     1575        }
     1576        if ((source->mode & PM_SOURCE_MODE_CR_LIMIT) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
     1577            continue;
     1578        }
     1579        xPSF->data.F32[nPSF] = source->moments->Mxx;
     1580        yPSF->data.F32[nPSF] = source->moments->Myy;
     1581        mPSF->data.F32[nPSF] = -2.5*log10(source->moments->Sum);
     1582        sPSF->data.F32[nPSF] = source->extNsigma;
     1583        nPSF++;
     1584    }
     1585    xSAT->n = nSAT;
     1586    ySAT->n = nSAT;
     1587    mSAT->n = nSAT;
     1588    sSAT->n = nSAT;
     1589
     1590    xPSF->n = nPSF;
     1591    yPSF->n = nPSF;
     1592    mPSF->n = nPSF;
     1593    sPSF->n = nPSF;
     1594
     1595    xEXT->n = nEXT;
     1596    yEXT->n = nEXT;
     1597    mEXT->n = nEXT;
     1598    sEXT->n = nEXT;
     1599
     1600    xCR->n = nCR;
     1601    yCR->n = nCR;
     1602    mCR->n = nCR;
     1603    sCR->n = nCR;
     1604
     1605    xDEF->n = nDEF;
     1606    yDEF->n = nDEF;
     1607    mDEF->n = nDEF;
     1608    sDEF->n = nDEF;
     1609
     1610    // four sections: MxxMyy, MagMxx, MagMyy, MagSigma
     1611
     1612    // first section: MxxMyy
     1613    section.dx = 0.75;
     1614    section.dy = 0.60;
     1615    section.x  = 0.00;
     1616    section.y  = 0.00;
     1617    section.name = psStringCopy ("MxxMyy");
     1618    KapaSetSection (myKapa, &section);
     1619    psFree (section.name);
     1620
     1621    graphdata.color = KapaColorByName ("black");
     1622    graphdata.xmin = Xmin;
     1623    graphdata.ymin = Ymin;
     1624    graphdata.xmax = Xmax;
     1625    graphdata.ymax = Ymax;
     1626    KapaSetLimits (myKapa, &graphdata);
     1627
     1628    KapaBox (myKapa, &graphdata);
     1629    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
     1630    KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM);
     1631
     1632    graphdata.color = KapaColorByName ("black");
     1633    graphdata.ptype = 0;
     1634    graphdata.size = 0.5;
     1635    graphdata.style = 2;
     1636    KapaPrepPlot   (myKapa, nPSF, &graphdata);
     1637    KapaPlotVector (myKapa, nPSF, xPSF->data.F32, "x");
     1638    KapaPlotVector (myKapa, nPSF, yPSF->data.F32, "y");
     1639
     1640    graphdata.color = KapaColorByName ("blue");
     1641    graphdata.ptype = 0;
     1642    graphdata.size = 0.5;
     1643    graphdata.style = 2;
     1644    KapaPrepPlot   (myKapa, nEXT, &graphdata);
     1645    KapaPlotVector (myKapa, nEXT, xEXT->data.F32, "x");
     1646    KapaPlotVector (myKapa, nEXT, yEXT->data.F32, "y");
     1647
     1648    graphdata.color = KapaColorByName ("red");
     1649    graphdata.ptype = 0;
     1650    graphdata.size = 0.5;
     1651    graphdata.style = 2;
     1652    KapaPrepPlot   (myKapa, nDEF, &graphdata);
     1653    KapaPlotVector (myKapa, nDEF, xDEF->data.F32, "x");
     1654    KapaPlotVector (myKapa, nDEF, yDEF->data.F32, "y");
     1655
     1656    graphdata.color = KapaColorByName ("red");
     1657    graphdata.ptype = 7;
     1658    graphdata.size = 1.0;
     1659    graphdata.style = 2;
     1660    KapaPrepPlot   (myKapa, nCR, &graphdata);
     1661    KapaPlotVector (myKapa, nCR, xCR->data.F32, "x");
     1662    KapaPlotVector (myKapa, nCR, yCR->data.F32, "y");
     1663
     1664    graphdata.color = KapaColorByName ("blue");
     1665    graphdata.ptype = 7;
     1666    graphdata.size = 1.0;
     1667    graphdata.style = 2;
     1668    KapaPrepPlot   (myKapa, nSAT, &graphdata);
     1669    KapaPlotVector (myKapa, nSAT, xSAT->data.F32, "x");
     1670    KapaPlotVector (myKapa, nSAT, ySAT->data.F32, "y");
     1671
     1672    // second section: MagMyy
     1673    section.dx = 0.75;
     1674    section.dy = 0.20;
     1675    section.x  = 0.00;
     1676    section.y  = 0.80;
     1677    section.name = psStringCopy ("MagMyy");
     1678    KapaSetSection (myKapa, &section);
     1679    psFree (section.name);
     1680
     1681    graphdata.color = KapaColorByName ("black");
     1682    graphdata.xmin = -17.1;
     1683    graphdata.xmax =  -7.9;
     1684    graphdata.ymin = Ymin;
     1685    graphdata.ymax = Ymax;
     1686    KapaSetLimits (myKapa, &graphdata);
     1687
     1688    strcpy (graphdata.labels, "0210");
     1689    KapaBox (myKapa, &graphdata);
     1690    KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_XP);
     1691    KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM);
     1692
     1693    graphdata.color = KapaColorByName ("black");
     1694    graphdata.ptype = 0;
     1695    graphdata.size = 0.5;
     1696    graphdata.style = 2;
     1697    KapaPrepPlot   (myKapa, nPSF, &graphdata);
     1698    KapaPlotVector (myKapa, nPSF, mPSF->data.F32, "x");
     1699    KapaPlotVector (myKapa, nPSF, yPSF->data.F32, "y");
     1700
     1701    graphdata.color = KapaColorByName ("blue");
     1702    graphdata.ptype = 0;
     1703    graphdata.size = 0.5;
     1704    graphdata.style = 2;
     1705    KapaPrepPlot   (myKapa, nEXT, &graphdata);
     1706    KapaPlotVector (myKapa, nEXT, mEXT->data.F32, "x");
     1707    KapaPlotVector (myKapa, nEXT, yEXT->data.F32, "y");
     1708
     1709    graphdata.color = KapaColorByName ("red");
     1710    graphdata.ptype = 0;
     1711    graphdata.size = 0.5;
     1712    graphdata.style = 2;
     1713    KapaPrepPlot   (myKapa, nDEF, &graphdata);
     1714    KapaPlotVector (myKapa, nDEF, mDEF->data.F32, "x");
     1715    KapaPlotVector (myKapa, nDEF, yDEF->data.F32, "y");
     1716
     1717    graphdata.color = KapaColorByName ("red");
     1718    graphdata.ptype = 7;
     1719    graphdata.size = 1.0;
     1720    graphdata.style = 2;
     1721    KapaPrepPlot   (myKapa, nCR, &graphdata);
     1722    KapaPlotVector (myKapa, nCR, mCR->data.F32, "x");
     1723    KapaPlotVector (myKapa, nCR, yCR->data.F32, "y");
     1724
     1725    graphdata.color = KapaColorByName ("blue");
     1726    graphdata.ptype = 7;
     1727    graphdata.size = 1.0;
     1728    graphdata.style = 2;
     1729    KapaPrepPlot   (myKapa, nSAT, &graphdata);
     1730    KapaPlotVector (myKapa, nSAT, mSAT->data.F32, "x");
     1731    KapaPlotVector (myKapa, nSAT, ySAT->data.F32, "y");
     1732
     1733    // third section: MagMxx
     1734    section.dx = 0.25;
     1735    section.dy = 0.60;
     1736    section.x  = 0.80;
     1737    section.y  = 0.00;
     1738    section.name = psStringCopy ("MagMxx");
     1739    KapaSetSection (myKapa, &section);
     1740    psFree (section.name);
     1741
     1742    graphdata.color = KapaColorByName ("black");
     1743    graphdata.xmin = Xmin;
     1744    graphdata.xmax = Xmax;
     1745    graphdata.ymin =  -7.9;
     1746    graphdata.ymax = -17.1;
     1747    KapaSetLimits (myKapa, &graphdata);
     1748
     1749    strcpy (graphdata.labels, "2001");
     1750    KapaBox (myKapa, &graphdata);
     1751    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
     1752    KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_YP);
     1753
     1754    graphdata.color = KapaColorByName ("black");
     1755    graphdata.ptype = 0;
     1756    graphdata.size = 0.5;
     1757    graphdata.style = 2;
     1758    KapaPrepPlot   (myKapa, nPSF, &graphdata);
     1759    KapaPlotVector (myKapa, nPSF, xPSF->data.F32, "x");
     1760    KapaPlotVector (myKapa, nPSF, mPSF->data.F32, "y");
     1761
     1762    graphdata.color = KapaColorByName ("blue");
     1763    graphdata.ptype = 0;
     1764    graphdata.size = 0.5;
     1765    graphdata.style = 2;
     1766    KapaPrepPlot   (myKapa, nEXT, &graphdata);
     1767    KapaPlotVector (myKapa, nEXT, xEXT->data.F32, "x");
     1768    KapaPlotVector (myKapa, nEXT, mEXT->data.F32, "y");
     1769
     1770    graphdata.color = KapaColorByName ("red");
     1771    graphdata.ptype = 0;
     1772    graphdata.size = 0.5;
     1773    graphdata.style = 2;
     1774    KapaPrepPlot   (myKapa, nDEF, &graphdata);
     1775    KapaPlotVector (myKapa, nDEF, xDEF->data.F32, "x");
     1776    KapaPlotVector (myKapa, nDEF, mDEF->data.F32, "y");
     1777
     1778    graphdata.color = KapaColorByName ("red");
     1779    graphdata.ptype = 7;
     1780    graphdata.size = 1.0;
     1781    graphdata.style = 2;
     1782    KapaPrepPlot   (myKapa, nCR, &graphdata);
     1783    KapaPlotVector (myKapa, nCR, xCR->data.F32, "x");
     1784    KapaPlotVector (myKapa, nCR, mCR->data.F32, "y");
     1785
     1786    graphdata.color = KapaColorByName ("blue");
     1787    graphdata.ptype = 7;
     1788    graphdata.size = 1.0;
     1789    graphdata.style = 2;
     1790    KapaPrepPlot   (myKapa, nSAT, &graphdata);
     1791    KapaPlotVector (myKapa, nSAT, xSAT->data.F32, "x");
     1792    KapaPlotVector (myKapa, nSAT, mSAT->data.F32, "y");
     1793
     1794    // fourth section: MagSigma
     1795    section.dx = 0.75;
     1796    section.dy = 0.15;
     1797    section.x  = 0.00;
     1798    section.y  = 0.65;
     1799    section.name = psStringCopy ("MagSigma");
     1800    KapaSetSection (myKapa, &section);
     1801    psFree (section.name);
     1802
     1803    graphdata.color = KapaColorByName ("black");
     1804    graphdata.xmax =  -7.9;
     1805    graphdata.xmin = -17.1;
     1806    graphdata.ymin = -20.1;
     1807    graphdata.ymax = +20.1;
     1808    KapaSetLimits (myKapa, &graphdata);
     1809
     1810    strcpy (graphdata.labels, "0100");
     1811    KapaBox (myKapa, &graphdata);
     1812    // KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_XM);
     1813    KapaSendLabel (myKapa, "EXT&ss&c", KAPA_LABEL_YM);
     1814
     1815    graphdata.color = KapaColorByName ("black");
     1816    graphdata.ptype = 0;
     1817    graphdata.size = 0.5;
     1818    graphdata.style = 2;
     1819    KapaPrepPlot   (myKapa, nPSF, &graphdata);
     1820    KapaPlotVector (myKapa, nPSF, mPSF->data.F32, "x");
     1821    KapaPlotVector (myKapa, nPSF, sPSF->data.F32, "y");
     1822
     1823    graphdata.color = KapaColorByName ("blue");
     1824    graphdata.ptype = 0;
     1825    graphdata.size = 0.5;
     1826    graphdata.style = 2;
     1827    KapaPrepPlot   (myKapa, nEXT, &graphdata);
     1828    KapaPlotVector (myKapa, nEXT, mEXT->data.F32, "x");
     1829    KapaPlotVector (myKapa, nEXT, sEXT->data.F32, "y");
     1830
     1831    graphdata.color = KapaColorByName ("red");
     1832    graphdata.ptype = 0;
     1833    graphdata.size = 0.5;
     1834    graphdata.style = 2;
     1835    KapaPrepPlot   (myKapa, nDEF, &graphdata);
     1836    KapaPlotVector (myKapa, nDEF, mDEF->data.F32, "x");
     1837    KapaPlotVector (myKapa, nDEF, sDEF->data.F32, "y");
     1838
     1839    graphdata.color = KapaColorByName ("red");
     1840    graphdata.ptype = 7;
     1841    graphdata.size = 1.0;
     1842    graphdata.style = 2;
     1843    KapaPrepPlot   (myKapa, nCR, &graphdata);
     1844    KapaPlotVector (myKapa, nCR, mCR->data.F32, "x");
     1845    KapaPlotVector (myKapa, nCR, sCR->data.F32, "y");
     1846
     1847    graphdata.color = KapaColorByName ("blue");
     1848    graphdata.ptype = 7;
     1849    graphdata.size = 1.0;
     1850    graphdata.style = 2;
     1851    KapaPrepPlot   (myKapa, nSAT, &graphdata);
     1852    KapaPlotVector (myKapa, nSAT, mSAT->data.F32, "x");
     1853    KapaPlotVector (myKapa, nSAT, sSAT->data.F32, "y");
     1854
     1855    // draw N circles to outline the clumps
     1856    {
     1857        KapaSelectSection (myKapa, "MxxMyy");
     1858
     1859        // draw a circle centered on psfX,Y with size of the psf limit
     1860        psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
     1861        psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
     1862
     1863        int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
     1864        float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
     1865
     1866        graphdata.color = KapaColorByName ("blue");
     1867        graphdata.style = 0;
     1868
     1869        graphdata.xmin = Xmin;
     1870        graphdata.ymin = Ymin;
     1871        graphdata.xmax = Xmax;
     1872        graphdata.ymax = Ymax;
     1873        KapaSetLimits (myKapa, &graphdata);
     1874
     1875        for (int n = 0; n < nRegions; n++) {
     1876
     1877            char regionName[64];
     1878            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
     1879            psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
     1880
     1881            float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
     1882            float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
     1883            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
     1884            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
     1885            float Rx = psfdX * PSF_CLUMP_NSIGMA;
     1886            float Ry = psfdY * PSF_CLUMP_NSIGMA;
     1887
     1888            for (int i = 0; i < xLimit->n; i++) {
     1889                xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX;
     1890                yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
     1891            }
     1892            KapaPrepPlot (myKapa, xLimit->n, &graphdata);
     1893            KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x");
     1894            KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y");
     1895        }
     1896        psFree (xLimit);
     1897        psFree (yLimit);
     1898    }
     1899
     1900    psFree (xSAT);
     1901    psFree (ySAT);
     1902    psFree (mSAT);
     1903
     1904    psFree (xEXT);
     1905    psFree (yEXT);
     1906    psFree (mEXT);
     1907
     1908    psFree (xPSF);
     1909    psFree (yPSF);
     1910    psFree (mPSF);
     1911
     1912    psFree (xDEF);
     1913    psFree (yDEF);
     1914    psFree (mDEF);
     1915
     1916    psFree (xCR);
     1917    psFree (yCR);
     1918    psFree (mCR);
     1919
     1920    pmVisualAskUser(NULL);
     1921    return true;
     1922}
     1923
     1924bool psphotVisualShowResidualImage (pmReadout *readout) {
     1925
     1926    if (!pmVisualIsVisual()) return true;
     1927
     1928    int myKapa = psphotKapaChannel (1);
     1929    if (myKapa == -1) return false;
     1930
     1931    psphotVisualScaleImage (myKapa, readout->image, "resid", 1);
     1932
     1933    pmVisualAskUser(NULL);
     1934    return true;
     1935}
     1936
     1937bool psphotVisualPlotApResid (psArray *sources, float mean, float error) {
     1938
     1939    Graphdata graphdata;
     1940    float lineX[2], lineY[2];
     1941
     1942    if (!pmVisualIsVisual()) return true;
     1943
     1944    int myKapa = psphotKapaChannel (2);
     1945    if (myKapa == -1) return false;
     1946
     1947    KapaClearPlots (myKapa);
    15991948    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);
    1609     psFree (section.name);
    16101949
    16111950    psVector *x = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    16121951    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);
     1952    psVector *dy = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    17701953
    17711954    graphdata.xmin = +32.0;
    … …  
    17851968        x->data.F32[n] = source->psfMag;
    17861969        y->data.F32[n] = source->apMag - source->psfMag;
     1970        dy->data.F32[n] = source->errMag;
    17871971        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
    17881972        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
    … …  
    17921976        n++;
    17931977    }
    1794     x->n = y->n = n;
     1978    x->n = y->n = dy->n = n;
    17951979
    17961980    float range;
    … …  
    18021986    graphdata.ymin -= 0.05*range;
    18031987
    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);
     1988    // XXX test
     1989    graphdata.xmin = -17.0;
     1990    graphdata.xmax =  -9.0;
     1991    graphdata.ymin = -0.31;
     1992    graphdata.ymax = +0.31;
     1993
     1994    KapaSetLimits (myKapa, &graphdata);
     1995
     1996    KapaSetFont (myKapa, "helvetica", 14);
     1997    KapaBox (myKapa, &graphdata);
     1998    KapaSendLabel (myKapa, "PSF Mag", KAPA_LABEL_XM);
     1999    KapaSendLabel (myKapa, "Ap Mag - PSF Mag", KAPA_LABEL_YM);
    18112000
    18122001    graphdata.color = KapaColorByName ("black");
    … …  
    18142003    graphdata.size = 0.5;
    18152004    graphdata.style = 2;
    1816     KapaPrepPlot (kapa3, n, &graphdata);
    1817     KapaPlotVector (kapa3, n, x->data.F32, "x");
    1818     KapaPlotVector (kapa3, n, y->data.F32, "y");
     2005    graphdata.etype |= 0x01;
     2006    KapaPrepPlot (myKapa, n, &graphdata);
     2007    KapaPlotVector (myKapa, n, x->data.F32, "x");
     2008    KapaPlotVector (myKapa, n, y->data.F32, "y");
     2009    KapaPlotVector (myKapa, n, dy->data.F32, "dym");
     2010    KapaPlotVector (myKapa, n, dy->data.F32, "dyp");
     2011
     2012    graphdata.color = KapaColorByName ("blue");
     2013    graphdata.ptype = 0;
     2014    graphdata.size = 0.5;
     2015    graphdata.style = 0;
     2016    graphdata.etype = 0;
     2017    lineX[0] = graphdata.xmin;
     2018    lineX[1] = graphdata.xmax;
     2019    lineY[0] = lineY[1] = mean;
     2020    KapaPrepPlot (myKapa, 2, &graphdata);
     2021    KapaPlotVector (myKapa, 2, lineX, "x");
     2022    KapaPlotVector (myKapa, 2, lineY, "y");
     2023
     2024    graphdata.color = KapaColorByName ("red");
     2025    graphdata.ptype = 0;
     2026    graphdata.size = 0.5;
     2027    graphdata.style = 0;
     2028    graphdata.etype = 0;
     2029    lineX[0] = graphdata.xmin;
     2030    lineX[1] = graphdata.xmax;
     2031    lineY[0] = lineY[1] = mean + error;
     2032    KapaPrepPlot (myKapa, 2, &graphdata);
     2033    KapaPlotVector (myKapa, 2, lineX, "x");
     2034    KapaPlotVector (myKapa, 2, lineY, "y");
     2035
     2036    graphdata.color = KapaColorByName ("red");
     2037    graphdata.ptype = 0;
     2038    graphdata.size = 0.5;
     2039    graphdata.style = 0;
     2040    graphdata.etype = 0;
     2041    lineX[0] = graphdata.xmin;
     2042    lineX[1] = graphdata.xmax;
     2043    lineY[0] = lineY[1] = mean - error;
     2044    KapaPrepPlot (myKapa, 2, &graphdata);
     2045    KapaPlotVector (myKapa, 2, lineX, "x");
     2046    KapaPlotVector (myKapa, 2, lineY, "y");
    18192047
    18202048    psFree (x);
    18212049    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     }
     2050    psFree (dy);
     2051
     2052    pmVisualAskUser(NULL);
     2053    return true;
     2054}
     2055
     2056bool psphotVisualShowPetrosians (psArray *sources) {
     2057
     2058    int Noverlay, NOVERLAY;
     2059    KiiOverlay *overlay;
     2060
     2061    if (!pmVisualIsVisual()) return true;
     2062
     2063    int kapa = psphotKapaChannel (1);
     2064    if (kapa == -1) return false;
     2065
     2066    Noverlay = 0;
     2067    NOVERLAY = 100;
     2068    ALLOCATE (overlay, KiiOverlay, NOVERLAY);
     2069
     2070    for (int i = 0; i < sources->n; i++) {
     2071        pmSource *source = sources->data[i];
     2072
     2073        if (!source) continue;
     2074        if (!source->extpars) continue;
     2075        if (!source->extpars->profile) continue;
     2076        if (!source->extpars->petrosian_80) continue;
     2077
     2078        pmSourceRadialProfile *profile = source->extpars->profile;
     2079        pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80;
     2080
     2081        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     2082        overlay[Noverlay].x = source->peak->xf;
     2083        overlay[Noverlay].y = source->peak->yf;
     2084        overlay[Noverlay].dx = 2.0*petrosian->radius;
     2085        overlay[Noverlay].dy = 2.0*petrosian->radius*profile->axes.minor/profile->axes.major;
     2086        overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD;
     2087        overlay[Noverlay].text = NULL;
     2088        Noverlay ++;
     2089        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     2090
     2091        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     2092        overlay[Noverlay].x = source->peak->xf;
     2093        overlay[Noverlay].y = source->peak->yf;
     2094        overlay[Noverlay].dx = 4.0*petrosian->radius;
     2095        overlay[Noverlay].dy = 4.0*petrosian->radius*profile->axes.minor/profile->axes.major;
     2096        overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD;
     2097        overlay[Noverlay].text = NULL;
     2098        Noverlay ++;
     2099        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     2100    }
     2101
     2102    KiiLoadOverlay (kapa, overlay, Noverlay, "red");
     2103    FREE (overlay);
     2104
     2105    pmVisualAskUser(NULL);
    18302106    return true;
    18312107}
    … …  
    18502126
    18512127# endif
     2128
     2129# if (0)
     2130// *** make a histogram of the source counts in the x and y directions
     2131psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0);
     2132psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0);
     2133psVectorHistogram (nX, xFaint, NULL, NULL, 0);
     2134psVectorHistogram (nY, yFaint, NULL, NULL, 0);
     2135psVector *dX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
     2136psVector *vX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
     2137psVector *dY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
     2138psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
     2139for (int i = 0; i < nX->nums->n; i++) {
     2140    dX->data.F32[i] = nX->nums->data.S32[i];
     2141    vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]);
     2142}
     2143for (int i = 0; i < nY->nums->n; i++) {
     2144    dY->data.F32[i] = nY->nums->data.S32[i];
     2145    vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]);
     2146}
     2147
     2148graphdata.color = KapaColorByName ("black");
     2149graphdata.ptype = 0;
     2150graphdata.size = 0.0;
     2151graphdata.style = 0;
     2152KapaPrepPlot (myKapa, dX->n, &graphdata);
     2153KapaPlotVector (myKapa, dX->n, dX->data.F32, "x");
     2154KapaPlotVector (myKapa, vX->n, vX->data.F32, "y");
     2155
     2156psFree (nX);
     2157psFree (dX);
     2158psFree (vX);
     2159
     2160psFree (nY);
     2161psFree (dY);
     2162psFree (vY);
     2163
     2164# endif
     2165
Note: See TracChangeset for help on using the changeset viewer.