Index: trunk/psModules/src/objects/pmPSFtryModel.c =================================================================== --- trunk/psModules/src/objects/pmPSFtryModel.c (revision 27782) +++ trunk/psModules/src/objects/pmPSFtryModel.c (revision 27818) @@ -107,6 +107,6 @@ // as we loop over orders, we need to refer to the initial selection, but we modify the // option values to match the current guess: save the max values here: - int Nx = options->psfTrendNx; - int Ny = options->psfTrendNy; + int Nx = (options->psfTrendMode == PM_TREND_MAP) ? options->psfTrendNx : options->psfTrendNx + 1; + int Ny = (options->psfTrendMode == PM_TREND_MAP) ? options->psfTrendNy : options->psfTrendNy + 1; for (int i = orderMin; i <= orderMax; i++) { Index: trunk/psModules/src/objects/pmSourceExtendedPars.c =================================================================== --- trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 27782) +++ trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 27818) @@ -36,23 +36,67 @@ #include "pmSourceExtendedPars.h" -// *** pmSourceRadialProfile describes the radial profile of a source in elliptical contours, and -// intermediate data used to measure the profile +// pmSourceRadialFlux carries the raw radial flux information, including angular bins +static void pmSourceRadialFluxFree(pmSourceRadialFlux *flux) +{ + if (!flux) return; + psFree(flux->radii); + psFree(flux->fluxes); + psFree(flux->theta); + psFree(flux->isophotalRadii); +} + +pmSourceRadialFlux *pmSourceRadialFluxAlloc() +{ + pmSourceRadialFlux *flux = (pmSourceRadialFlux *)psAlloc(sizeof(pmSourceRadialFlux)); + psMemSetDeallocator(flux, (psFreeFunc) pmSourceRadialFluxFree); + + flux->radii = NULL; + flux->fluxes = NULL; + flux->theta = NULL; + flux->isophotalRadii = NULL; + + return flux; +} + +bool psMemCheckSourceRadialFlux(psPtr ptr) +{ + PS_ASSERT_PTR(ptr, false); + return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceRadialFluxFree); +} + +// pmSourceEllipticalFlux carries the elliptical renormalized radial flux info +static void pmSourceEllipticalFluxFree(pmSourceEllipticalFlux *flux) +{ + if (!flux) return; + psFree(flux->radiusElliptical); + psFree(flux->fluxElliptical); +} + +pmSourceEllipticalFlux *pmSourceEllipticalFluxAlloc() +{ + pmSourceEllipticalFlux *flux = (pmSourceEllipticalFlux *)psAlloc(sizeof(pmSourceEllipticalFlux)); + psMemSetDeallocator(flux, (psFreeFunc) pmSourceEllipticalFluxFree); + + flux->radiusElliptical = NULL; + flux->fluxElliptical = NULL; + + return flux; +} + +bool psMemCheckSourceEllipticalFlux(psPtr ptr) +{ + PS_ASSERT_PTR(ptr, false); + return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceEllipticalFluxFree); +} + +// pmSourceRadialProfile defines flux information in radial bins static void pmSourceRadialProfileFree(pmSourceRadialProfile *profile) { - if (!profile) { - return; - } - psFree(profile->radii); - psFree(profile->fluxes); - psFree(profile->theta); - psFree(profile->isophotalRadii); - - psFree(profile->radiusElliptical); - psFree(profile->fluxElliptical); - + if (!profile) return; psFree(profile->binSB); psFree(profile->binSBstdev); psFree(profile->binSBerror); - + psFree(profile->binSum); + psFree(profile->binFill); psFree(profile->radialBins); psFree(profile->area); @@ -63,20 +107,12 @@ pmSourceRadialProfile *profile = (pmSourceRadialProfile *)psAlloc(sizeof(pmSourceRadialProfile)); psMemSetDeallocator(profile, (psFreeFunc) pmSourceRadialProfileFree); - - profile->radii = NULL; - profile->fluxes = NULL; - profile->theta = NULL; - profile->isophotalRadii = NULL; - - profile->radiusElliptical = NULL; - profile->fluxElliptical = NULL; profile->binSB = NULL; profile->binSBstdev = NULL; profile->binSBerror = NULL; - + profile->binSum = NULL; + profile->binFill = NULL; profile->radialBins = NULL; profile->area = NULL; - return profile; } @@ -88,6 +124,6 @@ } - -// *** pmSourceRadialProfileFreeVectors frees the intermediate data values +# if (0) +// pmSourceRadialProfileFreeVectors frees the intermediate data values bool pmSourceRadialProfileFreeVectors(pmSourceRadialProfile *profile) { @@ -124,4 +160,5 @@ return true; } +# endif // *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors @@ -150,7 +187,8 @@ if (!pars) return; - psFree(pars->profile); - psFree(pars->petrosian_50); - psFree(pars->petrosian_80); + psFree(pars->radFlux); + psFree(pars->ellipticalFlux); + psFree(pars->radProfile); + psFree(pars->petProfile); return; } @@ -160,8 +198,17 @@ psMemSetDeallocator(pars, (psFreeFunc) pmSourceExtendedParsFree); - pars->profile = NULL; - pars->petrosian_50 = NULL; - pars->petrosian_80 = NULL; - + pars->radFlux = NULL; + pars->ellipticalFlux = NULL; + pars->radProfile = NULL; + pars->petProfile = NULL; + + pars->petrosianFlux = NAN; + pars->petrosianFluxErr = NAN; + pars->petrosianRadius = NAN; + pars->petrosianRadiusErr = NAN; + pars->petrosianR90 = NAN; + pars->petrosianR90Err = NAN; + pars->petrosianR50 = NAN; + pars->petrosianR50Err = NAN; return pars; } Index: trunk/psModules/src/objects/pmSourceExtendedPars.h =================================================================== --- trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 27782) +++ trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 27818) @@ -19,30 +19,49 @@ psVector *theta; // angles corresponding to above radial profiles psVector *isophotalRadii; // isophotal radius for the above angles +} pmSourceRadialFlux; +typedef struct { psVector *radiusElliptical; // normalized radial coordinates for all relevant pixels psVector *fluxElliptical; // flux for the above radial coordinates +} pmSourceEllipticalFlux; +typedef struct { psVector *binSB; // mean surface brightness within radial bins psVector *binSBstdev; // scatter of mean surface brightness within radial bins psVector *binSBerror; // formal error on mean surface brightness within radial bins - - psVector *radialBins; // radii corresponding to above binnedBlux + psVector *binSum; // sum of flux within radial bins + psVector *binFill; // fraction of area actually lit + psVector *radialBins; // radii corresponding to above binnedFlux psVector *area; // differential area of the non-overlapping radial bins - - psEllipseAxes axes; // shape of elliptical contour } pmSourceRadialProfile; typedef struct { - float flux; - float fluxErr; - float radius; - float radiusErr; + float flux; + float fluxErr; + float radius; + float radiusErr; } pmSourceExtendedFlux; typedef struct { - pmSourceRadialProfile *profile; - pmSourceExtendedFlux *petrosian_50; - pmSourceExtendedFlux *petrosian_80; + pmSourceRadialFlux *radFlux; // raw radial flux information + pmSourceEllipticalFlux *ellipticalFlux; // flux for elliptically-renormalized radii + pmSourceRadialProfile *radProfile; // surface brightness profile in specified fixed bins + pmSourceRadialProfile *petProfile; // surface brightness profile in petrosian bins + psEllipseAxes axes; // shape of elliptical contour + float petrosianFlux; + float petrosianFluxErr; + float petrosianRadius; + float petrosianRadiusErr; + float petrosianR90; + float petrosianR90Err; + float petrosianR50; + float petrosianR50Err; } pmSourceExtendedPars; + +pmSourceRadialFlux *pmSourceRadialFluxAlloc(); +bool psMemCheckSourceRadialFlux(psPtr ptr); + +pmSourceEllipticalFlux *pmSourceEllipticalFluxAlloc(); +bool psMemCheckSourceEllipticalFlux(psPtr ptr); // *** pmSourceRadialProfile describes the radial profile of a source in elliptical contours, and @@ -50,7 +69,4 @@ pmSourceRadialProfile *pmSourceRadialProfileAlloc(); bool psMemCheckSourceRadialProfile(psPtr ptr); - -// *** pmSourceRadialProfileFreeVectors frees the intermediate data values -bool pmSourceRadialProfileFreeVectors(pmSourceRadialProfile *profile); // *** pmSourceExtendedPars describes the possible collection of extended flux measurements for a source @@ -65,4 +81,6 @@ bool pmSourceRadialProfileSortPair(psVector *index, psVector *extra); + + /// @} # endif /* PM_SOURCE_H */ Index: trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO.c (revision 27782) +++ trunk/psModules/src/objects/pmSourceIO.c (revision 27818) @@ -551,11 +551,11 @@ } if (!strcmp (exttype, "PS1_V1")) { - status &= pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe); + status &= pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); } if (!strcmp (exttype, "PS1_V2")) { - status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe); + status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); } if (!strcmp (exttype, "PS1_DV1")) { - status &= pmSourcesWrite_CMF_PS1_DV1_XSRC (file->fits, sources, xsrcname, recipe); + status &= pmSourcesWrite_CMF_PS1_DV1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); } } @@ -568,11 +568,11 @@ } if (!strcmp (exttype, "PS1_V1")) { - status &= pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname); + status &= pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, readout, sources, xfitname); } if (!strcmp (exttype, "PS1_V2")) { - status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname); + status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, readout, sources, xfitname); } if (!strcmp (exttype, "PS1_DV1")) { - status &= pmSourcesWrite_CMF_PS1_DV1_XFIT (file->fits, sources, xfitname); + status &= pmSourcesWrite_CMF_PS1_DV1_XFIT (file->fits, readout, sources, xfitname); } } Index: trunk/psModules/src/objects/pmSourceIO.h =================================================================== --- trunk/psModules/src/objects/pmSourceIO.h (revision 27782) +++ trunk/psModules/src/objects/pmSourceIO.h (revision 27818) @@ -36,14 +36,14 @@ bool pmSourcesWrite_CMF_PS1_V1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); -bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe); -bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, psArray *sources, char *extname); +bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname); bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); -bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe); -bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname); +bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname); bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); -bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe); -bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, psArray *sources, char *extname); +bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname); bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config); Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 27782) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 27818) @@ -364,5 +364,5 @@ // XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) { @@ -536,5 +536,5 @@ // XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, psArray *sources, char *extname) +bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname) { Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 27782) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 27818) @@ -339,5 +339,5 @@ // XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) { @@ -516,5 +516,5 @@ // XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, psArray *sources, char *extname) +bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname) { Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 27782) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 27818) @@ -45,6 +45,5 @@ // followed by a zero-size matrix, followed by the table data -bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, - psMetadata *imageHeader, psMetadata *tableHeader, char *extname) +bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -54,5 +53,4 @@ psArray *table; psMetadata *row; - int i; psF32 *PAR, *dPAR; psEllipseAxes axes; @@ -83,8 +81,8 @@ pmSource *source = (pmSource *) sources->data[0]; if (source->seq == -1) { - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); } else { - sources = psArraySort (sources, pmSourceSortBySeq); + sources = psArraySort (sources, pmSourceSortBySeq); } } @@ -94,5 +92,5 @@ // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. - for (i = 0; i < sources->n; i++) { + for (int i = 0; i < sources->n; i++) { pmSource *source = (pmSource *) sources->data[i]; @@ -102,5 +100,5 @@ // generated on Alloc, and would thus be wrong for read in sources. if (source->seq == -1) { - source->seq = i; + source->seq = i; } @@ -114,10 +112,10 @@ yPos = PAR[PM_PAR_YPOS]; if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; } else { - // in linear-fit mode, there is no error on the centroid - xErr = source->peak->dx; - yErr = source->peak->dy; + // in linear-fit mode, there is no error on the centroid + xErr = source->peak->dx; + yErr = source->peak->dy; } if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { @@ -212,4 +210,5 @@ } + // XXX why do we make a copy here to be supplemented with the masks? why not do this in the calling function? psMetadata *header = psMetadataCopy(NULL, tableHeader); pmSourceMasksHeader(header); @@ -344,6 +343,5 @@ } -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) { @@ -354,4 +352,5 @@ psF32 xPos, yPos; psF32 xErr, yErr; + int nRow = -1; // create a header to hold the output data @@ -361,4 +360,20 @@ psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); + pmChip *chip = readout->parent->parent; + pmFPA *fpa = chip->parent; + + // zero point corrections + bool status1 = false; + bool status2 = false; + float magOffset = 0.0; + float exptime = psMetadataLookupF32(&status1, fpa->concepts, "FPA.EXPOSURE"); + float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); + if (!isfinite(zeropt)) { + zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); + } + if (status1 && status2 && (exptime > 0.0)) { + magOffset = zeropt + 2.5*log10(exptime); + } + // let's write these out in S/N order sources = psArraySort (sources, pmSourceSortBySN); @@ -367,12 +382,24 @@ // which extended source analyses should we perform? - // bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); - // bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); - psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); - psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); - assert (radialBinsLower->n == radialBinsUpper->n); + psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + psAssert (radMin->n == radMax->n, "inconsistent annular bins"); + + // int nRadialBins = 0; + // if (doAnnuli) { + // // get the max count of radial bins + // for (int i = 0; i < sources->n; i++) { + // pmSource *source = sources->data[i]; + // if (!source->extpars) continue; + // if (!source->extpars->radProfile ) continue; + // if (!source->extpars->radProfile->binSB) continue; + // nRadialBins = PS_MAX(nRadialBins, source->extpars->radProfile->binSB->n); + // } + // } // we write out all sources, regardless of quality. the source flags tell us the state @@ -410,14 +437,31 @@ psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); -# if (0) + float AxialRatio = NAN; + float AxialTheta = NAN; + pmSourceExtendedPars *extpars = source->extpars; + if (extpars) { + AxialRatio = extpars->axes.minor / extpars->axes.major; + AxialTheta = extpars->axes.theta; + } + psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO", PS_DATA_F32, "Axial Ratio of radial profile", AxialRatio); + psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA", PS_DATA_F32, "Angle of radial profile ellipse", AxialTheta); + // Petrosian measurements // XXX insert header data: petrosian ref radius, flux ratio + // XXX check flags to see if Pet was measured if (doPetrosian) { - pmSourcePetrosianValues *petrosian = source->extpars->petrosian; - if (petrosian) { - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", petrosian->mag); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", petrosian->rad); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", petrosian->radErr); + pmSourceExtendedPars *extpars = source->extpars; + if (extpars) { + // XXX note that this mag is either calibrated or instrumental depending on existence of zero point + float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point + float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", mag); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", magErr); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", extpars->petrosianR50); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)", extpars->petrosianR50Err); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err); } else { psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); @@ -425,7 +469,12 @@ psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)",NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN); } } +# if (0) // Kron measurements if (doKron) { @@ -460,67 +509,74 @@ } } - - // Flux Annuli +# endif + + // Flux Annuli (if we have extended source measurements, we have these. only optionally save them) if (doAnnuli) { - pmSourceAnnuli *annuli = source->extpars->annuli; - if (annuli) { - psVector *fluxVal = annuli->flux; - psVector *fluxErr = annuli->fluxErr; - psVector *fluxVar = annuli->fluxVar; - - for (int j = 0; j < fluxVal->n; j++) { - char name[32]; - sprintf (name, "FLUX_VAL_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]); - sprintf (name, "FLUX_ERR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]); - sprintf (name, "FLUX_VAR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]); - } - } else { - for (int j = 0; j < radialBinsLower->n; j++) { - char name[32]; - sprintf (name, "FLUX_VAL_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN); - sprintf (name, "FLUX_ERR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN); - sprintf (name, "FLUX_VAR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN); - } - } - } - -# endif - psArrayAdd (table, 100, row); - psFree (row); - } - + psVector *radSB = psVectorAlloc(radMin->n, PS_TYPE_F32); + psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32); + psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32); + psVectorInit (radSB, NAN); + psVectorInit (radFlux, NAN); + psVectorInit (radFill, NAN); + if (!source->extpars) goto empty_annuli; + if (!source->extpars->radProfile) goto empty_annuli; + if (!source->extpars->radProfile->binSB) goto empty_annuli; + psAssert (source->extpars->radProfile->binSum, "programming error"); + psAssert (source->extpars->radProfile->binFill, "programming error"); + psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths"); + psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths"); + psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths"); + + // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux + for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) { + radSB->data.F32[j] = source->extpars->radProfile->binSB->data.F32[j]; + radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j]; + radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j]; + } + + empty_annuli: + psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB); + psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); + psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); + psFree (radSB); + psFree (radFlux); + psFree (radFill); + } + if (nRow < 0) { + nRow = row->list->n; + } else { + psAssert (nRow == row->list->n, "inconsistent row lengths"); + } + psArrayAdd (table, 100, row); + psFree (row); + } + if (table->n == 0) { - if (!psFitsWriteBlank (fits, outhead, extname)) { - psError(psErrorCodeLast(), false, "Unable to write empty sources file."); - psFree(outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; - } - + if (!psFitsWriteBlank (fits, outhead, extname)) { + psError(psErrorCodeLast(), false, "Unable to write empty sources file."); + psFree(outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; + } + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); - psFree (outhead); - psFree(table); - return false; + psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; } psFree (outhead); psFree (table); - + return true; } // XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname) +bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname) {