Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 17395) +++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 17396) @@ -237,6 +237,7 @@ if (!isfinite(shape.sxy)) return false; - PAR[PM_PAR_SKY] = moments->Sky; - PAR[PM_PAR_I0] = moments->Peak - moments->Sky; + // XXX turn this off here for now PAR[PM_PAR_SKY] = moments->Sky; + PAR[PM_PAR_SKY] = 0.0; + PAR[PM_PAR_I0] = moments->Peak; PAR[PM_PAR_XPOS] = peak->x; PAR[PM_PAR_YPOS] = peak->y; @@ -457,7 +458,10 @@ status = true; - status &= (dP < 0.5); +// status &= (dP < 0.5); status &= (PAR[PM_PAR_I0] > 0); status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + fprintf (stderr, "QGAUSS status pars: dP: %f, I0: %f, S/N: %f\n", + dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0])); return status; Index: trunk/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 17395) +++ trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 17396) @@ -19,5 +19,5 @@ note that a standard sersic model uses exp(-K*(z^(1/n) - 1). the additional elements (K, - the -1 offset) are absorbed in this model by the normalization, the exponenet, and the + the -1 offset) are absorbed in this model by the normalization, the exponent, and the radial scale. We fit the elements in this form, then re-normalize them on output. *****************************************************************************/ @@ -157,8 +157,8 @@ break; case PM_PAR_SXX: - params_min = 0.25; + params_min = 0.05; break; case PM_PAR_SYY: - params_min = 0.25; + params_min = 0.05; break; case PM_PAR_SXY: @@ -166,5 +166,5 @@ break; case PM_PAR_7: - params_min = 0.1; + params_min = 0.05; break; default: @@ -247,6 +247,7 @@ if (!isfinite(shape.sxy)) return false; - PAR[PM_PAR_SKY] = moments->Sky; - PAR[PM_PAR_I0] = moments->Peak - moments->Sky; + // XXX PAR[PM_PAR_SKY] = moments->Sky; + PAR[PM_PAR_SKY] = 0.0; + PAR[PM_PAR_I0] = moments->Peak; PAR[PM_PAR_XPOS] = peak->x; PAR[PM_PAR_YPOS] = peak->y; @@ -442,7 +443,10 @@ status = true; - status &= (dP < 0.5); +// status &= (dP < 0.5); status &= (PAR[PM_PAR_I0] > 0); status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + fprintf (stderr, "SERSIC status pars: dP: %f, I0: %f, S/N: %f\n", + dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0])); return status; Index: trunk/psModules/src/objects/pmPeaks.c =================================================================== --- trunk/psModules/src/objects/pmPeaks.c (revision 17395) +++ trunk/psModules/src/objects/pmPeaks.c (revision 17396) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-01-02 20:38:28 $ + * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -141,5 +141,5 @@ tmp->xf = x; tmp->yf = y; - + tmp->assigned = false; tmp->type = type; Index: trunk/psModules/src/objects/pmPeaks.h =================================================================== --- trunk/psModules/src/objects/pmPeaks.h (revision 17395) +++ trunk/psModules/src/objects/pmPeaks.h (revision 17396) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-01-02 20:42:46 $ + * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -58,5 +58,6 @@ float flux; ///< level in unsmoothed sci image float SN; ///< S/N implied by detection level - pmPeakType type; ///< Description of peak. + bool assigned; ///< is peak assigned to a source? + pmPeakType type; ///< Description of peak. } pmPeak; Index: trunk/psModules/src/objects/pmSource.c =================================================================== --- trunk/psModules/src/objects/pmSource.c (revision 17395) +++ trunk/psModules/src/objects/pmSource.c (revision 17396) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.51 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-03-19 00:51:09 $ + * @version $Revision: 1.52 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -50,5 +50,6 @@ psFree(tmp->modelPSF); psFree(tmp->modelEXT); - psFree(tmp->modelConv); + psFree(tmp->modelFits); + psFree(tmp->extpars); psFree(tmp->blends); psTrace("psModules.objects", 5, "---- end ----\n"); @@ -107,5 +108,5 @@ source->modelPSF = NULL; source->modelEXT = NULL; - source->modelConv = NULL; + source->modelFits = NULL; source->type = PM_SOURCE_TYPE_UNKNOWN; source->mode = PM_SOURCE_MODE_DEFAULT; @@ -947,11 +948,10 @@ return model; -// XXX when should I return the modelConv ?? + // the 'best' extended model is saved in source->modelEXT (may be a pointer to one of + // the elements of source->modelFits) case PM_SOURCE_TYPE_EXTENDED: - model = source->modelConv; - if (!model) { - model = source->modelEXT; - } + model = source->modelEXT; if (!model && source->modelPSF) { + // XXX raise an error or warning here? if (isPSF) { *isPSF = true; @@ -1015,5 +1015,6 @@ if (!strcasecmp (name, "DEFECT" )) return PM_SOURCE_MODE_DEFECT; if (!strcasecmp (name, "SATURATED" )) return PM_SOURCE_MODE_SATURATED; - if (!strcasecmp (name, "CRLIMIT" )) return PM_SOURCE_MODE_CRLIMIT; + if (!strcasecmp (name, "CRLIMIT" )) return PM_SOURCE_MODE_CR_LIMIT; + if (!strcasecmp (name, "EXTLIMIT" )) return PM_SOURCE_MODE_EXT_LIMIT; if (!strcasecmp (name, "SUBTRACTED")) return PM_SOURCE_MODE_SUBTRACTED; return PM_SOURCE_MODE_DEFAULT; @@ -1036,5 +1037,6 @@ case PM_SOURCE_MODE_DEFECT : return psStringCopy ("DEFECT" ); case PM_SOURCE_MODE_SATURATED : return psStringCopy ("SATURATED" ); - case PM_SOURCE_MODE_CRLIMIT : return psStringCopy ("CRLIMIT"); + case PM_SOURCE_MODE_CR_LIMIT : return psStringCopy ("CRLIMIT" ); + case PM_SOURCE_MODE_EXT_LIMIT : return psStringCopy ("EXTLIMIT" ); case PM_SOURCE_MODE_SUBTRACTED : return psStringCopy ("SUBTRACTED"); default: Index: trunk/psModules/src/objects/pmSource.h =================================================================== --- trunk/psModules/src/objects/pmSource.h (revision 17395) +++ trunk/psModules/src/objects/pmSource.h (revision 17396) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-03-05 01:08:08 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -22,14 +22,11 @@ * source. * - * XXX: The values given below are currently illustrative and will require - * some modification as the source classification code is developed. (TBD) - * */ typedef enum { - PM_SOURCE_TYPE_UNKNOWN, ///< a cosmic-ray - PM_SOURCE_TYPE_DEFECT, ///< a cosmic-ray - PM_SOURCE_TYPE_SATURATED, ///< random saturated pixels - PM_SOURCE_TYPE_STAR, ///< a good-quality star - PM_SOURCE_TYPE_EXTENDED, ///< an extended object (eg, galaxy) + PM_SOURCE_TYPE_UNKNOWN, ///< not yet classified + PM_SOURCE_TYPE_DEFECT, ///< a cosmic-ray + PM_SOURCE_TYPE_SATURATED, ///< random saturated pixels (eg, bleed trails) + PM_SOURCE_TYPE_STAR, ///< a good-quality star (subtracted model is PSF) + PM_SOURCE_TYPE_EXTENDED, ///< an extended object (eg, galaxy) (subtracted model is EXT) } pmSourceType; @@ -48,7 +45,8 @@ PM_SOURCE_MODE_BADPSF = 0x0400, ///< Failed to get good estimate of object's PSF PM_SOURCE_MODE_DEFECT = 0x0800, ///< Source is thought to be a defect - PM_SOURCE_MODE_SATURATED = 0x1000, ///< Source is thought to be saturation - PM_SOURCE_MODE_CRLIMIT = 0x2000, ///< Source has crNsigma above limit - PM_SOURCE_MODE_SUBTRACTED = 0x4000, ///< XXX this flag is actually only used internally (move) + PM_SOURCE_MODE_SATURATED = 0x1000, ///< Source is thought to be saturated pixels (bleed trail) + PM_SOURCE_MODE_CR_LIMIT = 0x2000, ///< Source has crNsigma above limit + PM_SOURCE_MODE_EXT_LIMIT = 0x4000, ///< Source has extNsigma above limit + PM_SOURCE_MODE_SUBTRACTED = 0x8000, ///< XXX this flag is actually only used internally (move) } pmSourceMode; @@ -60,4 +58,7 @@ * * XXX do I have to re-organize this (again!) to allow an arbitrary set of extended model fits?? + * XXX put the Mag and Err inside the pmModel? + * XXX keep the modelEXT or add to the psArray + * * */ @@ -65,5 +66,5 @@ const int id; ///< Unique ID for object int seq; ///< ID for output (generated on write) - pmPeak *peak; ///< Description of peak pixel. + pmPeak *peak; ///< Description of peak pixel. psImage *pixels; ///< Rectangular region including object pixels. psImage *weight; ///< Image variance. @@ -72,21 +73,21 @@ psImage *modelFlux; ///< cached copy of the best model for this source psImage *psfFlux; ///< cached copy of the psf model for this source - pmMoments *moments; ///< Basic moments measure for the object. + pmMoments *moments; ///< Basic moments measured for the object. pmModel *modelPSF; ///< PSF Model fit (parameters and type) - pmModel *modelEXT; ///< EXT (floating) Model fit (parameters and type). - pmModel *modelConv; ///< PSF-Convolved Model fit (parameters and type). + pmModel *modelEXT; ///< EXT Model fit used for subtraction (parameters and type) + psArray *modelFits; ///< collection of extended source models (best == modelEXT) pmSourceType type; ///< Best identification of object. - pmSourceMode mode; ///< Best identification of object. - psArray *blends; - float psfMag; ///< calculated from flux in modelPsf + pmSourceMode mode; ///< analysis flags set for object. + psArray *blends; ///< collection of sources thought to be confused with object + float psfMag; ///< calculated from flux in modelPSF float extMag; ///< calculated from flux in modelEXT float errMag; ///< error in psfMag OR extMag (depending on type) float apMag; ///< apMag corresponding to psfMag or extMag (depending on type) float pixWeight; ///< model-weighted coverage of valid pixels - float psfChisq; ///< probability of PSF + float psfChisq; ///< probability of PSF float crNsigma; ///< Nsigma deviation from PSF to CR float extNsigma; ///< Nsigma deviation from PSF to EXT + float sky, skyErr; ///< The sky and its error at the center of the object psRegion region; ///< area on image covered by selected pixels - float sky, skyErr; ///< The sky and its error at the center of the object pmSourceExtendedPars *extpars; ///< extended source parameters }; Index: trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO.c (revision 17395) +++ trunk/psModules/src/objects/pmSourceIO.c (revision 17396) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.55 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-03-17 22:04:27 $ + * @version $Revision: 1.56 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -327,6 +327,6 @@ // if this is not TRUE, the output files only contain the psf measurements. - bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "SAVE.XSRC"); - bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "SAVE.XFIT"); + bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_ANALYSIS"); + bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_FITS"); // define the EXTNAME values for the different data segments: @@ -445,5 +445,5 @@ if (xsrcname) { if (!strcmp (exttype, "PS1_DEV_1")) { - status = pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname); + status = pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe); } } @@ -457,4 +457,6 @@ psFree (headname); psFree (dataname); + psFree (xsrcname); + psFree (xfitname); psFree (outhead); return false; @@ -466,4 +468,6 @@ psFree (headname); psFree (dataname); + psFree (xsrcname); + psFree (xfitname); psFree (outhead); break; Index: trunk/psModules/src/objects/pmSourceIO.h =================================================================== --- trunk/psModules/src/objects/pmSourceIO.h (revision 17395) +++ trunk/psModules/src/objects/pmSourceIO.h (revision 17396) @@ -4,6 +4,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-03-05 01:08:08 $ + * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii * @@ -27,5 +27,5 @@ bool pmSourcesWrite_PS1_DEV_0 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, char *xsrcname); -bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname); +bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe); bool pmSourcesWrite_PS1_DEV_1_XFIT (psFits *fits, psArray *sources, char *extname); Index: trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 17395) +++ trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 17396) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-03-05 01:08:08 $ + * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -251,12 +251,11 @@ } -bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname) +bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) { + bool status; psArray *table; psMetadata *row; - int i; psF32 *PAR, *dPAR; - psEllipseAxes axes; psF32 xPos, yPos; psF32 xErr, yErr; @@ -273,6 +272,16 @@ table = psArrayAllocEmpty (sources->n); + // which extended source analyses should we perform? + 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); + // we write out all sources, regardless of quality. the source flags tell us the state - for (i = 0; i < sources->n; i++) { + for (int i = 0; i < sources->n; i++) { // skip source if it is not a ext sourc // XXX we have two places that extended source parameters are measured: @@ -283,13 +292,10 @@ pmSource *source = sources->data[i]; - // choose the convolved EXT model, if available, otherwise the simple one - pmModel *model = source->modelConv; - if (model == NULL) { - model = source->modelEXT; - } - if (model == NULL) continue; - - // XXX do we need to do this? + // skip sources without measurements if (source->extpars == NULL) continue; + + // we require a PSF model fit (ignore the real crud) + pmModel *model = source->modelPSF; + if (model == NULL) continue; // XXX I need to split the extended models from the extended aperture measurements @@ -301,8 +307,6 @@ yErr = dPAR[PM_PAR_YPOS]; - // XXX for the aperture values, I probably can / should just refer to the psf position and not write any of this - axes = pmPSF_ModelToAxes (PAR, 20.0); - row = psMetadataAlloc (); + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); @@ -311,19 +315,8 @@ psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG", PS_DATA_F32, "EXT fit instrumental magnitude", source->extMag); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); - - // XXX these should be major and minor, not 'x' and 'y' - psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", PS_DATA_F32, "EXT width in x coordinate", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", PS_DATA_F32, "EXT width in y coordinate", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA", PS_DATA_F32, "EXT orientation angle", axes.theta); - - // XXX at the moment, this will be a fixed sized table, but perhaps have multiple output formats depending on what is measured? - - // other values that I need to report: - // R, Mag Petrosian, ..... - if (source->extpars) { - - // Petrosian measurements + + // Petrosian measurements + // XXX insert header data: petrosian ref radius, flux ratio + if (doPetrosian) { pmSourcePetrosianValues *petrosian = source->extpars->petrosian; if (petrosian) { @@ -332,7 +325,14 @@ 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); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); + 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); } - - // Kron measurements + } + + // Kron measurements + if (doKron) { pmSourceKronValues *kron = source->extpars->kron; if (kron) { @@ -341,8 +341,15 @@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", kron->rad); psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", kron->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", NAN); } - - // Isophot measurements - // XXX insert header data: isophotal level + } + + // Isophot measurements + // XXX insert header data: isophotal level + if (doIsophotal) { pmSourceIsophotalValues *isophot = source->extpars->isophot; if (isophot) { @@ -351,6 +358,14 @@ psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", isophot->rad); psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", isophot->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", NAN); } - + } + + // Flux Annuli + if (doAnnuli) { pmSourceAnnuli *annuli = source->extpars->annuli; if (annuli) { @@ -367,24 +382,37 @@ 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); + } } } - psArrayAdd (table, 100, row); - psFree (row); + psArrayAdd (table, 100, row); + psFree (row); } if (table->n == 0) { - psFitsWriteBlank (fits, outhead, extname); - psFree (table); - return true; + psFitsWriteBlank (fits, outhead, extname); + psFree (outhead); + psFree (table); + return true; } psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(PS_ERR_IO, false, "writing ext data %s\n", extname); - psFree(table); - return false; - } + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); psFree (table); @@ -397,9 +425,9 @@ psArray *table; psMetadata *row; - int i; psF32 *PAR, *dPAR; psEllipseAxes axes; psF32 xPos, yPos; psF32 xErr, yErr; + char name[64]; // create a header to hold the output data @@ -412,72 +440,105 @@ sources = psArraySort (sources, pmSourceSortBySN); + // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams + int nParamMax = 0; + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (source->modelFits == NULL) continue; + for (int j = 0; j < source->modelFits->n; j++) { + pmModel *model = source->modelFits->data[j]; + assert (model); + nParamMax = PS_MAX (nParamMax, model->params->n); + } + } + table = psArrayAllocEmpty (sources->n); // we write out all sources, regardless of quality. the source flags tell us the state - for (i = 0; i < sources->n; i++) { - // skip source if it is not a ext sourc - // XXX we have two places that extended source parameters are measured: - // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models, - // psphotFitEXT, which does the simple extended source model fit (not psf-convolved) - // should we require both? + for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; - // XXX need to have an array of model fits; loop over the array and write out all - // which we calculated - - // choose the convolved EXT model, if available, otherwise the simple one - // XXX should not need to choose: write both out - pmModel *model = source->modelConv; - if (model == NULL) { - model = source->modelEXT; + // XXX if no model fits are saved, write out modelEXT? + if (source->modelFits == NULL) continue; + + // We have multiple sources : need to flag the one used to subtract the light (the 'best' model) + for (int j = 0; j < source->modelFits->n; j++) { + + // choose the convolved EXT model, if available, otherwise the simple one + pmModel *model = source->modelFits->data[j]; + assert (model); + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + + axes = pmPSF_ModelToAxes (PAR, 20.0); + + row = psMetadataAlloc (); + + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) + psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT", 0, "EXT model y coordinate", yPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); + + psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); + psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "name of model", pmModelClassGetName (model->type)); + + // XXX these should be major and minor, not 'x' and 'y' + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); + + // write out the other generic parameters + for (int k = 0; k < nParamMax; k++) { + if (k == PM_PAR_I0) continue; + if (k == PM_PAR_SKY) continue; + if (k == PM_PAR_XPOS) continue; + if (k == PM_PAR_YPOS) continue; + if (k == PM_PAR_SXX) continue; + if (k == PM_PAR_SXY) continue; + if (k == PM_PAR_SYY) continue; + + snprintf (name, 64, "EXT_PAR_%02d", k); + + if (k < model->params->n) { + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); + } else { + psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); + } + } + + // XXX other parameters which may be set. + // XXX flag / value to define the model + // XXX write out the model type, fit status flags + + psArrayAdd (table, 100, row); + psFree (row); } - if (model == NULL) continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - - axes = pmPSF_ModelToAxes (PAR, 20.0); - - row = psMetadataAlloc (); - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG", PS_DATA_F32, "EXT fit instrumental magnitude", source->extMag); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); - - // XXX these should be major and minor, not 'x' and 'y' - psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", PS_DATA_F32, "EXT width in x coordinate", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", PS_DATA_F32, "EXT width in y coordinate", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA", PS_DATA_F32, "EXT orientation angle", axes.theta); - - // XXX other parameters which may be set. - // XXX flag / value to define the model - - psArrayAdd (table, 100, row); - psFree (row); } if (table->n == 0) { - psFitsWriteBlank (fits, outhead, extname); - psFree (table); - return true; + psFitsWriteBlank (fits, outhead, extname); + psFree (outhead); + psFree (table); + return true; } psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(PS_ERR_IO, false, "writing ext data %s\n", extname); - psFree(table); - return false; - } + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); psFree (table); - return true; } Index: trunk/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- trunk/psModules/src/objects/pmSourcePhotometry.c (revision 17395) +++ trunk/psModules/src/objects/pmSourcePhotometry.c (revision 17396) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.40 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-04-02 22:40:36 $ + * @version $Revision: 1.41 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -119,6 +119,19 @@ } - // measure EXT model photometry (do both modelEXT and modelConv or just the one?) - status = pmSourcePhotometryModel (&source->extMag, source->modelEXT); + // if we have a collection of model fits, one of them is a pointer to modelEXT? + // XXX not guaranteed + if (source->modelFits) { + for (int i = 0; i < source->modelFits->n; i++) { + pmModel *model = source->modelFits->data[i]; + status = pmSourcePhotometryModel (&model->mag, model); + } + if (source->modelEXT) { + source->extMag = source->modelEXT->mag; + } + } else { + if (source->modelEXT) { + status = pmSourcePhotometryModel (&source->extMag, source->modelEXT); + } + } // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS