IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25452


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
Files:
24 edited
3 moved

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090715/ippconfig/mosaic2/camera.config

    r19019 r25452  
    9090  CMF.XSRC STR {CHIP.NAME}.xsrc # use .PSF and .EXT?
    9191  CMF.XFIT STR {CHIP.NAME}.xfit # use .PSF and .EXT?
     92  CMF.DETEFF STR {CHIP.NAME}.deteff
    9293
    9394  PSF.HEAD  STR {CHIP.NAME}.hdr

  • branches/eam_branches/20090715/ippconfig/mosaic2/psphot.config

    r25395 r25452  
    2424USE_FOOTPRINTS                      BOOL  TRUE            # use new pmFootprint peak packaging
    2525
     26PSF.TREND.MODE                      STR   MAP
     27PSF.TREND.NX                        S32   1
     28PSF.TREND.NY                        S32   1
     29
    2630PEAKS_NSIGMA_LIMIT   F32  15.0            # peak significance threshold
    2731MOMENTS_SN_MIN       F32  10.0           # min S/N to measure moments
    2832PSF_SN_LIM           F32  20.0            # minimum S/N for stars used for PSF model
    29 FULL_FIT_SN_LIM      F32  20.0
    30 EXT_MIN_SN           F32  20.0           # fit galaxies above this S/N limit
     33FULL_FIT_SN_LIM      F32  10.0
     34EXT_MIN_SN           F32  10.0           # fit galaxies above this S/N limit
    3135AP_MIN_SN            F32  10.0
    3236
    … …  
    4246PSPHOT.EXT.NSIGMA.LIMIT             F32   3.0  # sources with extNsigma greater that this get tagged as likely extended sources
    4347PSPHOT.EXT.NSIGMA.MOMENTS           F32   2.0  # sources with extNsigma greater that this get tagged as likely extended sources
     48
     49EXTENDED_SOURCE_SN_LIM              F32   3.0
     50EXTENDED_SOURCE_ANALYSIS            BOOL  TRUE  # perform any of the aperture-like measurements?
     51EXTENDED_SOURCE_PETROSIAN           BOOL  TRUE

  • branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c

    r25352 r25452  
    443443    float dSysBright = psVectorSystematicError (bright, brightErr, 0.1);
    444444    fprintf (stderr, "bright systematic error: %f\n", dSysBright);
     445    psFree(bright);
     446    psFree(brightErr);
    445447
    446448    // XXX test dump of fitted model (dump when tracing?)
    … …  
    11781180
    11791181    // free local allocations
     1182    psFree (mask);
     1183    psFree (chisq);
     1184    psFree (stats);
     1185    psFree (index);
     1186
    11801187    return (sqrt(S2guess));
    11811188}

  • branches/eam_branches/20090715/psModules/src/objects/pmPetrosian.c

    r25409 r25452  
    99 */
    1010
    11 # include "psphotInternal.h"
     11# include "pmPetrosian.h"
    1212
    1313static void pmPetrosianFree(pmPetrosian *petrosian)
    … …  
    5656}
    5757
    58 bool psphotPetrosianFreeVectors(pmPetrosian *petrosian) {
     58bool pmPetrosianFreeVectors(pmPetrosian *petrosian) {
    5959
    6060    psFree(petrosian->radii);
    … …  
    100100}
    101101
    102 bool psphotPetrosianSortPair (psVector *index, psVector *extra) {
     102bool pmPetrosianSortPair (psVector *index, psVector *extra) {
    103103
    104104    // sort the vector set by the radius

  • branches/eam_branches/20090715/psModules/src/objects/pmPetrosian.h

    r25433 r25452  
    1010#ifndef PM_PETROSIAN_H
    1111#define PM_PETROSIAN_H
     12
     13/// @addtogroup Objects Object Detection / Analysis Functions
     14/// @{
    1215
    1316typedef struct {
    … …  
    3336
    3437pmPetrosian *pmPetrosianAlloc();
    35 bool psphotPetrosianFreeVectors(pmPetrosian *petrosian);
    3638
    37 bool psphotPetrosianProfile (pmReadout *readout, pmSource *source, float skynoise);
    38 bool psphotRadialProfilesByAngles (pmSource *source, pmPetrosian *petro, int Nsec, float Rmax);
    39 float psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax);
    40 bool psphotRadiiFromProfiles (pmSource *source, pmPetrosian *petrosian, float fluxMin, float fluxMax);
    41 bool psphotEllipticalProfile (pmSource *source, pmPetrosian *petrosian);
    42 bool psphotEllipticalContour (pmSource *source, pmPetrosian *petrosian);
    43 bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax, float skynoise);
    44 bool psphotPetrosianStats (pmSource *source, pmPetrosian *petrosian);
    45 
    46 bool psphotPetrosianSortPair (psVector *index, psVector *extra);
    47 
    48 
    49 bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux);
    50 bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef);
    51 bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian);
    52 
    53 bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin,
    54                                  psVector *refRadius, psVector *meanSB,
    55                                  psVector *petRatio, psVector *petRatioErr, psVector *fluxSum,
    56                                  float petRadius, float ratioForRadius,
    57                                  float petFlux, float radiusForFlux);
     39bool pmPetrosianFreeVectors(pmPetrosian *petrosian);
     40bool pmPetrosianSortPair (psVector *index, psVector *extra);
    5841
    5942/// @}
    6043
    61 
    6244# endif /* PM_PETROSIAN_H */

  • branches/eam_branches/20090715/psModules/src/objects/pmSourceExtendedPars.c

    r23487 r25452  
    1717#endif
    1818
    19 #include <stdio.h>
    20 #include <math.h>
    21 #include <string.h>
     19// #include <stdio.h>
     20// #include <math.h>
     21// #include <string.h>
    2222#include <pslib.h>
    23 #include "pmHDU.h"
    24 #include "pmFPA.h"
    25 #include "pmFPAMaskWeight.h"
    26 #include "pmSpan.h"
    27 #include "pmFootprint.h"
    28 #include "pmPeaks.h"
    29 #include "pmMoments.h"
    30 #include "pmResiduals.h"
    31 #include "pmGrowthCurve.h"
    32 #include "pmTrend2D.h"
    33 #include "pmPSF.h"
    34 #include "pmModel.h"
    35 #include "pmSource.h"
    36 
     23// #include "pmHDU.h"
     24// #include "pmFPA.h"
     25// #include "pmFPAMaskWeight.h"
     26// #include "pmSpan.h"
     27// #include "pmFootprint.h"
     28// #include "pmPeaks.h"
     29// #include "pmMoments.h"
     30// #include "pmResiduals.h"
     31// #include "pmGrowthCurve.h"
     32// #include "pmTrend2D.h"
     33// #include "pmPSF.h"
     34// #include "pmModel.h"
     35// #include "pmSource.h"
     36#include "pmSourceExtendedPars.h"
     37
     38// *** pmSourceRadialProfile describes the radial profile of a source in elliptical contours, and
     39// intermediate data used to measure the profile
     40static void pmSourceRadialProfileFree(pmSourceRadialProfile *profile)
     41{
     42    if (!profile) {
     43        return;
     44    }
     45    psFree(profile->radii);
     46    psFree(profile->fluxes);
     47    psFree(profile->theta);
     48    psFree(profile->isophotalRadii);
     49
     50    psFree(profile->radiusElliptical);
     51    psFree(profile->fluxElliptical);
     52
     53    psFree(profile->binSB);
     54    psFree(profile->binSBstdev);
     55    psFree(profile->binSBerror);
     56
     57    psFree(profile->radialBins);
     58    psFree(profile->area);
     59}
     60
     61pmSourceRadialProfile *pmSourceRadialProfileAlloc()
     62{
     63    pmSourceRadialProfile *profile = (pmSourceRadialProfile *)psAlloc(sizeof(pmSourceRadialProfile));
     64    psMemSetDeallocator(profile, (psFreeFunc) pmSourceRadialProfileFree);
     65
     66    profile->radii = NULL;
     67    profile->fluxes = NULL;
     68    profile->theta = NULL;
     69    profile->isophotalRadii = NULL;
     70
     71    profile->radiusElliptical = NULL;
     72    profile->fluxElliptical = NULL;
     73
     74    profile->binSB = NULL;
     75    profile->binSBstdev = NULL;
     76    profile->binSBerror = NULL;
     77
     78    profile->radialBins = NULL;
     79    profile->area = NULL;
     80
     81    return profile;
     82}
     83
     84bool psMemCheckSourceRadialProfile(psPtr ptr)
     85{
     86    PS_ASSERT_PTR(ptr, false);
     87    return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceRadialProfileFree);
     88}
     89
     90
     91// *** pmSourceRadialProfileFreeVectors frees the intermediate data values
     92bool pmSourceRadialProfileFreeVectors(pmSourceRadialProfile *profile) {
     93
     94    psFree(profile->radii);
     95    psFree(profile->fluxes);
     96    psFree(profile->theta);
     97    psFree(profile->isophotalRadii);
     98
     99    psFree(profile->radiusElliptical);
     100    psFree(profile->fluxElliptical);
     101
     102    // psFree(profile->binSB);
     103    // psFree(profile->binSBstdev);
     104    // psFree(profile->binSBerror);
     105   
     106    // psFree(profile->radialBins);
     107    psFree(profile->area);
     108
     109    profile->radii = NULL;
     110    profile->fluxes = NULL;
     111    profile->theta = NULL;
     112    profile->isophotalRadii = NULL;
     113
     114    profile->radiusElliptical = NULL;
     115    profile->fluxElliptical = NULL;
     116
     117    // profile->binSB = NULL;
     118    // profile->binSBstdev = NULL;
     119    // profile->binSBerror = NULL;
     120   
     121    // profile->radialBins = NULL;
     122    profile->area = NULL;
     123
     124    return true;
     125}
     126
     127// *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors
     128# define COMPARE_INDEX(A,B) (index->data.F32[A] < index->data.F32[B])
     129# define SWAP_INDEX(TYPE,A,B) { \
     130  float tmp; \
     131  if (A != B) { \
     132    tmp = index->data.F32[A]; \
     133    index->data.F32[A] = index->data.F32[B]; \
     134    index->data.F32[B] = tmp; \
     135    tmp = extra->data.F32[A]; \
     136    extra->data.F32[A] = extra->data.F32[B]; \
     137    extra->data.F32[B] = tmp; \
     138  } \
     139}
     140
     141bool pmSourceRadialProfileSortPair (psVector *index, psVector *extra) {
     142
     143    // sort the vector set by the radius
     144    PSSORT (index->n, COMPARE_INDEX, SWAP_INDEX, NONE);
     145    return true;
     146}
     147
     148// *** pmSourceExtendedPars describes the possible collection of extended flux measurements for a source
    37149static void pmSourceExtendedParsFree (pmSourceExtendedPars *pars) {
    38150    if (!pars) return;
    39151
    40152    psFree(pars->profile);
    41     psFree(pars->annuli);
    42     psFree(pars->isophot);
    43     psFree(pars->petrosian);
    44     psFree(pars->kron);
     153    psFree(pars->petrosian_50);
     154    psFree(pars->petrosian_80);
    45155    return;
    46156}
    … …  
    51161
    52162    pars->profile = NULL;
    53     pars->annuli = NULL;
    54     pars->isophot = NULL;
    55     pars->petrosian = NULL;
    56     pars->kron = NULL;
     163    pars->petrosian_50 = NULL;
     164    pars->petrosian_80 = NULL;
    57165
    58166    return pars;
    … …  
    66174
    67175
    68 static void pmSourceRadialProfileFree (pmSourceRadialProfile *profile) {
    69     if (!profile) return;
    70 
    71     psFree(profile->radius);
    72     psFree(profile->flux);
    73     psFree(profile->variance);
     176// *** pmSourceExtendedFlux describes the flux within an elliptical aperture of some kind
     177static void pmSourceExtendedFluxFree (pmSourceExtendedFlux *flux) {
     178    if (!flux) return;
    74179    return;
    75180}
    76181
    77 pmSourceRadialProfile *pmSourceRadialProfileAlloc (void) {
    78 
    79     pmSourceRadialProfile *profile = (pmSourceRadialProfile *) psAlloc(sizeof(pmSourceRadialProfile));
    80     psMemSetDeallocator(profile, (psFreeFunc) pmSourceRadialProfileFree);
    81 
    82     profile->radius = NULL;
    83     profile->flux = NULL;
    84     profile->variance = NULL;
    85 
    86     return profile;
    87 }
    88 
    89 bool psMemCheckSourceRadialProfile(psPtr ptr)
     182pmSourceExtendedFlux *pmSourceExtendedFluxAlloc (void) {
     183
     184    pmSourceExtendedFlux *flux = (pmSourceExtendedFlux *) psAlloc(sizeof(pmSourceExtendedFlux));
     185    psMemSetDeallocator(flux, (psFreeFunc) pmSourceExtendedFluxFree);
     186
     187    flux->flux = 0.0;
     188    flux->fluxErr = 0.0;
     189    flux->radius = 0.0;
     190    flux->radiusErr = 0.0;
     191
     192    return flux;
     193}
     194
     195
     196bool psMemCheckSourceExtendedFlux(psPtr ptr)
    90197{
    91198    PS_ASSERT_PTR(ptr, false);
    92     return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceRadialProfileFree);
    93 }
    94 
    95 
    96 static void pmSourceIsophotalValuesFree (pmSourceIsophotalValues *isophot) {
    97     if (!isophot) return;
    98     return;
    99 }
    100 
    101 pmSourceIsophotalValues *pmSourceIsophotalValuesAlloc (void) {
    102 
    103     pmSourceIsophotalValues *isophot = (pmSourceIsophotalValues *) psAlloc(sizeof(pmSourceIsophotalValues));
    104     psMemSetDeallocator(isophot, (psFreeFunc) pmSourceIsophotalValuesFree);
    105 
    106     isophot->mag = 0.0;
    107     isophot->magErr = 0.0;
    108     isophot->rad = 0.0;
    109     isophot->radErr = 0.0;
    110 
    111     return isophot;
    112 }
    113 
    114 
    115 bool psMemCheckSourceIsophotalValues(psPtr ptr)
    116 {
    117     PS_ASSERT_PTR(ptr, false);
    118     return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceIsophotalValuesFree);
    119 }
    120 
    121 
    122 static void pmSourcePetrosianValuesFree (pmSourcePetrosianValues *petrosian) {
    123     if (!petrosian) return;
    124     return;
    125 }
    126 
    127 pmSourcePetrosianValues *pmSourcePetrosianValuesAlloc (void) {
    128 
    129     pmSourcePetrosianValues *petrosian = (pmSourcePetrosianValues *) psAlloc(sizeof(pmSourcePetrosianValues));
    130     psMemSetDeallocator(petrosian, (psFreeFunc) pmSourcePetrosianValuesFree);
    131 
    132     petrosian->mag = 0.0;
    133     petrosian->magErr = 0.0;
    134     petrosian->rad = 0.0;
    135     petrosian->radErr = 0.0;
    136 
    137     return petrosian;
    138 }
    139 
    140 
    141 bool psMemCheckSourcePetrosianValues(psPtr ptr)
    142 {
    143     PS_ASSERT_PTR(ptr, false);
    144     return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourcePetrosianValuesFree);
    145 }
    146 
    147 static void pmSourceKronValuesFree (pmSourceKronValues *kron) {
    148     if (!kron) return;
    149     return;
    150 }
    151 
    152 pmSourceKronValues *pmSourceKronValuesAlloc (void) {
    153 
    154     pmSourceKronValues *kron = (pmSourceKronValues *) psAlloc(sizeof(pmSourceKronValues));
    155     psMemSetDeallocator(kron, (psFreeFunc) pmSourceKronValuesFree);
    156 
    157     kron->mag = 0.0;
    158     kron->magErr = 0.0;
    159     kron->rad = 0.0;
    160     kron->radErr = 0.0;
    161 
    162     return kron;
    163 }
    164 
    165 
    166 bool psMemCheckSourceKronValues(psPtr ptr)
    167 {
    168     PS_ASSERT_PTR(ptr, false);
    169     return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceKronValuesFree);
    170 }
    171 
    172 
    173 static void pmSourceAnnuliFree (pmSourceAnnuli *annuli) {
    174     if (!annuli) return;
    175 
    176     psFree (annuli->flux);
    177     psFree (annuli->fluxErr);
    178     psFree (annuli->fluxVar);
    179 
    180     return;
    181 }
    182 
    183 pmSourceAnnuli *pmSourceAnnuliAlloc (void) {
    184 
    185     pmSourceAnnuli *annuli = (pmSourceAnnuli *) psAlloc(sizeof(pmSourceAnnuli));
    186     psMemSetDeallocator(annuli, (psFreeFunc) pmSourceAnnuliFree);
    187 
    188     annuli->flux = NULL;
    189     annuli->fluxErr = NULL;
    190     annuli->fluxVar = NULL;
    191 
    192     return annuli;
    193 }
    194 
    195 
    196 bool psMemCheckSourceAnnuli(psPtr ptr)
    197 {
    198     PS_ASSERT_PTR(ptr, false);
    199     return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceAnnuliFree);
    200 }
     199    return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceExtendedFluxFree);
     200}

  • branches/eam_branches/20090715/psModules/src/objects/pmSourceExtendedPars.h

    r23487 r25452  
    1515
    1616typedef struct {
    17   psVector *radius;
    18   psVector *flux;
    19   psVector *variance;
     17    psArray  *radii;                    // radii for raw radial profiles at evenly-spaced angles
     18    psArray  *fluxes;                   // fluxes measured at above radii
     19    psVector *theta;                    // angles corresponding to above radial profiles
     20    psVector *isophotalRadii;           // isophotal radius for the above angles
     21
     22    psVector *radiusElliptical;         // normalized radial coordinates for all relevant pixels
     23    psVector *fluxElliptical;           // flux for the above radial coordinates
     24
     25    psVector *binSB;                    // mean surface brightness within radial bins
     26    psVector *binSBstdev;               // scatter of mean surface brightness within radial bins
     27    psVector *binSBerror;               // formal error on mean surface brightness within radial bins
     28
     29    psVector *radialBins;               // radii corresponding to above binnedBlux
     30    psVector *area;                     // differential area of the non-overlapping radial bins
     31
     32    psEllipseAxes axes;                 // shape of elliptical contour
    2033} pmSourceRadialProfile;
    2134
    2235typedef struct {
    23   psVector *flux;
    24   psVector *fluxErr;
    25   psVector *fluxVar;
    26 } pmSourceAnnuli;
    27 
    28 typedef struct {
    29   float mag;
    30   float magErr;
    31   float rad;
    32   float radErr;
    33 } pmSourceIsophotalValues;
    34 
    35 typedef struct {
    36   float mag;
    37   float magErr;
    38   float rad;
    39   float radErr;
    40 } pmSourcePetrosianValues;
    41 
    42 typedef struct {
    43   float mag;
    44   float magErr;
    45   float rad;
    46   float radErr;
    47 } pmSourceKronValues;
     36  float flux;
     37  float fluxErr;
     38  float radius;
     39  float radiusErr;
     40} pmSourceExtendedFlux;
    4841
    4942typedef struct {
    5043  pmSourceRadialProfile   *profile;
    51   pmSourceAnnuli          *annuli;
    52   pmSourceIsophotalValues *isophot;
    53   pmSourcePetrosianValues *petrosian;
    54   pmSourceKronValues      *kron;
     44  pmSourceExtendedFlux    *petrosian_50;
     45  pmSourceExtendedFlux    *petrosian_80;
    5546} pmSourceExtendedPars;
    5647
    57 pmSourceExtendedPars *pmSourceExtendedParsAlloc(void);
     48// *** pmSourceRadialProfile describes the radial profile of a source in elliptical contours, and
     49// intermediate data used to measure the profile
     50pmSourceRadialProfile *pmSourceRadialProfileAlloc();
     51bool psMemCheckSourceRadialProfile(psPtr ptr);
     52
     53// *** pmSourceRadialProfileFreeVectors frees the intermediate data values
     54bool pmSourceRadialProfileFreeVectors(pmSourceRadialProfile *profile);
     55
     56// *** pmSourceExtendedPars describes the possible collection of extended flux measurements for a source
     57pmSourceExtendedPars *pmSourceExtendedParsAlloc (void);
    5858bool psMemCheckSourceExtendedPars(psPtr ptr);
    59 pmSourceRadialProfile *pmSourceRadialProfileAlloc(void);
    60 bool psMemCheckSourceRadialProfile(psPtr ptr);
    61 pmSourceIsophotalValues *pmSourceIsophotalValuesAlloc(void);
    62 bool psMemCheckSourceIsophotalValues(psPtr ptr);
    63 pmSourcePetrosianValues *pmSourcePetrosianValuesAlloc(void);
    64 bool psMemCheckSourcePetrosianValues(psPtr ptr);
    65 pmSourceKronValues *pmSourceKronValuesAlloc(void);
    66 bool psMemCheckSourceKronValues(psPtr ptr);
    67 pmSourceAnnuli *pmSourceAnnuliAlloc(void);
    68 bool psMemCheckSourceAnnuli(psPtr ptr);
     59
     60// *** pmSourceExtendedFlux describes the flux within an elliptical aperture of some kind
     61pmSourceExtendedFlux *pmSourceExtendedFluxAlloc(void);
     62bool psMemCheckSourceExtendedFlux(psPtr ptr);
     63
     64// *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors
     65bool pmSourceRadialProfileSortPair(psVector *index, psVector *extra);
    6966
    7067/// @}

  • branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c

    r24403 r25452  
    353353
    354354    // which extended source analyses should we perform?
    355     bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    356     bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    357     bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    358     bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
     355    // bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     356    // bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
     357    // bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
     358    // bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
    359359
    360360    psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    … …  
    396396        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    397397
     398# if (0)
    398399        // Petrosian measurements
    399400        // XXX insert header data: petrosian ref radius, flux ratio
    … …  
    476477        }
    477478
     479# endif
    478480        psArrayAdd (table, 100, row);
    479481        psFree (row);

  • branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_PS1_CAL_0.c

    r21516 r25452  
    349349
    350350    // which extended source analyses should we perform?
    351     bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    352     bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    353     bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    354     bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
     351    // bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     352    // bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
     353    // bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
     354    // bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
    355355
    356356    psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    … …  
    410410        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    411411
     412# if (0)
    412413        // Petrosian measurements
    413414        // XXX insert header data: petrosian ref radius, flux ratio
    … …  
    501502            }
    502503        }
     504# endif
    503505
    504506        psArrayAdd (table, 100, row);

  • branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_PS1_DEV_1.c

    r21516 r25452  
    300300
    301301    // which extended source analyses should we perform?
    302     bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    303     bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    304     bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    305     bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
     302    // bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     303    // bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
     304    // bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
     305    // bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
    306306
    307307    psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    … …  
    343343        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    344344
     345// XXX disable these outputs until we clean up the names
     346# if (0)
    345347        // Petrosian measurements
    346348        // XXX insert header data: petrosian ref radius, flux ratio
    … …  
    422424            }
    423425        }
     426# endif
    424427
    425428        psArrayAdd (table, 100, row);

  • 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.