IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27268


Ignore:
Timestamp:
Mar 12, 2010, 9:37:49 AM (16 years ago)
Author:
eugene
Message:

add psf model to radial flux plots

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20100225/psphot/src/psphotVisual.c

    r26894 r27268  
    180180bool psphotVisualShowImage (pmReadout *readout) {
    181181
    182     if (!pmVisualIsVisual()) return true;
     182    // if (!pmVisualIsVisual()) return true;
    183183
    184184    int kapa = psphotKapaChannel (1);
    … …  
    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);
    … …  
    14291551bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) {
    14301552
    1431     if (!pmVisualIsVisual()) return true;
     1553    // if (!pmVisualIsVisual()) return true;
    14321554
    14331555    int myKapa = psphotKapaChannel (1);
    … …  
    14561578    KapaSection section;
    14571579
    1458     if (!pmVisualIsVisual()) return true;
     1580    // if (!pmVisualIsVisual()) return true;
    14591581
    14601582    int myKapa = psphotKapaChannel (2);
    … …  
    15161638    psVector *sDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    15171639
     1640    psVector *xLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1641    psVector *yLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1642    psVector *mLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1643    psVector *sLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1644
    15181645    psVector *xCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    15191646    psVector *yCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    … …  
    15261653    int nPSF = 0;
    15271654    int nDEF = 0;
     1655    int nLOW = 0;
    15281656    int nCR  = 0;
    15291657    for (int i = 0; i < sources->n; i++) {
    15301658        pmSource *source = sources->data[i];
    15311659        if (source->moments == NULL) continue;
     1660
     1661        // only plot the measured sources...
     1662        if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) continue;
    15321663
    15331664        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) {
    … …  
    15611692            continue;
    15621693        }
    1563         if ((source->mode & PM_SOURCE_MODE_CR_LIMIT) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
     1694        if (source->errMag > 0.1) {
     1695            xLOW->data.F32[nLOW] = source->moments->Mxx;
     1696            yLOW->data.F32[nLOW] = source->moments->Myy;
     1697            mLOW->data.F32[nLOW] = -2.5*log10(source->moments->Sum);
     1698            sLOW->data.F32[nLOW] = source->extNsigma;
     1699            nLOW++;
    15641700            continue;
    15651701        }
    … …  
    15701706        nPSF++;
    15711707    }
     1708
    15721709    xSAT->n = nSAT;
    15731710    ySAT->n = nSAT;
    … …  
    15941731    mDEF->n = nDEF;
    15951732    sDEF->n = nDEF;
     1733
     1734    xLOW->n = nLOW;
     1735    yLOW->n = nLOW;
     1736    mLOW->n = nLOW;
     1737    sLOW->n = nLOW;
    15961738
    15971739    // four sections: MxxMyy, MagMxx, MagMyy, MagSigma
    … …  
    16571799    KapaPlotVector (myKapa, nSAT, ySAT->data.F32, "y");
    16581800
     1801    graphdata.color = KapaColorByName ("black");
     1802    graphdata.ptype = 7;
     1803    graphdata.size = 1.0;
     1804    graphdata.style = 2;
     1805    KapaPrepPlot   (myKapa, nLOW, &graphdata);
     1806    KapaPlotVector (myKapa, nLOW, xLOW->data.F32, "x");
     1807    KapaPlotVector (myKapa, nLOW, yLOW->data.F32, "y");
     1808
    16591809    // second section: MagMyy
    16601810    section.dx = 0.75;
    … …  
    16681818    graphdata.color = KapaColorByName ("black");
    16691819    graphdata.xmin = -17.1;
    1670     graphdata.xmax =  -7.9;
     1820    graphdata.xmax =  -6.9;
    16711821    graphdata.ymin = Ymin;
    16721822    graphdata.ymax = Ymax;
    … …  
    17301880    graphdata.xmin = Xmin;
    17311881    graphdata.xmax = Xmax;
    1732     graphdata.ymin =  -7.9;
     1882    graphdata.ymin =  -6.9;
    17331883    graphdata.ymax = -17.1;
    17341884    KapaSetLimits (myKapa, &graphdata);
    … …  
    17891939
    17901940    graphdata.color = KapaColorByName ("black");
    1791     graphdata.xmax =  -7.9;
     1941    graphdata.xmax =  -6.9;
    17921942    graphdata.xmin = -17.1;
    17931943    graphdata.ymin = -20.1;
Note: See TracChangeset for help on using the changeset viewer.