IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27819


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

Location:
trunk/psphot/src
Files:
1 added
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/Makefile.am

    r27660 r27819  
    185185        psphotEllipticalContour.c      \
    186186        psphotEllipticalProfile.c      \
     187        psphotRadialBins.c             \
    187188        psphotPetrosian.c              \
    188189        psphotPetrosianRadialBins.c    \
    189190        psphotPetrosianStats.c         \
     191        psphotPetrosianVisual.c        \
    190192        psphotEfficiency.c
    191193
    … …  
    197199#       psphotAnnuli.c                 \
    198200#       psphotKron.c                   \
    199 #       psphotPetrosianVisual.c        \
    200201#
    201202

  • trunk/psphot/src/psphot.h

    r27673 r27819  
    245245float           psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax);
    246246bool            psphotRadiiFromProfiles (pmSource *source, float fluxMin, float fluxMax);
    247 bool            psphotEllipticalProfile (pmSource *source);
     247bool            psphotEllipticalProfile (pmSource *source, bool RAW_RADIUS);
    248248bool            psphotEllipticalContour (pmSource *source);
    249249
    … …  
    281281
    282282// XXX visualization functions related to radial profiles (disabled)
    283 // bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux);
    284 // bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef);
    285 // bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian);
    286 // bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin,
    287 //                               psVector *refRadius, psVector *meanSB,
    288 //                               psVector *petRatio, psVector *petRatioErr, psVector *fluxSum,
    289 //                               float petRadius, float ratioForRadius,
    290 //                               float petFlux, float radiusForFlux);
     283bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux);
     284bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef);
     285bool psphotPetrosianVisualEllipticalContour (pmSourceRadialFlux *radFlux, pmSourceExtendedPars *extpars);
     286bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin,
     287                               psVector *refRadius, psVector *meanSB,
     288                               psVector *petRatio, psVector *petRatioErr, psVector *fluxSum,
     289                               float petRadius, float ratioForRadius,
     290                               float petFlux, float radiusForFlux);
     291
     292bool psphotRadialBins (psMetadata *recipe, pmSource *source, float radiusMax, float skynoise);
    291293
    292294// structures & functions to support psf-convolved model fitting

  • trunk/psphot/src/psphotEllipticalContour.c

    r25755 r27819  
    77bool psphotEllipticalContour (pmSource *source) {
    88
    9     pmSourceRadialProfile *profile = source->extpars->profile;
     9    psAssert (source, "missing source");
     10    psAssert (source->extpars, "missing extpars");
     11    psAssert (source->extpars->radFlux, "missing radFlux");
     12
     13    pmSourceRadialFlux *profile = source->extpars->radFlux;
     14    pmSourceExtendedPars *extpars = source->extpars;
    1015
    1116    // use LMM to fit theta vs radius to an ellipse
    … …  
    8590    /// XXX rationalize? if epsilon > 1, flip major and minor axes (rotate by 90 degrees)
    8691    if (params->data.F32[PAR_EPSILON] < 1.0) {
    87         profile->axes.major = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON];
    88         profile->axes.minor = params->data.F32[PAR_RMIN];
    89         profile->axes.theta = params->data.F32[PAR_PHI];
     92        extpars->axes.major = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON];
     93        extpars->axes.minor = params->data.F32[PAR_RMIN];
     94        extpars->axes.theta = params->data.F32[PAR_PHI];
    9095    } else {
    91         profile->axes.major = params->data.F32[PAR_RMIN];
    92         profile->axes.minor = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON];
    93         profile->axes.theta = params->data.F32[PAR_PHI] + 0.5*M_PI;
     96        extpars->axes.major = params->data.F32[PAR_RMIN];
     97        extpars->axes.minor = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON];
     98        extpars->axes.theta = params->data.F32[PAR_PHI] + 0.5*M_PI;
    9499    }
    95100
    96101    psTrace ("psphot", 4, "# fitted values:\n");
    97     psTrace ("psphot", 4, "Phi:   %f\n", profile->axes.theta*PS_DEG_RAD);
    98     psTrace ("psphot", 4, "Rmaj:  %f\n", profile->axes.major);
    99     psTrace ("psphot", 4, "Rmin:  %f\n", profile->axes.minor);
     102    psTrace ("psphot", 4, "Phi:   %f\n", extpars->axes.theta*PS_DEG_RAD);
     103    psTrace ("psphot", 4, "Rmaj:  %f\n", extpars->axes.major);
     104    psTrace ("psphot", 4, "Rmin:  %f\n", extpars->axes.minor);
    100105   
    101106    // show the results
    102     // psphotPetrosianVisualEllipticalContour (petrosian);
     107    // psphotPetrosianVisualEllipticalContour (profile, extpars);
    103108
    104109    psFree (x);

  • trunk/psphot/src/psphotEllipticalProfile.c

    r25755 r27819  
    11# include "psphotInternal.h"
    22
    3 bool psphotEllipticalProfile (pmSource *source) {
     3bool psphotEllipticalProfile (pmSource *source, bool RAW_RADIUS) {
    44
    5     pmSourceRadialProfile *profile = source->extpars->profile;
     5    psAssert (source, "missing source");
     6    psAssert (source->extpars, "missing extpars");
     7    psAssert (source->pixels, "missing pixels");
     8
     9    pmSourceExtendedPars *extpars = source->extpars;
     10
     11    if (!source->extpars->ellipticalFlux) {
     12        source->extpars->ellipticalFlux = pmSourceEllipticalFluxAlloc();
     13    }
     14    pmSourceEllipticalFlux *profile = source->extpars->ellipticalFlux;
    615
    716    profile->radiusElliptical = psVectorAllocEmpty(100, PS_TYPE_F32);
    … …  
    2130
    2231    psEllipseAxes axes;
    23     axes.major = M_SQRT1_2;
    24     axes.minor = M_SQRT1_2 * (profile->axes.minor / profile->axes.major);
     32    if (RAW_RADIUS) {
     33        // force circular profile
     34        axes.major = M_SQRT1_2;
     35        axes.minor = M_SQRT1_2;
     36    } else {
     37        axes.major = M_SQRT1_2;
     38        axes.minor = M_SQRT1_2 * (extpars->axes.minor / extpars->axes.major);
     39    }
    2540
    2641    // axes.major = 1.0;
    27     // axes.minor = profile->axes.minor / profile->axes.major;
     42    // axes.minor = extpars->axes.minor / extpars->axes.major;
    2843
    29     axes.theta = profile->axes.theta;
     44    axes.theta = extpars->axes.theta;
    3045    psEllipseShape shape = psEllipseAxesToShape (axes);
    3146
    … …  
    4661
    4762            float r2 = 0.5*PS_SQR(x/Sxx) + 0.5*PS_SQR(y/Syy) + x*y*Sxy;
     63            float Rraw = hypot(x, y);
    4864
    4965            psVectorAppend(radius, sqrt(r2));
    5066            psVectorAppend(flux, source->pixels->data.F32[iy][ix]);
    5167
    52             float Rraw = hypot(x, y);
    5368            psVectorAppend(radiusRaw, Rraw);
    5469            psVectorAppend(fluxRaw, source->pixels->data.F32[iy][ix]);
    … …  
    6782    // }
    6883
    69     // psphotPetrosianVisualProfileRadii (radius, flux, radiusRaw, fluxRaw, 0.0);
     84    psphotPetrosianVisualProfileRadii (radius, flux, radiusRaw, fluxRaw, source->peak->flux, 0.0);
    7085    // psphotPetrosianVisualProfileByAngle (radius, flux);
    7186

  • trunk/psphot/src/psphotExtendedSourceAnalysis.c

    r26894 r27819  
    3939    int Nkron = 0;
    4040
     41    psTimerStart ("psphot.extended");
     42
    4143    // find the currently selected readout
    4244    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
    … …  
    6668    float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
    6769
    68 # if (0)
    69     // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
    70     // XXX use this to set skynoise
    71     pmReadout *backModel = psphotSelectBackground (config, view);
    72     pmReadout *backStdev = psphotSelectBackgroundStdev (config, view);
    73 # endif
    74 
    7570    // S/N limit to perform full non-linear fits
    7671    float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
    … …  
    8176    bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    8277    bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
     78
     79# if (0)
     80    // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
     81    // XXX use this to set skynoise
     82    pmReadout *backModel = psphotSelectBackground (config, view);
     83    pmReadout *backStdev = psphotSelectBackgroundStdev (config, view);
     84# endif
    8385
    8486    // source analysis is done in S/N order (brightest first)
    … …  
    119121        Next ++;
    120122
     123        // force source image to be a bit larger...
     124        float radius = source->peak->xf - source->pixels->col0;
     125        radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
     126        radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
     127        radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
     128        pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);
     129
    121130        // if we request any of these measurements, we require the radial profile
    122131        if (doPetrosian || doIsophotal || doAnnuli || doKron) {
    … …  
    134143        if (doPetrosian) {
    135144            if (!psphotPetrosian (source, recipe, skynoise, maskVal)) {
    136                 psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
     145                psTrace ("psphot", 5, "FAILED petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    137146            } else {
    138147                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    … …  
    169178
    170179        if (source->extpars) {
    171             pmSourceRadialProfileFreeVectors(source->extpars->profile);
     180            psFree(source->extpars->radFlux);
     181            psFree(source->extpars->ellipticalFlux);
    172182        }
    173183    }
    174184
    175     psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), Next);
     185    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
    176186    psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
    177187    psLogMsg ("psphot", PS_LOG_INFO, "  %d isophotal\n", Nisophot);

  • trunk/psphot/src/psphotPetrosian.c

    r25755 r27819  
    77
    88    psAssert (source->extpars, "need to run psphotRadialProfile first");
    9     psAssert (source->extpars->profile, "need to run psphotRadialProfile first");
     9    psAssert (source->extpars->ellipticalFlux, "need to run psphotRadialProfile first");
    1010
    1111    // integrate the radial profile for radial bins defined for the petrosian measurement:
    … …  
    2424    psTrace ("psphot", 3, "source at %f,%f: petrosian radius: %f, flux: %f, axis ratio: %f, angle: %f",
    2525             source->peak->xf, source->peak->yf,
    26              source->extpars->petrosian_80->radius,
    27              source->extpars->petrosian_80->flux,
    28              source->extpars->profile->axes.minor/source->extpars->profile->axes.major,
    29              source->extpars->profile->axes.theta*PS_DEG_RAD);
     26             source->extpars->petrosianRadius,
     27             source->extpars->petrosianFlux,
     28             source->extpars->axes.minor/source->extpars->axes.major,
     29             source->extpars->axes.theta*PS_DEG_RAD);
    3030
    3131    return true;

  • trunk/psphot/src/psphotPetrosianRadialBins.c

    r25755 r27819  
    11# include "psphotInternal.h"
     2float InterpolateValues (float X0, float Y0, float X1, float Y1, float X);
    23
    34// convert the flux vs elliptical radius to annular bins
    … …  
    56// we are guaranteed to be limited by either the seeing (1 - few pixels) or by the pixels
    67// themselves.  this function does not attempt to measure the radial profiles accurately
    7 // for radii that are smaller than a minimum (currently 1.0 pixels). 
     8// for radii that are smaller than ~2 pixels
    89
    910// for small radii, we are measuring the mean surface brightness in non-overlapping radial
    … …  
    1314// track the non-overlapping radius values.
    1415
     16// Photo interpolates the image of interest to place the peak on the center of the central
     17// pixel, and then uses the exact fractions of the pixels in each of the first few annuli.
     18// Seems like a reasonable thing, but is there any significance to the difference?
     19
    1520// XXX move the resulting elements from profile to extpars->petrosian?
    1621bool psphotPetrosianRadialBins (pmSource *source, float radiusMax, float skynoise) {
    1722
    18     pmSourceRadialProfile *profile = source->extpars->profile;
    19 
    20     float skyModelErrorSQ = PS_SQR(skynoise);
    21 
    22     psVector *radius = profile->radiusElliptical;
    23     psVector *flux = profile->fluxElliptical;
     23    psAssert (source, "missing source");
     24    psAssert (source->extpars, "missing extpars");
     25    psAssert (source->extpars->ellipticalFlux, "missing ellipticalFlux");
     26
     27    psVector *radius = source->extpars->ellipticalFlux->radiusElliptical;
     28    psVector *flux = source->extpars->ellipticalFlux->fluxElliptical;
    2429
    2530    // sort incoming vectors by radius
    2631    pmSourceRadialProfileSortPair (radius, flux);
     32
     33    if (!source->extpars->petProfile) {
     34        source->extpars->petProfile = pmSourceRadialProfileAlloc();
     35    }
     36    pmSourceRadialProfile *profile = source->extpars->petProfile;
     37
     38    float skyModelErrorSQ = PS_SQR(skynoise);
    2739
    2840    int nMax = radiusMax;
    … …  
    107119    psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
    108120    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    109     // psStats *stats = psStatsAlloc(PS_STAT_FITTED_MEAN_V4 | PS_STAT_FITTED_STDEV_V4);
    110 
    111     // integrate flux, radius for each of these bins.  since flux is sorted by radius,
    112     // we can do this fairly quickly
    113121
    114122    bool done = false;
    … …  
    133141                dvalue = NAN;
    134142            }
    135             // binSB->data.F32[nOut] = stats->sampleMedian;
     143
    136144            binSB->data.F32[nOut] = value;
    137145            binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ);
    138             // binSB->data.F32[nOut] = stats->fittedMean;
    139             // binSBstdev->data.F32[nOut] = sqrt(PS_SQR(stats->fittedStdev) / values->n + skyModelErrorSQ);
    140146
    141147            // error in the SB is the stdev per bin / sqrt (number of pixels)
    … …  
    143149            // residual flux, but the sky from the sky model)
    144150
    145             psTrace ("psphot", 5, "%3d  %5.1f %5.1f : %5.1f  %5.2f\n",
    146                      nOut, radAlp->data.F32[nOut], radBet->data.F32[nOut], binSB->data.F32[nOut], binSBstdev->data.F32[nOut]);
     151            psTrace ("psphot", 5, "%3d  %5.1f %5.1f : %5.1f  %5.2f\n", nOut, radAlp->data.F32[nOut], radBet->data.F32[nOut], binSB->data.F32[nOut], binSBstdev->data.F32[nOut]);
    147152
    148153            nOut ++;
    … …  
    163168    // XXX I think this misses the last radial bin -- do we care?
    164169
     170    // interpolate any bins that were empty (extrapolate to center if needed)
     171    if (!isfinite(binSB->data.F32[0]) && !isfinite(binSB->data.F32[1])) {
     172        psWarning ("center 2 bins of source at %f, %f are NAN, skipping this source", source->peak->xf, source->peak->yf);
     173        // XXX raise a flag
     174        psFree(binSB);
     175        psFree(binSBstdev);
     176        psFree(binRad);
     177        psFree(binArea);
     178        psFree(radMin);
     179        psFree(radMax);
     180        psFree(radAlp);
     181        psFree(radBet);
     182        psFree(values);
     183        psFree(stats);
     184        return false;
     185    }
     186
     187    // if center bin is empty assume same SB as next radius (probably true due to PSF)
     188    if (!isfinite(binSB->data.F32[0])) {
     189        binSB->data.F32[0] = binSB->data.F32[1];
     190        binSBstdev->data.F32[0] = binSBstdev->data.F32[1];
     191    }
     192
     193    // interpolate any bins that were empty (if center if needed)
     194    for (int i = 1; i < binSB->n - 1; i++) {
     195        if (isfinite(binSB->data.F32[i])) continue;
     196        binSB->data.F32[i] = InterpolateValues (binRad->data.F32[i-1], binSB->data.F32[i-1], binRad->data.F32[i+1], binSB->data.F32[i+1], binRad->data.F32[i]);
     197        binSBstdev->data.F32[i] = InterpolateValues (binRad->data.F32[i-1], binSBstdev->data.F32[i-1], binRad->data.F32[i+1], binSBstdev->data.F32[i+1], binRad->data.F32[i]);
     198    }
     199
     200    psFree(profile->binSB);
     201    psFree(profile->binSBstdev);
     202    psFree(profile->radialBins);
     203    psFree(profile->area);
     204
    165205    // save the vectors
    166206    profile->radialBins = binRad;

  • trunk/psphot/src/psphotPetrosianStats.c

    r25755 r27819  
    1010bool psphotPetrosianStats (pmSource *source) {
    1111
    12     pmSourceRadialProfile *profile = source->extpars->profile;
    13 
    14     float petRadius = NAN;
    15     float petFlux = NAN;
     12    psAssert (source, "missing source");
     13    psAssert (source->extpars, "missing extpars");
     14    psAssert (source->extpars->petProfile, "missing petProfile");
     15
     16    pmSourceRadialProfile *profile = source->extpars->petProfile;
    1617
    1718    psVector *binSB      = profile->binSB;
    … …  
    2829    psVector *areaSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    2930
     31    float petRadius = NAN;
     32    float petFlux = NAN;
     33
    3034    bool anyPetro = false;
    3135    bool manyPetro = false;
    … …  
    3842    int lowestSignificantRadius = 0;
    3943    float lowestSignificantRatio = 1.0;
     44
     45    // find the Petrosian Radius and Petrosian Flux
    4046
    4147    int nOut = 0;
    … …  
    142148    }
    143149
    144     if (!source->extpars->petrosian_80) {
    145         source->extpars->petrosian_80 = pmSourceExtendedFluxAlloc ();
    146     }
    147     pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80;
     150    // now measure the radii R90 and R50 where flux = 0.9 (or 0.5) * petFlux;
     151    float flux90 = 0.9 * petFlux;
     152    float flux50 = 0.5 * petFlux;
     153    float R50 = NAN;
     154    float R90 = NAN;
     155    bool found50 = false;
     156    bool found90 = false;
     157    // XXX use bisection to do this faster:
     158    for (int i = 0; !(found50 && found90) && i < refRadius->n; i++) {
     159        if (!found50 && (fluxSum->data.F32[i] > flux50)) {
     160            if (i == 0) {
     161                psWarning ("does this case make any sense? (fluxSum[0] > flux50)");
     162                continue;
     163            } else {
     164                R50 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50);
     165                found50 = true;
     166            }
     167        }
     168        if (!found90 && (fluxSum->data.F32[i] > flux90)) {
     169            if (i == 0) {
     170                psWarning ("does this case make any sense? (fluxSum[0] > flux90)");
     171                continue;
     172            } else {
     173                R90 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90);
     174                found90 = true;
     175            }
     176        }
     177    }
     178
    148179
    149180    // XXX save flags (anyPetro, manyPetro)
    150     petrosian->radius = petRadius;
    151     petrosian->flux   = petFlux;
    152 
    153     // psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius);
     181    source->extpars->petrosianRadius = petRadius;
     182    source->extpars->petrosianFlux   = petFlux;
     183    source->extpars->petrosianR50    = R50;
     184    source->extpars->petrosianR90    = R90;
     185   
     186    // XXX add the errors
     187    source->extpars->petrosianRadiusErr = NAN;
     188    source->extpars->petrosianFluxErr   = NAN;
     189    source->extpars->petrosianR50Err    = NAN;
     190    source->extpars->petrosianR90Err    = NAN;
     191
     192    fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf);
     193    psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius);
    154194
    155195    psFree(fluxSum);

  • 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

  • trunk/psphot/src/psphotRadialProfile.c

    r25755 r27819  
    33bool psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal) {
    44
     5    bool status;
     6
    57    // allocate pmSourceExtendedParameters, if not already defined
    68    if (!source->extpars) {
    79        source->extpars = pmSourceExtendedParsAlloc ();
    8     }
    9 
    10     if (!source->extpars->profile) {
    11         source->extpars->profile = pmSourceRadialProfileAlloc ();
    1210    }
    1311
    … …  
    1715    float fluxMin = 0.0;
    1816    float fluxMax = source->peak->flux;
     17
     18    bool RAW_RADIUS = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_RAW_RADIUS");
    1919
    2020    // generate a series of radial profiles at Nsec evenly spaced angles.  the profile flux
    … …  
    2828    // use the radial profiles to determine the radius of a given isophote.  this isophote
    2929    // is used to determine the elliptical shape of the object, so it has a relatively high
    30     // value (nominally 50% of the peak)
     30    // value (nominally 25% of the peak)
    3131    if (!psphotRadiiFromProfiles (source, fluxMin, fluxMax)) {
    3232        psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles");
    … …  
    4242    // generate a single, normalized radial profile following the elliptical contours.
    4343    // the radius is normalized by the axis ratio so that on the major axis, 1 pixel = 1 pixel
    44     if (!psphotEllipticalProfile (source)) {
     44    if (!psphotEllipticalProfile (source, RAW_RADIUS)) {
    4545        psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
     46        return false;
     47    }
     48 
     49    // generated profile in averaged bins
     50    if (!psphotRadialBins (recipe, source, Rmax, skynoise)) {
     51        psError (PS_ERR_UNKNOWN, false, "failed to generate radial bins");
    4652        return false;
    4753    }

  • trunk/psphot/src/psphotRadialProfileByAngles.c

    r25755 r27819  
    1212bool psphotRadialProfilesByAngles (pmSource *source, int Nsec, float Rmax) {
    1313
     14    psAssert (source->extpars, "define extpars");
     15
    1416    // we want to have an even number of sectors so we can do 180 deg symmetrizing
    1517    Nsec = (Nsec % 2) ? Nsec + 1 : Nsec;
    1618    float dtheta = 2.0*M_PI / Nsec;
    1719
    18     pmSourceRadialProfile *profile = source->extpars->profile;
     20    if (!source->extpars->radFlux) {
     21        source->extpars->radFlux = pmSourceRadialFluxAlloc();
     22    }
     23    pmSourceRadialFlux *profile = source->extpars->radFlux;
    1924    psFree(profile->radii);
    2025    psFree(profile->fluxes);
    … …  
    2429    profile->fluxes = psArrayAllocEmpty(Nsec);
    2530    profile->theta = psVectorAllocEmpty(Nsec, PS_TYPE_F32);
    26 
    2731
    2832    for (int i = 0; i < Nsec; i++) {

  • trunk/psphot/src/psphotRadiiFromProfiles.c

    r25755 r27819  
    55bool psphotRadiiFromProfiles (pmSource *source, float fluxMin, float fluxMax) {
    66
    7     pmSourceRadialProfile *profile = source->extpars->profile;
     7    psAssert (source, "missing source");
     8    psAssert (source->extpars, "missing extpars");
     9    psAssert (source->extpars->radFlux, "missing radFlux");
     10
     11    pmSourceRadialFlux *profile = source->extpars->radFlux;
    812
    913    psFree(profile->isophotalRadii);
    … …  
    146150
    147151    // show the results
    148     // psphotPetrosianVisualProfileRadii (radius, flux, radiusBinned, fluxBinned, Ro);
     152    // psphotPetrosianVisualProfileRadii (radius, flux, radiusBinned, fluxBinned, fluxMax, Ro);
    149153
    150154    psFree(fluxBinned);
    … …  
    152156    return Ro;
    153157}
     158
     159
     160

  • trunk/psphot/src/psphotVisual.c

    r27532 r27819  
    446446    KapaSetFont (myKapa, "courier", 14);
    447447
     448    section.bg = KapaColorByName ("none"); // XXX probably should be 'none'
     449
    448450    float SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM");
    449451
    … …  
    538540    KapaSetLimits (myKapa, &graphdata);
    539541
     542    graphdata.padXm = NAN;
     543    graphdata.padYm = NAN;
     544    graphdata.padXp = 0.5;
     545    graphdata.padYp = 0.5;
    540546    KapaBox (myKapa, &graphdata);
     547
    541548    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
    542549    KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM);
    … …  
    547554    graphdata.style = 2;
    548555    KapaPrepPlot (myKapa, nF, &graphdata);
     556
    549557    KapaPlotVector (myKapa, nF, xFaint->data.F32, "x");
    550558    KapaPlotVector (myKapa, nF, yFaint->data.F32, "y");
    … …  
    562570    section.dy = 0.25;
    563571    section.x  = 0.00;
    564     section.y  = 0.80;
     572    section.y  = 0.75;
    565573    section.name = psStringCopy ("MagMyy");
    566574    KapaSetSection (myKapa, &section);
    … …  
    574582    KapaSetLimits (myKapa, &graphdata);
    575583
     584    graphdata.padXm = 0.5;
     585    graphdata.padYm = NAN;
     586    graphdata.padXp = NAN;
     587    graphdata.padYp = 0.5;
    576588    strcpy (graphdata.labels, "0210");
    577589    KapaBox (myKapa, &graphdata);
    … …  
    598610    section.dx = 0.25;
    599611    section.dy = 0.75;
    600     section.x  = 0.80;
     612    section.x  = 0.75;
    601613    section.y  = 0.00;
    602614    section.name = psStringCopy ("MagMxx");
    … …  
    611623    KapaSetLimits (myKapa, &graphdata);
    612624
     625    graphdata.padXm = NAN;
     626    graphdata.padYm = 0.5;
     627    graphdata.padXp = 0.5;
     628    graphdata.padYp = NAN;
    613629    strcpy (graphdata.labels, "2001");
    614630    KapaBox (myKapa, &graphdata);
    … …  
    13501366    assert (maskVal);
    13511367
     1368    section.bg  = KapaColorByName ("none"); // XXX probably should be 'none'
     1369
    13521370    KapaClearPlots (myKapa);
    13531371    // first section : mag vs CR nSigma
    … …  
    15861604    KapaSetFont (myKapa, "courier", 14);
    15871605
     1606    section.bg  = KapaColorByName ("none"); // XXX probably should be 'none'
     1607
    15881608    // select the max psfX,Y values for the plot limits
    15891609    float Xmin = 1000.0, Xmax = 0.0;
    … …  
    17541774    KapaSetLimits (myKapa, &graphdata);
    17551775
     1776    graphdata.padXm = NAN;
     1777    graphdata.padYm = NAN;
     1778    graphdata.padXp = 0.5;
     1779    graphdata.padYp = 0.5;
    17561780    KapaBox (myKapa, &graphdata);
    17571781    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
    … …  
    18221846    KapaSetLimits (myKapa, &graphdata);
    18231847
     1848    graphdata.padXm = 0.5;
     1849    graphdata.padYm = NAN;
     1850    graphdata.padXp = NAN;
     1851    graphdata.padYp = 0.5;
    18241852    strcpy (graphdata.labels, "0210");
    18251853    KapaBox (myKapa, &graphdata);
    … …  
    18701898    section.dx = 0.25;
    18711899    section.dy = 0.60;
    1872     section.x  = 0.80;
     1900    section.x  = 0.75;
    18731901    section.y  = 0.00;
    18741902    section.name = psStringCopy ("MagMxx");
    … …  
    18831911    KapaSetLimits (myKapa, &graphdata);
    18841912
     1913    graphdata.padXm = NAN;
     1914    graphdata.padYm = 0.5;
     1915    graphdata.padXp = 0.5;
     1916    graphdata.padYp = NAN;
    18851917    strcpy (graphdata.labels, "2001");
    18861918    KapaBox (myKapa, &graphdata);
    … …  
    19301962    // fourth section: MagSigma
    19311963    section.dx = 0.75;
    1932     section.dy = 0.15;
     1964    section.dy = 0.20;
    19331965    section.x  = 0.00;
    1934     section.y  = 0.65;
     1966    section.y  = 0.60;
    19351967    section.name = psStringCopy ("MagSigma");
    19361968    KapaSetSection (myKapa, &section);
    … …  
    19441976    KapaSetLimits (myKapa, &graphdata);
    19451977
     1978    graphdata.padXm = 0.5;
     1979    graphdata.padYm = NAN;
     1980    graphdata.padXp = 0.5;
     1981    graphdata.padYp = 0.5;
    19461982    strcpy (graphdata.labels, "0100");
    19471983    KapaBox (myKapa, &graphdata);
    … …  
    20532089    psFree (mDEF);
    20542090    psFree (sDEF);
     2091
     2092    psFree (xLOW);
     2093    psFree (yLOW);
     2094    psFree (mLOW);
     2095    psFree (sLOW);
    20552096
    20562097    psFree (xCR);
    … …  
    22972338        if (!source) continue;
    22982339        if (!source->extpars) continue;
    2299         if (!source->extpars->profile) continue;
    2300         if (!source->extpars->petrosian_80) continue;
    2301 
    2302         pmSourceRadialProfile *profile = source->extpars->profile;
    2303         pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80;
     2340        if (!source->extpars->petProfile) continue;
     2341
     2342        float petrosianRadius = source->extpars->petrosianRadius;
     2343        psEllipseAxes *axes = &source->extpars->axes;
    23042344
    23052345        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    23062346        overlay[Noverlay].x = source->peak->xf;
    23072347        overlay[Noverlay].y = source->peak->yf;
    2308         overlay[Noverlay].dx = 2.0*petrosian->radius;
    2309         overlay[Noverlay].dy = 2.0*petrosian->radius*profile->axes.minor/profile->axes.major;
    2310         overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD;
     2348        overlay[Noverlay].dx = 1.0*petrosianRadius;
     2349        overlay[Noverlay].dy = 1.0*petrosianRadius*axes->minor/axes->major;
     2350        overlay[Noverlay].angle = axes->theta * PS_DEG_RAD;
    23112351        overlay[Noverlay].text = NULL;
    23122352        Noverlay ++;
    23132353        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    23142354
    2315         overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    2316         overlay[Noverlay].x = source->peak->xf;
    2317         overlay[Noverlay].y = source->peak->yf;
    2318         overlay[Noverlay].dx = 4.0*petrosian->radius;
    2319         overlay[Noverlay].dy = 4.0*petrosian->radius*profile->axes.minor/profile->axes.major;
    2320         overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD;
    2321         overlay[Noverlay].text = NULL;
    2322         Noverlay ++;
    2323         CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     2355        // overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     2356        // overlay[Noverlay].x = source->peak->xf;
     2357        // overlay[Noverlay].y = source->peak->yf;
     2358        // overlay[Noverlay].dx = 2.0*petrosianRadius;
     2359        // overlay[Noverlay].dy = 2.0*petrosianRadius*axes->minor/axes->major;
     2360        // overlay[Noverlay].angle = axes->theta * PS_DEG_RAD;
     2361        // overlay[Noverlay].text = NULL;
     2362        // Noverlay ++;
     2363        // CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    23242364    }
    23252365
Note: See TracChangeset for help on using the changeset viewer.