Index: branches/eam_branches/ipp-dev-20210817/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- branches/eam_branches/ipp-dev-20210817/psModules/src/objects/pmSourceIO_CMF.c.in (revision 41820) +++ branches/eam_branches/ipp-dev-20210817/psModules/src/objects/pmSourceIO_CMF.c.in (revision 41821) @@ -57,12 +57,30 @@ // followed by a zero-size matrix, followed by the table data -bool pmSourcesWrite_CMF_@CMFMODE@ (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) +/** this file is used to generate c-code for a number of different formats. the perl program + ** 'mksource.pl' reads this template file and generates an output file for one of the named + ** formats, e.g., PS1_V1 or PS1_DV4 (see list in mksource.pl). + + ** Any instances of the word @CMFMODE@ are replaced by the named version of the format. + + ** Any line which starts with a command of the form: @COMMAND[,COMMAND...]@ is kept or removed + ** from the output c-code depending on the commands. These may be: ALL (lines are kept for + ** all formats) or the name of a format. A line with just a plain format will be kept only + ** for that format. The format name may be optionally proceded by >, <, <=, >=, ! or =. In + ** these cases, the line is kept if the format is greater / lesser / not the same, etc as the + ** listed format. The sequence of formats is defined for series (V, SV, DV) in mksource.pl so + ** that, e.g., >PS1_V3 would include e.g., PS1_V4 and PS1_V5, but not PS1_DV4. The number in + ** a series may be replaced with ? to stand for the full set. The format + ** commands may be grouped on a line separated by commas. Note that only ! (not) is + ** restrictive: only this rule will reduce the number of matches. So for example, + ** @ALL,= PS1_DV3 will be matched by + ** ALL. But @ALL,!PS1_DV3@ will generate a line for any format except PS1_DV3. + + **/ + +bool pmSourcesWrite_CMF_@CMFMODE@_New (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) { PS_ASSERT_PTR_NON_NULL(fits, false); PS_ASSERT_PTR_NON_NULL(sources, false); PS_ASSERT_PTR_NON_NULL(extname, false); - - psArray *table; - psMetadata *row; pmChip *chip = readout->parent->parent; @@ -73,12 +91,10 @@ pmSource *source = (pmSource *) sources->data[0]; if (source->seq == -1) { - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortByFlux); } else { - sources = psArraySort (sources, pmSourceSortBySeq); - } - } - - table = psArrayAllocEmpty (sources->n); + sources = psArraySort (sources, pmSourceSortBySeq); + } + } float magOffset; @@ -88,12 +104,206 @@ pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); + // before we can generate the table structure, we need to know if certain fields were measured. + // this information is not available a priori, so we need to check all sources + + bool haveLensingOBJsmear = false; + bool haveLensingOBJshear = false; + bool haveLensingPSFsmear = false; + bool haveLensingPSFshear = false; + bool haveLensingPSF = false; + + for (int i = 0; (i < sources->n) && !haveLensingOBJsmear && !haveLensingOBJshear && !haveLensingPSFsmear && !haveLensingPSFshear && !haveLensingPSF; 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->lensingOBJ && source->lensingOBJ->smear) haveLensingOBJsmear = true; + if (source->lensingOBJ && source->lensingOBJ->shear) haveLensingOBJshear = true; + if (source->lensingPSF && source->lensingPSF->smear) haveLensingPSFsmear = true; + if (source->lensingPSF && source->lensingPSF->shear) haveLensingPSFshear = true; + if (source->lensingPSF ) haveLensingPSF = true; + } + + /************ generate the table columns *****************/ + + // before we allocate the table, generate an empty array of table columns and generate them + psArray *tableColumns = psArrayAllocEmpty (100); + + // add the named / typed columns to the collection of columns + @ALL@ psFitsTableColumnAdd (tableColumns, "IPP_IDET", PS_DATA_U32); // "IPP detection identifier index" + @ALL@ psFitsTableColumnAdd (tableColumns, "X_PSF", PS_DATA_F32); // "PSF x coordinate" + @ALL@ psFitsTableColumnAdd (tableColumns, "Y_PSF", PS_DATA_F32); // "PSF y coordinate" + @ALL@ psFitsTableColumnAdd (tableColumns, "X_PSF_SIG", PS_DATA_F32); // "Sigma in PSF x coordinate" + @ALL@ psFitsTableColumnAdd (tableColumns, "Y_PSF_SIG", PS_DATA_F32); // "Sigma in PSF y coordinate" + + // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision + @PS1_V1@ psFitsTableColumnAdd (tableColumns, "RA_PSF", PS_DATA_F32); // "PSF RA coordinate (degrees)" + @PS1_V1@ psFitsTableColumnAdd (tableColumns, "DEC_PSF", PS_DATA_F32); // "PSF DEC coordinate (degrees)" + + @ALL@ psFitsTableColumnAdd (tableColumns, "POSANGLE", PS_DATA_F32); // "position angle at source (degrees)" + @ALL@ psFitsTableColumnAdd (tableColumns, "PLTSCALE", PS_DATA_F32); // "plate scale at source (arcsec/pixel)" + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_INST_MAG", PS_DATA_F32); // "PSF fit instrumental magnitude" + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_INST_MAG_SIG", PS_DATA_F32); // "Sigma of PSF instrumental magnitude" + + @ALL,!PS1_V1,!PS1_V2@ psFitsTableColumnAdd (tableColumns, "PSF_INST_FLUX", PS_DATA_F32); // "PSF fit instrumental flux (counts)" + @ALL,!PS1_V1,!PS1_V2@ psFitsTableColumnAdd (tableColumns, "PSF_INST_FLUX_SIG", PS_DATA_F32); // "Sigma of PSF instrumental flux" + + @ALL@ psFitsTableColumnAdd (tableColumns, "AP_MAG", PS_DATA_F32); // "magnitude in standard aperture" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "AP_MAG_RAW", PS_DATA_F32); // "magnitude in reported aperture" + @ALL@ psFitsTableColumnAdd (tableColumns, "AP_MAG_RADIUS", PS_DATA_F32); // "radius used for aperture mags" + @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableColumnAdd (tableColumns, "AP_FLUX", PS_DATA_F32); // "instrumental flux in standard aperture" + @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableColumnAdd (tableColumns, "AP_FLUX_SIG", PS_DATA_F32); // "aperture flux error" + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "AP_NPIX", PS_DATA_S32); // "aperture unmasked pixels" + + @=PS1_V3,>PS1_SV1@ psFitsTableColumnAdd (tableColumns, "PEAK_FLUX_AS_MAG", PS_DATA_F32); // "Peak flux expressed as magnitude" + @ALL@ psFitsTableColumnAdd (tableColumns, "SKY", PS_DATA_F32); // "Sky level" + @ALL@ psFitsTableColumnAdd (tableColumns, "SKY_SIGMA", PS_DATA_F32); // "Sigma of sky level" + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_CHISQ", PS_DATA_F32); // "Chisq of PSF-fit" + @ALL@ psFitsTableColumnAdd (tableColumns, "CR_NSIGMA", PS_DATA_F32); // "Nsigma deviations from PSF to CF" + @ALL@ psFitsTableColumnAdd (tableColumns, "EXT_NSIGMA", PS_DATA_F32); // "Nsigma deviations from PSF to EXT" + + // PSF shape parameters: + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_MAJOR", PS_DATA_F32); // "PSF width (major axis)" + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_MINOR", PS_DATA_F32); // "PSF width (minor axis)" + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_THETA", PS_DATA_F32); // "PSF orientation angle" + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "PSF_CORE", PS_DATA_F32); // "k term if defined" + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "PSF_FWHM_MAJ", PS_DATA_F32); // "PSF FWHM (major axis)" + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "PSF_FWHM_MIN", PS_DATA_F32); // "PSF FWHM (minor axis)" + + // psf data quality + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_QF", PS_DATA_F32); // "PSF coverage/quality factor (bad)" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "PSF_QF_PERFECT", PS_DATA_F32); // "PSF coverage/quality factor (poor)" + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_NDOF", PS_DATA_S32); // "degrees of freedom" + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_NPIX", PS_DATA_S32); // "number of pixels in fit" + + @ALL@ psFitsTableColumnAdd (tableColumns, "MOMENTS_XX", PS_DATA_F32); // "second moments (X^2)" + @ALL@ psFitsTableColumnAdd (tableColumns, "MOMENTS_XY", PS_DATA_F32); // "second moments (X*Y)" + @ALL@ psFitsTableColumnAdd (tableColumns, "MOMENTS_YY", PS_DATA_F32); // "second moments (Y*Y)" + @>PS1_V2,PS1_SV?@ psFitsTableColumnAdd (tableColumns, "MOMENTS_M3C", PS_DATA_F32); // "third momemt cos theta" + @>PS1_V2,PS1_SV?@ psFitsTableColumnAdd (tableColumns, "MOMENTS_M3S", PS_DATA_F32); // "third momemt sin theta" + @>PS1_V2,PS1_SV?@ psFitsTableColumnAdd (tableColumns, "MOMENTS_M4C", PS_DATA_F32); // "fourth momemt cos theta" + @>PS1_V2,PS1_SV?@ psFitsTableColumnAdd (tableColumns, "MOMENTS_M4S", PS_DATA_F32); // "fourth momemt sin theta" + + // Lensing parameters are only written if they are measured + if (haveLensingOBJsmear) { + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X11_SM_OBJ", PS_DATA_F32); // "smear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X12_SM_OBJ", PS_DATA_F32); // "smear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X22_SM_OBJ", PS_DATA_F32); // "smear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_SM_OBJ", PS_DATA_F32); // "smear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_SM_OBJ", PS_DATA_F32); // "smear polarizability element (objects)" + } + if (haveLensingOBJshear) { + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X11_SH_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X12_SH_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X22_SH_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_SH_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_SH_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" + } + if (false) { + // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" + } + + if (haveLensingPSFsmear) { + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X11_SM_PSF", PS_DATA_F32); // "smear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X12_SM_PSF", PS_DATA_F32); // "smear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X22_SM_PSF", PS_DATA_F32); // "smear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_SM_PSF", PS_DATA_F32); // "smear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_SM_PSF", PS_DATA_F32); // "smear polarizability element (PSFs)" + } + if (haveLensingPSFshear) { + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X11_SH_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X12_SH_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X22_SH_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_SH_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_SH_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" + } + if (haveLensingPSF) { + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" + } + + // if lensing params exist also include the backmapped chipID and chip coordinates + if (haveLensingPSF) { + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_NUM", PS_DATA_S16); // "id of warp input chip" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_X", PS_DATA_S16); // "x coord in warp input chip" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_Y", PS_DATA_S16); // "y coord in warp input chip" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "PADDING3", PS_DATA_S16); // "more padding" + } + + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "MOMENTS_R1", PS_DATA_F32); // "first radial moment" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "MOMENTS_RH", PS_DATA_F32); // "half radial moment" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "KRON_FLUX", PS_DATA_F32); // "Kron Flux (in 2.5 R1)" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "KRON_FLUX_ERR", PS_DATA_F32); // "Kron Flux Error" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "KRON_FLUX_INNER", PS_DATA_F32); // "Kron Flux (in 1.0 R1)" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "KRON_FLUX_OUTER", PS_DATA_F32); // "Kron Flux (in 4.0 R1)" + + @>PS1_V3@ psFitsTableColumnAdd (tableColumns, "SKY_LIMIT_RAD", PS_DATA_F32); // "Radius where object hits sky" + @>PS1_V3@ psFitsTableColumnAdd (tableColumns, "SKY_LIMIT_FLUX", PS_DATA_F32); // "Flux / pix where object hits sky" + @>PS1_V3@ psFitsTableColumnAdd (tableColumns, "SKY_LIMIT_SLOPE", PS_DATA_F32); // "d(Flux/pix)/dRadius where object hits sky" + + @>PS1_DV4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_NUM", PS_DATA_S16); // "id of warp input chip" + @>PS1_DV4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_X", PS_DATA_S16); // "x coord in warp input chip" + @>PS1_DV4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_Y", PS_DATA_S16); // "y coord in warp input chip" + @>PS1_DV4@ psFitsTableColumnAdd (tableColumns, "PADDING3", PS_DATA_S16); // "more padding" + + @PS1_DV?@ psFitsTableColumnAdd (tableColumns, "DIFF_NPOS", PS_DATA_S32); // "nPos (n pix > 3 sigma)" + @PS1_DV?@ psFitsTableColumnAdd (tableColumns, "DIFF_FRATIO", PS_DATA_F32); // "fPos / (fPos + fNeg)" + @PS1_DV?@ psFitsTableColumnAdd (tableColumns, "DIFF_NRATIO_BAD", PS_DATA_F32); // "nPos / (nPos + nNeg)" + @PS1_DV?@ psFitsTableColumnAdd (tableColumns, "DIFF_NRATIO_MASK", PS_DATA_F32); // "nPos / (nPos + nMask)" + @PS1_DV?@ psFitsTableColumnAdd (tableColumns, "DIFF_NRATIO_ALL", PS_DATA_F32); // "nPos / (nGood + nMask + nBad)" + + @>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "DIFF_R_P", PS_DATA_F32); // "distance to positive match source" + @>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "DIFF_SN_P", PS_DATA_F32); // "signal-to-noise of pos match src" + @>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "DIFF_R_M", PS_DATA_F32); // "distance to negative match source" + @>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "DIFF_SN_M", PS_DATA_F32); // "signal-to-noise of neg match src" + + @ALL@ psFitsTableColumnAdd (tableColumns, "FLAGS", PS_DATA_U32); // "psphot analysis flags" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "FLAGS2", PS_DATA_U32); // "psphot analysis flags" + @PS1_V3,PS1_SV2@ psFitsTableColumnAdd (tableColumns, "PADDING2", PS_DATA_S32); // "more padding" + @ALL@ psFitsTableColumnAdd (tableColumns, "N_FRAMES", PS_DATA_U16); // "Number of frames overlapping source center" + + @ALL@ psFitsTableColumnAdd (tableColumns, "PADDING", PS_DATA_S16); // "more padding" + + // 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. + + // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not + // subtracted + + // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image + // edge; 3) any pixels in the 3x3 peak region are masked; + + // generate an FITS table using the array of columns defined above, with space for all rows + psFitsTable *table = psFitsTableCreate (tableColumns, sources->n); + + /************ write the data to the table *****************/ + // 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 (int i = 0; i < sources->n; i++) { - // this is the source associated with this image + // 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; + // 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 @@ -102,133 +312,360 @@ // generated on Alloc, and would thus be wrong for read in sources. if (source->seq == -1) { - source->seq = i; - } - - // set the 'best' values for various output fields: - pmSourceOutputs outputs; - pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); - - pmSourceOutputsMoments moments; - 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 (); + source->seq = i; + } + + // set the 'best' values for various output fields: + pmSourceOutputs outputs; + pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); + + pmSourceOutputsMoments moments; + pmSourceOutputsSetMoments (&moments, source); + + @PS1_DV?@ pmSourceDiffStats diffStats; + @PS1_DV?@ pmSourceDiffStatsInit(&diffStats); + @PS1_DV?@ if (source->diffStats) { diffStats = *source->diffStats;} // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation // This set of psMetadataAdd Entries marks the "----" "Start of the PSF segment" - @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); - - @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_SV?,>PS1_DV1@ 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", source->apRadius); - @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", source->apFlux); - @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", source->apFluxErr); - @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "AP_NPIX", PS_DATA_S32, "aperture unmasked pixels", source->apNpixels); - - @=PS1_V3,>PS1_SV1@ 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@ psFitsTableSetU32 (table, i, "IPP_IDET", source->seq); // "IPP detection identifier index", + @ALL@ psFitsTableSetF32 (table, i, "X_PSF", outputs.xPos); // "PSF x coordinate", + @ALL@ psFitsTableSetF32 (table, i, "Y_PSF", outputs.yPos); // "PSF y coordinate", + @ALL@ psFitsTableSetF32 (table, i, "X_PSF_SIG", outputs.xErr); // "Sigma in PSF x coordinate", + @ALL@ psFitsTableSetF32 (table, i, "Y_PSF_SIG", outputs.yErr); // "Sigma in PSF y coordinate", + + // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision + @PS1_V1@ psFitsTableSetF32 (table, i, "RA_PSF", outputs.ra); // "PSF RA coordinate (degrees)", + @PS1_V1@ psFitsTableSetF32 (table, i, "DEC_PSF", outputs.dec); // "PSF DEC coordinate (degrees)", + + @ALL@ psFitsTableSetF32 (table, i, "POSANGLE", outputs.posAngle); // "position angle at source (degrees)", + @ALL@ psFitsTableSetF32 (table, i, "PLTSCALE", outputs.pltScale); // "plate scale at source (arcsec/pixel)", + @ALL@ psFitsTableSetF32 (table, i, "PSF_INST_MAG", source->psfMag); // "PSF fit instrumental magnitude", + @ALL@ psFitsTableSetF32 (table, i, "PSF_INST_MAG_SIG", source->psfMagErr); // "Sigma of PSF instrumental magnitude", + + @ALL,!PS1_V1,!PS1_V2@ psFitsTableSetF32 (table, i, "PSF_INST_FLUX", source->psfFlux); // "PSF fit instrumental flux (counts)", + @ALL,!PS1_V1,!PS1_V2@ psFitsTableSetF32 (table, i, "PSF_INST_FLUX_SIG", source->psfFluxErr); // "Sigma of PSF instrumental flux", + + @ALL@ psFitsTableSetF32 (table, i, "AP_MAG", source->apMag); // "magnitude in standard aperture", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "AP_MAG_RAW", source->apMagRaw); // "magnitude in reported aperture", + @ALL@ psFitsTableSetF32 (table, i, "AP_MAG_RADIUS", source->apRadius); // "radius used for aperture mags", + @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableSetF32 (table, i, "AP_FLUX", source->apFlux); // "instrumental flux in standard aperture", + @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableSetF32 (table, i, "AP_FLUX_SIG", source->apFluxErr); // "aperture flux error", + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetS32 (table, i, "AP_NPIX", source->apNpixels); // "aperture unmasked pixels", + + @=PS1_V3,>PS1_SV1@ psFitsTableSetF32 (table, i, "PEAK_FLUX_AS_MAG", outputs.peakMag); // "Peak flux expressed as magnitude", + @ALL@ psFitsTableSetF32 (table, i, "SKY", source->sky); // "Sky level", + @ALL@ psFitsTableSetF32 (table, i, "SKY_SIGMA", source->skyErr); // "Sigma of sky level", + + @ALL@ psFitsTableSetF32 (table, i, "PSF_CHISQ", outputs.chisq); // "Chisq of PSF-fit", + @ALL@ psFitsTableSetF32 (table, i, "CR_NSIGMA", source->crNsigma); // "Nsigma deviations from PSF to CF", + @ALL@ psFitsTableSetF32 (table, i, "EXT_NSIGMA", source->extNsigma); // "Nsigma deviations from PSF to EXT", // PSF shape parameters: - @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@ psFitsTableSetF32 (table, i, "PSF_MAJOR", outputs.psfMajor); // "PSF width (major axis)", + @ALL@ psFitsTableSetF32 (table, i, "PSF_MINOR", outputs.psfMinor); // "PSF width (minor axis)", + @ALL@ psFitsTableSetF32 (table, i, "PSF_THETA", outputs.psfTheta); // "PSF orientation angle", + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetF32 (table, i, "PSF_CORE", outputs.psfCore); // "k term if defined", + + // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetF32 (table, i, "PSF_FWHM_MAJ", outputs.psfMajorFWHM); // "PSF FWHM (major axis)", + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetF32 (table, i, "PSF_FWHM_MIN", outputs.psfMinorFWHM); // "PSF FWHM (minor axis)", + + // psf data quality + @ALL@ psFitsTableSetF32 (table, i, "PSF_QF", source->pixWeightNotBad); // "PSF coverage/quality factor (bad)", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "PSF_QF_PERFECT", source->pixWeightNotPoor); // "PSF coverage/quality factor (poor)", + @ALL@ psFitsTableSetS32 (table, i, "PSF_NDOF", outputs.nDOF); // "degrees of freedom", + @ALL@ psFitsTableSetS32 (table, i, "PSF_NPIX", outputs.nPix); // "number of pixels in fit", + + @ALL@ psFitsTableSetF32 (table, i, "MOMENTS_XX", moments.Mxx); // "second moments (X^2)", + @ALL@ psFitsTableSetF32 (table, i, "MOMENTS_XY", moments.Mxy); // "second moments (X*Y)", + @ALL@ psFitsTableSetF32 (table, i, "MOMENTS_YY", moments.Myy); // "second moments (Y*Y)", + + @>PS1_V2,PS1_SV?@ psFitsTableSetF32 (table, i, "MOMENTS_M3C", moments.M_c3); // "third momemt cos theta", + @>PS1_V2,PS1_SV?@ psFitsTableSetF32 (table, i, "MOMENTS_M3S", moments.M_s3); // "third momemt sin theta", + @>PS1_V2,PS1_SV?@ psFitsTableSetF32 (table, i, "MOMENTS_M4C", moments.M_c4); // "fourth momemt cos theta", + @>PS1_V2,PS1_SV?@ psFitsTableSetF32 (table, i, "MOMENTS_M4S", moments.M_s4); // "fourth momemt sin theta", + + // Lensing parameters: + if (source->lensingOBJ && source->lensingOBJ->smear) { + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X11_SM_OBJ", source->lensingOBJ->smear->X11); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X12_SM_OBJ", source->lensingOBJ->smear->X12); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X22_SM_OBJ", source->lensingOBJ->smear->X22); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_SM_OBJ", source->lensingOBJ->smear->e1); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_SM_OBJ", source->lensingOBJ->smear->e2); // "smear polarizability element (objects)", + } + if (source->lensingOBJ && source->lensingOBJ->shear) { + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X11_SH_OBJ", source->lensingOBJ->shear->X11); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X12_SH_OBJ", source->lensingOBJ->shear->X12); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X22_SH_OBJ", source->lensingOBJ->shear->X22); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_SH_OBJ", source->lensingOBJ->shear->e1); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_SH_OBJ", source->lensingOBJ->shear->e2); // "shear polarizability element (objects)", + } + if (false && source->lensingOBJ) { + // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_PSF", source->lensingOBJ->e1); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_PSF", source->lensingOBJ->e2); // "shear polarizability element (objects)", + } + + if (source->lensingPSF && source->lensingPSF->smear) { + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X11_SM_PSF", source->lensingPSF->smear->X11); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X12_SM_PSF", source->lensingPSF->smear->X12); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X22_SM_PSF", source->lensingPSF->smear->X22); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_SM_PSF", source->lensingPSF->smear->e1); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_SM_PSF", source->lensingPSF->smear->e2); // "smear polarizability element (objects)", + } + if (source->lensingPSF && source->lensingPSF->shear) { + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X11_SH_PSF", source->lensingPSF->shear->X11); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X12_SH_PSF", source->lensingPSF->shear->X12); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X22_SH_PSF", source->lensingPSF->shear->X22); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_SH_PSF", source->lensingPSF->shear->e1); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_SH_PSF", source->lensingPSF->shear->e2); // "shear polarizability element (objects)", + } + if (source->lensingPSF) { + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_PSF", source->lensingPSF->e1); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_PSF", source->lensingPSF->e2); // "shear polarizability element (objects)", + } + + // if lensing params exist also include the backmapped chipID and chip coordinates + if (source->lensingPSF && source->lensingPSF->shear) { + @>PS1_V4@ psFitsTableSetS16 (table, i, "SRC_CHIP_NUM", source->chipNum); // "id of warp input chip", + @>PS1_V4@ psFitsTableSetS16 (table, i, "SRC_CHIP_X", source->chipX); // "x coord in warp input chip", + @>PS1_V4@ psFitsTableSetS16 (table, i, "SRC_CHIP_Y", source->chipY); // "y coord in warp input chip", + @>PS1_V4@ psFitsTableSetS16 (table, i, "PADDING3", 0); // "more padding", + } + + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "MOMENTS_R1", moments.Mrf); // "first radial moment", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "MOMENTS_RH", moments.Mrh); // "half radial moment", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "KRON_FLUX", moments.Krf); // "Kron Flux (in 2.5 R1)", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "KRON_FLUX_ERR", moments.dKrf); // "Kron Flux Error", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "KRON_FLUX_INNER", moments.Kinner); // "Kron Flux (in 1.0 R1)", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "KRON_FLUX_OUTER", moments.Kouter); // "Kron Flux (in 4.0 R1)", + @>PS1_V3@ psFitsTableSetF32 (table, i, "SKY_LIMIT_RAD", source->skyRadius); // "Radius where object hits sky", + @>PS1_V3@ psFitsTableSetF32 (table, i, "SKY_LIMIT_FLUX", source->skyFlux); // "Flux / pix where object hits sky", + @>PS1_V3@ psFitsTableSetF32 (table, i, "SKY_LIMIT_SLOPE", source->skySlope); // "d(Flux/pix)/dRadius where object hits sky", + + @>PS1_DV4@ psFitsTableSetS16 (table, i, "SRC_CHIP_NUM", source->chipNum); // "id of warp input chip", + @>PS1_DV4@ psFitsTableSetS16 (table, i, "SRC_CHIP_X", source->chipX); // "x coord in warp input chip", + @>PS1_DV4@ psFitsTableSetS16 (table, i, "SRC_CHIP_Y", source->chipY); // "y coord in warp input chip", + @>PS1_DV4@ psFitsTableSetS16 (table, i, "PADDING3", 0); // "more padding", + + @PS1_DV?@ psFitsTableSetS32 (table, i, "DIFF_NPOS", diffStats.nGood); // "nPos (n pix > 3 sigma)", + @PS1_DV?@ psFitsTableSetF32 (table, i, "DIFF_FRATIO", diffStats.fRatio); // "fPos / (fPos + fNeg)", + @PS1_DV?@ psFitsTableSetF32 (table, i, "DIFF_NRATIO_BAD", diffStats.nRatioBad); // "nPos / (nPos + nNeg)", + @PS1_DV?@ psFitsTableSetF32 (table, i, "DIFF_NRATIO_MASK", diffStats.nRatioMask); // "nPos / (nPos + nMask)", + @PS1_DV?@ psFitsTableSetF32 (table, i, "DIFF_NRATIO_ALL", diffStats.nRatioAll); // "nPos / (nGood + nMask + nBad)", + + @>PS1_DV1@ psFitsTableSetF32 (table, i, "DIFF_R_P", diffStats.Rp); // "distance to positive match source", + @>PS1_DV1@ psFitsTableSetF32 (table, i, "DIFF_SN_P", diffStats.SNp); // "signal-to-noise of pos match src", + @>PS1_DV1@ psFitsTableSetF32 (table, i, "DIFF_R_M", diffStats.Rm); // "distance to negative match source", + @>PS1_DV1@ psFitsTableSetF32 (table, i, "DIFF_SN_M", diffStats.SNm); // "signal-to-noise of neg match src", + + @ALL@ psFitsTableSetU32 (table, i, "FLAGS", source->mode); // "psphot analysis flags", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetU32 (table, i, "FLAGS2", source->mode2); // "psphot analysis flags", + @PS1_V3,PS1_SV2@ psFitsTableSetS32 (table, i, "PADDING2", 0); // "padding", + + @ALL@ psFitsTableSetU16 (table, i, "N_FRAMES", source->nFrames); // "Number of frames overlapping source center", + @ALL@ psFitsTableSetS16 (table, i, "PADDING", 0); // "padding", + } + + // 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); + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTableNew(fits, header, table, extname)) { + psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); + psFree(table); + psFree(header); + return false; + } + psFree(table); + psFree(header); + + return true; +} + +bool pmSourcesWrite_CMF_@CMFMODE@_Old (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + PS_ASSERT_PTR_NON_NULL(extname, false); + + psArray *table; + psMetadata *row; + + pmChip *chip = readout->parent->parent; + + // if the sequence is defined, write these in seq order; otherwise + // write them in S/N order: + if (sources->n > 0) { + pmSource *source = (pmSource *) sources->data[0]; + if (source->seq == -1) { + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortByFlux); + } else { + sources = psArraySort (sources, pmSourceSortBySeq); + } + } + + table = psArrayAllocEmpty (sources->n); + + float magOffset; + float zeroptErr; + float fwhmMajor; + float fwhmMinor; + pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); + + // 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 (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; + + // If source->seq is -1, source was generated in this analysis. If source->seq is + // not -1, source was read from elsewhere: in the latter case, preserve the source + // ID. source.seq is used instead of source.id since the latter is a const + // generated on Alloc, and would thus be wrong for read in sources. + if (source->seq == -1) { + source->seq = i; + } + + // set the 'best' values for various output fields: + pmSourceOutputs outputs; + pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); + + pmSourceOutputsMoments moments; + 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 (); + + // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation + // This set of psMetadataAdd Entries marks the "----" "Start of the PSF segment" + @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); + + @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_SV?,>PS1_DV1@ 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", source->apRadius); + @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", source->apFlux); + @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", source->apFluxErr); + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "AP_NPIX", PS_DATA_S32, "aperture unmasked pixels", source->apNpixels); + + @=PS1_V3,>PS1_SV1@ 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); + + // PSF shape parameters: + @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); @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_CORE", PS_DATA_F32, "k term if defined", outputs.psfCore); - // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values + // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_FWHM_MAJ", PS_DATA_F32, "PSF FWHM (major axis)", outputs.psfMajorFWHM); @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_FWHM_MIN", PS_DATA_F32, "PSF FWHM (minor axis)", outputs.psfMinorFWHM); // psf data quality - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); - @>PS1_V2,PS1_SV?,>PS1_DV1@ 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_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); - @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); - @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); - @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); - - // Lensing parameters: - if (source->lensingOBJ && source->lensingOBJ->smear) { - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->X11); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->X12); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->X22); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->e1); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->e2); - } - if (source->lensingOBJ && source->lensingOBJ->shear) { - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->X11); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->X12); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->X22); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->e1); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->e2); - } - if (false && source->lensingOBJ) { - // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->e1); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->e2); - } - - if (source->lensingPSF && source->lensingPSF->smear) { - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X11); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X12); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X22); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->e1); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->e2); - } - if (source->lensingPSF && source->lensingPSF->shear) { - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X11); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X12); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X22); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->e1); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->e2); - } - if (source->lensingPSF) { - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->e1); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->e2); - } + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); + @>PS1_V2,PS1_SV?,>PS1_DV1@ 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_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); + @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); + @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); + @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); + + // Lensing parameters: + if (source->lensingOBJ && source->lensingOBJ->smear) { + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->X11); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->X12); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->X22); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->e1); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->e2); + } + if (source->lensingOBJ && source->lensingOBJ->shear) { + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->X11); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->X12); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->X22); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->e1); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->e2); + } + if (false && source->lensingOBJ) { + // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->e1); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->e2); + } + + if (source->lensingPSF && source->lensingPSF->smear) { + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X11); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X12); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X22); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->e1); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->e2); + } + if (source->lensingPSF && source->lensingPSF->shear) { + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X11); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X12); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X22); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->e1); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->e2); + } + if (source->lensingPSF) { + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->e1); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->e2); + } // if lensing params exist also include the backmapped chipID and chip coordinates - if (source->lensingPSF && source->lensingPSF->shear) { + if (source->lensingPSF && source->lensingPSF->shear) { @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM", PS_DATA_S16, "id of warp input chip", source->chipNum); @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X", PS_DATA_S16, "x coord in warp input chip", source->chipX); @@ -244,40 +681,40 @@ @>PS1_V2,PS1_SV?,>PS1_DV1@ 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_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM", PS_DATA_S16, "id of warp input chip", source->chipNum); - @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X", PS_DATA_S16, "x coord in warp input chip", source->chipX); - @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_Y", PS_DATA_S16, "y coord in warp input chip", source->chipY); - @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING3", PS_DATA_S16, "more padding", 0); - - @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_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P", PS_DATA_F32, "distance to positive match source", diffStats.Rp); - @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P", PS_DATA_F32, "signal-to-noise of pos match src", diffStats.SNp); - @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M", PS_DATA_F32, "distance to negative match source", diffStats.Rm); - @>PS1_DV1@ 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_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); - @PS1_V3,PS1_SV2@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2", PS_DATA_S32, "more padding", 0); - @PS1_SV?@ - - // 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); + @>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_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM", PS_DATA_S16, "id of warp input chip", source->chipNum); + @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X", PS_DATA_S16, "x coord in warp input chip", source->chipX); + @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_Y", PS_DATA_S16, "y coord in warp input chip", source->chipY); + @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING3", PS_DATA_S16, "more padding", 0); + + @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_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P", PS_DATA_F32, "distance to positive match source", diffStats.Rp); + @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P", PS_DATA_F32, "signal-to-noise of pos match src", diffStats.SNp); + @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M", PS_DATA_F32, "distance to negative match source", diffStats.Rm); + @>PS1_DV1@ 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_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); + @PS1_V3,PS1_SV2@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2", PS_DATA_S32, "more padding", 0); + @PS1_SV?@ + + // 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); @@ -320,6 +757,13 @@ } +bool pmSourcesWrite_CMF_@CMFMODE@ (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) +{ + // bool status = pmSourcesWrite_CMF_@CMFMODE@_Old (fits, readout, sources, imageHeader, tableHeader, extname, recipe); + bool status = pmSourcesWrite_CMF_@CMFMODE@_New (fits, readout, sources, imageHeader, tableHeader, extname, recipe); + return status; +} + // read in a readout from the fits file -psArray *pmSourcesRead_CMF_@CMFMODE@ (psFits *fits, psMetadata *header) +psArray *pmSourcesRead_CMF_@CMFMODE@_New (psFits *fits, psMetadata *header) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -355,5 +799,214 @@ psArray *sources = psArrayAlloc(numSources); // Array of sources, to return - // convert the table to the pmSource entriesa + // reads full table into memory + psFitsTable *table = psFitsReadTableNew (fits); + + // convert the table to the pmSource entries + for (int i = 0; i < numSources; i++) { + pmSource *source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags + + // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + @ALL@ source->seq = psFitsTableGetU32 (&status, table, i, "IPP_IDET"); + @ALL@ PAR[PM_PAR_XPOS] = psFitsTableGetF32 (&status, table, i, "X_PSF"); + @ALL@ PAR[PM_PAR_YPOS] = psFitsTableGetF32 (&status, table, i, "Y_PSF"); + @ALL@ dPAR[PM_PAR_XPOS] = psFitsTableGetF32 (&status, table, i, "X_PSF_SIG"); + @ALL@ dPAR[PM_PAR_YPOS] = psFitsTableGetF32 (&status, table, i, "Y_PSF_SIG"); + @ALL@ axes.major = psFitsTableGetF32 (&status, table, i, "PSF_MAJOR"); + @ALL@ axes.minor = psFitsTableGetF32 (&status, table, i, "PSF_MINOR"); + @ALL@ axes.theta = psFitsTableGetF32 (&status, table, i, "PSF_THETA"); + @ALL@ axes.theta = axes.theta * PS_RAD_DEG; + + @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) { + @>PS1_V4,>PS1_SV2,>PS1_DV3@ PAR[PM_PAR_7] = psFitsTableGetF32 (&status, table, i, "PSF_CORE"); + @>PS1_V4,>PS1_SV2,>PS1_DV3@ } + + @ALL@ PAR[PM_PAR_SKY] = psFitsTableGetF32 (&status, table, i, "SKY"); + @ALL@ dPAR[PM_PAR_SKY] = psFitsTableGetF32 (&status, table, i, "SKY_SIGMA"); + @ALL@ source->sky = PAR[PM_PAR_SKY]; + @ALL@ source->skyErr = dPAR[PM_PAR_SKY]; + + // XXX use these to determine PAR[PM_PAR_I0]? + @ALL@ source->psfMag = psFitsTableGetF32 (&status, table, i, "PSF_INST_MAG"); + @ALL@ source->psfMagErr = psFitsTableGetF32 (&status, table, i, "PSF_INST_MAG_SIG"); + @ALL@ source->apMag = psFitsTableGetF32 (&status, table, i, "AP_MAG"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->apMagRaw = psFitsTableGetF32 (&status, table, i, "AP_MAG_RAW"); + @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux = psFitsTableGetF32 (&status, table, i, "AP_FLUX"); + @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr = psFitsTableGetF32 (&status, table, i, "AP_FLUX_SIG"); + + // XXX use these to determine PAR[PM_PAR_I0] if they exist? + // XXX add these to PS1_SV1? + @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFlux = psFitsTableGetF32 (&status, table, i, "PSF_INST_FLUX"); + @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFluxErr = psFitsTableGetF32 (&status, table, i, "PSF_INST_FLUX_SIG"); + + // XXX this scaling is incorrect: does not include the 2 \pi AREA factor + @ALL@ PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; + @ALL@ dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; + + pmPSF_AxesToModel (PAR, axes, model->class->useReff); + + @ALL@ float peakMag = psFitsTableGetF32 (&status, table, i, "PEAK_FLUX_AS_MAG"); + @ALL@ float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; + + // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) + @ALL@ source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); + @ALL@ source->peak->rawFlux = peakFlux; + @ALL@ source->peak->smoothFlux = peakFlux; + @ALL@ source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position + @ALL@ source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position + @ALL@ source->peak->dx = dPAR[PM_PAR_XPOS]; + @ALL@ source->peak->dy = dPAR[PM_PAR_YPOS]; + + @ALL@ source->pixWeightNotBad = psFitsTableGetF32 (&status, table, i, "PSF_QF"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->pixWeightNotPoor = psFitsTableGetF32 (&status, table, i, "PSF_QF_PERFECT"); + @ALL@ source->crNsigma = psFitsTableGetF32 (&status, table, i, "CR_NSIGMA"); + @ALL@ source->extNsigma = psFitsTableGetF32 (&status, table, i, "EXT_NSIGMA"); + @ALL@ source->apRadius = psFitsTableGetF32 (&status, table, i, "AP_MAG_RADIUS"); + @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels = psFitsTableGetS32 (&status, table, i, "AP_NPIX"); + + // note that some older versions used PSF_PROBABILITY: this was not well defined. + @ALL@ model->chisq = psFitsTableGetF32 (&status, table, i, "PSF_CHISQ"); + @ALL@ model->nDOF = psFitsTableGetS32 (&status, table, i, "PSF_NDOF"); + @ALL@ model->nPix = psFitsTableGetS32 (&status, table, i, "PSF_NPIX"); + + @ALL@ source->moments = pmMomentsAlloc (); + @ALL@ source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf + @ALL@ source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf + + @ALL@ source->moments->Mxx = psFitsTableGetF32 (&status, table, i, "MOMENTS_XX"); + @ALL@ source->moments->Mxy = psFitsTableGetF32 (&status, table, i, "MOMENTS_XY"); + @ALL@ source->moments->Myy = psFitsTableGetF32 (&status, table, i, "MOMENTS_YY"); + + // 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_SV?@ source->moments->Mxxx = +1.00 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M3C"); + @>PS1_V2,PS1_SV?@ source->moments->Mxxy = 0.00; + @>PS1_V2,PS1_SV?@ source->moments->Mxyy = 0.00; + @>PS1_V2,PS1_SV?@ source->moments->Myyy = -1.00 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M3S"); + @>PS1_V2,PS1_SV?@ source->moments->Mxxxx = +1.00 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M4C"); + @>PS1_V2,PS1_SV?@ source->moments->Mxxxy = 0.00; + @>PS1_V2,PS1_SV?@ source->moments->Mxxyy = 0.00; + @>PS1_V2,PS1_SV?@ source->moments->Mxyyy = -0.25 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M4S"); + @>PS1_V2,PS1_SV?@ source->moments->Myyyy = 0.00; + + // Lensing parameters (on read if PS1_V5+) + if (haveLensOBJ) { + source->lensingOBJ = pmSourceLensingAlloc (); + source->lensingOBJ->smear = pmLensingParsAlloc(); + source->lensingOBJ->shear = pmLensingParsAlloc(); + + @>PS1_V4@ source->lensingOBJ->smear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->e1 = psFitsTableGetF32 (&status, table, i, "LENS_E1_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->e2 = psFitsTableGetF32 (&status, table, i, "LENS_E2_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->e1 = psFitsTableGetF32 (&status, table, i, "LENS_E1_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->e2 = psFitsTableGetF32 (&status, table, i, "LENS_E2_SH_OBJ"); + } + + @>PS1_V4@ source->chipNum = psFitsTableGetS16 (&status, table, i, "SRC_CHIP_NUM"); + @>PS1_V4@ source->chipX = psFitsTableGetS16 (&status, table, i, "SRC_CHIP_X"); + @>PS1_V4@ source->chipY = psFitsTableGetS16 (&status, table, i, "SRC_CHIP_Y"); + + if (haveLensPSF) { + source->lensingPSF = pmSourceLensingAlloc (); + source->lensingPSF->smear = pmLensingParsAlloc(); + source->lensingPSF->shear = pmLensingParsAlloc(); + + @>PS1_V4@ source->lensingPSF->smear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->e1 = psFitsTableGetF32 (&status, table, i, "LENS_E1_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->e2 = psFitsTableGetF32 (&status, table, i, "LENS_E2_SM_PSF"); + @>PS1_V4@ source->lensingPSF->shear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->e1 = psFitsTableGetF32 (&status, table, i, "LENS_E1_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->e2 = psFitsTableGetF32 (&status, table, i, "LENS_E2_SH_PSF"); + @>PS1_V4@ source->lensingPSF->e1 = psFitsTableGetF32 (&status, table, i, "LENS_E1_PSF"); + @>PS1_V4@ source->lensingPSF->e2 = psFitsTableGetF32 (&status, table, i, "LENS_E2_PSF"); + } + + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrf = psFitsTableGetF32 (&status, table, i, "MOMENTS_R1"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrh = psFitsTableGetF32 (&status, table, i, "MOMENTS_RH"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFlux = psFitsTableGetF32 (&status, table, i, "KRON_FLUX"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFluxErr = psFitsTableGetF32 (&status, table, i, "KRON_FLUX_ERR"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFinner = psFitsTableGetF32 (&status, table, i, "KRON_FLUX_INNER"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFouter = psFitsTableGetF32 (&status, table, i, "KRON_FLUX_OUTER"); + + @>PS1_V3@ source->skyRadius = psFitsTableGetF32 (&status, table, i, "SKY_LIMIT_RAD"); + @>PS1_V3@ source->skyFlux = psFitsTableGetF32 (&status, table, i, "SKY_LIMIT_FLUX"); + @>PS1_V3@ source->skySlope = psFitsTableGetF32 (&status, table, i, "SKY_LIMIT_SLOPE"); + + @PS1_DV?@ int nPos = psFitsTableGetS32 (&status, table, i, "DIFF_NPOS"); + @PS1_DV?@ if (nPos) { + @PS1_DV?@ source->diffStats = pmSourceDiffStatsAlloc(); + @PS1_DV?@ source->diffStats->nGood = nPos; + @PS1_DV?@ source->diffStats->fRatio = psFitsTableGetF32 (&status, table, i, "DIFF_FRATIO"); + @PS1_DV?@ source->diffStats->nRatioBad = psFitsTableGetF32 (&status, table, i, "DIFF_NRATIO_BAD"); + @PS1_DV?@ source->diffStats->nRatioMask = psFitsTableGetF32 (&status, table, i, "DIFF_NRATIO_MASK"); + @PS1_DV?@ source->diffStats->nRatioAll = psFitsTableGetF32 (&status, table, i, "DIFF_NRATIO_ALL"); + + @>PS1_DV1@ source->diffStats->Rp = psFitsTableGetF32 (&status, table, i, "DIFF_R_P"); + @>PS1_DV1@ source->diffStats->SNp = psFitsTableGetF32 (&status, table, i, "DIFF_SN_P"); + @>PS1_DV1@ source->diffStats->Rm = psFitsTableGetF32 (&status, table, i, "DIFF_R_M"); + @>PS1_DV1@ source->diffStats->SNm = psFitsTableGetF32 (&status, table, i, "DIFF_SN_M"); + @PS1_DV?@ } + + @ALL@ source->mode = psFitsTableGetU32 (&status, table, i, "FLAGS"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->mode2 = psFitsTableGetU32 (&status, table, i, "FLAGS2"); + @ALL@ source->nFrames = psFitsTableGetU16 (&status, table, i, "N_FRAMES"); + assert (status); + + sources->data[i] = source; + } + psFree (table); + return sources; +} + +// read in a readout from the fits file +psArray *pmSourcesRead_CMF_@CMFMODE@_Old (psFits *fits, psMetadata *header) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(header, false); + + bool status; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + + // define PSF model type + int defaultModelType = pmModelClassGetType ("PS_MODEL_GAUSS"); + int modelType = -1; + + // if header does not define the model, default to a gaussian + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSFMODEL"); + if (PSF_NAME != NULL) { + modelType = pmModelClassGetType (PSF_NAME); + } + // work around bug in psphotFullForce + if (modelType < 0) { + modelType = defaultModelType; + } + // assert (modelType > -1); + + // do we expect to find lensing parameters? + bool haveLensOBJ = psMetadataLookupBool (&status, header, "LENS_OBJ"); + bool haveLensPSF = psMetadataLookupBool (&status, header, "LENS_PSF"); + + // We get the size of the table, and allocate the array of sources first because the table + // is large and ephemeral --- when the table gets blown away, whatever is allocated after + // the table is read blocks the free. In fact, it's better to read the table row by row. + long numSources = psFitsTableSize(fits); // Number of sources in table + psArray *sources = psArrayAlloc(numSources); // Array of sources, to return + + // convert the table to the pmSource entries for (int i = 0; i < numSources; i++) { psMetadata *row = psFitsReadTableRow(fits, i); // Table row @@ -382,8 +1035,8 @@ @ALL@ axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); @ALL@ axes.theta = axes.theta * PS_RAD_DEG; - - @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) { - @>PS1_V4,>PS1_SV2,>PS1_DV3@ PAR[PM_PAR_7] = psMetadataLookupF32 (&status, row, "PSF_CORE"); - @>PS1_V4,>PS1_SV2,>PS1_DV3@ } + + @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) { + @>PS1_V4,>PS1_SV2,>PS1_DV3@ PAR[PM_PAR_7] = psMetadataLookupF32 (&status, row, "PSF_CORE"); + @>PS1_V4,>PS1_SV2,>PS1_DV3@ } @ALL@ PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); @@ -396,12 +1049,12 @@ @ALL@ source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); @ALL@ source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); - @>PS1_V2,PS1_SV?,>PS1_DV1@ source->apMagRaw = psMetadataLookupF32 (&status, row, "AP_MAG_RAW"); - @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux = psMetadataLookupF32 (&status, row, "AP_FLUX"); - @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr = psMetadataLookupF32 (&status, row, "AP_FLUX_SIG"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->apMagRaw = psMetadataLookupF32 (&status, row, "AP_MAG_RAW"); + @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux = psMetadataLookupF32 (&status, row, "AP_FLUX"); + @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr = psMetadataLookupF32 (&status, row, "AP_FLUX_SIG"); // XXX use these to determine PAR[PM_PAR_I0] if they exist? - // XXX add these to PS1_SV1? - @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); - @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); + // XXX add these to PS1_SV1? + @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); + @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); // XXX this scaling is incorrect: does not include the 2 \pi AREA factor @@ -424,9 +1077,9 @@ @ALL@ source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); - @>PS1_V2,PS1_SV?,>PS1_DV1@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT"); @ALL@ source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); @ALL@ source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); @ALL@ source->apRadius = psMetadataLookupF32 (&status, row, "AP_MAG_RADIUS"); - @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels = psMetadataLookupS32 (&status, row, "AP_NPIX"); + @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels = psMetadataLookupS32 (&status, row, "AP_NPIX"); // note that some older versions used PSF_PROBABILITY: this was not well defined. @@ -443,6 +1096,6 @@ @ALL@ source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); - // 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 + // 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_SV?@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C"); @>PS1_V2,PS1_SV?@ source->moments->Mxxy = 0.0; @@ -456,44 +1109,44 @@ @>PS1_V2,PS1_SV?@ source->moments->Myyyy = 0.0; - // Lensing parameters (on read if PS1_V5+) - if (haveLensOBJ) { - source->lensingOBJ = pmSourceLensingAlloc (); - source->lensingOBJ->smear = pmLensingParsAlloc(); - source->lensingOBJ->shear = pmLensingParsAlloc(); - - @>PS1_V4@ source->lensingOBJ->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_OBJ"); - @>PS1_V4@ source->lensingOBJ->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_OBJ"); - @>PS1_V4@ source->lensingOBJ->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_OBJ"); - @>PS1_V4@ source->lensingOBJ->smear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SM_OBJ"); - @>PS1_V4@ source->lensingOBJ->smear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SM_OBJ"); - @>PS1_V4@ source->lensingOBJ->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_OBJ"); - @>PS1_V4@ source->lensingOBJ->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_OBJ"); - @>PS1_V4@ source->lensingOBJ->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_OBJ"); - @>PS1_V4@ source->lensingOBJ->shear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SH_OBJ"); - @>PS1_V4@ source->lensingOBJ->shear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SH_OBJ"); - } - - @>PS1_V4@ source->chipNum = psMetadataLookupS16 (&status, row, "SRC_CHIP_NUM"); - @>PS1_V4@ source->chipX = psMetadataLookupS16 (&status, row, "SRC_CHIP_X"); - @>PS1_V4@ source->chipY = psMetadataLookupS16 (&status, row, "SRC_CHIP_Y"); - - if (haveLensPSF) { - source->lensingPSF = pmSourceLensingAlloc (); - source->lensingPSF->smear = pmLensingParsAlloc(); - source->lensingPSF->shear = pmLensingParsAlloc(); - - @>PS1_V4@ source->lensingPSF->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_PSF"); - @>PS1_V4@ source->lensingPSF->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_PSF"); - @>PS1_V4@ source->lensingPSF->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_PSF"); - @>PS1_V4@ source->lensingPSF->smear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SM_PSF"); - @>PS1_V4@ source->lensingPSF->smear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SM_PSF"); - @>PS1_V4@ source->lensingPSF->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_PSF"); - @>PS1_V4@ source->lensingPSF->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_PSF"); - @>PS1_V4@ source->lensingPSF->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_PSF"); - @>PS1_V4@ source->lensingPSF->shear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SH_PSF"); - @>PS1_V4@ source->lensingPSF->shear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SH_PSF"); - @>PS1_V4@ source->lensingPSF->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_PSF"); - @>PS1_V4@ source->lensingPSF->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_PSF"); - } + // Lensing parameters (on read if PS1_V5+) + if (haveLensOBJ) { + source->lensingOBJ = pmSourceLensingAlloc (); + source->lensingOBJ->smear = pmLensingParsAlloc(); + source->lensingOBJ->shear = pmLensingParsAlloc(); + + @>PS1_V4@ source->lensingOBJ->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SH_OBJ"); + } + + @>PS1_V4@ source->chipNum = psMetadataLookupS16 (&status, row, "SRC_CHIP_NUM"); + @>PS1_V4@ source->chipX = psMetadataLookupS16 (&status, row, "SRC_CHIP_X"); + @>PS1_V4@ source->chipY = psMetadataLookupS16 (&status, row, "SRC_CHIP_Y"); + + if (haveLensPSF) { + source->lensingPSF = pmSourceLensingAlloc (); + source->lensingPSF->smear = pmLensingParsAlloc(); + source->lensingPSF->shear = pmLensingParsAlloc(); + + @>PS1_V4@ source->lensingPSF->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SM_PSF"); + @>PS1_V4@ source->lensingPSF->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SH_PSF"); + @>PS1_V4@ source->lensingPSF->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_PSF"); + @>PS1_V4@ source->lensingPSF->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_PSF"); + } @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrf = psMetadataLookupF32 (&status, row, "MOMENTS_R1"); @@ -508,18 +1161,18 @@ @>PS1_V3@ source->skySlope = psMetadataLookupF32 (&status, row, "SKY_LIMIT_SLOPE"); - @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_DV1@ source->diffStats->Rp = psMetadataLookupF32 (&status, row, "DIFF_R_P"); - @>PS1_DV1@ source->diffStats->SNp = psMetadataLookupF32 (&status, row, "DIFF_SN_P"); - @>PS1_DV1@ source->diffStats->Rm = psMetadataLookupF32 (&status, row, "DIFF_R_M"); - @>PS1_DV1@ source->diffStats->SNm = psMetadataLookupF32 (&status, row, "DIFF_SN_M"); - @PS1_DV?@ } + @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_DV1@ source->diffStats->Rp = psMetadataLookupF32 (&status, row, "DIFF_R_P"); + @>PS1_DV1@ source->diffStats->SNp = psMetadataLookupF32 (&status, row, "DIFF_SN_P"); + @>PS1_DV1@ source->diffStats->Rm = psMetadataLookupF32 (&status, row, "DIFF_R_M"); + @>PS1_DV1@ source->diffStats->SNm = psMetadataLookupF32 (&status, row, "DIFF_SN_M"); + @PS1_DV?@ } @ALL@ source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); @@ -533,4 +1186,10 @@ return sources; +} + +psArray *pmSourcesRead_CMF_@CMFMODE@ (psFits *fits, psMetadata *header) { +// psArray *array = pmSourcesRead_CMF_@CMFMODE@_Old (fits, header); + psArray *array = pmSourcesRead_CMF_@CMFMODE@_New (fits, header); + return array; } @@ -599,17 +1258,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, "max 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, "max 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++) { - // this is the source associated with this image + // 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; + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; // skip sources without measurements @@ -639,11 +1298,11 @@ psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); - float AxialRatio = NAN; - float AxialTheta = NAN; - pmSourceExtendedPars *extpars = source->extpars; - if (extpars) { - AxialRatio = extpars->axes.minor / extpars->axes.major; - AxialTheta = extpars->axes.theta; - } + 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); @@ -651,22 +1310,22 @@ // Petrosian measurements // XXX insert header data: petrosian ref radius, flux ratio - // XXX check flags to see if Pet was measured + // XXX check flags to see if Pet was measured if (doPetrosian) { - pmSourceExtendedPars *extpars = source->extpars; + 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 - // NOTE EAM 20140806 : PETRO_MAG_ERR was inverted!! this allows for it to be repaired - float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point - if (repairMagErrors) { - // I need to add the kron error in quadrature becasue pet_error ignores the object flux - float Krf = source->moments ? source->moments->KronFlux : NAN; - float dKrf = source->moments ? source->moments->KronFluxErr : NAN; - if (isfinite (Krf) && isfinite (dKrf)) { - magErr = sqrt(PS_SQR(1.0 / magErr) + PS_SQR(dKrf / Krf)); - } else { - magErr = 1.0 / magErr; - } - } + // 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 + // NOTE EAM 20140806 : PETRO_MAG_ERR was inverted!! this allows for it to be repaired + float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point + if (repairMagErrors) { + // I need to add the kron error in quadrature becasue pet_error ignores the object flux + float Krf = source->moments ? source->moments->KronFlux : NAN; + float dKrf = source->moments ? source->moments->KronFluxErr : NAN; + if (isfinite (Krf) && isfinite (dKrf)) { + magErr = sqrt(PS_SQR(1.0 / magErr) + PS_SQR(dKrf / Krf)); + } else { + magErr = 1.0 / magErr; + } + } 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); @@ -705,61 +1364,61 @@ // Flux Annuli (if we have extended source measurements, we have these. only optionally save them) if (doAnnuli) { - 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); + 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); @@ -885,25 +1544,25 @@ extpars->axes.theta = psMetadataLookupF32(&status, row, "F25_THETA"); - // magErr may have been saved in inverted form + // magErr may have been saved in inverted form 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)) { - if (repairMagErrors) { - // I need to add the kron error in quadrature becasue pet_error ignores the object flux - float Krf = source->moments ? source->moments->KronFlux : NAN; - float dKrf = source->moments ? source->moments->KronFluxErr : NAN; - if (isfinite (Krf) && isfinite (dKrf)) { - magErr = 1.0 / sqrt(PS_SQR(magErr) - PS_SQR(dKrf / Krf)); - } else { - magErr = 1.0 / magErr; - } - extpars->petrosianFluxErr = extpars->petrosianFlux * magErr; - } else { - extpars->petrosianFluxErr = extpars->petrosianFlux / magErr; - } - } - } + extpars->petrosianFlux = pow(10., (magOffset - mag) / 2.5); + if (isfinite(magErr)) { + if (repairMagErrors) { + // I need to add the kron error in quadrature becasue pet_error ignores the object flux + float Krf = source->moments ? source->moments->KronFlux : NAN; + float dKrf = source->moments ? source->moments->KronFluxErr : NAN; + if (isfinite (Krf) && isfinite (dKrf)) { + magErr = 1.0 / sqrt(PS_SQR(magErr) - PS_SQR(dKrf / Krf)); + } else { + magErr = 1.0 / magErr; + } + extpars->petrosianFluxErr = extpars->petrosianFlux * magErr; + } else { + extpars->petrosianFluxErr = extpars->petrosianFlux / magErr; + } + } + } extpars->petrosianRadius = psMetadataLookupF32(&status, row, "PETRO_RADIUS"); @@ -973,9 +1632,9 @@ int nParamMax = 0; for (int i = 0; i < sources->n; i++) { - // this is the source associated with this image + // 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; + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; if (source->modelFits == NULL) continue; @@ -1004,6 +1663,6 @@ pmSource *thisSource = sources->data[i]; - // this is the "real" version of this source - pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + // 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? @@ -1018,9 +1677,9 @@ // pmSourceExtFitPars *extPars = source->extFitPars->data[j]; - // assert (extPars); - - // skip models which were not actually fitted - // XXX - if (model->flags & badModel) continue; + // assert (extPars); + + // skip models which were not actually fitted + // XXX + if (model->flags & badModel) continue; PAR = model->params->data.F32; @@ -1029,37 +1688,37 @@ yPos = PAR[PM_PAR_YPOS]; - // for the extended source models, we do not always fit the centroid in the non-linear fitting process - // current situation (hard-wired into psphotSourceFits.c:psphotFitPCM, - // SERSIC, DEV, EXP : X,Y not fitted (PCM and not PCM) - // TRAIL : X,Y are fitted - // - - // XXX this should be based on what happened, not on the model type - if (model->type == modelTypeTrail) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - } else { - // this is definitely an underestimate since it does not - // account for the extent of the source - xErr = fwhmMajor * model->magErr / 2.35; - yErr = fwhmMinor * model->magErr / 2.35; - } - - @>PS1_DV2,>PS1_SV3@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; - @>PS1_DV2,>PS1_SV3@ float posAngle = 0.0; - @>PS1_DV2,>PS1_SV3@ float pltScale = 0.0; - @>PS1_DV2,>PS1_SV3@ pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); - @>PS1_DV2,>PS1_SV3@ double raPos = ptSky.r*PS_DEG_RAD; - @>PS1_DV2,>PS1_SV3@ double decPos = ptSky.d*PS_DEG_RAD; - @>PS1_DV2,>PS1_SV3@ posAngle *= PS_DEG_RAD; - @>PS1_DV2,>PS1_SV3@ pltScale *= PS_DEG_RAD*3600.0; - - float kronFlux = source->moments ? source->moments->KronFlux : NAN; - float kronMag = isfinite(kronFlux) ? -2.5*log10(kronFlux) : NAN; + // for the extended source models, we do not always fit the centroid in the non-linear fitting process + // current situation (hard-wired into psphotSourceFits.c:psphotFitPCM, + // SERSIC, DEV, EXP : X,Y not fitted (PCM and not PCM) + // TRAIL : X,Y are fitted + // + + // XXX this should be based on what happened, not on the model type + if (model->type == modelTypeTrail) { + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + } else { + // this is definitely an underestimate since it does not + // account for the extent of the source + xErr = fwhmMajor * model->magErr / 2.35; + yErr = fwhmMinor * model->magErr / 2.35; + } + + @>PS1_DV2,>PS1_SV3@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; + @>PS1_DV2,>PS1_SV3@ float posAngle = 0.0; + @>PS1_DV2,>PS1_SV3@ float pltScale = 0.0; + @>PS1_DV2,>PS1_SV3@ pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + @>PS1_DV2,>PS1_SV3@ double raPos = ptSky.r*PS_DEG_RAD; + @>PS1_DV2,>PS1_SV3@ double decPos = ptSky.d*PS_DEG_RAD; + @>PS1_DV2,>PS1_SV3@ posAngle *= PS_DEG_RAD; + @>PS1_DV2,>PS1_SV3@ pltScale *= PS_DEG_RAD*3600.0; + + float kronFlux = source->moments ? source->moments->KronFlux : NAN; + float kronMag = isfinite(kronFlux) ? -2.5*log10(kronFlux) : NAN; row = psMetadataAlloc (); - // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation - // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment" + // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation + // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment" 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); @@ -1070,21 +1729,21 @@ @>PS1_DV2,>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "RA_EXT", 0, "EXT model ra coordinate", raPos); @>PS1_DV2,>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "DEC_EXT", 0, "EXT model dec coordinate", decPos); - @>PS1_DV2@ float instFlux = isfinite(model->mag) ? pow(10.0, -0.4*model->mag) : NAN; - @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_FLUX", 0, "EXT fit instrumental counts", instFlux); + @>PS1_DV2@ float instFlux = isfinite(model->mag) ? pow(10.0, -0.4*model->mag) : NAN; + @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_FLUX", 0, "EXT fit instrumental counts", instFlux); 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); - @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN; - @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration", calMag); - @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ", PS_DATA_F32, "EXT Model Chisq", model->chisq); - @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF", PS_DATA_S32, "EXT Model num degrees of freedom", model->nDOF); - @>PS1_SV1,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_MODEL_TYPE", PS_DATA_S32, "type for chosen EXT_MODEL", source->modelEXT ? source->modelEXT->type : -1); - - // EAM : adding for PV2 outputs: - @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", model->flags); + @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN; + @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration", calMag); + @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ", PS_DATA_F32, "EXT Model Chisq", model->chisq); + @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF", PS_DATA_S32, "EXT Model num degrees of freedom", model->nDOF); + @>PS1_SV1,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_MODEL_TYPE", PS_DATA_S32, "type for chosen EXT_MODEL", source->modelEXT ? source->modelEXT->type : -1); + + // EAM : adding for PV2 outputs: + @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", model->flags); @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE", 0, "position angle at source (degrees)", posAngle); - @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "PLTSCALE", 0, "plate scale at source (arcsec/pixel)", pltScale); + @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "PLTSCALE", 0, "plate scale at source (arcsec/pixel)", pltScale); // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_XX", 0, "second moment in x", extPars->Mxx); @@ -1103,29 +1762,29 @@ // XXX these should be major and minor, not 'x' and 'y' - 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]); // this is not fitted - 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)", NAN); // this is not fitted, so error is NAN - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_THETA]); - } else { - if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY])) { - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (SXX, isnan)", PAR[PM_PAR_SXX]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (SYY, isnan)", PAR[PM_PAR_SYY]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT angle (SXY, isnan)", PAR[PM_PAR_SXY]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR", 0, "EXT width err (SXX, isnan)", dPAR[PM_PAR_SXX]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR", 0, "EXT width err (SYY, isnan)", dPAR[PM_PAR_SYY]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]); - } else { - psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->class->useReff); - 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_SXX]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)", dPAR[PM_PAR_SYY]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_SXY]); - } - } + 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]); // this is not fitted + 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)", NAN); // this is not fitted, so error is NAN + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_THETA]); + } else { + if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY])) { + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (SXX, isnan)", PAR[PM_PAR_SXX]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (SYY, isnan)", PAR[PM_PAR_SYY]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT angle (SXY, isnan)", PAR[PM_PAR_SXY]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR", 0, "EXT width err (SXX, isnan)", dPAR[PM_PAR_SXX]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR", 0, "EXT width err (SYY, isnan)", dPAR[PM_PAR_SYY]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]); + } else { + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->class->useReff); + 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_SXX]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)", dPAR[PM_PAR_SYY]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_SXY]); + } + } // write out the other generic parameters @@ -1140,5 +1799,5 @@ snprintf (name, 64, "EXT_PAR_%02d", k); - + if (k < model->params->n) { psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]); @@ -1148,15 +1807,15 @@ } - // optionally, write out the covariance matrix values - // XXX do I need to pad this to match the biggest covar matrix? - if (false && 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]); - - } - } - } + // optionally, write out the covariance matrix values + // XXX do I need to pad this to match the biggest covar matrix? + if (false && 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); @@ -1280,28 +1939,28 @@ PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07"); // XXX add an error: - // dPAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07_"); - } - - // NOTE: we no longer write out the covariance matrix - if (false) { - // 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; - } - - // we are only saving the values stored in dPAR[SXX,etc] + // dPAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07_"); + } + + // NOTE: we no longer write out the covariance matrix + if (false) { + // 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; + } + + // we are only saving the values stored in dPAR[SXX,etc] dPAR[PM_PAR_SXX] = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ_ERR"); dPAR[PM_PAR_SYY] = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN_ERR"); dPAR[PM_PAR_SXY] = psMetadataLookupF32(&status, row, "EXT_THETA_ERR"); - // other parameters that we need to read + // other parameters that we need to read PAR[PM_PAR_SKY] = psMetadataLookupF32(&status, row, "SKY_EXT"); @@ -1341,6 +2000,6 @@ // 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; + psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n"); + return true; } @@ -1385,77 +2044,77 @@ for (int i = 0; i < sources->n; i++) { - // this is the source associated with this image + // 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; + // 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; + 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) - // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation - // This set of psMetadataAdd Entries marks the "----" "Start of the XRAD segment" - 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: - psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX", PS_META_REPLACE, "flux within annuli", radFlux); - psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_META_REPLACE, "flux error in annuli", radFluxErr); - psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_META_REPLACE, "flux standard deviation", radFluxStdev); - psMetadataAddVector (row, PS_LIST_TAIL, "APER_FILL", PS_META_REPLACE, "fill factor of annuli", radFill); - psFree (radFlux); - psFree (radFluxErr); - psFree (radFluxStdev); - psFree (radFill); - - psArrayAdd (table, 100, row); - psFree (row); - } + 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) + // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation + // This set of psMetadataAdd Entries marks the "----" "Start of the XRAD segment" + 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: + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX", PS_META_REPLACE, "flux within annuli", radFlux); + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_META_REPLACE, "flux error in annuli", radFluxErr); + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_META_REPLACE, "flux standard deviation", radFluxStdev); + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FILL", PS_META_REPLACE, "fill factor of annuli", radFill); + psFree (radFlux); + psFree (radFluxErr); + psFree (radFluxStdev); + psFree (radFill); + + psArrayAdd (table, 100, row); + psFree (row); + } } @@ -1581,6 +2240,6 @@ // perform full non-linear fits / extended source analysis? if (!psMetadataLookupBool (&status, recipe, "GALAXY_SHAPES")) { - psLogMsg ("psphot", PS_LOG_INFO, "galaxy shapes were not measured, skipping\n"); - return true; + psLogMsg ("psphot", PS_LOG_INFO, "galaxy shapes were not measured, skipping\n"); + return true; } @@ -1609,11 +2268,11 @@ pmSource *thisSource = sources->data[i]; - // this is the "real" version of this source - pmSource *source = thisSource->parent ? thisSource->parent : thisSource; - - // if we did not fit the galaxy model, modelFits will be NULL + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + + // if we did not fit the galaxy model, modelFits will be NULL if (source->modelFits == NULL) continue; - // if we did not fit the galaxy model, galaxyFits will also be NULL + // if we did not fit the galaxy model, galaxyFits will also be NULL if (source->galaxyFits == NULL) continue;