IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 30, 2010, 1:33:47 PM (16 years ago)
Author:
eugene
Message:

some API re-work to support psphotForced and detection efficiency measurements; better visualization; CR masking (optional); better source size analysis; some motion towards stack photometry

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotVisual.c

    r26894 r27532  
    10991099    psFree (outsat);
    11001100    return true;
     1101}
     1102
     1103static void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1)
     1104{
     1105    float x[2], y[2];
     1106    x[0] = x0;
     1107    x[1] = x1;
     1108    y[0] = y0;
     1109    y[1] = y1;
     1110    KapaPrepPlot   (myKapa, 2, graphdata);
     1111    KapaPlotVector (myKapa, 2, x, "x");
     1112    KapaPlotVector (myKapa, 2, y, "y");
    11011113}
    11021114
    … …  
    11381150    }
    11391151
     1152    // generate model profiles (major and minor axis):
     1153    // create a model with theta = 0.0 so major and minor axes are equiv to x and y:
     1154    psEllipseShape rawShape, rotShape;
     1155
     1156    rawShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
     1157    rawShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
     1158    rawShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
     1159
     1160    psEllipseAxes axes = psEllipseShapeToAxes (rawShape, 20.0);
     1161
     1162    axes.theta = 0.0;
     1163
     1164    rotShape = psEllipseAxesToShape (axes);
     1165
     1166    psVector *params = psVectorAlloc(source->modelPSF->params->n, PS_TYPE_F32);
     1167    for (int i = 0; i < source->modelPSF->params->n; i++) {
     1168        params->data.F32[i] = source->modelPSF->params->data.F32[i];
     1169    }
     1170    params->data.F32[PM_PAR_SXX] = rotShape.sx * M_SQRT2;
     1171    params->data.F32[PM_PAR_SYY] = rotShape.sy * M_SQRT2;
     1172    params->data.F32[PM_PAR_SXY] = rotShape.sxy;
     1173    params->data.F32[PM_PAR_XPOS] = 0.0;
     1174    params->data.F32[PM_PAR_YPOS] = 0.0;
     1175
     1176    psVector *rmod = psVectorAlloc(300, PS_TYPE_F32);
     1177    psVector *fmaj = psVectorAlloc(300, PS_TYPE_F32);
     1178    psVector *fmin = psVectorAlloc(300, PS_TYPE_F32);
     1179
     1180    psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     1181
     1182    float r = 0.0;
     1183    for (int i = 0; i < rmod->n; i++) {
     1184        r = i*0.1;
     1185        rmod->data.F32[i] = r;
     1186
     1187        coord->data.F32[1] = r;
     1188        coord->data.F32[0] = 0.0;
     1189        fmaj->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
     1190
     1191        coord->data.F32[0] = r;
     1192        coord->data.F32[1] = 0.0;
     1193        fmin->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
     1194    }
     1195    psFree (coord);
     1196    psFree (params);
     1197
     1198    float FWHM_MAJOR = 2.0*source->modelPSF->modelRadius (source->modelPSF->params, 0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
     1199    float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
     1200    if (FWHM_MAJOR < FWHM_MINOR) PS_SWAP (FWHM_MAJOR, FWHM_MINOR);
     1201
     1202    psEllipseMoments emoments;
     1203    emoments.x2 = source->moments->Mxx;
     1204    emoments.xy = source->moments->Mxy;
     1205    emoments.y2 = source->moments->Myy;
     1206    axes = psEllipseMomentsToAxes (emoments, 20.0);
     1207    float MOMENTS_MAJOR = 2.355*axes.major;
     1208    float MOMENTS_MINOR = 2.355*axes.minor;
     1209
     1210    float logHM = log10(0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
     1211
    11401212    // reset source Add/Sub state to recorded
    11411213    if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    … …  
    11741246    KapaPlotVector (myKapa, nb, fb->data.F32, "y");
    11751247
     1248    graphdata.color = KapaColorByName ("blue");
     1249    graphdata.ptype = 0;
     1250    graphdata.size = 0.0;
     1251    graphdata.style = 0;
     1252    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1253    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1254    KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
     1255    plotline (myKapa, &graphdata, 0.0, logHM, 30.0, logHM);
     1256    plotline (myKapa, &graphdata, 0.5*FWHM_MINOR, 0.0, 0.5*FWHM_MINOR, 5.0);
     1257    graphdata.ltype = 1;
     1258    plotline (myKapa, &graphdata, 0.5*MOMENTS_MINOR, 0.0, 0.5*MOMENTS_MINOR, 5.0);
     1259    graphdata.ltype = 0;
     1260       
     1261    graphdata.color = KapaColorByName ("green");
     1262    graphdata.ptype = 0;
     1263    graphdata.size = 0.0;
     1264    graphdata.style = 0;
     1265    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1266    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1267    KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
     1268    plotline (myKapa, &graphdata, 0.5*FWHM_MAJOR, 0.0, 0.5*FWHM_MAJOR, 5.0);
     1269    graphdata.ltype = 1;
     1270    plotline (myKapa, &graphdata, 0.5*MOMENTS_MAJOR, 0.0, 0.5*MOMENTS_MAJOR, 5.0);
     1271    graphdata.ltype = 0;
     1272       
     1273    for (int i = 0; i < rmod->n; i++) {
     1274        rmod->data.F32[i] = log10(rmod->data.F32[i]);
     1275    }
     1276
    11761277    // ** loglog **
    11771278    KapaSelectSection (myKapa, "loglog");
    … …  
    11821283    graphdata.ymin = -0.05;
    11831284    graphdata.ymax = +5.05;
     1285    graphdata.color = KapaColorByName ("black");
    11841286    KapaSetLimits (myKapa, &graphdata);
    11851287
    … …  
    12041306    KapaPlotVector (myKapa, nb, Rb->data.F32, "x");
    12051307    KapaPlotVector (myKapa, nb, fb->data.F32, "y");
     1308
     1309    graphdata.color = KapaColorByName ("blue");
     1310    graphdata.ptype = 0;
     1311    graphdata.size = 0.0;
     1312    graphdata.style = 0;
     1313    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1314    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1315    KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
     1316
     1317    graphdata.color = KapaColorByName ("green");
     1318    graphdata.ptype = 0;
     1319    graphdata.size = 0.0;
     1320    graphdata.style = 0;
     1321    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1322    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1323    KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
     1324
     1325    psFree (rmod);
     1326    psFree (fmin);
     1327    psFree (fmaj);
    12061328
    12071329    psFree (rg);
    … …  
    13921514        if (source == NULL) continue;
    13931515
    1394         // if (source->type != type) continue;
    13951516        if (mode) {
    13961517            if (keep) {
    … …  
    15161637    psVector *sDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    15171638
     1639    psVector *xLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1640    psVector *yLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1641    psVector *mLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1642    psVector *sLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1643
    15181644    psVector *xCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    15191645    psVector *yCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    … …  
    15261652    int nPSF = 0;
    15271653    int nDEF = 0;
     1654    int nLOW = 0;
    15281655    int nCR  = 0;
    15291656    for (int i = 0; i < sources->n; i++) {
    15301657        pmSource *source = sources->data[i];
    15311658        if (source->moments == NULL) continue;
     1659
     1660        // only plot the measured sources...
     1661        if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) continue;
    15321662
    15331663        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) {
    … …  
    15611691            continue;
    15621692        }
    1563         if ((source->mode & PM_SOURCE_MODE_CR_LIMIT) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
     1693        if (source->errMag > 0.1) {
     1694            xLOW->data.F32[nLOW] = source->moments->Mxx;
     1695            yLOW->data.F32[nLOW] = source->moments->Myy;
     1696            mLOW->data.F32[nLOW] = -2.5*log10(source->moments->Sum);
     1697            sLOW->data.F32[nLOW] = source->extNsigma;
     1698            nLOW++;
    15641699            continue;
    15651700        }
    … …  
    15701705        nPSF++;
    15711706    }
     1707
    15721708    xSAT->n = nSAT;
    15731709    ySAT->n = nSAT;
    … …  
    15941730    mDEF->n = nDEF;
    15951731    sDEF->n = nDEF;
     1732
     1733    xLOW->n = nLOW;
     1734    yLOW->n = nLOW;
     1735    mLOW->n = nLOW;
     1736    sLOW->n = nLOW;
    15961737
    15971738    // four sections: MxxMyy, MagMxx, MagMyy, MagSigma
    … …  
    16571798    KapaPlotVector (myKapa, nSAT, ySAT->data.F32, "y");
    16581799
     1800    graphdata.color = KapaColorByName ("black");
     1801    graphdata.ptype = 7;
     1802    graphdata.size = 1.0;
     1803    graphdata.style = 2;
     1804    KapaPrepPlot   (myKapa, nLOW, &graphdata);
     1805    KapaPlotVector (myKapa, nLOW, xLOW->data.F32, "x");
     1806    KapaPlotVector (myKapa, nLOW, yLOW->data.F32, "y");
     1807
    16591808    // second section: MagMyy
    16601809    section.dx = 0.75;
    … …  
    16681817    graphdata.color = KapaColorByName ("black");
    16691818    graphdata.xmin = -17.1;
    1670     graphdata.xmax =  -7.9;
     1819    graphdata.xmax =  -6.9;
    16711820    graphdata.ymin = Ymin;
    16721821    graphdata.ymax = Ymax;
    … …  
    17301879    graphdata.xmin = Xmin;
    17311880    graphdata.xmax = Xmax;
    1732     graphdata.ymin =  -7.9;
     1881    graphdata.ymin =  -6.9;
    17331882    graphdata.ymax = -17.1;
    17341883    KapaSetLimits (myKapa, &graphdata);
    … …  
    17891938
    17901939    graphdata.color = KapaColorByName ("black");
    1791     graphdata.xmax =  -7.9;
     1940    graphdata.xmax =  -6.9;
    17921941    graphdata.xmin = -17.1;
    17931942    graphdata.ymin = -20.1;
Note: See TracChangeset for help on using the changeset viewer.