IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27818


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

APIs mods to pmSourcesWrite_CMF_ XSRC functions; added more (nearly) complete extended source parameters; fix bug in psf for polynomial of order 0; updates to the pmSourceExtendedPars structure and related

Location:
trunk/psModules/src/objects
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/pmPSFtryModel.c

    r27330 r27818  
    107107    // as we loop over orders, we need to refer to the initial selection, but we modify the
    108108    // option values to match the current guess: save the max values here:
    109     int Nx = options->psfTrendNx;
    110     int Ny = options->psfTrendNy;
     109    int Nx = (options->psfTrendMode == PM_TREND_MAP) ? options->psfTrendNx : options->psfTrendNx + 1;
     110    int Ny = (options->psfTrendMode == PM_TREND_MAP) ? options->psfTrendNy : options->psfTrendNy + 1;
    111111    for (int i = orderMin; i <= orderMax; i++) {
    112112

  • trunk/psModules/src/objects/pmSourceExtendedPars.c

    r25754 r27818  
    3636#include "pmSourceExtendedPars.h"
    3737
    38 // *** pmSourceRadialProfile describes the radial profile of a source in elliptical contours, and
    39 // intermediate data used to measure the profile
     38// pmSourceRadialFlux carries the raw radial flux information, including angular bins
     39static void pmSourceRadialFluxFree(pmSourceRadialFlux *flux)
     40{
     41    if (!flux) return;
     42    psFree(flux->radii);
     43    psFree(flux->fluxes);
     44    psFree(flux->theta);
     45    psFree(flux->isophotalRadii);
     46}
     47
     48pmSourceRadialFlux *pmSourceRadialFluxAlloc()
     49{
     50    pmSourceRadialFlux *flux = (pmSourceRadialFlux *)psAlloc(sizeof(pmSourceRadialFlux));
     51    psMemSetDeallocator(flux, (psFreeFunc) pmSourceRadialFluxFree);
     52
     53    flux->radii = NULL;
     54    flux->fluxes = NULL;
     55    flux->theta = NULL;
     56    flux->isophotalRadii = NULL;
     57
     58    return flux;
     59}
     60
     61bool psMemCheckSourceRadialFlux(psPtr ptr)
     62{
     63    PS_ASSERT_PTR(ptr, false);
     64    return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceRadialFluxFree);
     65}
     66
     67// pmSourceEllipticalFlux carries the elliptical renormalized radial flux info
     68static void pmSourceEllipticalFluxFree(pmSourceEllipticalFlux *flux)
     69{
     70    if (!flux) return;
     71    psFree(flux->radiusElliptical);
     72    psFree(flux->fluxElliptical);
     73}
     74
     75pmSourceEllipticalFlux *pmSourceEllipticalFluxAlloc()
     76{
     77    pmSourceEllipticalFlux *flux = (pmSourceEllipticalFlux *)psAlloc(sizeof(pmSourceEllipticalFlux));
     78    psMemSetDeallocator(flux, (psFreeFunc) pmSourceEllipticalFluxFree);
     79
     80    flux->radiusElliptical = NULL;
     81    flux->fluxElliptical = NULL;
     82
     83    return flux;
     84}
     85
     86bool psMemCheckSourceEllipticalFlux(psPtr ptr)
     87{
     88    PS_ASSERT_PTR(ptr, false);
     89    return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceEllipticalFluxFree);
     90}
     91
     92// pmSourceRadialProfile defines flux information in radial bins
    4093static void pmSourceRadialProfileFree(pmSourceRadialProfile *profile)
    4194{
    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 
     95    if (!profile) return;
    5396    psFree(profile->binSB);
    5497    psFree(profile->binSBstdev);
    5598    psFree(profile->binSBerror);
    56 
     99    psFree(profile->binSum);
     100    psFree(profile->binFill);
    57101    psFree(profile->radialBins);
    58102    psFree(profile->area);
    … …  
    63107    pmSourceRadialProfile *profile = (pmSourceRadialProfile *)psAlloc(sizeof(pmSourceRadialProfile));
    64108    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;
    73109
    74110    profile->binSB = NULL;
    75111    profile->binSBstdev = NULL;
    76112    profile->binSBerror = NULL;
    77 
     113    profile->binSum = NULL;
     114    profile->binFill = NULL;
    78115    profile->radialBins = NULL;
    79116    profile->area = NULL;
    80 
    81117    return profile;
    82118}
    … …  
    88124}
    89125
    90 
    91 // *** pmSourceRadialProfileFreeVectors frees the intermediate data values
     126# if (0)
     127// pmSourceRadialProfileFreeVectors frees the intermediate data values
    92128bool pmSourceRadialProfileFreeVectors(pmSourceRadialProfile *profile) {
    93129
    … …  
    124160    return true;
    125161}
     162# endif
    126163
    127164// *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors
    … …  
    150187    if (!pars) return;
    151188
    152     psFree(pars->profile);
    153     psFree(pars->petrosian_50);
    154     psFree(pars->petrosian_80);
     189    psFree(pars->radFlux);
     190    psFree(pars->ellipticalFlux);
     191    psFree(pars->radProfile);
     192    psFree(pars->petProfile);
    155193    return;
    156194}
    … …  
    160198    psMemSetDeallocator(pars, (psFreeFunc) pmSourceExtendedParsFree);
    161199
    162     pars->profile = NULL;
    163     pars->petrosian_50 = NULL;
    164     pars->petrosian_80 = NULL;
    165 
     200    pars->radFlux = NULL;
     201    pars->ellipticalFlux = NULL;
     202    pars->radProfile = NULL;
     203    pars->petProfile = NULL;
     204
     205    pars->petrosianFlux = NAN;
     206    pars->petrosianFluxErr = NAN;
     207    pars->petrosianRadius = NAN;
     208    pars->petrosianRadiusErr = NAN;
     209    pars->petrosianR90 = NAN;
     210    pars->petrosianR90Err = NAN;
     211    pars->petrosianR50 = NAN;
     212    pars->petrosianR50Err = NAN;
    166213    return pars;
    167214}

  • trunk/psModules/src/objects/pmSourceExtendedPars.h

    r25754 r27818  
    1919    psVector *theta;                    // angles corresponding to above radial profiles
    2020    psVector *isophotalRadii;           // isophotal radius for the above angles
     21} pmSourceRadialFlux;
    2122
     23typedef struct {
    2224    psVector *radiusElliptical;         // normalized radial coordinates for all relevant pixels
    2325    psVector *fluxElliptical;           // flux for the above radial coordinates
     26} pmSourceEllipticalFlux;
    2427
     28typedef struct {
    2529    psVector *binSB;                    // mean surface brightness within radial bins
    2630    psVector *binSBstdev;               // scatter of mean surface brightness within radial bins
    2731    psVector *binSBerror;               // formal error on mean surface brightness within radial bins
    28 
    29     psVector *radialBins;               // radii corresponding to above binnedBlux
     32    psVector *binSum;                   // sum of flux within radial bins
     33    psVector *binFill;                  // fraction of area actually lit
     34    psVector *radialBins;               // radii corresponding to above binnedFlux
    3035    psVector *area;                     // differential area of the non-overlapping radial bins
    31 
    32     psEllipseAxes axes;                 // shape of elliptical contour
    3336} pmSourceRadialProfile;
    3437
    3538typedef struct {
    36   float flux;
    37   float fluxErr;
    38   float radius;
    39   float radiusErr;
     39    float flux;
     40    float fluxErr;
     41    float radius;
     42    float radiusErr;
    4043} pmSourceExtendedFlux;
    4144
    4245typedef struct {
    43   pmSourceRadialProfile   *profile;
    44   pmSourceExtendedFlux    *petrosian_50;
    45   pmSourceExtendedFlux    *petrosian_80;
     46    pmSourceRadialFlux     *radFlux;        // raw radial flux information
     47    pmSourceEllipticalFlux *ellipticalFlux; // flux for elliptically-renormalized radii
     48    pmSourceRadialProfile  *radProfile;     // surface brightness profile in specified fixed bins
     49    pmSourceRadialProfile  *petProfile;     // surface brightness profile in petrosian bins
     50    psEllipseAxes axes;                     // shape of elliptical contour
     51    float petrosianFlux;
     52    float petrosianFluxErr;
     53    float petrosianRadius;
     54    float petrosianRadiusErr;
     55    float petrosianR90;
     56    float petrosianR90Err;
     57    float petrosianR50;
     58    float petrosianR50Err;
    4659} pmSourceExtendedPars;
     60
     61pmSourceRadialFlux *pmSourceRadialFluxAlloc();
     62bool psMemCheckSourceRadialFlux(psPtr ptr);
     63
     64pmSourceEllipticalFlux *pmSourceEllipticalFluxAlloc();
     65bool psMemCheckSourceEllipticalFlux(psPtr ptr);
    4766
    4867// *** pmSourceRadialProfile describes the radial profile of a source in elliptical contours, and
    … …  
    5069pmSourceRadialProfile *pmSourceRadialProfileAlloc();
    5170bool psMemCheckSourceRadialProfile(psPtr ptr);
    52 
    53 // *** pmSourceRadialProfileFreeVectors frees the intermediate data values
    54 bool pmSourceRadialProfileFreeVectors(pmSourceRadialProfile *profile);
    5571
    5672// *** pmSourceExtendedPars describes the possible collection of extended flux measurements for a source
    … …  
    6581bool pmSourceRadialProfileSortPair(psVector *index, psVector *extra);
    6682
     83
     84
    6785/// @}
    6886# endif /* PM_SOURCE_H */

  • trunk/psModules/src/objects/pmSourceIO.c

    r27531 r27818  
    551551                }
    552552                if (!strcmp (exttype, "PS1_V1")) {
    553                     status &= pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe);
     553                    status &= pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
    554554                }
    555555                if (!strcmp (exttype, "PS1_V2")) {
    556                     status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe);
     556                    status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
    557557                }
    558558                if (!strcmp (exttype, "PS1_DV1")) {
    559                     status &= pmSourcesWrite_CMF_PS1_DV1_XSRC (file->fits, sources, xsrcname, recipe);
     559                    status &= pmSourcesWrite_CMF_PS1_DV1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
    560560                }
    561561            }
    … …  
    568568                }
    569569                if (!strcmp (exttype, "PS1_V1")) {
    570                     status &= pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname);
     570                    status &= pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, readout, sources, xfitname);
    571571                }
    572572                if (!strcmp (exttype, "PS1_V2")) {
    573                     status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname);
     573                    status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, readout, sources, xfitname);
    574574                }
    575575                if (!strcmp (exttype, "PS1_DV1")) {
    576                     status &= pmSourcesWrite_CMF_PS1_DV1_XFIT (file->fits, sources, xfitname);
     576                    status &= pmSourcesWrite_CMF_PS1_DV1_XFIT (file->fits, readout, sources, xfitname);
    577577                }
    578578            }

  • trunk/psModules/src/objects/pmSourceIO.h

    r27531 r27818  
    3636
    3737bool pmSourcesWrite_CMF_PS1_V1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
    38 bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe);
    39 bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, psArray *sources, char *extname);
     38bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
     39bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname);
    4040
    4141bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
    42 bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe);
    43 bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname);
     42bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
     43bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname);
    4444
    4545bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
    46 bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe);
    47 bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, psArray *sources, char *extname);
     46bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
     47bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname);
    4848
    4949bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config);

  • trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c

    r27531 r27818  
    364364
    365365// XXX this layout is still the same as PS1_DEV_1
    366 bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe)
     366bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
    367367{
    368368
    … …  
    536536
    537537// XXX this layout is still the same as PS1_DEV_1
    538 bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, psArray *sources, char *extname)
     538bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname)
    539539{
    540540

  • trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c

    r27177 r27818  
    339339
    340340// XXX this layout is still the same as PS1_DEV_1
    341 bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe)
     341bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
    342342{
    343343
    … …  
    516516
    517517// XXX this layout is still the same as PS1_DEV_1
    518 bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, psArray *sources, char *extname)
     518bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname)
    519519{
    520520

  • trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c

    r27177 r27818  
    4545// followed by a zero-size matrix, followed by the table data
    4646
    47 bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources,
    48                                 psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
     47bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
    4948{
    5049    PS_ASSERT_PTR_NON_NULL(fits, false);
    … …  
    5453    psArray *table;
    5554    psMetadata *row;
    56     int i;
    5755    psF32 *PAR, *dPAR;
    5856    psEllipseAxes axes;
    … …  
    8381        pmSource *source = (pmSource *) sources->data[0];
    8482        if (source->seq == -1) {
    85           // let's write these out in S/N order
    86           sources = psArraySort (sources, pmSourceSortBySN);
     83            // let's write these out in S/N order
     84            sources = psArraySort (sources, pmSourceSortBySN);
    8785        } else {
    88           sources = psArraySort (sources, pmSourceSortBySeq);
     86            sources = psArraySort (sources, pmSourceSortBySeq);
    8987        }
    9088    }
    … …  
    9492    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    9593    // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
    96     for (i = 0; i < sources->n; i++) {
     94    for (int i = 0; i < sources->n; i++) {
    9795        pmSource *source = (pmSource *) sources->data[i];
    9896
    … …  
    102100        // generated on Alloc, and would thus be wrong for read in sources.
    103101        if (source->seq == -1) {
    104           source->seq = i;
     102            source->seq = i;
    105103        }
    106104
    … …  
    114112            yPos = PAR[PM_PAR_YPOS];
    115113            if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
    116               xErr = dPAR[PM_PAR_XPOS];
    117               yErr = dPAR[PM_PAR_YPOS];
     114                xErr = dPAR[PM_PAR_XPOS];
     115                yErr = dPAR[PM_PAR_YPOS];
    118116            } else {
    119               // in linear-fit mode, there is no error on the centroid
    120               xErr = source->peak->dx;
    121               yErr = source->peak->dy;
     117                // in linear-fit mode, there is no error on the centroid
     118                xErr = source->peak->dx;
     119                yErr = source->peak->dy;
    122120            }
    123121            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
    … …  
    212210    }
    213211
     212    // XXX why do we make a copy here to be supplemented with the masks?  why not do this in the calling function?
    214213    psMetadata *header = psMetadataCopy(NULL, tableHeader);
    215214    pmSourceMasksHeader(header);
    … …  
    344343}
    345344
    346 // XXX this layout is still the same as PS1_DEV_1
    347 bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe)
     345bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
    348346{
    349347
    … …  
    354352    psF32 xPos, yPos;
    355353    psF32 xErr, yErr;
     354    int nRow = -1;
    356355
    357356    // create a header to hold the output data
    … …  
    361360    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname);
    362361
     362    pmChip *chip = readout->parent->parent;
     363    pmFPA  *fpa  = chip->parent;
     364
     365    // zero point corrections
     366    bool status1 = false;
     367    bool status2 = false;
     368    float magOffset = 0.0;
     369    float exptime   = psMetadataLookupF32(&status1, fpa->concepts, "FPA.EXPOSURE");
     370    float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
     371    if (!isfinite(zeropt)) {
     372        zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
     373    }
     374    if (status1 && status2 && (exptime > 0.0)) {
     375        magOffset = zeropt + 2.5*log10(exptime);
     376    }
     377
    363378    // let's write these out in S/N order
    364379    sources = psArraySort (sources, pmSourceSortBySN);
    … …  
    367382
    368383    // which extended source analyses should we perform?
    369     // bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     384    bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
     385    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    370386    // bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    371     // bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    372387    // bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
    373388
    374     psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    375     psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
    376     assert (radialBinsLower->n == radialBinsUpper->n);
     389    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
     390    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
     391    psAssert (radMin->n == radMax->n, "inconsistent annular bins");
     392
     393    // int nRadialBins = 0;
     394    // if (doAnnuli) {
     395    //  // get the max count of radial bins
     396    //  for (int i = 0; i < sources->n; i++) {
     397    //      pmSource *source = sources->data[i];
     398    //      if (!source->extpars) continue;
     399    //      if (!source->extpars->radProfile ) continue;
     400    //         if (!source->extpars->radProfile->binSB) continue;
     401    //      nRadialBins = PS_MAX(nRadialBins, source->extpars->radProfile->binSB->n);
     402    //  }
     403    // }
    377404
    378405    // we write out all sources, regardless of quality.  the source flags tell us the state
    … …  
    410437        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    411438
    412 # if (0)
     439        float AxialRatio = NAN;
     440        float AxialTheta = NAN;
     441        pmSourceExtendedPars *extpars = source->extpars;
     442        if (extpars) {
     443            AxialRatio = extpars->axes.minor / extpars->axes.major;
     444            AxialTheta = extpars->axes.theta;
     445        }
     446        psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO",       PS_DATA_F32, "Axial Ratio of radial profile",              AxialRatio);
     447        psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA",        PS_DATA_F32, "Angle of radial profile ellipse",                  AxialTheta);
     448
    413449        // Petrosian measurements
    414450        // XXX insert header data: petrosian ref radius, flux ratio
     451        // XXX check flags to see if Pet was measured
    415452        if (doPetrosian) {
    416             pmSourcePetrosianValues *petrosian = source->extpars->petrosian;
    417             if (petrosian) {
    418                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       petrosian->mag);
    419                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr);
    420                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          petrosian->rad);
    421                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    petrosian->radErr);
     453            pmSourceExtendedPars *extpars = source->extpars;
     454            if (extpars) {
     455                // XXX note that this mag is either calibrated or instrumental depending on existence of zero point
     456                float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point
     457                float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point
     458                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude", mag);
     459                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", magErr);
     460                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius);
     461                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr);
     462                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50",     PS_DATA_F32, "Petrosian R50 (pix)", extpars->petrosianR50);
     463                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)", extpars->petrosianR50Err);
     464                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90);
     465                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err);
    422466            } else {
    423467                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
    … …  
    425469                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          NAN);
    426470                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    NAN);
     471                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50",     PS_DATA_F32, "Petrosian R50 (pix)", NAN);
     472                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)",NAN);
     473                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", NAN);
     474                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN);
    427475            }
    428476        }
    429477
     478# if (0)
    430479        // Kron measurements
    431480        if (doKron) {
    … …  
    460509            }
    461510        }
    462 
    463         // Flux Annuli
     511# endif
     512
     513        // Flux Annuli (if we have extended source measurements, we have these.  only optionally save them)
    464514        if (doAnnuli) {
    465             pmSourceAnnuli *annuli = source->extpars->annuli;
    466             if (annuli) {
    467                 psVector *fluxVal = annuli->flux;
    468                 psVector *fluxErr = annuli->fluxErr;
    469                 psVector *fluxVar = annuli->fluxVar;
    470 
    471                 for (int j = 0; j < fluxVal->n; j++) {
    472                     char name[32];
    473                     sprintf (name, "FLUX_VAL_R_%02d", j);
    474                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]);
    475                     sprintf (name, "FLUX_ERR_R_%02d", j);
    476                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]);
    477                     sprintf (name, "FLUX_VAR_R_%02d", j);
    478                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]);
    479                 }
    480             } else {
    481                 for (int j = 0; j < radialBinsLower->n; j++) {
    482                     char name[32];
    483                     sprintf (name, "FLUX_VAL_R_%02d", j);
    484                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN);
    485                     sprintf (name, "FLUX_ERR_R_%02d", j);
    486                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN);
    487                     sprintf (name, "FLUX_VAR_R_%02d", j);
    488                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN);
    489                 }
    490             }
    491         }
    492 
    493 # endif
    494         psArrayAdd (table, 100, row);
    495         psFree (row);
    496     }
    497 
     515            psVector *radSB   = psVectorAlloc(radMin->n, PS_TYPE_F32);
     516            psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32);
     517            psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32);
     518            psVectorInit (radSB, NAN);
     519            psVectorInit (radFlux, NAN);
     520            psVectorInit (radFill, NAN);
     521            if (!source->extpars) goto empty_annuli;
     522            if (!source->extpars->radProfile) goto empty_annuli;
     523            if (!source->extpars->radProfile->binSB) goto empty_annuli;
     524            psAssert (source->extpars->radProfile->binSum, "programming error");
     525            psAssert (source->extpars->radProfile->binFill, "programming error");
     526            psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths");
     527            psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths");
     528            psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths");
     529
     530            // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
     531            for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) {
     532                radSB->data.F32[j]   = source->extpars->radProfile->binSB->data.F32[j];
     533                radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j];
     534                radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j];
     535            }
     536
     537        empty_annuli:
     538            psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB);
     539            psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux);
     540            psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill);
     541            psFree (radSB);
     542            psFree (radFlux);
     543            psFree (radFill);
     544        }
     545        if (nRow < 0) {
     546            nRow = row->list->n;
     547        } else {
     548            psAssert (nRow == row->list->n, "inconsistent row lengths");
     549        }
     550        psArrayAdd (table, 100, row);
     551        psFree (row);
     552    }
     553   
    498554    if (table->n == 0) {
    499         if (!psFitsWriteBlank (fits, outhead, extname)) {
    500             psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
    501             psFree(outhead);
    502             psFree(table);
    503             return false;
    504         }
    505         psFree (outhead);
    506         psFree (table);
    507         return true;
    508     }
    509 
     555        if (!psFitsWriteBlank (fits, outhead, extname)) {
     556            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
     557            psFree(outhead);
     558            psFree(table);
     559            return false;
     560        }
     561        psFree (outhead);
     562        psFree (table);
     563        return true;
     564    }
     565   
    510566    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
    511567    if (!psFitsWriteTable (fits, outhead, table, extname)) {
    512         psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
    513         psFree (outhead);
    514         psFree(table);
    515         return false;
     568        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
     569        psFree (outhead);
     570    psFree(table);
     571    return false;
    516572    }
    517573    psFree (outhead);
    518574    psFree (table);
    519 
     575   
    520576    return true;
    521577}
    522578
    523579// XXX this layout is still the same as PS1_DEV_1
    524 bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname)
     580bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname)
    525581{
    526582
Note: See TracChangeset for help on using the changeset viewer.