IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 2, 2010, 11:32:52 AM (16 years ago)
Author:
eugene
Message:

updates to psphotVisual mandated by kapa API changes; added psphotRadialBins; fixed up petrosian calculations (now correctly generate petRad, petMag, petR50, and petR90); updates to psphotPetrosianVisual to handle structure changes; option to use elliptical or circular profiles for radial flux and petrosians

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotPetrosianVisual.c

    r25755 r27819  
    11# include "psphotInternal.h"
     2# define FORCE_VISUAL 0
    23
    34// this function displays representative images as the psphot analysis progresses:
    … …  
    5455
    5556    // return true;
    56     if (!pmVisualIsVisual()) return true;
     57    if (!FORCE_VISUAL && !pmVisualIsVisual()) return true;
    5758
    5859    if (kapa2 == -1) {
    … …  
    101102
    102103    // return true;
    103     if (!pmVisualIsVisual()) return true;
     104    if (!FORCE_VISUAL && !pmVisualIsVisual()) return true;
    104105
    105106    if (kapa == -1) {
    … …  
    172173    KapaSection section;
    173174
    174     if (!pmVisualIsVisual()) return true;
     175    if (!FORCE_VISUAL && !pmVisualIsVisual()) return true;
     176
     177    if (kapa2 == -1) {
     178        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:stats");
     179        if (kapa2 == -1) {
     180            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     181            pmVisualSetVisual(false);
     182            return false;
     183        }
     184    }
     185
     186    KapaClearPlots (kapa2);
     187    KapaInitGraph (&graphdata);
     188    KapaSetFont (kapa2, "courier", 14);
     189
     190    // radius vs flux
     191    // radius vs mean SB
     192    // radius vs petRatio
     193
     194    // *** section 1: radius vs mean SB
     195    section.dx = 1.00;
     196    section.dy = 0.33;
     197    section.x  = 0.00;
     198    section.y  = 0.00;
     199    section.name = psStringCopy ("meanSB");
     200    KapaSetSection (kapa2, &section);
     201    psFree (section.name);
     202
     203    graphdata.color = KapaColorByName ("black");
     204    pmVisualLimitsFromVectors (&graphdata, radBin, fluxBin);
     205    KapaSetLimits (kapa2, &graphdata);
     206
     207    KapaBox (kapa2, &graphdata);
     208    KapaSendLabel (kapa2, "radius", KAPA_LABEL_XM);
     209    KapaSendLabel (kapa2, "mean SB", KAPA_LABEL_YM);
     210
     211    graphdata.color = KapaColorByName ("black");
     212    graphdata.style = 2;
     213    graphdata.ptype = 0;
     214    graphdata.size = 1.0;
     215    KapaPrepPlot (kapa2, radBin->n, &graphdata);
     216    KapaPlotVector (kapa2, radBin->n, radBin->data.F32, "x");
     217    KapaPlotVector (kapa2, radBin->n, fluxBin->data.F32, "y");
     218
     219    graphdata.color = KapaColorByName ("red");
     220    graphdata.style = 2;
     221    graphdata.ptype = 1;
     222    graphdata.size = 2.0;
     223    KapaPrepPlot (kapa2, refRadius->n, &graphdata);
     224    KapaPlotVector (kapa2, refRadius->n, refRadius->data.F32, "x");
     225    KapaPlotVector (kapa2, refRadius->n, meanSB->data.F32, "y");
     226
     227    // *** section 2: radius vs petrosian ratio
     228    section.dx = 1.00;
     229    section.dy = 0.33;
     230    section.x  = 0.00;
     231    section.y  = 0.33;
     232    section.name = psStringCopy ("ratio");
     233    KapaSetSection (kapa2, &section);
     234    psFree (section.name);
     235
     236    graphdata.color = KapaColorByName ("black");
     237    graphdata.ymax = +1.05;
     238    graphdata.ymin = -0.05;
     239    pmVisualLimitsFromVectors (&graphdata, radBin, NULL);
     240    KapaSetLimits (kapa2, &graphdata);
     241
     242    KapaBox (kapa2, &graphdata);
     243    KapaSendLabel (kapa2, "radius", KAPA_LABEL_XM);
     244    KapaSendLabel (kapa2, "ratio", KAPA_LABEL_YM);
     245
     246    graphdata.color = KapaColorByName ("black");
     247    graphdata.style = 2;
     248    graphdata.ptype = 0;
     249    graphdata.size = 1.0;
     250    graphdata.etype = 0x01;
     251    KapaPrepPlot (kapa2, refRadius->n, &graphdata);
     252    KapaPlotVector (kapa2, refRadius->n, refRadius->data.F32, "x");
     253    KapaPlotVector (kapa2, refRadius->n, petRatio->data.F32, "y");
     254    KapaPlotVector (kapa2, refRadius->n, petRatioErr->data.F32, "dym");
     255    KapaPlotVector (kapa2, refRadius->n, petRatioErr->data.F32, "dyp");
     256    graphdata.etype = 0;
     257
     258    graphdata.color = KapaColorByName ("red");
     259    graphdata.style = 2;
     260    graphdata.ptype = 2;
     261    graphdata.size = 2.0;
     262    KapaPrepPlot   (kapa2, 1, &graphdata);
     263    KapaPlotVector (kapa2, 1, &petRadius, "x");
     264    KapaPlotVector (kapa2, 1, &ratioForRadius, "y");
     265
     266    // *** section 3: radius vs integrated flux
     267    section.dx = 1.00;
     268    section.dy = 0.33;
     269    section.x  = 0.00;
     270    section.y  = 0.66;
     271    section.name = psStringCopy ("flux");
     272    KapaSetSection (kapa2, &section);
     273    psFree (section.name);
     274
     275    graphdata.color = KapaColorByName ("black");
     276    pmVisualLimitsFromVectors (&graphdata, radBin, fluxSum);
     277    KapaSetLimits (kapa2, &graphdata);
     278
     279    KapaBox (kapa2, &graphdata);
     280    KapaSendLabel (kapa2, "radius", KAPA_LABEL_XM);
     281    KapaSendLabel (kapa2, "integrated flux", KAPA_LABEL_YM);
     282
     283    graphdata.color = KapaColorByName ("black");
     284    graphdata.style = 2;
     285    graphdata.ptype = 0;
     286    graphdata.size = 1.0;
     287    KapaPrepPlot   (kapa2, refRadius->n, &graphdata);
     288    KapaPlotVector (kapa2, refRadius->n, refRadius->data.F32, "x");
     289    KapaPlotVector (kapa2, refRadius->n, fluxSum->data.F32, "y");
     290
     291    graphdata.color = KapaColorByName ("red");
     292    graphdata.ptype = 2;
     293    graphdata.style = 2;
     294    graphdata.size = 2.0;
     295    KapaPrepPlot   (kapa2, 1, &graphdata);
     296    KapaPlotVector (kapa2, 1, &radiusForFlux, "x");
     297    KapaPlotVector (kapa2, 1, &petFlux, "y");
     298
     299    // pause and wait for user input:
     300    // continue, save (provide name), ??
     301    char key[10];
     302    fprintf (stdout, "[c]ontinue? ");
     303    if (!fgets(key, 8, stdin)) {
     304        psWarning("Unable to read option");
     305    }
     306    return true;
     307}
     308
     309bool psphotPetrosianVisualEllipticalContour (pmSourceRadialFlux *radFlux, pmSourceExtendedPars *extpars) {
     310
     311    Graphdata graphdata;
     312
     313    if (!FORCE_VISUAL && !pmVisualIsVisual()) return true;
    175314
    176315    if (kapa == -1) {
    … …  
    187326    KapaSetFont (kapa, "courier", 14);
    188327
    189     // radius vs flux
    190     // radius vs mean SB
    191     // radius vs petRatio
    192 
    193     // *** section 1: radius vs mean SB
    194     section.dx = 1.00;
    195     section.dy = 0.33;
    196     section.x  = 0.00;
    197     section.y  = 0.00;
    198     section.name = psStringCopy ("meanSB");
    199     KapaSetSection (kapa, &section);
    200     psFree (section.name);
    201 
    202     graphdata.color = KapaColorByName ("black");
    203     pmVisualLimitsFromVectors (&graphdata, radBin, fluxBin);
     328    psVector *theta = radFlux->theta;
     329    psVector *radius = radFlux->isophotalRadii;
     330
     331    // find Rmin and Rmax for the initial guess
     332    float Rmin = radius->data.F32[0];
     333    float Rmax = radius->data.F32[0];
     334
     335    psVector *Rx = psVectorAlloc(radius->n, PS_TYPE_F32);
     336    psVector *Ry = psVectorAlloc(radius->n, PS_TYPE_F32);
     337
     338    for (int i = 0; i < theta->n; i++) {
     339        Rx->data.F32[i] = radius->data.F32[i]*cos(theta->data.F32[i]);
     340        Ry->data.F32[i] = radius->data.F32[i]*sin(theta->data.F32[i]);
     341
     342        // check the radius range
     343        Rmin = MIN (Rmin, radius->data.F32[i]);
     344        Rmax = MAX (Rmax, radius->data.F32[i]);
     345    }   
     346
     347    psVector *rx = psVectorAlloc(361, PS_TYPE_F32);
     348    psVector *ry = psVectorAlloc(361, PS_TYPE_F32);
     349
     350    float epsilon = extpars->axes.minor / extpars->axes.major;
     351
     352    for (int i = 0; i < 361; i++) {
     353
     354        float alpha = PS_RAD_DEG * i;
     355
     356        float cs_alpha = cos(alpha);
     357        float sn_alpha = sin(alpha);
     358
     359        float cs_phi = cos(alpha - extpars->axes.theta);
     360        float sn_phi = sin(alpha - extpars->axes.theta);
     361
     362        float r = 1.0 / sqrt(SQ(sn_phi) + SQ(epsilon*cs_phi));
     363
     364        // generate the model fit here
     365        rx->data.F32[i] = extpars->axes.minor * cs_alpha * r;
     366        ry->data.F32[i] = extpars->axes.minor * sn_alpha * r;
     367    }   
     368
     369    graphdata.color = KapaColorByName ("black");
     370    graphdata.xmin = -1.1*Rmax;
     371    graphdata.ymin = -1.1*Rmax;
     372    graphdata.xmax = +1.1*Rmax;
     373    graphdata.ymax = +1.1*Rmax;
    204374    KapaSetLimits (kapa, &graphdata);
    205375
    206376    KapaBox (kapa, &graphdata);
    207     KapaSendLabel (kapa, "radius", KAPA_LABEL_XM);
    208     KapaSendLabel (kapa, "mean SB", KAPA_LABEL_YM);
    209 
    210     graphdata.color = KapaColorByName ("black");
    211     graphdata.style = 2;
    212     graphdata.ptype = 0;
    213     graphdata.size = 1.0;
    214     KapaPrepPlot (kapa, radBin->n, &graphdata);
    215     KapaPlotVector (kapa, radBin->n, radBin->data.F32, "x");
    216     KapaPlotVector (kapa, radBin->n, fluxBin->data.F32, "y");
     377    KapaSendLabel (kapa, "R_x", KAPA_LABEL_XM);
     378    KapaSendLabel (kapa, "R_y", KAPA_LABEL_YM);
    217379
    218380    graphdata.color = KapaColorByName ("red");
    219381    graphdata.style = 2;
    220     graphdata.ptype = 1;
    221     graphdata.size = 2.0;
    222     KapaPrepPlot (kapa, refRadius->n, &graphdata);
    223     KapaPlotVector (kapa, refRadius->n, refRadius->data.F32, "x");
    224     KapaPlotVector (kapa, refRadius->n, meanSB->data.F32, "y");
    225 
    226     // *** section 2: radius vs petrosian ratio
    227     section.dx = 1.00;
    228     section.dy = 0.33;
    229     section.x  = 0.00;
    230     section.y  = 0.33;
    231     section.name = psStringCopy ("ratio");
    232     KapaSetSection (kapa, &section);
    233     psFree (section.name);
    234 
    235     graphdata.color = KapaColorByName ("black");
    236     graphdata.ymax = +1.05;
    237     graphdata.ymin = -0.05;
    238     pmVisualLimitsFromVectors (&graphdata, radBin, NULL);
    239     KapaSetLimits (kapa, &graphdata);
    240 
    241     KapaBox (kapa, &graphdata);
    242     KapaSendLabel (kapa, "radius", KAPA_LABEL_XM);
    243     KapaSendLabel (kapa, "ratio", KAPA_LABEL_YM);
    244 
    245     graphdata.color = KapaColorByName ("black");
    246     graphdata.style = 2;
    247     graphdata.ptype = 0;
    248     graphdata.size = 1.0;
    249     graphdata.etype = 0x01;
    250     KapaPrepPlot (kapa, refRadius->n, &graphdata);
    251     KapaPlotVector (kapa, refRadius->n, refRadius->data.F32, "x");
    252     KapaPlotVector (kapa, refRadius->n, petRatio->data.F32, "y");
    253     KapaPlotVector (kapa, refRadius->n, petRatioErr->data.F32, "dym");
    254     KapaPlotVector (kapa, refRadius->n, petRatioErr->data.F32, "dyp");
    255     graphdata.etype = 0;
    256 
    257     graphdata.color = KapaColorByName ("red");
    258     graphdata.style = 2;
    259382    graphdata.ptype = 2;
    260     graphdata.size = 2.0;
    261     KapaPrepPlot   (kapa, 1, &graphdata);
    262     KapaPlotVector (kapa, 1, &petRadius, "x");
    263     KapaPlotVector (kapa, 1, &ratioForRadius, "y");
    264 
    265     // *** section 3: radius vs integrated flux
    266     section.dx = 1.00;
    267     section.dy = 0.33;
    268     section.x  = 0.00;
    269     section.y  = 0.66;
    270     section.name = psStringCopy ("flux");
    271     KapaSetSection (kapa, &section);
    272     psFree (section.name);
    273 
    274     graphdata.color = KapaColorByName ("black");
    275     pmVisualLimitsFromVectors (&graphdata, radBin, fluxSum);
    276     KapaSetLimits (kapa, &graphdata);
    277 
    278     KapaBox (kapa, &graphdata);
    279     KapaSendLabel (kapa, "radius", KAPA_LABEL_XM);
    280     KapaSendLabel (kapa, "integrated flux", KAPA_LABEL_YM);
    281 
    282     graphdata.color = KapaColorByName ("black");
    283     graphdata.style = 2;
    284     graphdata.ptype = 0;
    285     graphdata.size = 1.0;
    286     KapaPrepPlot   (kapa, refRadius->n, &graphdata);
    287     KapaPlotVector (kapa, refRadius->n, refRadius->data.F32, "x");
    288     KapaPlotVector (kapa, refRadius->n, fluxSum->data.F32, "y");
    289 
    290     graphdata.color = KapaColorByName ("red");
    291     graphdata.ptype = 2;
    292     graphdata.style = 2;
    293     graphdata.size = 2.0;
    294     KapaPrepPlot   (kapa, 1, &graphdata);
    295     KapaPlotVector (kapa, 1, &radiusForFlux, "x");
    296     KapaPlotVector (kapa, 1, &petFlux, "y");
     383    graphdata.size = 1.0;
     384    KapaPrepPlot (kapa, Rx->n, &graphdata);
     385    KapaPlotVector (kapa, Rx->n, Rx->data.F32, "x");
     386    KapaPlotVector (kapa, Rx->n, Ry->data.F32, "y");
     387
     388    graphdata.color = KapaColorByName ("black");
     389    graphdata.style = 0;
     390    graphdata.ptype = 0;
     391    graphdata.size = 1.0;
     392    KapaPrepPlot (kapa, rx->n, &graphdata);
     393    KapaPlotVector (kapa, rx->n, rx->data.F32, "x");
     394    KapaPlotVector (kapa, rx->n, ry->data.F32, "y");
    297395
    298396    // pause and wait for user input:
    … …  
    306404}
    307405
    308 bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian) {
    309 
    310     Graphdata graphdata;
    311 
    312     if (!pmVisualIsVisual()) return true;
    313 
    314     if (kapa == -1) {
    315         kapa = KapaOpenNamedSocket ("kapa", "psphot:plots");
    316         if (kapa == -1) {
    317             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    318             pmVisualSetVisual(false);
    319             return false;
    320         }
    321     }
    322 
    323     KapaClearPlots (kapa);
    324     KapaInitGraph (&graphdata);
    325     KapaSetFont (kapa, "courier", 14);
    326 
    327     psVector *theta = petrosian->theta;
    328     psVector *radius = petrosian->isophotalRadii;
    329 
    330     // find Rmin and Rmax for the initial guess
    331     float Rmin = radius->data.F32[0];
    332     float Rmax = radius->data.F32[0];
    333 
    334     psVector *Rx = psVectorAlloc(radius->n, PS_TYPE_F32);
    335     psVector *Ry = psVectorAlloc(radius->n, PS_TYPE_F32);
    336 
    337     for (int i = 0; i < theta->n; i++) {
    338         Rx->data.F32[i] = radius->data.F32[i]*cos(theta->data.F32[i]);
    339         Ry->data.F32[i] = radius->data.F32[i]*sin(theta->data.F32[i]);
    340 
    341         // check the radius range
    342         Rmin = MIN (Rmin, radius->data.F32[i]);
    343         Rmax = MAX (Rmax, radius->data.F32[i]);
    344     }   
    345 
    346     psVector *rx = psVectorAlloc(361, PS_TYPE_F32);
    347     psVector *ry = psVectorAlloc(361, PS_TYPE_F32);
    348 
    349     float epsilon = petrosian->axes.minor / petrosian->axes.major;
    350 
    351     for (int i = 0; i < 361; i++) {
    352 
    353         float alpha = PS_RAD_DEG * i;
    354 
    355         float cs_alpha = cos(alpha);
    356         float sn_alpha = sin(alpha);
    357 
    358         float cs_phi = cos(alpha - petrosian->axes.theta);
    359         float sn_phi = sin(alpha - petrosian->axes.theta);
    360 
    361         float r = 1.0 / sqrt(SQ(sn_phi) + SQ(epsilon*cs_phi));
    362 
    363         // generate the model fit here
    364         rx->data.F32[i] = petrosian->axes.minor * cs_alpha * r;
    365         ry->data.F32[i] = petrosian->axes.minor * sn_alpha * r;
    366     }   
    367 
    368     graphdata.color = KapaColorByName ("black");
    369     graphdata.xmin = -1.1*Rmax;
    370     graphdata.ymin = -1.1*Rmax;
    371     graphdata.xmax = +1.1*Rmax;
    372     graphdata.ymax = +1.1*Rmax;
    373     KapaSetLimits (kapa, &graphdata);
    374 
    375     KapaBox (kapa, &graphdata);
    376     KapaSendLabel (kapa, "R_x", KAPA_LABEL_XM);
    377     KapaSendLabel (kapa, "R_y", KAPA_LABEL_YM);
    378 
    379     graphdata.color = KapaColorByName ("red");
    380     graphdata.style = 2;
    381     graphdata.ptype = 2;
    382     graphdata.size = 1.0;
    383     KapaPrepPlot (kapa, Rx->n, &graphdata);
    384     KapaPlotVector (kapa, Rx->n, Rx->data.F32, "x");
    385     KapaPlotVector (kapa, Rx->n, Ry->data.F32, "y");
    386 
    387     graphdata.color = KapaColorByName ("black");
    388     graphdata.style = 0;
    389     graphdata.ptype = 0;
    390     graphdata.size = 1.0;
    391     KapaPrepPlot (kapa, rx->n, &graphdata);
    392     KapaPlotVector (kapa, rx->n, rx->data.F32, "x");
    393     KapaPlotVector (kapa, rx->n, ry->data.F32, "y");
    394 
    395     // pause and wait for user input:
    396     // continue, save (provide name), ??
    397     char key[10];
    398     fprintf (stdout, "[c]ontinue? ");
    399     if (!fgets(key, 8, stdin)) {
    400         psWarning("Unable to read option");
    401     }
    402     return true;
    403 }
    404 
    405406# endif
Note: See TracChangeset for help on using the changeset viewer.