IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 19, 2009, 8:49:41 AM (17 years ago)
Author:
eugene
Message:

further work on the petrosian analysis & extended sources

Location:
branches/eam_branches/20090715/psphot/src
Files:
16 edited
1 moved

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090715/psphot/src/Makefile.am

    r25433 r25452  
    2525libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    2626
    27 bin_PROGRAMS = psphot psphotTest psphotMomentsStudy psphotPetrosianStudy
     27bin_PROGRAMS = psphot psphotTest psphotMomentsStudy
    2828# bin_PROGRAMS = psphotPetrosianStudy
    2929# bin_PROGRAMS = psphot
    … …  
    4141psphotMomentsStudy_LDADD = libpsphot.la
    4242
    43 psphotPetrosianStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    44 psphotPetrosianStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    45 psphotPetrosianStudy_LDADD = libpsphot.la
     43# psphotPetrosianStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     44# psphotPetrosianStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
     45# psphotPetrosianStudy_LDADD = libpsphot.la
    4646
    4747psphot_SOURCES = \
    … …  
    6868        psphotMomentsStudy.c
    6969
    70 psphotPetrosianStudy_SOURCES = \
    71         psphotPetrosianStudy.c
     70# psphotPetrosianStudy_SOURCES = \
     71#         psphotPetrosianStudy.c
    7272
    7373libpsphot_la_SOURCES = \
    … …  
    113113        psphotExtendedSourceAnalysis.c \
    114114        psphotExtendedSourceFits.c     \
    115         psphotRadialProfile.c          \
    116         psphotPetrosian.c              \
    117         psphotIsophotal.c              \
    118         psphotAnnuli.c                 \
    119         psphotKron.c                   \
    120115        psphotKernelFromPSF.c          \
    121116        psphotPSFConvModel.c           \
    … …  
    138133        psphotThreadTools.c            \
    139134        psphotAddNoise.c               \
    140         psphotPetrosianProfile.c       \
    141         psphotRadialProfileByAngle.c   \
     135        psphotRadialProfile.c          \
     136        psphotRadialProfileByAngles.c  \
    142137        psphotRadiiFromProfiles.c      \
    143138        psphotEllipticalContour.c      \
    144139        psphotEllipticalProfile.c      \
     140        psphotPetrosian.c              \
    145141        psphotPetrosianRadialBins.c    \
    146         psphotPetrosianVisual.c        \
    147142        psphotPetrosianStats.c         \
    148         psphotPetrosianAnalysis.c      \
    149         pmPetrosian.c                  \
    150143        psphotEfficiency.c
     144
     145# re-instate these
     146#       psphotIsophotal.c              \
     147#       psphotAnnuli.c                 \
     148#       psphotKron.c                   \
     149#       psphotPetrosianVisual.c        \
     150#
     151
     152# test versions
     153#       psphotPetrosianProfile.c       \
     154#       psphotPetrosianAnalysis.c      \
     155#
    151156
    152157include_HEADERS = \
    153158        psphot.h \
    154         pmPetrosian.h \
    155159        psphotErrorCodes.h
    156160

  • branches/eam_branches/20090715/psphot/src/psphot.h

    r25433 r25452  
    88#include <psmodules.h>
    99#include "psphotErrorCodes.h"
    10 #include "pmPetrosian.h"
    1110
    1211#define PSPHOT_RECIPE "PSPHOT" // Name of the recipe to use
    … …  
    169168psKernel       *psphotKernelFromPSF (pmSource *source, int nPix);
    170169
    171 bool            psphotRadialProfile (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
    172 bool            psphotPetrosian (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
    173 bool            psphotIsophotal (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
    174 bool            psphotAnnuli (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
    175 bool            psphotKron (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
     170// functions related to extended source analysis
     171bool  psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal);
     172bool  psphotRadialProfilesByAngles (pmSource *source, int Nsec, float Rmax);
     173float psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax);
     174bool  psphotRadiiFromProfiles (pmSource *source, float fluxMin, float fluxMax);
     175bool  psphotEllipticalProfile (pmSource *source);
     176bool  psphotEllipticalContour (pmSource *source);
    176177
    177178// psphotVisual functions
    … …  
    194195bool psphotVisualPlotApResid (psArray *sources);
    195196bool psphotVisualPlotSourceSize (psMetadata *recipe, psArray *sources);
    196 
    197 bool psphotVisualShowPetrosian (pmSource *source, pmPetrosian *petrosian);
    198 bool psphotPetrosianAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe);
     197bool psphotVisualShowPetrosians (psArray *sources);
     198
     199// bool psphotPetrosianAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe);
     200// bool psphotPetrosianProfile (pmReadout *readout, pmSource *source, float skynoise);
     201
     202bool psphotPetrosian (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal);
     203bool psphotPetrosianRadialBins (pmSource *source, float radiusMax, float skynoise);
     204bool psphotPetrosianStats (pmSource *source);
     205
     206// XXX old versions, currently disabled
     207// bool            psphotIsophotal (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
     208// bool            psphotAnnuli (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
     209// bool            psphotKron (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
     210
     211// XXX visualization functions related to radial profiles (disabled)
     212// bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux);
     213// bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef);
     214// bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian);
     215// bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin,
     216//                               psVector *refRadius, psVector *meanSB,
     217//                               psVector *petRatio, psVector *petRatioErr, psVector *fluxSum,
     218//                               float petRadius, float ratioForRadius,
     219//                               float petFlux, float radiusForFlux);
    199220
    200221bool psphotImageQuality (psMetadata *recipe, psArray *sources);

  • branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c

    r25433 r25452  
    55psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord);
    66
    7 bool psphotEllipticalContour (pmSource *source, pmPetrosian *petrosian) {
     7bool psphotEllipticalContour (pmSource *source) {
     8
     9    pmSourceRadialProfile *profile = source->extpars->profile;
    810
    911    // use LMM to fit theta vs radius to an ellipse
    10     psVector *theta = petrosian->theta;
    11     psVector *radius = petrosian->isophotalRadii;
     12    psVector *theta = profile->theta;
     13    psVector *radius = profile->isophotalRadii;
    1214
    1315    // find Rmin and Rmax for the initial guess
    … …  
    8385    /// XXX rationalize? if epsilon > 1, flip major and minor axes (rotate by 90 degrees)
    8486    if (params->data.F32[PAR_EPSILON] < 1.0) {
    85         petrosian->axes.major = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON];
    86         petrosian->axes.minor = params->data.F32[PAR_RMIN];
    87         petrosian->axes.theta = params->data.F32[PAR_PHI];
     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];
    8890    } else {
    89         petrosian->axes.major = params->data.F32[PAR_RMIN];
    90         petrosian->axes.minor = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON];
    91         petrosian->axes.theta = params->data.F32[PAR_PHI] + 0.5*M_PI;
     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;
    9294    }
    9395
    9496    psTrace ("psphot", 4, "# fitted values:\n");
    95     psTrace ("psphot", 4, "Phi:   %f\n", petrosian->axes.theta*PS_DEG_RAD);
    96     psTrace ("psphot", 4, "Rmaj:  %f\n", petrosian->axes.major);
    97     psTrace ("psphot", 4, "Rmin:  %f\n", petrosian->axes.minor);
     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);
    98100   
    99101    // show the results

  • branches/eam_branches/20090715/psphot/src/psphotEllipticalProfile.c

    r25433 r25452  
    11# include "psphotInternal.h"
    22
    3 bool psphotEllipticalProfile (pmSource *source, pmPetrosian *petrosian) {
     3bool psphotEllipticalProfile (pmSource *source) {
    44
    5     petrosian->radiusElliptical = psVectorAllocEmpty(100, PS_TYPE_F32);
    6     petrosian->fluxElliptical = psVectorAllocEmpty(100, PS_TYPE_F32);
     5    pmSourceRadialProfile *profile = source->extpars->profile;
    76
    8     psVector *radius = petrosian->radiusElliptical;
    9     psVector *flux = petrosian->fluxElliptical;
     7    profile->radiusElliptical = psVectorAllocEmpty(100, PS_TYPE_F32);
     8    profile->fluxElliptical = psVectorAllocEmpty(100, PS_TYPE_F32);
     9
     10    psVector *radius = profile->radiusElliptical;
     11    psVector *flux = profile->fluxElliptical;
    1012
    1113    // the psEllipse functions use z = 0.5(x/Sxx)^2 + 0.5(y/Syy)^2 + x y Sxy
    … …  
    2022    psEllipseAxes axes;
    2123    axes.major = M_SQRT1_2;
    22     axes.minor = M_SQRT1_2 * (petrosian->axes.minor / petrosian->axes.major);
     24    axes.minor = M_SQRT1_2 * (profile->axes.minor / profile->axes.major);
    2325
    2426    // axes.major = 1.0;
    25     // axes.minor = petrosian->axes.minor / petrosian->axes.major;
     27    // axes.minor = profile->axes.minor / profile->axes.major;
    2628
    27     axes.theta = petrosian->axes.theta;
     29    axes.theta = profile->axes.theta;
    2830    psEllipseShape shape = psEllipseAxesToShape (axes);
    2931
    … …  
    5658    // psVector *radiusRaw = psVectorAllocEmpty(100, PS_TYPE_F32);
    5759    // psVector *fluxRaw = psVectorAllocEmpty(100, PS_TYPE_F32);
    58     // for (int i = 0; i < petrosian->radii->n; i++) {
    59     //   psVector *r = petrosian->radii->data[i];
    60     //   psVector *f = petrosian->fluxes->data[i];
     60    // for (int i = 0; i < profile->radii->n; i++) {
     61    //   psVector *r = profile->radii->data[i];
     62    //   psVector *f = profile->fluxes->data[i];
    6163    //   for (int j = 0; j < r->n; j++) {
    6264    //  psVectorAppend(radiusRaw, r->data.F32[j]);

  • branches/eam_branches/20090715/psphot/src/psphotExtendedSourceAnalysis.c

    r21519 r25452  
    2121    }
    2222
    23     // option to limit analysis to a specific region
    24     char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
    25     psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
    26     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
     23    // XXX require petrosian analysis for non-linear fits?
     24
     25    // XXX temporary user-supplied systematic sky noise measurement (derive from background model)
     26    float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
     27
     28# if (0)
     29    // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
     30    // XXX use this to set skynoise
     31    pmReadout *backModel = psphotSelectBackground (config, view);
     32    pmReadout *backStdev = psphotSelectBackgroundStdev (config, view);
     33# endif
    2734
    2835    // S/N limit to perform full non-linear fits
    … …  
    3845    sources = psArraySort (sources, pmSourceSortBySN);
    3946
    40     // XXX some init functions for the extended source recipe options?
     47    // option to limit analysis to a specific region
     48    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
     49    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
     50    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
    4151
    4252    // choose the sources of interest
    … …  
    4959        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
    5060        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
     61        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
     62        if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
     63        if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
    5164
    5265        // limit selection to some SN limit
    … …  
    6881        // if we request any of these measurements, we require the radial profile
    6982        if (doPetrosian || doIsophotal || doAnnuli || doKron) {
    70             if (!psphotRadialProfile (source, recipe, maskVal)) {
     83            if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) {
    7184                // all measurements below require the radial profile; skip them all
    7285                // re-subtract the object, leave local sky
    … …  
    7992        }
    8093
     94        // Petrosian Mags
     95        if (doPetrosian) {
     96            if (!psphotPetrosian (source, recipe, skynoise, maskVal)) {
     97                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
     98            } else {
     99                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
     100                Npetro ++;
     101                source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
     102            }
     103        }
     104
     105# if (0)
    81106        // Isophotal Mags
    82107        if (doIsophotal) {
    … …  
    89114            }
    90115        }
    91 
    92         // Petrosian Mags
    93         if (doPetrosian) {
    94             if (!psphotPetrosian (source, recipe, maskVal)) {
    95                 psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    96             } else {
    97                 psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    98                 Npetro ++;
    99                 source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
    100             }
    101         }
    102 
    103116        // Kron Mags
    104117        if (doKron) {
    … …  
    111124            }
    112125        }
    113 
    114         // Radial Annuli
    115         if (doAnnuli) {
    116             if (!psphotAnnuli (source, recipe, maskVal)) {
    117                 psError(PSPHOT_ERR_UNKNOWN, false, "failure in Annuli analysis");
    118                 return false;
    119             }
    120             psTrace ("psphot", 5, "measured annuli for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    121             Nannuli ++;
    122             source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
    123         }
     126# endif
    124127
    125128        // re-subtract the object, leave local sky
    126129        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    127130        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     131
     132        if (source->extpars) {
     133            pmSourceRadialProfileFreeVectors(source->extpars->profile);
     134        }
    128135    }
    129136
    … …  
    133140    psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
    134141    psLogMsg ("psphot", PS_LOG_INFO, "  %d kron\n", Nkron);
     142
     143    psphotVisualShowResidualImage (readout);
     144
     145    if (doPetrosian) {
     146        psphotVisualShowPetrosians (sources);
     147    }
     148
    135149    return true;
    136150}

  • branches/eam_branches/20090715/psphot/src/psphotFitSourcesLinear.c

    r25433 r25452  
    7676            if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
    7777        } else {
    78             if (source->mode & PM_SOURCE_MODE_BLEND) continue;
     78            // if (source->mode & PM_SOURCE_MODE_BLEND) continue;
    7979        }
    8080

  • branches/eam_branches/20090715/psphot/src/psphotPetrosian.c

    r21183 r25452  
    11# include "psphotInternal.h"
    22
    3 bool psphotPetrosian (pmSource *source, psMetadata *recipe, psImageMaskType maskVal) {
     3bool psphotPetrosian (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal) {
    44
    5   bool status;
     5    // XXX these need to go into recipe values
     6    float Rmax = 200;
    67
    7   assert (source->extpars);
    8   assert (source->extpars->profile);
    9   assert (source->extpars->profile->radius);
    10   assert (source->extpars->profile->flux);
     8    psAssert (source->extpars, "need to run psphotRadialProfile first");
     9    psAssert (source->extpars->profile, "need to run psphotRadialProfile first");
    1110
    12   psVector *radius = source->extpars->profile->radius;
    13   psVector *flux = source->extpars->profile->flux;
     11    // integrate the radial profile for radial bins defined for the petrosian measurement:
     12    // SB_i (r_i) where \alpha r_i < r < \beta r_i
     13    if (!psphotPetrosianRadialBins (source, Rmax, skynoise)) {
     14        psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
     15        return false;
     16    }
     17 
     18    // use the SB_i from above to calculate the petrosian radius and the flux within that radius
     19    if (!psphotPetrosianStats (source)) {
     20        psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
     21        return false;
     22    }
     23 
     24    psTrace ("psphot", 3, "source at %f,%f: petrosian radius: %f, flux: %f, axis ratio: %f, angle: %f",
     25             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);
    1430
    15   // flux at which to measure isophotal parameters
    16   float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0");
    17   float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO");
    18   assert (status);
    19 
    20   // first find flux at R0
    21   int firstAbove = -1;
    22   int lastBelow = -1;
    23   for (int i = 0; i < radius->n; i++) {
    24     if (radius->data.F32[i] < PETROSIAN_R0) lastBelow = i;
    25     if ((firstAbove < 0) && (radius->data.F32[i] > PETROSIAN_R0)) firstAbove = i;
    26   }
    27   // if we don't go out far enough, we have a problem...
    28   if (lastBelow == radius->n - 1) {
    29     psTrace ("psphot", 5, "did not go out far enough to reach petrosian reference radius...");
    30     // XXX skip object? raise a flag ?
    31     return false;
    32   }
    33   if (firstAbove < 0) {
    34     psTrace ("psphot", 5, "did not go out far enough to bound petrosian reference radius");
    35     // XXX raise a flag ?
    36     return false;
    37   }
    38 
    39   // average flux in this range
    40   float fluxR0 = 0.0;
    41   int fluxRn = 0;
    42   for (int i = PS_MIN(firstAbove, lastBelow); i <= PS_MAX(firstAbove, lastBelow); i++) {
    43     fluxR0 += flux->data.F32[i];
    44     fluxRn ++;
    45   }
    46   fluxR0 /= (float)(fluxRn);
    47 
    48   // target flux for petrosian radius
    49   float fluxRP = fluxR0 * PETROSIAN_RF;
    50 
    51   // find the first bin below the flux level and the last above the level
    52   // XXX can this be done faster with bisection?
    53   // XXX do I need to worry about crazy outliers?
    54   // XXX should i be smoothing or fitting the curve?
    55   int firstBelow = -1;
    56   int lastAbove = -1;
    57   for (int i = 0; i < flux->n; i++) {
    58     if (flux->data.F32[i] > fluxRP) lastAbove = i;
    59     if ((firstBelow < 0) && (flux->data.F32[i] < fluxRP)) firstBelow = i;
    60   }
    61   // if we don't go out far enough, we have a problem...
    62   if (lastAbove == radius->n - 1) {
    63     psTrace ("psphot", 5, "did not go out far enough to reach petrosian radius...");
    64     // XXX skip object? raise a flag ?
    65     return false;
    66   }
    67   if (firstBelow < 0) {
    68     psTrace ("psphot", 5, "did not go out far enough to bound petrosian radius");
    69     // XXX raise a flag ?
    70     return false;
    71   }
    72 
    73   // need to examine pixels in this vicinity
    74   float fluxFirst = 0;
    75   float fluxLast = 0;
    76   for (int i = 0; i <= PS_MAX(firstBelow, lastAbove); i++) {
    77     if (i <= firstBelow) {
    78       fluxFirst += flux->data.F32[i];
    79     }
    80     if (i <= lastAbove) {
    81       fluxLast += flux->data.F32[i];
    82     }
    83   }
    84   float fluxRPSum    = 0.5*(fluxLast + fluxFirst);
    85   float fluxRPSumErr = 0.5*fabs(fluxLast - fluxFirst);
    86   // XXX need to use the weight appropriately here...
    87 
    88   float rad     = 0.5*(radius->data.F32[firstBelow] + radius->data.F32[lastAbove]);
    89   float radErr  = 0.5*fabs(radius->data.F32[firstBelow] - radius->data.F32[lastAbove]);
    90 
    91   if (!source->extpars->petrosian) {
    92     source->extpars->petrosian = pmSourcePetrosianValuesAlloc ();
    93   }
    94 
    95   // these are uncalibrated: instrumental mags and pixel units
    96   source->extpars->petrosian->mag    = -2.5*log10(fluxRPSum);
    97   source->extpars->petrosian->magErr = fluxRPSumErr / fluxRPSum;
    98 
    99   source->extpars->petrosian->rad    = rad;
    100   source->extpars->petrosian->radErr = radErr;
    101 
    102   psTrace ("psphot", 5, "Petrosian flux:%f +/- %f @ %f +/- %f for %f, %f\n",
    103            source->extpars->petrosian->mag, source->extpars->petrosian->magErr,
    104            source->extpars->petrosian->rad, source->extpars->petrosian->radErr,
    105            source->peak->xf, source->peak->yf);
    106 
    107   return true;
    108 
     31    return true;
    10932}

  • branches/eam_branches/20090715/psphot/src/psphotPetrosianAnalysis.c

    r25433 r25452  
    6666
    6767    psphotVisualShowResidualImage (readout);
    68 
    69     // pause and wait for user input:
    70     // continue, save (provide name), ??
    71     char key[10];
    72     fprintf (stdout, "[c]ontinue? ");
    73     if (!fgets(key, 8, stdin)) {
    74         psWarning("Unable to read option");
    75     }
    76 
    7768    return true;
    7869}

  • branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c

    r25433 r25452  
    1212    pmPetrosian *petrosian = pmPetrosianAlloc();
    1313
     14    // XXX these need to go into recipe values
    1415    int Nsec = 24;
    1516    float Rmax = 200;

  • branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c

    r25433 r25452  
    1313// track the non-overlapping radius values.
    1414
    15 bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax, float skynoise) {
     15// XXX move the resulting elements from profile to extpars->petrosian?
     16bool psphotPetrosianRadialBins (pmSource *source, float radiusMax, float skynoise) {
    1617
    17     // XXX for testing, let's just set this to a value
    18    
     18    pmSourceRadialProfile *profile = source->extpars->profile;
     19
    1920    float skyModelErrorSQ = PS_SQR(skynoise);
    2021
    21     psVector *radius = petrosian->radiusElliptical;
    22     psVector *flux = petrosian->fluxElliptical;
     22    psVector *radius = profile->radiusElliptical;
     23    psVector *flux = profile->fluxElliptical;
    2324
    2425    // sort incoming vectors by radius
    25     psphotPetrosianSortPair (radius, flux);
     26    pmSourceRadialProfileSortPair (radius, flux);
    2627
    2728    int nMax = radiusMax;
    … …  
    163164
    164165    // save the vectors
    165     petrosian->radialBins = binRad;
    166     petrosian->area       = binArea;
    167     petrosian->binSB      = binSB;
    168     petrosian->binSBstdev = binSBstdev;
     166    profile->radialBins = binRad;
     167    profile->area       = binArea;
     168    profile->binSB      = binSB;
     169    profile->binSBstdev = binSBstdev;
    169170
    170171    // psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, source->peak->flux, 0.0);

  • branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c

    r25433 r25452  
    88float InterpolateValues (float X0, float Y0, float X1, float Y1, float X);
    99
    10 bool psphotPetrosianStats (pmSource *source, pmPetrosian *petrosian) {
     10bool psphotPetrosianStats (pmSource *source) {
     11
     12    pmSourceRadialProfile *profile = source->extpars->profile;
    1113
    1214    float petRadius, petFlux;
    1315
    14     psVector *binSB      = petrosian->binSB;
    15     psVector *binSBstdev = petrosian->binSBstdev;
    16     psVector *binRad     = petrosian->radialBins;
    17     psVector *area       = petrosian->area;
     16    psVector *binSB      = profile->binSB;
     17    psVector *binSBstdev = profile->binSBstdev;
     18    psVector *binRad     = profile->radialBins;
     19    psVector *area       = profile->area;
    1820
    1921    psVector *fluxSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    … …  
    139141    }
    140142
    141     // save petRadius, petFlux
     143    if (!source->extpars->petrosian_80) {
     144        source->extpars->petrosian_80 = pmSourceExtendedFluxAlloc ();
     145    }
     146    pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80;
     147
    142148    // XXX save flags (anyPetro, manyPetro)
    143     petrosian->petrosianRadius = petRadius;
    144     petrosian->petrosianFlux   = petFlux;
     149    petrosian->radius = petRadius;
     150    petrosian->flux   = petFlux;
    145151
    146     psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius);
     152    // psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius);
    147153
    148154    psFree(fluxSum);

  • branches/eam_branches/20090715/psphot/src/psphotRadialProfile.c

    r21366 r25452  
    11# include "psphotInternal.h"
    22
    3 # define COMPARE_RADIUS(A,B) (radius->data.F32[A] < radius->data.F32[B])
    4 # define SWAP_RADIUS(TYPE,A,B) { \
    5   float tmp; \
    6   if (A != B) { \
    7     tmp = radius->data.F32[A]; \
    8     radius->data.F32[A] = radius->data.F32[B]; \
    9     radius->data.F32[B] = tmp; \
    10     tmp = flux->data.F32[A]; \
    11     flux->data.F32[A] = flux->data.F32[B]; \
    12     flux->data.F32[B] = tmp; \
    13     tmp = variance->data.F32[A]; \
    14     variance->data.F32[A] = variance->data.F32[B]; \
    15     variance->data.F32[B] = tmp; \
    16   } \
    17 }
    18 
    19 bool psphotRadialProfile (pmSource *source, psMetadata *recipe, psImageMaskType maskVal) {
     3bool psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal) {
    204
    215    // allocate pmSourceExtendedParameters, if not already defined
    … …  
    2812    }
    2913
    30     int nPts = source->pixels->numRows * source->pixels->numCols;
    31     source->extpars->profile->radius = psVectorAllocEmpty (nPts, PS_TYPE_F32);
    32     source->extpars->profile->flux   = psVectorAllocEmpty (nPts, PS_TYPE_F32);
    33     source->extpars->profile->variance = psVectorAllocEmpty (nPts, PS_TYPE_F32);
     14    // XXX these need to go into recipe values
     15    int Nsec = 24;
     16    float Rmax = 200;
     17    float fluxMin = 0.0;
     18    float fluxMax = source->peak->flux;
    3419
    35     psVector *radius = source->extpars->profile->radius;
    36     psVector *flux   = source->extpars->profile->flux;
    37     psVector *variance = source->extpars->profile->variance;
     20    // generate a series of radial profiles at Nsec evenly spaced angles.  the profile flux
     21    // is measured by interpolation for small radii; for large radii, the pixels in a box
     22    // are averaged to increase the S/N
     23    if (!psphotRadialProfilesByAngles (source, Nsec, Rmax)) {
     24        psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian");
     25        return false;
     26    }
    3827
    39     // XXX use the extended source model here for Xo, Yo?
    40     // XXX define a radius scaled to the elliptical contour?
     28    // use the radial profiles to determine the radius of a given isophote.  this isophote
     29    // is used to determine the elliptical shape of the object, so it has a relatively high
     30    // value (nominally 50% of the peak)
     31    if (!psphotRadiiFromProfiles (source, fluxMin, fluxMax)) {
     32        psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles");
     33        return false;
     34    }
    4135
    42     int n = 0;
    43 
    44     float Xo = 0.0;
    45     float Yo = 0.0;
    46 
    47     if (source->modelEXT) {
    48       Xo = source->modelEXT->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
    49       Yo = source->modelEXT->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
    50     } else {
    51       Xo = source->peak->xf - source->pixels->col0;
    52       Yo = source->peak->yf - source->pixels->row0;
     36    // convert the isophotal radius vs angle measurements to an elliptical contour
     37    if (!psphotEllipticalContour (source)) {
     38        psLogMsg ("psphot", 3, "failed to measure elliptical contour");
     39        return false;
    5340    }
    54     for (int iy = 0; iy < source->pixels->numRows; iy++) {
    55         for (int ix = 0; ix < source->pixels->numCols; ix++) {
    56             if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) continue;
    57             radius->data.F32[n] = hypot (ix - Xo, iy - Yo) ;
    58             flux->data.F32[n]   = source->pixels->data.F32[iy][ix];
    59             variance->data.F32[n] = source->variance->data.F32[iy][ix];
    60             n++;
    61         }
     41 
     42    // generate a single, normalized radial profile following the elliptical contours.
     43    // the radius is normalized by the axis ratio so that on the major axis, 1 pixel = 1 pixel
     44    if (!psphotEllipticalProfile (source)) {
     45        psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
     46        return false;
    6247    }
    63     radius->n = n;
    64     variance->n = n;
    65     flux->n = n;
    66 
    67     // sort the vector set by the radius
    68     PSSORT (radius->n, COMPARE_RADIUS, SWAP_RADIUS, NONE);
    69 
     48 
    7049    return true;
    7150}

  • branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngles.c

    r25433 r25452  
    1010psVector *psphotLineValuesBresen (psImage *image, int X1, int Y1, int X2, int Y2, int dW, int swapcoords);
    1111
    12 bool psphotRadialProfilesByAngles (pmSource *source, pmPetrosian *petrosian, int Nsec, float Rmax) {
     12bool psphotRadialProfilesByAngles (pmSource *source, int Nsec, float Rmax) {
    1313
    1414    // we want to have an even number of sectors so we can do 180 deg symmetrizing
    … …  
    1616    float dtheta = 2.0*M_PI / Nsec;
    1717
    18     psFree(petrosian->radii);
    19     psFree(petrosian->fluxes);
    20     psFree(petrosian->theta);
    21 
    22     petrosian->radii = psArrayAllocEmpty(Nsec);
    23     petrosian->fluxes = psArrayAllocEmpty(Nsec);
    24     petrosian->theta = psVectorAllocEmpty(Nsec, PS_TYPE_F32);
     18    pmSourceRadialProfile *profile = source->extpars->profile;
     19    psFree(profile->radii);
     20    psFree(profile->fluxes);
     21    psFree(profile->theta);
     22
     23    profile->radii = psArrayAllocEmpty(Nsec);
     24    profile->fluxes = psArrayAllocEmpty(Nsec);
     25    profile->theta = psVectorAllocEmpty(Nsec, PS_TYPE_F32);
    2526
    2627
    … …  
    7273        }
    7374
    74         psArrayAdd (petrosian->radii, 100, radius);
    75         psArrayAdd (petrosian->fluxes, 100, flux);
    76         psVectorAppend (petrosian->theta, theta);
     75        psArrayAdd (profile->radii, 100, radius);
     76        psArrayAdd (profile->fluxes, 100, flux);
     77        psVectorAppend (profile->theta, theta);
    7778
    7879        // psphotPetrosianVisualProfileByAngle (radius, flux);
    … …  
    8485    for (int i = 0; i < Nsec / 2; i++) {
    8586
    86         psVector *r1 = petrosian->radii->data[i];
    87         psVector *r2 = petrosian->radii->data[i+Nsec/2];
    88 
    89         psVector *f1 = petrosian->fluxes->data[i];
    90         psVector *f2 = petrosian->fluxes->data[i+Nsec/2];
     87        psVector *r1 = profile->radii->data[i];
     88        psVector *r2 = profile->radii->data[i+Nsec/2];
     89
     90        psVector *f1 = profile->fluxes->data[i];
     91        psVector *f2 = profile->fluxes->data[i+Nsec/2];
    9192
    9293        psAssert (r1->n == r2->n, "mis-matched vectors");

  • branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c

    r25433 r25452  
    11# include "psphotInternal.h"
    22
    3 // Given the Petrosian data (radii, fluxes) determine the radius for each profile at the desisred isophote
     3// Given the Radial Profiles (radii, fluxes) determine the radius for each profile at the desired isophote
    44
    5 bool psphotRadiiFromProfiles (pmSource *source, pmPetrosian *petrosian, float fluxMin, float fluxMax) {
     5bool psphotRadiiFromProfiles (pmSource *source, float fluxMin, float fluxMax) {
    66
    7   petrosian->isophotalRadii = psVectorAlloc(petrosian->theta->n, PS_TYPE_F32);
     7    pmSourceRadialProfile *profile = source->extpars->profile;
    88
    9   for (int i = 0; i < petrosian->theta->n; i++) {
    10       psVector *radii = petrosian->radii->data[i];
    11       psVector *fluxes = petrosian->fluxes->data[i];
    12       float radius = psphotRadiusFromProfile (source, radii, fluxes, fluxMin, fluxMax);
     9    psFree(profile->isophotalRadii);
     10    profile->isophotalRadii = psVectorAlloc(profile->theta->n, PS_TYPE_F32);
    1311
    14       // psphotPetrosianVisualProfileByAngle (radii, fluxes, radius);
     12    for (int i = 0; i < profile->theta->n; i++) {
     13        psVector *radii = profile->radii->data[i];
     14        psVector *fluxes = profile->fluxes->data[i];
     15        float radius = psphotRadiusFromProfile (source, radii, fluxes, fluxMin, fluxMax);
    1516
    16       // warn on NAN?
    17       petrosian->isophotalRadii->data.F32[i] = radius;
    18   }
    19   return true;
     17        // psphotPetrosianVisualProfileByAngle (radii, fluxes, radius);
     18
     19        // warn on NAN?
     20        profile->isophotalRadii->data.F32[i] = radius;
     21    }
     22    return true;
    2023}
    2124
    … …  
    7982
    8083        // sort the flux by the radius
    81         psphotPetrosianSortPair (radius, flux);
     84        pmSourceRadialProfileSortPair (radius, flux);
    8285
    8386        int nOut = 0;
    … …  
    128131            // XXX is there a macro in psLib that does this interpolation?
    129132            if (i == 0) {
    130                 psLogMsg ("psphot", 3, "bogus radial profile for ..., skipping");
     133                psTrace ("psphot", 4, "bogus radial profile for source at %f, %f, skipping", source->peak->xf, source->peak->yf);
    131134                psFree (fluxBinned);
    132135                psFree (radiusBinned);

  • branches/eam_branches/20090715/psphot/src/psphotReadout.c

    r25409 r25452  
    212212finish:
    213213
    214     psphotPetrosianAnalysis (readout, sources, recipe);
     214    // XXX drop this : test dev version
     215    // psphotPetrosianAnalysis (readout, sources, recipe);
    215216
    216217    // plot positive sources

  • branches/eam_branches/20090715/psphot/src/psphotSourceSize.c

    r25433 r25452  
    334334        bool isPSF = (fabs(nSigma) < options->nSigmaApResid) && (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY);
    335335        if (isPSF) {
     336            if (Mxx > 3.0) {
     337                fprintf (stderr, "!");
     338            }
    336339            Npsf ++;
    337340            continue;
    … …  
    354357
    355358        // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)?
    356         bool isEXT = (nSigma > options->nSigmaApResid) && (Mxx > psfClump->X) && (Myy > psfClump->Y);
     359        bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y));
    357360        if (isEXT) {
    358361            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    … …  
    361364        }
    362365
     366        fprintf (stderr, "miss %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma);
    363367        Nmiss ++;
    364368    }

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

    r25433 r25452  
    20942094}
    20952095
    2096 bool psphotVisualShowPetrosian (pmSource *source, pmPetrosian *petrosian) {
    2097 
    2098     KiiOverlay overlay[2];
    2099 
    2100     if (source == NULL) return true;
     2096bool psphotVisualShowPetrosians (psArray *sources) {
     2097
     2098    int Noverlay, NOVERLAY;
     2099    KiiOverlay *overlay;
    21012100
    21022101    // if (!pmVisualIsVisual()) return true;
    … …  
    21052104    if (kapa == -1) return false;
    21062105
    2107     if (kapa == -1) {
    2108         fprintf (stderr, "kapa not opened, skipping\n");
    2109         return false;
    2110     }
    2111 
    2112     overlay[0].type = KII_OVERLAY_CIRCLE;
    2113     overlay[0].x = source->peak->xf;
    2114     overlay[0].y = source->peak->yf;
    2115     overlay[0].dx = 2.0*petrosian->petrosianRadius;
    2116     overlay[0].dy = 2.0*petrosian->petrosianRadius*petrosian->axes.minor/petrosian->axes.major;
    2117     overlay[0].angle = petrosian->axes.theta * PS_DEG_RAD;
    2118     overlay[0].text = NULL;
    2119 
    2120     overlay[1].type = KII_OVERLAY_CIRCLE;
    2121     overlay[1].x = source->peak->xf;
    2122     overlay[1].y = source->peak->yf;
    2123     overlay[1].dx = 4.0*petrosian->petrosianRadius;
    2124     overlay[1].dy = 4.0*petrosian->petrosianRadius*petrosian->axes.minor/petrosian->axes.major;
    2125     overlay[1].angle = petrosian->axes.theta * PS_DEG_RAD;
    2126     overlay[1].text = NULL;
    2127 
    2128     KiiLoadOverlay (kapa, overlay, 2, "red");
     2106    Noverlay = 0;
     2107    NOVERLAY = 100;
     2108    ALLOCATE (overlay, KiiOverlay, NOVERLAY);
     2109
     2110    for (int i = 0; i < sources->n; i++) {
     2111        pmSource *source = sources->data[i];
     2112
     2113        if (!source) continue;
     2114        if (!source->extpars) continue;
     2115        if (!source->extpars->profile) continue;
     2116        if (!source->extpars->petrosian_80) continue;
     2117
     2118        pmSourceRadialProfile *profile = source->extpars->profile;
     2119        pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80;
     2120
     2121        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     2122        overlay[Noverlay].x = source->peak->xf;
     2123        overlay[Noverlay].y = source->peak->yf;
     2124        overlay[Noverlay].dx = 2.0*petrosian->radius;
     2125        overlay[Noverlay].dy = 2.0*petrosian->radius*profile->axes.minor/profile->axes.major;
     2126        overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD;
     2127        overlay[Noverlay].text = NULL;
     2128        Noverlay ++;
     2129        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     2130
     2131        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     2132        overlay[Noverlay].x = source->peak->xf;
     2133        overlay[Noverlay].y = source->peak->yf;
     2134        overlay[Noverlay].dx = 4.0*petrosian->radius;
     2135        overlay[Noverlay].dy = 4.0*petrosian->radius*profile->axes.minor/profile->axes.major;
     2136        overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD;
     2137        overlay[Noverlay].text = NULL;
     2138        Noverlay ++;
     2139        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     2140    }
     2141
     2142    KiiLoadOverlay (kapa, overlay, Noverlay, "red");
     2143    FREE (overlay);
     2144
     2145    // pause and wait for user input:
     2146    // continue, save (provide name), ??
     2147    char key[10];
     2148    fprintf (stdout, "[c]ontinue? ");
     2149    if (!fgets(key, 8, stdin)) {
     2150        psWarning("Unable to read option");
     2151    }
    21292152
    21302153    return true;
    … …  
    21522175
    21532176# if (0)
    2154     // *** make a histogram of the source counts in the x and y directions
    2155     psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0);
    2156     psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0);
    2157     psVectorHistogram (nX, xFaint, NULL, NULL, 0);
    2158     psVectorHistogram (nY, yFaint, NULL, NULL, 0);
    2159     psVector *dX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
    2160     psVector *vX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
    2161     psVector *dY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
    2162     psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
    2163     for (int i = 0; i < nX->nums->n; i++) {
    2164         dX->data.F32[i] = nX->nums->data.S32[i];
    2165         vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]);
    2166     }
    2167     for (int i = 0; i < nY->nums->n; i++) {
    2168         dY->data.F32[i] = nY->nums->data.S32[i];
    2169         vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]);
    2170     }
    2171 
    2172     graphdata.color = KapaColorByName ("black");
    2173     graphdata.ptype = 0;
    2174     graphdata.size = 0.0;
    2175     graphdata.style = 0;
    2176     KapaPrepPlot (myKapa, dX->n, &graphdata);
    2177     KapaPlotVector (myKapa, dX->n, dX->data.F32, "x");
    2178     KapaPlotVector (myKapa, vX->n, vX->data.F32, "y");
    2179 
    2180     psFree (nX);
    2181     psFree (dX);
    2182     psFree (vX);
    2183 
    2184     psFree (nY);
    2185     psFree (dY);
    2186     psFree (vY);
     2177// *** make a histogram of the source counts in the x and y directions
     2178psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0);
     2179psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0);
     2180psVectorHistogram (nX, xFaint, NULL, NULL, 0);
     2181psVectorHistogram (nY, yFaint, NULL, NULL, 0);
     2182psVector *dX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
     2183psVector *vX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
     2184psVector *dY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
     2185psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
     2186for (int i = 0; i < nX->nums->n; i++) {
     2187    dX->data.F32[i] = nX->nums->data.S32[i];
     2188    vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]);
     2189}
     2190for (int i = 0; i < nY->nums->n; i++) {
     2191    dY->data.F32[i] = nY->nums->data.S32[i];
     2192    vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]);
     2193}
     2194
     2195graphdata.color = KapaColorByName ("black");
     2196graphdata.ptype = 0;
     2197graphdata.size = 0.0;
     2198graphdata.style = 0;
     2199KapaPrepPlot (myKapa, dX->n, &graphdata);
     2200KapaPlotVector (myKapa, dX->n, dX->data.F32, "x");
     2201KapaPlotVector (myKapa, vX->n, vX->data.F32, "y");
     2202
     2203psFree (nX);
     2204psFree (dX);
     2205psFree (vX);
     2206
     2207psFree (nY);
     2208psFree (dY);
     2209psFree (vY);
    21872210
    21882211
Note: See TracChangeset for help on using the changeset viewer.