Index: /branches/eam_branches/ipp-20101205/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 30093) +++ /branches/eam_branches/ipp-20101205/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 30094) @@ -20,4 +20,5 @@ #include "pmConfig.h" +#include "pmErrorCodes.h" #include "pmDetrendDB.h" @@ -71,5 +72,4 @@ psF32 errMag, chisq, apRadius; psS32 nPix, nDOF; - char keyword1[80], keyword2[80]; pmChip *chip = readout->parent->parent; @@ -103,4 +103,5 @@ table = psArrayAllocEmpty (sources->n); +# if (0) // we use this just to define the output vectors (which must be present for all objects) bool status = false; @@ -118,5 +119,5 @@ psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]); } - +# endif // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state @@ -258,4 +259,6 @@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); +# if (0) + // XXX if we have raw radial apertures, write them out here psVector *radFlux = psVectorAlloc(radMax->n, PS_TYPE_F32); psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32); @@ -284,4 +287,5 @@ psFree (radFluxErr); psFree (radFill); +# endif // XXX not sure how to get this : need to load Nimages with weight? @@ -781,2 +785,137 @@ return true; } + +// **** write out the radial flux values for the sources for a given matched-PSF image +// **** how do we distinguish the matched-PSF images from the non-matched version +bool pmSourcesWrite_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) +{ + + psArray *table; + psMetadata *row; + psF32 xPos, yPos; + char keyword1[80], keyword2[80]; + + // create a header to hold the output data + psMetadata *outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "radial flux table extension", extname); + + // we use this just to define the output vectors (which must be present for all objects) + bool status = false; + psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + psAssert (radMax, "this must have been defined and tested earlier!"); + psAssert (radMax->n, "this must have been defined and tested earlier!"); + psAssert (radMin->n == radMax->n, "inconsistent annular bins"); + + // write the radial profile apertures to header + for (int i = 0; i < radMax->n; i++) { + sprintf (keyword1, "RMIN_%02d", i); + sprintf (keyword2, "RMAX_%02d", i); + psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]); + psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]); + } + + psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES"); + if (!fwhmValues) { + psError (PM_ERR_CONFIG, true, "convolved or measured FWHM is not defined for this readout"); + return false; + } + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + + table = psArrayAllocEmpty (sources->n); + + // we write out all sources, regardless of quality. the source flags tell us the state + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // skip sources without radial aper measurements (or insufficient) + if (source->radialAper == NULL) continue; + psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set"); + + for (int entry = 0; entry < fwhmValues->n; entry++) { + + // choose the convolved EXT model, if available, otherwise the simple one + pmSourceRadialApertures *radialAper = source->radialAper->data[entry]; + assert (radialAper); + + bool useMoments = true; + useMoments = (useMoments && source->moments); // can't if there are no moments + useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured + useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed... + + if (useMoments) { + xPos = source->moments->Mx; + yPos = source->moments->My; + } else { + xPos = source->peak->xf; + yPos = source->peak->yf; + } + + 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_APER", 0, "Center of aperture measurements", xPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER", 0, "Center of aperture measurements", yPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM", 0, "FWHM of matched PSF", fwhmValues->data.F32[entry]); + + // XXX if we have raw radial apertures, write them out here + psVector *radFlux = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVector *radFill = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVectorInit (radFlux, NAN); + psVectorInit (radFluxErr, NAN); + psVectorInit (radFill, NAN); + if (!radialAper->flux) goto write_annuli; + if (!radialAper->fill) goto write_annuli; + psAssert (radialAper->flux->n <= radFlux->n, "inconsistent vector lengths"); + psAssert (radialAper->fill->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 < radialAper->flux->n; j++) { + radFlux->data.F32[j] = radialAper->flux->data.F32[j]; + radFluxErr->data.F32[j] = radialAper->fluxErr->data.F32[j]; + radFill->data.F32[j] = radialAper->fill->data.F32[j]; + } + + write_annuli: + psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); + psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli", radFluxErr); + psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); + psFree (radFlux); + psFree (radFluxErr); + psFree (radFill); + + 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; + } + + 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; + } + psFree (outhead); + psFree (table); + return true; +}