Index: trunk/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 33693) +++ trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 34259) @@ -89,5 +89,9 @@ // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. for (int i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; + // this is the source associated with this image + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; // If source->seq is -1, source was generated in this analysis. If source->seq is @@ -106,78 +110,107 @@ pmSourceOutputsSetMoments (&moments, source); + @PS1_DV?@ pmSourceDiffStats diffStats; + @PS1_DV?@ pmSourceDiffStatsInit(&diffStats); + @PS1_DV?@ if (source->diffStats) { + @PS1_DV?@ diffStats = *source->diffStats; + @PS1_DV?@ } + row = psMetadataAlloc (); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision - @=PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F32, "PSF RA coordinate (degrees)", outputs.ra); - @=PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F32, "PSF DEC coordinate (degrees)", outputs.dec); - - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); - - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux", source->psfFluxErr); - - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in reported aperture", source->apMagRaw); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); - @psfMag); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); + + @ALL,!PS1_V1,!PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); + @ALL,!PS1_V1,!PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux", source->psfFluxErr); + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); + @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in reported aperture", source->apMagRaw); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); + @PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", source->apFlux); + @PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", source->apFluxErr); + + @PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); - @>PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); - - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); - - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", source->pixWeightNotPoor); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); - - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); - - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); - - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kinner); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kouter); - - @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_RAD", PS_DATA_F32, "Radius where object hits sky", source->skyRadius); - @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_FLUX", PS_DATA_F32, "Flux / pix where object hits sky", source->skyFlux); - @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_SLOPE", PS_DATA_F32, "d(Flux/pix)/dRadius where object hits sky", source->skySlope); - - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2", PS_DATA_S32, "more padding", 0); - - // XXX not sure how to get this : need to load Nimages with weight? - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", source->nFrames); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); + @ALL,!PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); + @ALL,!PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); + + @>=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); + @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", source->pixWeightNotPoor); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); + + @>PS1_V2,PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); + @>PS1_V2,PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); + @>PS1_V2,PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); + @>PS1_V2,PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); + + @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); + @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); + @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); + @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); + @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 1.0 R1)", moments.Kinner); + @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 4.0 R1)", moments.Kouter); + + @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_RAD", PS_DATA_F32, "Radius where object hits sky", source->skyRadius); + @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_FLUX", PS_DATA_F32, "Flux / pix where object hits sky", source->skyFlux); + @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_SLOPE", PS_DATA_F32, "d(Flux/pix)/dRadius where object hits sky", source->skySlope); + + @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); + @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); + @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); + @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); + @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); + + @PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P", PS_DATA_F32, "distance to positive match source", diffStats.Rp); + @PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P", PS_DATA_F32, "signal-to-noise of pos match src", diffStats.SNp); + @PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M", PS_DATA_F32, "distance to negative match source", diffStats.Rm); + @PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M", PS_DATA_F32, "signal-to-noise of neg match src", diffStats.SNm); + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); + @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); + @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2", PS_DATA_S32, "more padding", 0); + @PS1_SV1@ + + // note that this definition is inconsistent with the definition in + // Ohana/src/libautocode/def. This version creates a table with data not + // properly aligned with the 8-byte boundaries. The structure defined by + // libautocode does this, but has a different order of elements (and adds + // padding2 to fix things). We have generated may files with PS1_SV1 as is, so + // I'll leave it. But in future a PS1_SV2 should be forced to match + // libautocode. Note that addstar knows to detect the alternate version of + // PS1_SV1 and correctly interpret its fields. + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", source->nFrames); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); psArrayAdd (table, 100, row); @@ -282,9 +315,10 @@ @ALL@ source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); @ALL@ source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); - @>PS1_V2@ source->apMagRaw = psMetadataLookupF32 (&status, row, "AP_MAG_RAW"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->apMagRaw = psMetadataLookupF32 (&status, row, "AP_MAG_RAW"); // XXX use these to determine PAR[PM_PAR_I0] if they exist? - @>PS1_V2@ source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); - @>PS1_V2@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); + // XXX add these to PS1_SV1? + @>PS1_V2,PS1_SV1,PS1_DV?@ source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); + @>PS1_V2,PS1_SV1,PS1_DV?@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); // XXX this scaling is incorrect: does not include the 2 \pi AREA factor @@ -307,5 +341,5 @@ @ALL@ source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); - @>PS1_V2@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT"); @ALL@ source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); @ALL@ source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); @@ -325,11 +359,23 @@ @ALL@ source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); - @>PS1_V2@ source->moments->Mrf = psMetadataLookupF32 (&status, row, "MOMENTS_R1"); - @>PS1_V2@ source->moments->Mrh = psMetadataLookupF32 (&status, row, "MOMENTS_RH"); - @>PS1_V2@ source->moments->KronFlux = psMetadataLookupF32 (&status, row, "KRON_FLUX"); - @>PS1_V2@ source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR"); - - @>PS1_V2@ source->moments->KronFinner = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER"); - @>PS1_V2@ source->moments->KronFouter = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER"); + // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data, + // we are storing enough information so the output will be consistent with the input + @>PS1_V2,PS1_SV1@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C"); + @>PS1_V2,PS1_SV1@ source->moments->Mxxy = 0.0; + @>PS1_V2,PS1_SV1@ source->moments->Mxyy = 0.0; + @>PS1_V2,PS1_SV1@ source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S"); + + @>PS1_V2,PS1_SV1@ source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C"); + @>PS1_V2,PS1_SV1@ source->moments->Mxxxy = 0.0; + @>PS1_V2,PS1_SV1@ source->moments->Mxxyy = 0.0; + @>PS1_V2,PS1_SV1@ source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S"); + @>PS1_V2,PS1_SV1@ source->moments->Myyyy = 0.0; + + @>PS1_V2,PS1_SV1,PS1_DV2@ source->moments->Mrf = psMetadataLookupF32 (&status, row, "MOMENTS_R1"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->moments->Mrh = psMetadataLookupF32 (&status, row, "MOMENTS_RH"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->moments->KronFlux = psMetadataLookupF32 (&status, row, "KRON_FLUX"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->moments->KronFinner = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->moments->KronFouter = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER"); @>PS1_V3@ source->skyRadius = psMetadataLookupF32 (&status, row, "SKY_LIMIT_RAD"); @@ -337,19 +383,21 @@ @>PS1_V3@ source->skySlope = psMetadataLookupF32 (&status, row, "SKY_LIMIT_SLOPE"); - // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data, - // we are storing enough information so the output will be consistent with the input - @>PS1_V2@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C"); - @>PS1_V2@ source->moments->Mxxy = 0.0; - @>PS1_V2@ source->moments->Mxyy = 0.0; - @>PS1_V2@ source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S"); - - @>PS1_V2@ source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C"); - @>PS1_V2@ source->moments->Mxxxy = 0.0; - @>PS1_V2@ source->moments->Mxxyy = 0.0; - @>PS1_V2@ source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S"); - @>PS1_V2@ source->moments->Myyyy = 0.0; - - @ALL@ source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); - @>PS1_V2@ source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2"); + @PS1_DV?@ int nPos = psMetadataLookupS32 (&status, row, "DIFF_NPOS"); + @PS1_DV?@ if (nPos) { + @PS1_DV?@ source->diffStats = pmSourceDiffStatsAlloc(); + @PS1_DV?@ source->diffStats->nGood = nPos; + @PS1_DV?@ source->diffStats->fRatio = psMetadataLookupF32 (&status, row, "DIFF_FRATIO"); + @PS1_DV?@ source->diffStats->nRatioBad = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_BAD"); + @PS1_DV?@ source->diffStats->nRatioMask = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_MASK"); + @PS1_DV?@ source->diffStats->nRatioAll = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_ALL"); + + @PS1_DV2@ source->diffStats->Rp = psMetadataLookupF32 (&status, row, "DIFF_R_P"); + @PS1_DV2@ source->diffStats->SNp = psMetadataLookupF32 (&status, row, "DIFF_SN_P"); + @PS1_DV2@ source->diffStats->Rm = psMetadataLookupF32 (&status, row, "DIFF_R_M"); + @PS1_DV2@ source->diffStats->SNm = psMetadataLookupF32 (&status, row, "DIFF_SN_M"); + @PS1_DV?@ } + + @ALL@ source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2"); assert (status); @@ -409,13 +457,17 @@ // 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 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]); - psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]); + sprintf (keyword1, "RMIN_%02d", i); + sprintf (keyword2, "RMAX_%02d", i); + psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]); + psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]); } // 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]; + // this is the source associated with this image + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; // skip sources without measurements @@ -470,4 +522,5 @@ 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); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL", PS_DATA_F32, "Petrosian Fill Factor", extpars->petrosianFill); } else { psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); @@ -479,4 +532,5 @@ 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); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL", PS_DATA_F32, "Petrosian Fill Factor", NAN); } } @@ -539,10 +593,104 @@ psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); psFree (outhead); - psFree(table); - return false; + psFree(table); + return false; } psFree (outhead); psFree (table); + return true; +} + +bool pmSourcesRead_CMF_@CMFMODE@_XSRC(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + + bool status; + long numSources = psFitsTableSize(fits); // Number of sources in table + if (numSources == 0) { + psError(psErrorCodeLast(), false, "XSRC Table contains no entries\n"); + return false; + } + + // petrosian mags are not saved, we need to calculate fluxes. For this we need exptime and zero point + float zeropt = psMetadataLookupF32(&status, hduHeader, "FPA.ZP"); + float exptime = psMetadataLookupF32(&status, hduHeader, "EXPTIME"); + float magOffset = zeropt + 2.5*log10(exptime); + + for (long i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + if (!row) { + psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); + psFree(row); + return false; + } + // Find the source with this sequence number. + // XXX: I am assuming that sources is sorted in order of seq + long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + pmSource *source = NULL; +#ifndef ASSUME_SORTED + long j = seq < sources->n ? seq : sources->n - 1; + for (; j >= 0; j--) { + source = sources->data[j]; + if (source->seq == seq) { + break; + } + } +#else + long j = sourceIndex[seq]; + psAssert(j >= 0 && j < sources->n, "invalid sourceIndex"); + source = sources->data[j]; +#endif + if (!source) { + psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); + psFree(row); + return false; + } + + if (!source->extpars) { + source->extpars = pmSourceExtendedParsAlloc (); + } + pmSourceExtendedPars *extpars = source->extpars; + + // Assume that X_EXT Y_EXT and sigmas match the psf src so skip + + // We don't have enough information to calculate the major and minor axis. Set major to 1. Should we scale this by + // psf size or something? + extpars->axes.major = 1.0; + extpars->axes.minor = extpars->axes.major * psMetadataLookupF32(&status, row, "F25_ARATIO"); + extpars->axes.theta = psMetadataLookupF32(&status, row, "F25_THETA"); + + float mag = psMetadataLookupF32(&status, row, "PETRO_MAG"); + float magErr = psMetadataLookupF32(&status, row, "PETRO_MAG_ERR"); + if (isfinite(mag)) { + extpars->petrosianFlux = pow(10., (magOffset - mag) / 2.5); + if (isfinite(magErr)) { + extpars->petrosianFluxErr = extpars->petrosianFlux / magErr; + } + } + + extpars->petrosianRadius = psMetadataLookupF32(&status, row, "PETRO_RADIUS"); + extpars->petrosianRadiusErr= psMetadataLookupF32(&status, row, "PETRO_RADIUS_ERR"); + extpars->petrosianR50 = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50"); + extpars->petrosianR50Err = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50_ERR"); + extpars->petrosianR90 = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90"); + extpars->petrosianR90Err = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90_ERR"); + extpars->petrosianFill = psMetadataLookupF32(&status, row, "PETRO_FILL"); + + psVector *radSB = psMetadataLookupVector(&status, row, "PROF_SB"); + psVector *radFlux = psMetadataLookupVector(&status, row, "PROF_FLUX"); + psVector *radFill = psMetadataLookupVector(&status, row, "PROF_FILL"); + + if (radSB && radSB->n > 0) { + extpars->radProfile = pmSourceRadialProfileAlloc(); + extpars->radProfile->binSB = psMemIncrRefCounter(radSB); + extpars->radProfile->binSum = psMemIncrRefCounter(radFlux); + extpars->radProfile->binFill = psMemIncrRefCounter(radFill); + } + + psFree(row); + } + return true; } @@ -572,5 +720,10 @@ int nParamMax = 0; for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; + // this is the source associated with this image + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + if (source->modelFits == NULL) continue; for (int j = 0; j < source->modelFits->n; j++) { @@ -586,5 +739,8 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; // XXX if no model fits are saved, write out modelEXT? @@ -641,7 +797,19 @@ // 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); + if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) { + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (major axis), length for trail", PAR[PM_PAR_LENGTH]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (minor axis), sigma for trail", PAR[PM_PAR_SIGMA]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", PAR[PM_PAR_THETA]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)", dPAR[PM_PAR_LENGTH]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)", dPAR[PM_PAR_SIGMA]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_THETA]); + } else { + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (major axis), length for trail", axes.major); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (minor axis), sigma for trail", axes.minor); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)", dPAR[PM_PAR_LENGTH]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)", dPAR[PM_PAR_SIGMA]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_THETA]); + } // write out the other generic parameters @@ -658,14 +826,21 @@ if (k < model->params->n) { - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); + psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]); } else { - psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); + psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", NAN); } } - // XXX other parameters which may be set. - // XXX flag / value to define the model - // XXX write out the model type, fit status flags - + // optionally, write out the covariance matrix values + // XXX do I need to pad this to match the biggest covar matrix? + if (model->covar) { + for (int iy = 0; iy < model->covar->numCols; iy++) { + for (int ix = iy; ix < model->covar->numCols; ix++) { + snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix); + psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]); + + } + } + } psArrayAdd (table, 100, row); psFree (row); @@ -697,6 +872,338 @@ } +bool pmSourcesRead_CMF_@CMFMODE@_XFIT(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + + bool status; + long numSources = psFitsTableSize(fits); // Number of sources in table + if (numSources == 0) { + psError(psErrorCodeLast(), false, "XFIT Table contains no entries\n"); + return false; + } + + for (long i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + if (!row) { + psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); + psFree(row); + return false; + } + // Find the source with this sequence number. + // XXX: I am assuming that sources is sorted in order of seq. + long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + long j = seq < sources->n ? seq : sources->n - 1; + pmSource *source = NULL; + for (; j >= 0; j--) { + source = sources->data[j]; + if (source->seq == seq) { + break; + } + } + if (!source) { + psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); + psFree(row); + return false; + } + if (!source->modelFits) { + // XXX: where to find the number of models to expect? + source->modelFits = psArrayAllocEmpty(5); + } + psString modelName = psMetadataLookupStr(&status, row, "MODEL_TYPE"); + if (!modelName) { + psError(PS_ERR_UNKNOWN, true, "Failed to find model name for row %ld\n", i); + psFree(row); + return false; + } + pmModelType modelType = pmModelClassGetType(modelName); + if (modelType < 0) { + psError(PS_ERR_UNKNOWN, true, "Failed to find model type for %s\n", modelName); + psFree(row); + return false; + } + pmModel *model = pmModelAlloc(modelType); + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + PAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT"); + PAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT"); + dPAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT_SIG"); + dPAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT_SIG"); + + model->mag = psMetadataLookupF32(&status, row, "EXT_INST_MAG"); + model->magErr = psMetadataLookupF32(&status, row, "EXT_INST_MAG_SIG"); + + psEllipseAxes axes; + axes.major = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ"); + axes.minor = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN"); + axes.theta = psMetadataLookupF32(&status, row, "EXT_THETA"); + if (!pmPSF_AxesToModel(PAR, axes, modelType)) { + // Do we need to fail here or can this happen? + psError(PS_ERR_UNKNOWN, false, "Failed to convert psf axes to model"); + psFree(model); + psFree(row); + return false; + } + // XXX: clean this up + if (model->params->n > 7) { + PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07"); + } + // read the covariance matrix + int nparams = model->params->n; + psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32); + for (int y = 0; y < nparams; y++) { + for (int x = 0; x < nparams; x++) { + char name[64]; + snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x); + covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name); + } + } + model->covar = covar; + + psArrayAdd(source->modelFits, 1, model); + psFree(model); + + psFree(row); + } + + 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_@CMFMODE@_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) { + bool status = false; + psArray *table; + psMetadata *row; + psF32 xPos, yPos; + char keyword1[80], keyword2[80]; + + // perform full non-linear fits / extended source analysis? + if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) { + psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n"); + return true; + } + + // 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) + 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 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]); + psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]); + } + + // the FWHM values are available if we measured a psf-matched convolved set + psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES"); + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortByFlux); + + 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++) { + + // this is the source associated with this image + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + + // 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 < source->radialAper->n; entry++) { + + // choose the convolved EXT model, if available, otherwise the simple one + pmSourceRadialApertures *radialAper = source->radialAper->data[entry]; + assert (radialAper); + + if (pmSourcePositionUseMoments(source)) { + 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); + if (fwhmValues) { + psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM", 0, "FWHM of matched PSF", fwhmValues->data.F32[entry]); + } else { + psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM", 0, "image is not FWHM-matched", NAN); + } + + // 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); + psVector *radFluxStdev = 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]; + radFluxStdev->data.F32[j] = radialAper->fluxStdev->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_FLUX_STDEV", PS_DATA_VECTOR, "flux standard deviation", radFluxStdev); + psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); + psFree (radFlux); + psFree (radFluxErr); + psFree (radFluxStdev); + 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; } + +bool pmSourcesRead_CMF_@CMFMODE@_XRAD(psFits *fits, pmReadout *readout, psMetadata *hduHeader, psArray *sources, long *sourceIndex) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + + bool status; + long numSources = psFitsTableSize(fits); // Number of sources in table + if (numSources == 0) { + psError(psErrorCodeLast(), false, "XRAD Table contains no entries\n"); + return false; + } + + long seq_first = -1; + long seq_last = -1; + psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32); + long max_entries = -1; + long num_entries = -1; + + for (long i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + if (!row) { + psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); + psFree(row); + return false; + } + // Find the source with this sequence number. + // XXX: I am assuming that sources is sorted in order of seq. + long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + long j = seq < sources->n ? seq : sources->n - 1; + pmSource *source = NULL; + for (; j >= 0; j--) { + source = sources->data[j]; + if (source->seq == seq) { + break; + } + } + if (!source) { + psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); + psFree(row); + return false; + } + if (seq_first == -1) { + seq_first = seq; + } + if (seq == seq_first) { + psF32 value = psMetadataLookupF32(&status, row, "PSF_FWHM"); + psVectorAppend(fwhmValues, value); + } + if (seq == seq_last) { + num_entries++; + } else { + num_entries = 1; + seq_last = seq; + } + if (num_entries > max_entries) { + max_entries = num_entries; + } + + if (!source->radialAper) { + // XXX: where to find the number of models to expect? + source->radialAper = psArrayAllocEmpty(5); + } + pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc(); + + radialAper->flux = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX")); + radialAper->fluxStdev = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_STDEV")); + radialAper->fluxErr = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_ERR")); + radialAper->fill = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FILL")); + + psArrayAdd(source->radialAper, 1, radialAper); + + psFree(radialAper); + psFree(row); + } + + // check for consistency between the length of fwhmValues and the maximum number of entries for each row + if (fwhmValues->n != max_entries) { + psError(PS_ERR_PROGRAMMING, true, "number of PSF_FWHM values found %ld does not match expected number: %ld\n", + fwhmValues->n, max_entries); + psAssert(0, "fixme"); + } + + if (!readout->analysis) { + readout->analysis = psMetadataAlloc(); + } + + psMetadataAddVector(readout->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues); + psFree(fwhmValues); + + return true; +}