Changeset 41821 for branches/eam_branches/ipp-dev-20210817
- Timestamp:
- Sep 28, 2021, 5:06:06 PM (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branches/ipp-dev-20210817/psModules/src/objects/pmSourceIO_CMF.c.in
r38114 r41821 57 57 // followed by a zero-size matrix, followed by the table data 58 58 59 bool pmSourcesWrite_CMF_@CMFMODE@ (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) 59 /** this file is used to generate c-code for a number of different formats. the perl program 60 ** 'mksource.pl' reads this template file and generates an output file for one of the named 61 ** formats, e.g., PS1_V1 or PS1_DV4 (see list in mksource.pl). 62 63 ** Any instances of the word @CMFMODE@ are replaced by the named version of the format. 64 65 ** Any line which starts with a command of the form: @COMMAND[,COMMAND...]@ is kept or removed 66 ** from the output c-code depending on the commands. These may be: ALL (lines are kept for 67 ** all formats) or the name of a format. A line with just a plain format will be kept only 68 ** for that format. The format name may be optionally proceded by >, <, <=, >=, ! or =. In 69 ** these cases, the line is kept if the format is greater / lesser / not the same, etc as the 70 ** listed format. The sequence of formats is defined for series (V, SV, DV) in mksource.pl so 71 ** that, e.g., >PS1_V3 would include e.g., PS1_V4 and PS1_V5, but not PS1_DV4. The number in 72 ** a series may be replaced with ? to stand for the full set. The format 73 ** commands may be grouped on a line separated by commas. Note that only ! (not) is 74 ** restrictive: only this rule will reduce the number of matches. So for example, 75 ** @ALL,<PS1_DV3@ will have the effect of @ALL@ since any format >= PS1_DV3 will be matched by 76 ** ALL. But @ALL,!PS1_DV3@ will generate a line for any format except PS1_DV3. 77 78 **/ 79 80 bool pmSourcesWrite_CMF_@CMFMODE@_New (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) 60 81 { 61 82 PS_ASSERT_PTR_NON_NULL(fits, false); 62 83 PS_ASSERT_PTR_NON_NULL(sources, false); 63 84 PS_ASSERT_PTR_NON_NULL(extname, false); 64 65 psArray *table;66 psMetadata *row;67 85 68 86 pmChip *chip = readout->parent->parent; … … 73 91 pmSource *source = (pmSource *) sources->data[0]; 74 92 if (source->seq == -1) { 75 76 93 // let's write these out in S/N order 94 sources = psArraySort (sources, pmSourceSortByFlux); 77 95 } else { 78 sources = psArraySort (sources, pmSourceSortBySeq); 79 } 80 } 81 82 table = psArrayAllocEmpty (sources->n); 96 sources = psArraySort (sources, pmSourceSortBySeq); 97 } 98 } 83 99 84 100 float magOffset; … … 88 104 pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); 89 105 106 // before we can generate the table structure, we need to know if certain fields were measured. 107 // this information is not available a priori, so we need to check all sources 108 109 bool haveLensingOBJsmear = false; 110 bool haveLensingOBJshear = false; 111 bool haveLensingPSFsmear = false; 112 bool haveLensingPSFshear = false; 113 bool haveLensingPSF = false; 114 115 for (int i = 0; (i < sources->n) && !haveLensingOBJsmear && !haveLensingOBJshear && !haveLensingPSFsmear && !haveLensingPSFshear && !haveLensingPSF; i++) { 116 // this is the source associated with this image 117 pmSource *thisSource = sources->data[i]; 118 119 // this is the "real" version of this source 120 pmSource *source = thisSource->parent ? thisSource->parent : thisSource; 121 122 if (source->lensingOBJ && source->lensingOBJ->smear) haveLensingOBJsmear = true; 123 if (source->lensingOBJ && source->lensingOBJ->shear) haveLensingOBJshear = true; 124 if (source->lensingPSF && source->lensingPSF->smear) haveLensingPSFsmear = true; 125 if (source->lensingPSF && source->lensingPSF->shear) haveLensingPSFshear = true; 126 if (source->lensingPSF ) haveLensingPSF = true; 127 } 128 129 /************ generate the table columns *****************/ 130 131 // before we allocate the table, generate an empty array of table columns and generate them 132 psArray *tableColumns = psArrayAllocEmpty (100); 133 134 // add the named / typed columns to the collection of columns 135 @ALL@ psFitsTableColumnAdd (tableColumns, "IPP_IDET", PS_DATA_U32); // "IPP detection identifier index" 136 @ALL@ psFitsTableColumnAdd (tableColumns, "X_PSF", PS_DATA_F32); // "PSF x coordinate" 137 @ALL@ psFitsTableColumnAdd (tableColumns, "Y_PSF", PS_DATA_F32); // "PSF y coordinate" 138 @ALL@ psFitsTableColumnAdd (tableColumns, "X_PSF_SIG", PS_DATA_F32); // "Sigma in PSF x coordinate" 139 @ALL@ psFitsTableColumnAdd (tableColumns, "Y_PSF_SIG", PS_DATA_F32); // "Sigma in PSF y coordinate" 140 141 // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision 142 @PS1_V1@ psFitsTableColumnAdd (tableColumns, "RA_PSF", PS_DATA_F32); // "PSF RA coordinate (degrees)" 143 @PS1_V1@ psFitsTableColumnAdd (tableColumns, "DEC_PSF", PS_DATA_F32); // "PSF DEC coordinate (degrees)" 144 145 @ALL@ psFitsTableColumnAdd (tableColumns, "POSANGLE", PS_DATA_F32); // "position angle at source (degrees)" 146 @ALL@ psFitsTableColumnAdd (tableColumns, "PLTSCALE", PS_DATA_F32); // "plate scale at source (arcsec/pixel)" 147 @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_INST_MAG", PS_DATA_F32); // "PSF fit instrumental magnitude" 148 @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_INST_MAG_SIG", PS_DATA_F32); // "Sigma of PSF instrumental magnitude" 149 150 @ALL,!PS1_V1,!PS1_V2@ psFitsTableColumnAdd (tableColumns, "PSF_INST_FLUX", PS_DATA_F32); // "PSF fit instrumental flux (counts)" 151 @ALL,!PS1_V1,!PS1_V2@ psFitsTableColumnAdd (tableColumns, "PSF_INST_FLUX_SIG", PS_DATA_F32); // "Sigma of PSF instrumental flux" 152 153 @ALL@ psFitsTableColumnAdd (tableColumns, "AP_MAG", PS_DATA_F32); // "magnitude in standard aperture" 154 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "AP_MAG_RAW", PS_DATA_F32); // "magnitude in reported aperture" 155 @ALL@ psFitsTableColumnAdd (tableColumns, "AP_MAG_RADIUS", PS_DATA_F32); // "radius used for aperture mags" 156 @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableColumnAdd (tableColumns, "AP_FLUX", PS_DATA_F32); // "instrumental flux in standard aperture" 157 @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableColumnAdd (tableColumns, "AP_FLUX_SIG", PS_DATA_F32); // "aperture flux error" 158 @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "AP_NPIX", PS_DATA_S32); // "aperture unmasked pixels" 159 160 @<PS1_V3,PS1_SV1,PS1_DV?@ psFitsTableColumnAdd (tableColumns, "PEAK_FLUX_AS_MAG", PS_DATA_F32); // "Peak flux expressed as magnitude" 161 162 @ALL@ psFitsTableColumnAdd (tableColumns, "CAL_PSF_MAG", PS_DATA_F32); // "PSF Magnitude using supplied calibration" 163 @ALL@ psFitsTableColumnAdd (tableColumns, "CAL_PSF_MAG_SIG", PS_DATA_F32); // "measured scatter of zero point calibration" 164 165 // NOTE: RA & DEC (both double) need to be on an 8-byte boundary... 166 @ALL,!PS1_V1@ psFitsTableColumnAdd (tableColumns, "RA_PSF", PS_DATA_F64); // "PSF RA coordinate (degrees)" 167 @ALL,!PS1_V1@ psFitsTableColumnAdd (tableColumns, "DEC_PSF", PS_DATA_F64); // "PSF DEC coordinate (degrees)" 168 169 @>=PS1_V3,>PS1_SV1@ psFitsTableColumnAdd (tableColumns, "PEAK_FLUX_AS_MAG", PS_DATA_F32); // "Peak flux expressed as magnitude" 170 @ALL@ psFitsTableColumnAdd (tableColumns, "SKY", PS_DATA_F32); // "Sky level" 171 @ALL@ psFitsTableColumnAdd (tableColumns, "SKY_SIGMA", PS_DATA_F32); // "Sigma of sky level" 172 @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_CHISQ", PS_DATA_F32); // "Chisq of PSF-fit" 173 @ALL@ psFitsTableColumnAdd (tableColumns, "CR_NSIGMA", PS_DATA_F32); // "Nsigma deviations from PSF to CF" 174 @ALL@ psFitsTableColumnAdd (tableColumns, "EXT_NSIGMA", PS_DATA_F32); // "Nsigma deviations from PSF to EXT" 175 176 // PSF shape parameters: 177 @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_MAJOR", PS_DATA_F32); // "PSF width (major axis)" 178 @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_MINOR", PS_DATA_F32); // "PSF width (minor axis)" 179 @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_THETA", PS_DATA_F32); // "PSF orientation angle" 180 @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "PSF_CORE", PS_DATA_F32); // "k term if defined" 181 @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "PSF_FWHM_MAJ", PS_DATA_F32); // "PSF FWHM (major axis)" 182 @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "PSF_FWHM_MIN", PS_DATA_F32); // "PSF FWHM (minor axis)" 183 184 // psf data quality 185 @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_QF", PS_DATA_F32); // "PSF coverage/quality factor (bad)" 186 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "PSF_QF_PERFECT", PS_DATA_F32); // "PSF coverage/quality factor (poor)" 187 @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_NDOF", PS_DATA_S32); // "degrees of freedom" 188 @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_NPIX", PS_DATA_S32); // "number of pixels in fit" 189 190 @ALL@ psFitsTableColumnAdd (tableColumns, "MOMENTS_XX", PS_DATA_F32); // "second moments (X^2)" 191 @ALL@ psFitsTableColumnAdd (tableColumns, "MOMENTS_XY", PS_DATA_F32); // "second moments (X*Y)" 192 @ALL@ psFitsTableColumnAdd (tableColumns, "MOMENTS_YY", PS_DATA_F32); // "second moments (Y*Y)" 193 @>PS1_V2,PS1_SV?@ psFitsTableColumnAdd (tableColumns, "MOMENTS_M3C", PS_DATA_F32); // "third momemt cos theta" 194 @>PS1_V2,PS1_SV?@ psFitsTableColumnAdd (tableColumns, "MOMENTS_M3S", PS_DATA_F32); // "third momemt sin theta" 195 @>PS1_V2,PS1_SV?@ psFitsTableColumnAdd (tableColumns, "MOMENTS_M4C", PS_DATA_F32); // "fourth momemt cos theta" 196 @>PS1_V2,PS1_SV?@ psFitsTableColumnAdd (tableColumns, "MOMENTS_M4S", PS_DATA_F32); // "fourth momemt sin theta" 197 198 // Lensing parameters are only written if they are measured 199 if (haveLensingOBJsmear) { 200 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X11_SM_OBJ", PS_DATA_F32); // "smear polarizability element (objects)" 201 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X12_SM_OBJ", PS_DATA_F32); // "smear polarizability element (objects)" 202 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X22_SM_OBJ", PS_DATA_F32); // "smear polarizability element (objects)" 203 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_SM_OBJ", PS_DATA_F32); // "smear polarizability element (objects)" 204 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_SM_OBJ", PS_DATA_F32); // "smear polarizability element (objects)" 205 } 206 if (haveLensingOBJshear) { 207 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X11_SH_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" 208 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X12_SH_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" 209 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X22_SH_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" 210 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_SH_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" 211 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_SH_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" 212 } 213 if (false) { 214 // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above 215 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" 216 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" 217 } 218 219 if (haveLensingPSFsmear) { 220 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X11_SM_PSF", PS_DATA_F32); // "smear polarizability element (PSFs)" 221 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X12_SM_PSF", PS_DATA_F32); // "smear polarizability element (PSFs)" 222 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X22_SM_PSF", PS_DATA_F32); // "smear polarizability element (PSFs)" 223 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_SM_PSF", PS_DATA_F32); // "smear polarizability element (PSFs)" 224 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_SM_PSF", PS_DATA_F32); // "smear polarizability element (PSFs)" 225 } 226 if (haveLensingPSFshear) { 227 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X11_SH_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" 228 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X12_SH_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" 229 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X22_SH_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" 230 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_SH_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" 231 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_SH_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" 232 } 233 if (haveLensingPSF) { 234 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" 235 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" 236 } 237 238 // if lensing params exist also include the backmapped chipID and chip coordinates 239 if (haveLensingPSF) { 240 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_NUM", PS_DATA_S16); // "id of warp input chip" 241 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_X", PS_DATA_S16); // "x coord in warp input chip" 242 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_Y", PS_DATA_S16); // "y coord in warp input chip" 243 @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "PADDING3", PS_DATA_S16); // "more padding" 244 } 245 246 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "MOMENTS_R1", PS_DATA_F32); // "first radial moment" 247 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "MOMENTS_RH", PS_DATA_F32); // "half radial moment" 248 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "KRON_FLUX", PS_DATA_F32); // "Kron Flux (in 2.5 R1)" 249 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "KRON_FLUX_ERR", PS_DATA_F32); // "Kron Flux Error" 250 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "KRON_FLUX_INNER", PS_DATA_F32); // "Kron Flux (in 1.0 R1)" 251 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "KRON_FLUX_OUTER", PS_DATA_F32); // "Kron Flux (in 4.0 R1)" 252 253 @>PS1_V3@ psFitsTableColumnAdd (tableColumns, "SKY_LIMIT_RAD", PS_DATA_F32); // "Radius where object hits sky" 254 @>PS1_V3@ psFitsTableColumnAdd (tableColumns, "SKY_LIMIT_FLUX", PS_DATA_F32); // "Flux / pix where object hits sky" 255 @>PS1_V3@ psFitsTableColumnAdd (tableColumns, "SKY_LIMIT_SLOPE", PS_DATA_F32); // "d(Flux/pix)/dRadius where object hits sky" 256 257 @>PS1_DV4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_NUM", PS_DATA_S16); // "id of warp input chip" 258 @>PS1_DV4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_X", PS_DATA_S16); // "x coord in warp input chip" 259 @>PS1_DV4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_Y", PS_DATA_S16); // "y coord in warp input chip" 260 @>PS1_DV4@ psFitsTableColumnAdd (tableColumns, "PADDING3", PS_DATA_S16); // "more padding" 261 262 @PS1_DV?@ psFitsTableColumnAdd (tableColumns, "DIFF_NPOS", PS_DATA_S32); // "nPos (n pix > 3 sigma)" 263 @PS1_DV?@ psFitsTableColumnAdd (tableColumns, "DIFF_FRATIO", PS_DATA_F32); // "fPos / (fPos + fNeg)" 264 @PS1_DV?@ psFitsTableColumnAdd (tableColumns, "DIFF_NRATIO_BAD", PS_DATA_F32); // "nPos / (nPos + nNeg)" 265 @PS1_DV?@ psFitsTableColumnAdd (tableColumns, "DIFF_NRATIO_MASK", PS_DATA_F32); // "nPos / (nPos + nMask)" 266 @PS1_DV?@ psFitsTableColumnAdd (tableColumns, "DIFF_NRATIO_ALL", PS_DATA_F32); // "nPos / (nGood + nMask + nBad)" 267 268 @>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "DIFF_R_P", PS_DATA_F32); // "distance to positive match source" 269 @>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "DIFF_SN_P", PS_DATA_F32); // "signal-to-noise of pos match src" 270 @>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "DIFF_R_M", PS_DATA_F32); // "distance to negative match source" 271 @>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "DIFF_SN_M", PS_DATA_F32); // "signal-to-noise of neg match src" 272 273 @ALL@ psFitsTableColumnAdd (tableColumns, "FLAGS", PS_DATA_U32); // "psphot analysis flags" 274 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "FLAGS2", PS_DATA_U32); // "psphot analysis flags" 275 @PS1_V3,PS1_SV2@ psFitsTableColumnAdd (tableColumns, "PADDING2", PS_DATA_S32); // "more padding" 276 @ALL@ psFitsTableColumnAdd (tableColumns, "N_FRAMES", PS_DATA_U16); // "Number of frames overlapping source center" 277 278 @ALL@ psFitsTableColumnAdd (tableColumns, "PADDING", PS_DATA_S16); // "more padding" 279 280 // note that this definition is inconsistent with the definition in 281 // Ohana/src/libautocode/def. This version creates a table with data not 282 // properly aligned with the 8-byte boundaries. The structure defined by 283 // libautocode does this, but has a different order of elements (and adds 284 // padding2 to fix things). We have generated may files with PS1_SV1 as is, so 285 // I'll leave it. But in future a PS1_SV2 should be forced to match 286 // libautocode. Note that addstar knows to detect the alternate version of 287 // PS1_SV1 and correctly interpret its fields. 288 289 // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not 290 // subtracted 291 292 // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image 293 // edge; 3) any pixels in the 3x3 peak region are masked; 294 295 // generate an FITS table using the array of columns defined above, with space for all rows 296 psFitsTable *table = psFitsTableCreate (tableColumns, sources->n); 297 298 /************ write the data to the table *****************/ 299 90 300 // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state 91 301 // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. 92 302 for (int i = 0; i < sources->n; i++) { 93 303 // this is the source associated with this image 94 304 pmSource *thisSource = sources->data[i]; 95 305 96 97 306 // this is the "real" version of this source 307 pmSource *source = thisSource->parent ? thisSource->parent : thisSource; 98 308 99 309 // If source->seq is -1, source was generated in this analysis. If source->seq is … … 102 312 // generated on Alloc, and would thus be wrong for read in sources. 103 313 if (source->seq == -1) { 104 source->seq = i; 105 } 106 107 // set the 'best' values for various output fields: 108 pmSourceOutputs outputs; 109 pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); 110 111 pmSourceOutputsMoments moments; 112 pmSourceOutputsSetMoments (&moments, source); 113 114 @PS1_DV?@ pmSourceDiffStats diffStats; 115 @PS1_DV?@ pmSourceDiffStatsInit(&diffStats); 116 @PS1_DV?@ if (source->diffStats) { 117 @PS1_DV?@ diffStats = *source->diffStats; 118 @PS1_DV?@ } 119 120 row = psMetadataAlloc (); 314 source->seq = i; 315 } 316 317 // set the 'best' values for various output fields: 318 pmSourceOutputs outputs; 319 pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); 320 321 pmSourceOutputsMoments moments; 322 pmSourceOutputsSetMoments (&moments, source); 323 324 @PS1_DV?@ pmSourceDiffStats diffStats; 325 @PS1_DV?@ pmSourceDiffStatsInit(&diffStats); 326 @PS1_DV?@ if (source->diffStats) { diffStats = *source->diffStats;} 121 327 122 328 // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation 123 329 // This set of psMetadataAdd Entries marks the "----" "Start of the PSF segment" 124 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq);125 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos);126 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos);127 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr);128 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr);129 130 // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision 131 @PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F32, "PSF RA coordinate (degrees)", outputs.ra); 132 @PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F32, "PSF DEC coordinate (degrees)", outputs.dec); 133 134 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle);135 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale);136 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag);137 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr);138 139 @ALL,!PS1_V1,!PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); 140 @ALL,!PS1_V1,!PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux", source->psfFluxErr); 141 142 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag);143 @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in reported aperture", source->apMagRaw); 144 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", source->apRadius);145 @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", source->apFlux); 146 @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", source->apFluxErr); 147 @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "AP_NPIX", PS_DATA_S32, "aperture unmasked pixels", source->apNpixels); 148 149 @<PS1_V3,PS1_SV1,PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); 150 151 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag);152 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);153 154 // NOTE: RA & DEC (both double) need to be on an 8-byte boundary... 155 @ALL,!PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); 156 @ALL,!PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); 157 158 @>=PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); 159 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky);160 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr);161 162 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq);163 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma);164 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma);330 @ALL@ psFitsTableSetU32 (table, i, "IPP_IDET", source->seq); // "IPP detection identifier index", 331 @ALL@ psFitsTableSetF32 (table, i, "X_PSF", outputs.xPos); // "PSF x coordinate", 332 @ALL@ psFitsTableSetF32 (table, i, "Y_PSF", outputs.yPos); // "PSF y coordinate", 333 @ALL@ psFitsTableSetF32 (table, i, "X_PSF_SIG", outputs.xErr); // "Sigma in PSF x coordinate", 334 @ALL@ psFitsTableSetF32 (table, i, "Y_PSF_SIG", outputs.yErr); // "Sigma in PSF y coordinate", 335 336 // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision 337 @PS1_V1@ psFitsTableSetF32 (table, i, "RA_PSF", outputs.ra); // "PSF RA coordinate (degrees)", 338 @PS1_V1@ psFitsTableSetF32 (table, i, "DEC_PSF", outputs.dec); // "PSF DEC coordinate (degrees)", 339 340 @ALL@ psFitsTableSetF32 (table, i, "POSANGLE", outputs.posAngle); // "position angle at source (degrees)", 341 @ALL@ psFitsTableSetF32 (table, i, "PLTSCALE", outputs.pltScale); // "plate scale at source (arcsec/pixel)", 342 @ALL@ psFitsTableSetF32 (table, i, "PSF_INST_MAG", source->psfMag); // "PSF fit instrumental magnitude", 343 @ALL@ psFitsTableSetF32 (table, i, "PSF_INST_MAG_SIG", source->psfMagErr); // "Sigma of PSF instrumental magnitude", 344 345 @ALL,!PS1_V1,!PS1_V2@ psFitsTableSetF32 (table, i, "PSF_INST_FLUX", source->psfFlux); // "PSF fit instrumental flux (counts)", 346 @ALL,!PS1_V1,!PS1_V2@ psFitsTableSetF32 (table, i, "PSF_INST_FLUX_SIG", source->psfFluxErr); // "Sigma of PSF instrumental flux", 347 348 @ALL@ psFitsTableSetF32 (table, i, "AP_MAG", source->apMag); // "magnitude in standard aperture", 349 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "AP_MAG_RAW", source->apMagRaw); // "magnitude in reported aperture", 350 @ALL@ psFitsTableSetF32 (table, i, "AP_MAG_RADIUS", source->apRadius); // "radius used for aperture mags", 351 @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableSetF32 (table, i, "AP_FLUX", source->apFlux); // "instrumental flux in standard aperture", 352 @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableSetF32 (table, i, "AP_FLUX_SIG", source->apFluxErr); // "aperture flux error", 353 @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetS32 (table, i, "AP_NPIX", source->apNpixels); // "aperture unmasked pixels", 354 355 @<PS1_V3,PS1_SV1,PS1_DV?@ psFitsTableSetF32 (table, i, "PEAK_FLUX_AS_MAG", outputs.peakMag); // "Peak flux expressed as magnitude" 356 357 @ALL@ psFitsTableSetF32 (table, i, "CAL_PSF_MAG", outputs.calMag); // "PSF Magnitude using supplied calibration", 358 @ALL@ psFitsTableSetF32 (table, i, "CAL_PSF_MAG_SIG", zeroptErr); // "measured scatter of zero point calibration", 359 360 // NOTE: RA & DEC (both double) need to be on an 8-byte boundary... 361 @ALL,!PS1_V1@ psFitsTableSetF64 (table, i, "RA_PSF", outputs.ra); // "PSF RA coordinate (degrees)", 362 @ALL,!PS1_V1@ psFitsTableSetF64 (table, i, "DEC_PSF", outputs.dec); // "PSF DEC coordinate (degrees)", 363 364 @>=PS1_V3,>PS1_SV1@ psFitsTableSetF32 (table, i, "PEAK_FLUX_AS_MAG", outputs.peakMag); // "Peak flux expressed as magnitude", 365 @ALL@ psFitsTableSetF32 (table, i, "SKY", source->sky); // "Sky level", 366 @ALL@ psFitsTableSetF32 (table, i, "SKY_SIGMA", source->skyErr); // "Sigma of sky level", 367 368 @ALL@ psFitsTableSetF32 (table, i, "PSF_CHISQ", outputs.chisq); // "Chisq of PSF-fit", 369 @ALL@ psFitsTableSetF32 (table, i, "CR_NSIGMA", source->crNsigma); // "Nsigma deviations from PSF to CF", 370 @ALL@ psFitsTableSetF32 (table, i, "EXT_NSIGMA", source->extNsigma); // "Nsigma deviations from PSF to EXT", 165 371 166 372 // PSF shape parameters: 167 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); 168 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); 169 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); 373 @ALL@ psFitsTableSetF32 (table, i, "PSF_MAJOR", outputs.psfMajor); // "PSF width (major axis)", 374 @ALL@ psFitsTableSetF32 (table, i, "PSF_MINOR", outputs.psfMinor); // "PSF width (minor axis)", 375 @ALL@ psFitsTableSetF32 (table, i, "PSF_THETA", outputs.psfTheta); // "PSF orientation angle", 376 @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetF32 (table, i, "PSF_CORE", outputs.psfCore); // "k term if defined", 377 378 // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values 379 @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetF32 (table, i, "PSF_FWHM_MAJ", outputs.psfMajorFWHM); // "PSF FWHM (major axis)", 380 @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetF32 (table, i, "PSF_FWHM_MIN", outputs.psfMinorFWHM); // "PSF FWHM (minor axis)", 381 382 // psf data quality 383 @ALL@ psFitsTableSetF32 (table, i, "PSF_QF", source->pixWeightNotBad); // "PSF coverage/quality factor (bad)", 384 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "PSF_QF_PERFECT", source->pixWeightNotPoor); // "PSF coverage/quality factor (poor)", 385 @ALL@ psFitsTableSetS32 (table, i, "PSF_NDOF", outputs.nDOF); // "degrees of freedom", 386 @ALL@ psFitsTableSetS32 (table, i, "PSF_NPIX", outputs.nPix); // "number of pixels in fit", 387 388 @ALL@ psFitsTableSetF32 (table, i, "MOMENTS_XX", moments.Mxx); // "second moments (X^2)", 389 @ALL@ psFitsTableSetF32 (table, i, "MOMENTS_XY", moments.Mxy); // "second moments (X*Y)", 390 @ALL@ psFitsTableSetF32 (table, i, "MOMENTS_YY", moments.Myy); // "second moments (Y*Y)", 391 392 @>PS1_V2,PS1_SV?@ psFitsTableSetF32 (table, i, "MOMENTS_M3C", moments.M_c3); // "third momemt cos theta", 393 @>PS1_V2,PS1_SV?@ psFitsTableSetF32 (table, i, "MOMENTS_M3S", moments.M_s3); // "third momemt sin theta", 394 @>PS1_V2,PS1_SV?@ psFitsTableSetF32 (table, i, "MOMENTS_M4C", moments.M_c4); // "fourth momemt cos theta", 395 @>PS1_V2,PS1_SV?@ psFitsTableSetF32 (table, i, "MOMENTS_M4S", moments.M_s4); // "fourth momemt sin theta", 396 397 // Lensing parameters: 398 if (source->lensingOBJ && source->lensingOBJ->smear) { 399 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X11_SM_OBJ", source->lensingOBJ->smear->X11); // "smear polarizability element (objects)", 400 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X12_SM_OBJ", source->lensingOBJ->smear->X12); // "smear polarizability element (objects)", 401 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X22_SM_OBJ", source->lensingOBJ->smear->X22); // "smear polarizability element (objects)", 402 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_SM_OBJ", source->lensingOBJ->smear->e1); // "smear polarizability element (objects)", 403 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_SM_OBJ", source->lensingOBJ->smear->e2); // "smear polarizability element (objects)", 404 } 405 if (source->lensingOBJ && source->lensingOBJ->shear) { 406 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X11_SH_OBJ", source->lensingOBJ->shear->X11); // "shear polarizability element (objects)", 407 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X12_SH_OBJ", source->lensingOBJ->shear->X12); // "shear polarizability element (objects)", 408 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X22_SH_OBJ", source->lensingOBJ->shear->X22); // "shear polarizability element (objects)", 409 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_SH_OBJ", source->lensingOBJ->shear->e1); // "shear polarizability element (objects)", 410 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_SH_OBJ", source->lensingOBJ->shear->e2); // "shear polarizability element (objects)", 411 } 412 if (false && source->lensingOBJ) { 413 // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above 414 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_PSF", source->lensingOBJ->e1); // "shear polarizability element (objects)", 415 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_PSF", source->lensingOBJ->e2); // "shear polarizability element (objects)", 416 } 417 418 if (source->lensingPSF && source->lensingPSF->smear) { 419 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X11_SM_PSF", source->lensingPSF->smear->X11); // "smear polarizability element (objects)", 420 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X12_SM_PSF", source->lensingPSF->smear->X12); // "smear polarizability element (objects)", 421 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X22_SM_PSF", source->lensingPSF->smear->X22); // "smear polarizability element (objects)", 422 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_SM_PSF", source->lensingPSF->smear->e1); // "smear polarizability element (objects)", 423 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_SM_PSF", source->lensingPSF->smear->e2); // "smear polarizability element (objects)", 424 } 425 if (source->lensingPSF && source->lensingPSF->shear) { 426 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X11_SH_PSF", source->lensingPSF->shear->X11); // "shear polarizability element (objects)", 427 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X12_SH_PSF", source->lensingPSF->shear->X12); // "shear polarizability element (objects)", 428 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X22_SH_PSF", source->lensingPSF->shear->X22); // "shear polarizability element (objects)", 429 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_SH_PSF", source->lensingPSF->shear->e1); // "shear polarizability element (objects)", 430 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_SH_PSF", source->lensingPSF->shear->e2); // "shear polarizability element (objects)", 431 } 432 if (source->lensingPSF) { 433 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_PSF", source->lensingPSF->e1); // "shear polarizability element (objects)", 434 @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_PSF", source->lensingPSF->e2); // "shear polarizability element (objects)", 435 } 436 437 // if lensing params exist also include the backmapped chipID and chip coordinates 438 if (source->lensingPSF && source->lensingPSF->shear) { 439 @>PS1_V4@ psFitsTableSetS16 (table, i, "SRC_CHIP_NUM", source->chipNum); // "id of warp input chip", 440 @>PS1_V4@ psFitsTableSetS16 (table, i, "SRC_CHIP_X", source->chipX); // "x coord in warp input chip", 441 @>PS1_V4@ psFitsTableSetS16 (table, i, "SRC_CHIP_Y", source->chipY); // "y coord in warp input chip", 442 @>PS1_V4@ psFitsTableSetS16 (table, i, "PADDING3", 0); // "more padding", 443 } 444 445 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "MOMENTS_R1", moments.Mrf); // "first radial moment", 446 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "MOMENTS_RH", moments.Mrh); // "half radial moment", 447 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "KRON_FLUX", moments.Krf); // "Kron Flux (in 2.5 R1)", 448 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "KRON_FLUX_ERR", moments.dKrf); // "Kron Flux Error", 449 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "KRON_FLUX_INNER", moments.Kinner); // "Kron Flux (in 1.0 R1)", 450 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "KRON_FLUX_OUTER", moments.Kouter); // "Kron Flux (in 4.0 R1)", 451 @>PS1_V3@ psFitsTableSetF32 (table, i, "SKY_LIMIT_RAD", source->skyRadius); // "Radius where object hits sky", 452 @>PS1_V3@ psFitsTableSetF32 (table, i, "SKY_LIMIT_FLUX", source->skyFlux); // "Flux / pix where object hits sky", 453 @>PS1_V3@ psFitsTableSetF32 (table, i, "SKY_LIMIT_SLOPE", source->skySlope); // "d(Flux/pix)/dRadius where object hits sky", 454 455 @>PS1_DV4@ psFitsTableSetS16 (table, i, "SRC_CHIP_NUM", source->chipNum); // "id of warp input chip", 456 @>PS1_DV4@ psFitsTableSetS16 (table, i, "SRC_CHIP_X", source->chipX); // "x coord in warp input chip", 457 @>PS1_DV4@ psFitsTableSetS16 (table, i, "SRC_CHIP_Y", source->chipY); // "y coord in warp input chip", 458 @>PS1_DV4@ psFitsTableSetS16 (table, i, "PADDING3", 0); // "more padding", 459 460 @PS1_DV?@ psFitsTableSetS32 (table, i, "DIFF_NPOS", diffStats.nGood); // "nPos (n pix > 3 sigma)", 461 @PS1_DV?@ psFitsTableSetF32 (table, i, "DIFF_FRATIO", diffStats.fRatio); // "fPos / (fPos + fNeg)", 462 @PS1_DV?@ psFitsTableSetF32 (table, i, "DIFF_NRATIO_BAD", diffStats.nRatioBad); // "nPos / (nPos + nNeg)", 463 @PS1_DV?@ psFitsTableSetF32 (table, i, "DIFF_NRATIO_MASK", diffStats.nRatioMask); // "nPos / (nPos + nMask)", 464 @PS1_DV?@ psFitsTableSetF32 (table, i, "DIFF_NRATIO_ALL", diffStats.nRatioAll); // "nPos / (nGood + nMask + nBad)", 465 466 @>PS1_DV1@ psFitsTableSetF32 (table, i, "DIFF_R_P", diffStats.Rp); // "distance to positive match source", 467 @>PS1_DV1@ psFitsTableSetF32 (table, i, "DIFF_SN_P", diffStats.SNp); // "signal-to-noise of pos match src", 468 @>PS1_DV1@ psFitsTableSetF32 (table, i, "DIFF_R_M", diffStats.Rm); // "distance to negative match source", 469 @>PS1_DV1@ psFitsTableSetF32 (table, i, "DIFF_SN_M", diffStats.SNm); // "signal-to-noise of neg match src", 470 471 @ALL@ psFitsTableSetU32 (table, i, "FLAGS", source->mode); // "psphot analysis flags", 472 @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetU32 (table, i, "FLAGS2", source->mode2); // "psphot analysis flags", 473 @PS1_V3,PS1_SV2@ psFitsTableSetS32 (table, i, "PADDING2", 0); // "padding", 474 475 @ALL@ psFitsTableSetU16 (table, i, "N_FRAMES", source->nFrames); // "Number of frames overlapping source center", 476 @ALL@ psFitsTableSetS16 (table, i, "PADDING", 0); // "padding", 477 } 478 479 // XXX why do we make a copy here to be supplemented with the masks? why not do this in the calling function? 480 psMetadata *header = psMetadataCopy(NULL, tableHeader); 481 pmSourceMasksHeader(header); 482 483 psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); 484 if (!psFitsWriteTableNew(fits, header, table, extname)) { 485 psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); 486 psFree(table); 487 psFree(header); 488 return false; 489 } 490 psFree(table); 491 psFree(header); 492 493 return true; 494 } 495 496 bool pmSourcesWrite_CMF_@CMFMODE@_Old (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) 497 { 498 PS_ASSERT_PTR_NON_NULL(fits, false); 499 PS_ASSERT_PTR_NON_NULL(sources, false); 500 PS_ASSERT_PTR_NON_NULL(extname, false); 501 502 psArray *table; 503 psMetadata *row; 504 505 pmChip *chip = readout->parent->parent; 506 507 // if the sequence is defined, write these in seq order; otherwise 508 // write them in S/N order: 509 if (sources->n > 0) { 510 pmSource *source = (pmSource *) sources->data[0]; 511 if (source->seq == -1) { 512 // let's write these out in S/N order 513 sources = psArraySort (sources, pmSourceSortByFlux); 514 } else { 515 sources = psArraySort (sources, pmSourceSortBySeq); 516 } 517 } 518 519 table = psArrayAllocEmpty (sources->n); 520 521 float magOffset; 522 float zeroptErr; 523 float fwhmMajor; 524 float fwhmMinor; 525 pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); 526 527 // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state 528 // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. 529 for (int i = 0; i < sources->n; i++) { 530 // this is the source associated with this image 531 pmSource *thisSource = sources->data[i]; 532 533 // this is the "real" version of this source 534 pmSource *source = thisSource->parent ? thisSource->parent : thisSource; 535 536 // If source->seq is -1, source was generated in this analysis. If source->seq is 537 // not -1, source was read from elsewhere: in the latter case, preserve the source 538 // ID. source.seq is used instead of source.id since the latter is a const 539 // generated on Alloc, and would thus be wrong for read in sources. 540 if (source->seq == -1) { 541 source->seq = i; 542 } 543 544 // set the 'best' values for various output fields: 545 pmSourceOutputs outputs; 546 pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); 547 548 pmSourceOutputsMoments moments; 549 pmSourceOutputsSetMoments (&moments, source); 550 551 @PS1_DV?@ pmSourceDiffStats diffStats; 552 @PS1_DV?@ pmSourceDiffStatsInit(&diffStats); 553 @PS1_DV?@ if (source->diffStats) { 554 @PS1_DV?@ diffStats = *source->diffStats; 555 @PS1_DV?@ } 556 557 row = psMetadataAlloc (); 558 559 // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation 560 // This set of psMetadataAdd Entries marks the "----" "Start of the PSF segment" 561 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); 562 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); 563 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); 564 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); 565 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); 566 567 // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision 568 @PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F32, "PSF RA coordinate (degrees)", outputs.ra); 569 @PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F32, "PSF DEC coordinate (degrees)", outputs.dec); 570 571 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); 572 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); 573 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); 574 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); 575 576 @ALL,!PS1_V1,!PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); 577 @ALL,!PS1_V1,!PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux", source->psfFluxErr); 578 579 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); 580 @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in reported aperture", source->apMagRaw); 581 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", source->apRadius); 582 @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", source->apFlux); 583 @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", source->apFluxErr); 584 @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "AP_NPIX", PS_DATA_S32, "aperture unmasked pixels", source->apNpixels); 585 586 @<PS1_V3,PS1_SV1,PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); 587 588 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); 589 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); 590 591 // NOTE: RA & DEC (both double) need to be on an 8-byte boundary... 592 @ALL,!PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); 593 @ALL,!PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); 594 595 @>=PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); 596 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); 597 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); 598 599 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); 600 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); 601 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); 602 603 // PSF shape parameters: 604 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); 605 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); 606 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); 170 607 @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_CORE", PS_DATA_F32, "k term if defined", outputs.psfCore); 171 608 172 609 // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values 173 610 @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_FWHM_MAJ", PS_DATA_F32, "PSF FWHM (major axis)", outputs.psfMajorFWHM); 174 611 @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_FWHM_MIN", PS_DATA_F32, "PSF FWHM (minor axis)", outputs.psfMinorFWHM); 175 612 176 613 // psf data quality 177 @ALL@ 178 179 @ALL@ 180 @ALL@ 181 182 @ALL@ 183 @ALL@ 184 @ALL@ 185 186 @>PS1_V2,PS1_SV?@psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3);187 @>PS1_V2,PS1_SV?@psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3);188 @>PS1_V2,PS1_SV?@psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4);189 @>PS1_V2,PS1_SV?@psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4);190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 614 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); 615 @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", source->pixWeightNotPoor); 616 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); 617 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); 618 619 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); 620 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); 621 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); 622 623 @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); 624 @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); 625 @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); 626 @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); 627 628 // Lensing parameters: 629 if (source->lensingOBJ && source->lensingOBJ->smear) { 630 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->X11); 631 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->X12); 632 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->X22); 633 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->e1); 634 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->e2); 635 } 636 if (source->lensingOBJ && source->lensingOBJ->shear) { 637 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->X11); 638 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->X12); 639 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->X22); 640 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->e1); 641 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->e2); 642 } 643 if (false && source->lensingOBJ) { 644 // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above 645 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->e1); 646 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->e2); 647 } 648 649 if (source->lensingPSF && source->lensingPSF->smear) { 650 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X11); 651 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X12); 652 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X22); 653 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->e1); 654 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->e2); 655 } 656 if (source->lensingPSF && source->lensingPSF->shear) { 657 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X11); 658 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X12); 659 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X22); 660 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->e1); 661 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->e2); 662 } 663 if (source->lensingPSF) { 664 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->e1); 665 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->e2); 666 } 230 667 231 668 // if lensing params exist also include the backmapped chipID and chip coordinates 232 669 if (source->lensingPSF && source->lensingPSF->shear) { 233 670 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM", PS_DATA_S16, "id of warp input chip", source->chipNum); 234 671 @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X", PS_DATA_S16, "x coord in warp input chip", source->chipX); … … 244 681 @>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); 245 682 246 @>PS1_V3@ 247 @>PS1_V3@ 248 @>PS1_V3@ 249 250 251 @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X", PS_DATA_S16, "x coord in warp input chip",source->chipX);252 @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_Y", PS_DATA_S16, "y coord in warp input chip",source->chipY);253 254 255 @PS1_DV?@ 256 @PS1_DV?@ 257 @PS1_DV?@ 258 @PS1_DV?@ 259 @PS1_DV?@ 260 261 @>PS1_DV1@ 262 @>PS1_DV1@ 263 @>PS1_DV1@ 264 @>PS1_DV1@ 265 266 @ALL@ 267 268 269 270 271 272 273 274 275 276 277 278 279 280 @ALL@ 281 @ALL@ 683 @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_RAD", PS_DATA_F32, "Radius where object hits sky", source->skyRadius); 684 @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_FLUX", PS_DATA_F32, "Flux / pix where object hits sky", source->skyFlux); 685 @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_SLOPE", PS_DATA_F32, "d(Flux/pix)/dRadius where object hits sky", source->skySlope); 686 687 @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM", PS_DATA_S16, "id of warp input chip", source->chipNum); 688 @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X", PS_DATA_S16, "x coord in warp input chip", source->chipX); 689 @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_Y", PS_DATA_S16, "y coord in warp input chip", source->chipY); 690 @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING3", PS_DATA_S16, "more padding", 0); 691 692 @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); 693 @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); 694 @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); 695 @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); 696 @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); 697 698 @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P", PS_DATA_F32, "distance to positive match source", diffStats.Rp); 699 @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P", PS_DATA_F32, "signal-to-noise of pos match src", diffStats.SNp); 700 @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M", PS_DATA_F32, "distance to negative match source", diffStats.Rm); 701 @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M", PS_DATA_F32, "signal-to-noise of neg match src", diffStats.SNm); 702 703 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); 704 @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); 705 @PS1_V3,PS1_SV2@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2", PS_DATA_S32, "more padding", 0); 706 @PS1_SV?@ 707 708 // note that this definition is inconsistent with the definition in 709 // Ohana/src/libautocode/def. This version creates a table with data not 710 // properly aligned with the 8-byte boundaries. The structure defined by 711 // libautocode does this, but has a different order of elements (and adds 712 // padding2 to fix things). We have generated may files with PS1_SV1 as is, so 713 // I'll leave it. But in future a PS1_SV2 should be forced to match 714 // libautocode. Note that addstar knows to detect the alternate version of 715 // PS1_SV1 and correctly interpret its fields. 716 717 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", source->nFrames); 718 @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); 282 719 283 720 psArrayAdd (table, 100, row); … … 320 757 } 321 758 759 bool pmSourcesWrite_CMF_@CMFMODE@ (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) 760 { 761 // bool status = pmSourcesWrite_CMF_@CMFMODE@_Old (fits, readout, sources, imageHeader, tableHeader, extname, recipe); 762 bool status = pmSourcesWrite_CMF_@CMFMODE@_New (fits, readout, sources, imageHeader, tableHeader, extname, recipe); 763 return status; 764 } 765 322 766 // read in a readout from the fits file 323 psArray *pmSourcesRead_CMF_@CMFMODE@ (psFits *fits, psMetadata *header)767 psArray *pmSourcesRead_CMF_@CMFMODE@_New (psFits *fits, psMetadata *header) 324 768 { 325 769 PS_ASSERT_PTR_NON_NULL(fits, false); … … 355 799 psArray *sources = psArrayAlloc(numSources); // Array of sources, to return 356 800 357 // convert the table to the pmSource entriesa 801 // reads full table into memory 802 psFitsTable *table = psFitsReadTableNew (fits); 803 804 // convert the table to the pmSource entries 805 for (int i = 0; i < numSources; i++) { 806 pmSource *source = pmSourceAlloc (); 807 pmModel *model = pmModelAlloc (modelType); 808 source->modelPSF = model; 809 source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags 810 811 // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. 812 PAR = model->params->data.F32; 813 dPAR = model->dparams->data.F32; 814 815 @ALL@ source->seq = psFitsTableGetU32 (&status, table, i, "IPP_IDET"); 816 @ALL@ PAR[PM_PAR_XPOS] = psFitsTableGetF32 (&status, table, i, "X_PSF"); 817 @ALL@ PAR[PM_PAR_YPOS] = psFitsTableGetF32 (&status, table, i, "Y_PSF"); 818 @ALL@ dPAR[PM_PAR_XPOS] = psFitsTableGetF32 (&status, table, i, "X_PSF_SIG"); 819 @ALL@ dPAR[PM_PAR_YPOS] = psFitsTableGetF32 (&status, table, i, "Y_PSF_SIG"); 820 @ALL@ axes.major = psFitsTableGetF32 (&status, table, i, "PSF_MAJOR"); 821 @ALL@ axes.minor = psFitsTableGetF32 (&status, table, i, "PSF_MINOR"); 822 @ALL@ axes.theta = psFitsTableGetF32 (&status, table, i, "PSF_THETA"); 823 @ALL@ axes.theta = axes.theta * PS_RAD_DEG; 824 825 @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) { 826 @>PS1_V4,>PS1_SV2,>PS1_DV3@ PAR[PM_PAR_7] = psFitsTableGetF32 (&status, table, i, "PSF_CORE"); 827 @>PS1_V4,>PS1_SV2,>PS1_DV3@ } 828 829 @ALL@ PAR[PM_PAR_SKY] = psFitsTableGetF32 (&status, table, i, "SKY"); 830 @ALL@ dPAR[PM_PAR_SKY] = psFitsTableGetF32 (&status, table, i, "SKY_SIGMA"); 831 @ALL@ source->sky = PAR[PM_PAR_SKY]; 832 @ALL@ source->skyErr = dPAR[PM_PAR_SKY]; 833 834 // XXX use these to determine PAR[PM_PAR_I0]? 835 @ALL@ source->psfMag = psFitsTableGetF32 (&status, table, i, "PSF_INST_MAG"); 836 @ALL@ source->psfMagErr = psFitsTableGetF32 (&status, table, i, "PSF_INST_MAG_SIG"); 837 @ALL@ source->apMag = psFitsTableGetF32 (&status, table, i, "AP_MAG"); 838 @>PS1_V2,PS1_SV?,>PS1_DV1@ source->apMagRaw = psFitsTableGetF32 (&status, table, i, "AP_MAG_RAW"); 839 @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux = psFitsTableGetF32 (&status, table, i, "AP_FLUX"); 840 @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr = psFitsTableGetF32 (&status, table, i, "AP_FLUX_SIG"); 841 842 // XXX use these to determine PAR[PM_PAR_I0] if they exist? 843 // XXX add these to PS1_SV1? 844 @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFlux = psFitsTableGetF32 (&status, table, i, "PSF_INST_FLUX"); 845 @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFluxErr = psFitsTableGetF32 (&status, table, i, "PSF_INST_FLUX_SIG"); 846 847 // XXX this scaling is incorrect: does not include the 2 \pi AREA factor 848 @ALL@ PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; 849 @ALL@ dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; 850 851 pmPSF_AxesToModel (PAR, axes, model->class->useReff); 852 853 @ALL@ float peakMag = psFitsTableGetF32 (&status, table, i, "PEAK_FLUX_AS_MAG"); 854 @ALL@ float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; 855 856 // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) 857 @ALL@ source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); 858 @ALL@ source->peak->rawFlux = peakFlux; 859 @ALL@ source->peak->smoothFlux = peakFlux; 860 @ALL@ source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position 861 @ALL@ source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position 862 @ALL@ source->peak->dx = dPAR[PM_PAR_XPOS]; 863 @ALL@ source->peak->dy = dPAR[PM_PAR_YPOS]; 864 865 @ALL@ source->pixWeightNotBad = psFitsTableGetF32 (&status, table, i, "PSF_QF"); 866 @>PS1_V2,PS1_SV?,>PS1_DV1@ source->pixWeightNotPoor = psFitsTableGetF32 (&status, table, i, "PSF_QF_PERFECT"); 867 @ALL@ source->crNsigma = psFitsTableGetF32 (&status, table, i, "CR_NSIGMA"); 868 @ALL@ source->extNsigma = psFitsTableGetF32 (&status, table, i, "EXT_NSIGMA"); 869 @ALL@ source->apRadius = psFitsTableGetF32 (&status, table, i, "AP_MAG_RADIUS"); 870 @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels = psFitsTableGetS32 (&status, table, i, "AP_NPIX"); 871 872 // note that some older versions used PSF_PROBABILITY: this was not well defined. 873 @ALL@ model->chisq = psFitsTableGetF32 (&status, table, i, "PSF_CHISQ"); 874 @ALL@ model->nDOF = psFitsTableGetS32 (&status, table, i, "PSF_NDOF"); 875 @ALL@ model->nPix = psFitsTableGetS32 (&status, table, i, "PSF_NPIX"); 876 877 @ALL@ source->moments = pmMomentsAlloc (); 878 @ALL@ source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf 879 @ALL@ source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf 880 881 @ALL@ source->moments->Mxx = psFitsTableGetF32 (&status, table, i, "MOMENTS_XX"); 882 @ALL@ source->moments->Mxy = psFitsTableGetF32 (&status, table, i, "MOMENTS_XY"); 883 @ALL@ source->moments->Myy = psFitsTableGetF32 (&status, table, i, "MOMENTS_YY"); 884 885 // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data, 886 // we are storing enough information so the output will be consistent with the input 887 @>PS1_V2,PS1_SV?@ source->moments->Mxxx = +1.00 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M3C"); 888 @>PS1_V2,PS1_SV?@ source->moments->Mxxy = 0.00; 889 @>PS1_V2,PS1_SV?@ source->moments->Mxyy = 0.00; 890 @>PS1_V2,PS1_SV?@ source->moments->Myyy = -1.00 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M3S"); 891 @>PS1_V2,PS1_SV?@ source->moments->Mxxxx = +1.00 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M4C"); 892 @>PS1_V2,PS1_SV?@ source->moments->Mxxxy = 0.00; 893 @>PS1_V2,PS1_SV?@ source->moments->Mxxyy = 0.00; 894 @>PS1_V2,PS1_SV?@ source->moments->Mxyyy = -0.25 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M4S"); 895 @>PS1_V2,PS1_SV?@ source->moments->Myyyy = 0.00; 896 897 // Lensing parameters (on read if PS1_V5+) 898 if (haveLensOBJ) { 899 source->lensingOBJ = pmSourceLensingAlloc (); 900 source->lensingOBJ->smear = pmLensingParsAlloc(); 901 source->lensingOBJ->shear = pmLensingParsAlloc(); 902 903 @>PS1_V4@ source->lensingOBJ->smear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SM_OBJ"); 904 @>PS1_V4@ source->lensingOBJ->smear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SM_OBJ"); 905 @>PS1_V4@ source->lensingOBJ->smear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SM_OBJ"); 906 @>PS1_V4@ source->lensingOBJ->smear->e1 = psFitsTableGetF32 (&status, table, i, "LENS_E1_SM_OBJ"); 907 @>PS1_V4@ source->lensingOBJ->smear->e2 = psFitsTableGetF32 (&status, table, i, "LENS_E2_SM_OBJ"); 908 @>PS1_V4@ source->lensingOBJ->shear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SH_OBJ"); 909 @>PS1_V4@ source->lensingOBJ->shear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SH_OBJ"); 910 @>PS1_V4@ source->lensingOBJ->shear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SH_OBJ"); 911 @>PS1_V4@ source->lensingOBJ->shear->e1 = psFitsTableGetF32 (&status, table, i, "LENS_E1_SH_OBJ"); 912 @>PS1_V4@ source->lensingOBJ->shear->e2 = psFitsTableGetF32 (&status, table, i, "LENS_E2_SH_OBJ"); 913 } 914 915 @>PS1_V4@ source->chipNum = psFitsTableGetS16 (&status, table, i, "SRC_CHIP_NUM"); 916 @>PS1_V4@ source->chipX = psFitsTableGetS16 (&status, table, i, "SRC_CHIP_X"); 917 @>PS1_V4@ source->chipY = psFitsTableGetS16 (&status, table, i, "SRC_CHIP_Y"); 918 919 if (haveLensPSF) { 920 source->lensingPSF = pmSourceLensingAlloc (); 921 source->lensingPSF->smear = pmLensingParsAlloc(); 922 source->lensingPSF->shear = pmLensingParsAlloc(); 923 924 @>PS1_V4@ source->lensingPSF->smear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SM_PSF"); 925 @>PS1_V4@ source->lensingPSF->smear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SM_PSF"); 926 @>PS1_V4@ source->lensingPSF->smear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SM_PSF"); 927 @>PS1_V4@ source->lensingPSF->smear->e1 = psFitsTableGetF32 (&status, table, i, "LENS_E1_SM_PSF"); 928 @>PS1_V4@ source->lensingPSF->smear->e2 = psFitsTableGetF32 (&status, table, i, "LENS_E2_SM_PSF"); 929 @>PS1_V4@ source->lensingPSF->shear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SH_PSF"); 930 @>PS1_V4@ source->lensingPSF->shear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SH_PSF"); 931 @>PS1_V4@ source->lensingPSF->shear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SH_PSF"); 932 @>PS1_V4@ source->lensingPSF->shear->e1 = psFitsTableGetF32 (&status, table, i, "LENS_E1_SH_PSF"); 933 @>PS1_V4@ source->lensingPSF->shear->e2 = psFitsTableGetF32 (&status, table, i, "LENS_E2_SH_PSF"); 934 @>PS1_V4@ source->lensingPSF->e1 = psFitsTableGetF32 (&status, table, i, "LENS_E1_PSF"); 935 @>PS1_V4@ source->lensingPSF->e2 = psFitsTableGetF32 (&status, table, i, "LENS_E2_PSF"); 936 } 937 938 @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrf = psFitsTableGetF32 (&status, table, i, "MOMENTS_R1"); 939 @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrh = psFitsTableGetF32 (&status, table, i, "MOMENTS_RH"); 940 @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFlux = psFitsTableGetF32 (&status, table, i, "KRON_FLUX"); 941 @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFluxErr = psFitsTableGetF32 (&status, table, i, "KRON_FLUX_ERR"); 942 @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFinner = psFitsTableGetF32 (&status, table, i, "KRON_FLUX_INNER"); 943 @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFouter = psFitsTableGetF32 (&status, table, i, "KRON_FLUX_OUTER"); 944 945 @>PS1_V3@ source->skyRadius = psFitsTableGetF32 (&status, table, i, "SKY_LIMIT_RAD"); 946 @>PS1_V3@ source->skyFlux = psFitsTableGetF32 (&status, table, i, "SKY_LIMIT_FLUX"); 947 @>PS1_V3@ source->skySlope = psFitsTableGetF32 (&status, table, i, "SKY_LIMIT_SLOPE"); 948 949 @PS1_DV?@ int nPos = psFitsTableGetS32 (&status, table, i, "DIFF_NPOS"); 950 @PS1_DV?@ if (nPos) { 951 @PS1_DV?@ source->diffStats = pmSourceDiffStatsAlloc(); 952 @PS1_DV?@ source->diffStats->nGood = nPos; 953 @PS1_DV?@ source->diffStats->fRatio = psFitsTableGetF32 (&status, table, i, "DIFF_FRATIO"); 954 @PS1_DV?@ source->diffStats->nRatioBad = psFitsTableGetF32 (&status, table, i, "DIFF_NRATIO_BAD"); 955 @PS1_DV?@ source->diffStats->nRatioMask = psFitsTableGetF32 (&status, table, i, "DIFF_NRATIO_MASK"); 956 @PS1_DV?@ source->diffStats->nRatioAll = psFitsTableGetF32 (&status, table, i, "DIFF_NRATIO_ALL"); 957 958 @>PS1_DV1@ source->diffStats->Rp = psFitsTableGetF32 (&status, table, i, "DIFF_R_P"); 959 @>PS1_DV1@ source->diffStats->SNp = psFitsTableGetF32 (&status, table, i, "DIFF_SN_P"); 960 @>PS1_DV1@ source->diffStats->Rm = psFitsTableGetF32 (&status, table, i, "DIFF_R_M"); 961 @>PS1_DV1@ source->diffStats->SNm = psFitsTableGetF32 (&status, table, i, "DIFF_SN_M"); 962 @PS1_DV?@ } 963 964 @ALL@ source->mode = psFitsTableGetU32 (&status, table, i, "FLAGS"); 965 @>PS1_V2,PS1_SV?,>PS1_DV1@ source->mode2 = psFitsTableGetU32 (&status, table, i, "FLAGS2"); 966 @ALL@ source->nFrames = psFitsTableGetU16 (&status, table, i, "N_FRAMES"); 967 assert (status); 968 969 sources->data[i] = source; 970 } 971 psFree (table); 972 return sources; 973 } 974 975 // read in a readout from the fits file 976 psArray *pmSourcesRead_CMF_@CMFMODE@_Old (psFits *fits, psMetadata *header) 977 { 978 PS_ASSERT_PTR_NON_NULL(fits, false); 979 PS_ASSERT_PTR_NON_NULL(header, false); 980 981 bool status; 982 psF32 *PAR, *dPAR; 983 psEllipseAxes axes; 984 985 // define PSF model type 986 int defaultModelType = pmModelClassGetType ("PS_MODEL_GAUSS"); 987 int modelType = -1; 988 989 // if header does not define the model, default to a gaussian 990 char *PSF_NAME = psMetadataLookupStr (&status, header, "PSFMODEL"); 991 if (PSF_NAME != NULL) { 992 modelType = pmModelClassGetType (PSF_NAME); 993 } 994 // work around bug in psphotFullForce 995 if (modelType < 0) { 996 modelType = defaultModelType; 997 } 998 // assert (modelType > -1); 999 1000 // do we expect to find lensing parameters? 1001 bool haveLensOBJ = psMetadataLookupBool (&status, header, "LENS_OBJ"); 1002 bool haveLensPSF = psMetadataLookupBool (&status, header, "LENS_PSF"); 1003 1004 // We get the size of the table, and allocate the array of sources first because the table 1005 // is large and ephemeral --- when the table gets blown away, whatever is allocated after 1006 // the table is read blocks the free. In fact, it's better to read the table row by row. 1007 long numSources = psFitsTableSize(fits); // Number of sources in table 1008 psArray *sources = psArrayAlloc(numSources); // Array of sources, to return 1009 1010 // convert the table to the pmSource entries 358 1011 for (int i = 0; i < numSources; i++) { 359 1012 psMetadata *row = psFitsReadTableRow(fits, i); // Table row … … 382 1035 @ALL@ axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); 383 1036 @ALL@ axes.theta = axes.theta * PS_RAD_DEG; 384 385 386 387 1037 1038 @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) { 1039 @>PS1_V4,>PS1_SV2,>PS1_DV3@ PAR[PM_PAR_7] = psMetadataLookupF32 (&status, row, "PSF_CORE"); 1040 @>PS1_V4,>PS1_SV2,>PS1_DV3@ } 388 1041 389 1042 @ALL@ PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); … … 396 1049 @ALL@ source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); 397 1050 @ALL@ source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); 398 @>PS1_V2,PS1_SV?,>PS1_DV1@ source->apMagRaw = psMetadataLookupF32 (&status, row, "AP_MAG_RAW");399 @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux= psMetadataLookupF32 (&status, row, "AP_FLUX");400 1051 @>PS1_V2,PS1_SV?,>PS1_DV1@ source->apMagRaw = psMetadataLookupF32 (&status, row, "AP_MAG_RAW"); 1052 @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux = psMetadataLookupF32 (&status, row, "AP_FLUX"); 1053 @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr = psMetadataLookupF32 (&status, row, "AP_FLUX_SIG"); 401 1054 402 1055 // XXX use these to determine PAR[PM_PAR_I0] if they exist? 403 404 405 1056 // XXX add these to PS1_SV1? 1057 @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); 1058 @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); 406 1059 407 1060 // XXX this scaling is incorrect: does not include the 2 \pi AREA factor … … 424 1077 425 1078 @ALL@ source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); 426 1079 @>PS1_V2,PS1_SV?,>PS1_DV1@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT"); 427 1080 @ALL@ source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); 428 1081 @ALL@ source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); 429 1082 @ALL@ source->apRadius = psMetadataLookupF32 (&status, row, "AP_MAG_RADIUS"); 430 1083 @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels = psMetadataLookupS32 (&status, row, "AP_NPIX"); 431 1084 432 1085 // note that some older versions used PSF_PROBABILITY: this was not well defined. … … 443 1096 @ALL@ source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); 444 1097 445 446 1098 // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data, 1099 // we are storing enough information so the output will be consistent with the input 447 1100 @>PS1_V2,PS1_SV?@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C"); 448 1101 @>PS1_V2,PS1_SV?@ source->moments->Mxxy = 0.0; … … 456 1109 @>PS1_V2,PS1_SV?@ source->moments->Myyyy = 0.0; 457 1110 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 1111 // Lensing parameters (on read if PS1_V5+) 1112 if (haveLensOBJ) { 1113 source->lensingOBJ = pmSourceLensingAlloc (); 1114 source->lensingOBJ->smear = pmLensingParsAlloc(); 1115 source->lensingOBJ->shear = pmLensingParsAlloc(); 1116 1117 @>PS1_V4@ source->lensingOBJ->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_OBJ"); 1118 @>PS1_V4@ source->lensingOBJ->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_OBJ"); 1119 @>PS1_V4@ source->lensingOBJ->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_OBJ"); 1120 @>PS1_V4@ source->lensingOBJ->smear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SM_OBJ"); 1121 @>PS1_V4@ source->lensingOBJ->smear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SM_OBJ"); 1122 @>PS1_V4@ source->lensingOBJ->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_OBJ"); 1123 @>PS1_V4@ source->lensingOBJ->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_OBJ"); 1124 @>PS1_V4@ source->lensingOBJ->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_OBJ"); 1125 @>PS1_V4@ source->lensingOBJ->shear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SH_OBJ"); 1126 @>PS1_V4@ source->lensingOBJ->shear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SH_OBJ"); 1127 } 1128 1129 @>PS1_V4@ source->chipNum = psMetadataLookupS16 (&status, row, "SRC_CHIP_NUM"); 1130 @>PS1_V4@ source->chipX = psMetadataLookupS16 (&status, row, "SRC_CHIP_X"); 1131 @>PS1_V4@ source->chipY = psMetadataLookupS16 (&status, row, "SRC_CHIP_Y"); 1132 1133 if (haveLensPSF) { 1134 source->lensingPSF = pmSourceLensingAlloc (); 1135 source->lensingPSF->smear = pmLensingParsAlloc(); 1136 source->lensingPSF->shear = pmLensingParsAlloc(); 1137 1138 @>PS1_V4@ source->lensingPSF->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_PSF"); 1139 @>PS1_V4@ source->lensingPSF->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_PSF"); 1140 @>PS1_V4@ source->lensingPSF->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_PSF"); 1141 @>PS1_V4@ source->lensingPSF->smear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SM_PSF"); 1142 @>PS1_V4@ source->lensingPSF->smear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SM_PSF"); 1143 @>PS1_V4@ source->lensingPSF->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_PSF"); 1144 @>PS1_V4@ source->lensingPSF->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_PSF"); 1145 @>PS1_V4@ source->lensingPSF->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_PSF"); 1146 @>PS1_V4@ source->lensingPSF->shear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SH_PSF"); 1147 @>PS1_V4@ source->lensingPSF->shear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SH_PSF"); 1148 @>PS1_V4@ source->lensingPSF->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_PSF"); 1149 @>PS1_V4@ source->lensingPSF->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_PSF"); 1150 } 498 1151 499 1152 @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrf = psMetadataLookupF32 (&status, row, "MOMENTS_R1"); … … 508 1161 @>PS1_V3@ source->skySlope = psMetadataLookupF32 (&status, row, "SKY_LIMIT_SLOPE"); 509 1162 510 511 512 513 514 515 516 517 518 519 520 521 522 523 1163 @PS1_DV?@ int nPos = psMetadataLookupS32 (&status, row, "DIFF_NPOS"); 1164 @PS1_DV?@ if (nPos) { 1165 @PS1_DV?@ source->diffStats = pmSourceDiffStatsAlloc(); 1166 @PS1_DV?@ source->diffStats->nGood = nPos; 1167 @PS1_DV?@ source->diffStats->fRatio = psMetadataLookupF32 (&status, row, "DIFF_FRATIO"); 1168 @PS1_DV?@ source->diffStats->nRatioBad = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_BAD"); 1169 @PS1_DV?@ source->diffStats->nRatioMask = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_MASK"); 1170 @PS1_DV?@ source->diffStats->nRatioAll = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_ALL"); 1171 1172 @>PS1_DV1@ source->diffStats->Rp = psMetadataLookupF32 (&status, row, "DIFF_R_P"); 1173 @>PS1_DV1@ source->diffStats->SNp = psMetadataLookupF32 (&status, row, "DIFF_SN_P"); 1174 @>PS1_DV1@ source->diffStats->Rm = psMetadataLookupF32 (&status, row, "DIFF_R_M"); 1175 @>PS1_DV1@ source->diffStats->SNm = psMetadataLookupF32 (&status, row, "DIFF_SN_M"); 1176 @PS1_DV?@ } 524 1177 525 1178 @ALL@ source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); … … 533 1186 534 1187 return sources; 1188 } 1189 1190 psArray *pmSourcesRead_CMF_@CMFMODE@ (psFits *fits, psMetadata *header) { 1191 // psArray *array = pmSourcesRead_CMF_@CMFMODE@_Old (fits, header); 1192 psArray *array = pmSourcesRead_CMF_@CMFMODE@_New (fits, header); 1193 return array; 535 1194 } 536 1195 … … 599 1258 // write the radial profile apertures to header 600 1259 for (int i = 0; i < radMax->n; i++) { 601 602 603 604 1260 sprintf (keyword1, "RMIN_%02d", i); 1261 sprintf (keyword2, "RMAX_%02d", i); 1262 psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]); 1263 psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "max radius for SB profile", radMax->data.F32[i]); 605 1264 } 606 1265 607 1266 // we write out all sources, regardless of quality. the source flags tell us the state 608 1267 for (int i = 0; i < sources->n; i++) { 609 1268 // this is the source associated with this image 610 1269 pmSource *thisSource = sources->data[i]; 611 1270 612 613 1271 // this is the "real" version of this source 1272 pmSource *source = thisSource->parent ? thisSource->parent : thisSource; 614 1273 615 1274 // skip sources without measurements … … 639 1298 psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); 640 1299 641 642 643 644 645 646 647 1300 float AxialRatio = NAN; 1301 float AxialTheta = NAN; 1302 pmSourceExtendedPars *extpars = source->extpars; 1303 if (extpars) { 1304 AxialRatio = extpars->axes.minor / extpars->axes.major; 1305 AxialTheta = extpars->axes.theta; 1306 } 648 1307 psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO", PS_DATA_F32, "Axial Ratio of radial profile", AxialRatio); 649 1308 psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA", PS_DATA_F32, "Angle of radial profile ellipse", AxialTheta); … … 651 1310 // Petrosian measurements 652 1311 // XXX insert header data: petrosian ref radius, flux ratio 653 1312 // XXX check flags to see if Pet was measured 654 1313 if (doPetrosian) { 655 1314 pmSourceExtendedPars *extpars = source->extpars; 656 1315 if (extpars) { 657 658 659 660 661 662 663 664 665 666 667 668 669 670 1316 // XXX note that this mag is either calibrated or instrumental depending on existence of zero point 1317 float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point 1318 // NOTE EAM 20140806 : PETRO_MAG_ERR was inverted!! this allows for it to be repaired 1319 float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point 1320 if (repairMagErrors) { 1321 // I need to add the kron error in quadrature becasue pet_error ignores the object flux 1322 float Krf = source->moments ? source->moments->KronFlux : NAN; 1323 float dKrf = source->moments ? source->moments->KronFluxErr : NAN; 1324 if (isfinite (Krf) && isfinite (dKrf)) { 1325 magErr = sqrt(PS_SQR(1.0 / magErr) + PS_SQR(dKrf / Krf)); 1326 } else { 1327 magErr = 1.0 / magErr; 1328 } 1329 } 671 1330 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", mag); 672 1331 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", magErr); … … 705 1364 // Flux Annuli (if we have extended source measurements, we have these. only optionally save them) 706 1365 if (doAnnuli) { 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 1366 psVector *radSB = psVectorAlloc(radMin->n, PS_TYPE_F32); 1367 psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32); 1368 psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32); 1369 psVectorInit (radSB, NAN); 1370 psVectorInit (radFlux, NAN); 1371 psVectorInit (radFill, NAN); 1372 if (!source->extpars) goto empty_annuli; 1373 if (!source->extpars->radProfile) goto empty_annuli; 1374 if (!source->extpars->radProfile->binSB) goto empty_annuli; 1375 psAssert (source->extpars->radProfile->binSum, "programming error"); 1376 psAssert (source->extpars->radProfile->binFill, "programming error"); 1377 psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths"); 1378 psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths"); 1379 psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths"); 1380 1381 // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux 1382 for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) { 1383 radSB->data.F32[j] = source->extpars->radProfile->binSB->data.F32[j]; 1384 radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j]; 1385 radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j]; 1386 } 1387 1388 empty_annuli: 1389 psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB); 1390 psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); 1391 psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); 1392 psFree (radSB); 1393 psFree (radFlux); 1394 psFree (radFill); 1395 } 1396 if (nRow < 0) { 1397 nRow = row->list->n; 1398 } else { 1399 psAssert (nRow == row->list->n, "inconsistent row lengths"); 1400 } 1401 psArrayAdd (table, 100, row); 1402 psFree (row); 744 1403 } 745 1404 746 1405 if (table->n == 0) { 747 748 749 750 751 752 753 754 755 1406 if (!psFitsWriteBlank (fits, outhead, extname)) { 1407 psError(psErrorCodeLast(), false, "Unable to write empty sources file."); 1408 psFree(outhead); 1409 psFree(table); 1410 return false; 1411 } 1412 psFree (outhead); 1413 psFree (table); 1414 return true; 756 1415 } 757 1416 758 1417 psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); 759 1418 if (!psFitsWriteTable (fits, outhead, table, extname)) { 760 761 762 763 1419 psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); 1420 psFree (outhead); 1421 psFree(table); 1422 return false; 764 1423 } 765 1424 psFree (outhead); … … 885 1544 extpars->axes.theta = psMetadataLookupF32(&status, row, "F25_THETA"); 886 1545 887 1546 // magErr may have been saved in inverted form 888 1547 float mag = psMetadataLookupF32(&status, row, "PETRO_MAG"); 889 1548 float magErr = psMetadataLookupF32(&status, row, "PETRO_MAG_ERR"); 890 1549 if (isfinite(mag)) { 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 1550 extpars->petrosianFlux = pow(10., (magOffset - mag) / 2.5); 1551 if (isfinite(magErr)) { 1552 if (repairMagErrors) { 1553 // I need to add the kron error in quadrature becasue pet_error ignores the object flux 1554 float Krf = source->moments ? source->moments->KronFlux : NAN; 1555 float dKrf = source->moments ? source->moments->KronFluxErr : NAN; 1556 if (isfinite (Krf) && isfinite (dKrf)) { 1557 magErr = 1.0 / sqrt(PS_SQR(magErr) - PS_SQR(dKrf / Krf)); 1558 } else { 1559 magErr = 1.0 / magErr; 1560 } 1561 extpars->petrosianFluxErr = extpars->petrosianFlux * magErr; 1562 } else { 1563 extpars->petrosianFluxErr = extpars->petrosianFlux / magErr; 1564 } 1565 } 1566 } 908 1567 909 1568 extpars->petrosianRadius = psMetadataLookupF32(&status, row, "PETRO_RADIUS"); … … 973 1632 int nParamMax = 0; 974 1633 for (int i = 0; i < sources->n; i++) { 975 1634 // this is the source associated with this image 976 1635 pmSource *thisSource = sources->data[i]; 977 1636 978 979 1637 // this is the "real" version of this source 1638 pmSource *source = thisSource->parent ? thisSource->parent : thisSource; 980 1639 981 1640 if (source->modelFits == NULL) continue; … … 1004 1663 pmSource *thisSource = sources->data[i]; 1005 1664 1006 1007 1665 // this is the "real" version of this source 1666 pmSource *source = thisSource->parent ? thisSource->parent : thisSource; 1008 1667 1009 1668 // XXX if no model fits are saved, write out modelEXT? … … 1018 1677 1019 1678 // pmSourceExtFitPars *extPars = source->extFitPars->data[j]; 1020 1021 1022 1023 1024 1679 // assert (extPars); 1680 1681 // skip models which were not actually fitted 1682 // XXX 1683 if (model->flags & badModel) continue; 1025 1684 1026 1685 PAR = model->params->data.F32; … … 1029 1688 yPos = PAR[PM_PAR_YPOS]; 1030 1689 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1690 // for the extended source models, we do not always fit the centroid in the non-linear fitting process 1691 // current situation (hard-wired into psphotSourceFits.c:psphotFitPCM, 1692 // SERSIC, DEV, EXP : X,Y not fitted (PCM and not PCM) 1693 // TRAIL : X,Y are fitted 1694 // 1695 1696 // XXX this should be based on what happened, not on the model type 1697 if (model->type == modelTypeTrail) { 1698 xErr = dPAR[PM_PAR_XPOS]; 1699 yErr = dPAR[PM_PAR_YPOS]; 1700 } else { 1701 // this is definitely an underestimate since it does not 1702 // account for the extent of the source 1703 xErr = fwhmMajor * model->magErr / 2.35; 1704 yErr = fwhmMinor * model->magErr / 2.35; 1705 } 1706 1707 @>PS1_DV2,>PS1_SV3@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; 1708 @>PS1_DV2,>PS1_SV3@ float posAngle = 0.0; 1709 @>PS1_DV2,>PS1_SV3@ float pltScale = 0.0; 1710 @>PS1_DV2,>PS1_SV3@ pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); 1711 @>PS1_DV2,>PS1_SV3@ double raPos = ptSky.r*PS_DEG_RAD; 1712 @>PS1_DV2,>PS1_SV3@ double decPos = ptSky.d*PS_DEG_RAD; 1713 @>PS1_DV2,>PS1_SV3@ posAngle *= PS_DEG_RAD; 1714 @>PS1_DV2,>PS1_SV3@ pltScale *= PS_DEG_RAD*3600.0; 1715 1716 float kronFlux = source->moments ? source->moments->KronFlux : NAN; 1717 float kronMag = isfinite(kronFlux) ? -2.5*log10(kronFlux) : NAN; 1059 1718 1060 1719 row = psMetadataAlloc (); 1061 1720 1062 1063 1721 // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation 1722 // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment" 1064 1723 psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); 1065 1724 psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); … … 1070 1729 @>PS1_DV2,>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "RA_EXT", 0, "EXT model ra coordinate", raPos); 1071 1730 @>PS1_DV2,>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "DEC_EXT", 0, "EXT model dec coordinate", decPos); 1072 1073 1731 @>PS1_DV2@ float instFlux = isfinite(model->mag) ? pow(10.0, -0.4*model->mag) : NAN; 1732 @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_FLUX", 0, "EXT fit instrumental counts", instFlux); 1074 1733 1075 1734 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); 1076 1735 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); 1077 1736 1078 1079 1080 1081 1082 1083 1084 1085 1737 @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN; 1738 @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration", calMag); 1739 @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ", PS_DATA_F32, "EXT Model Chisq", model->chisq); 1740 @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF", PS_DATA_S32, "EXT Model num degrees of freedom", model->nDOF); 1741 @>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); 1742 1743 // EAM : adding for PV2 outputs: 1744 @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", model->flags); 1086 1745 1087 1746 @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE", 0, "position angle at source (degrees)", posAngle); 1088 1747 @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "PLTSCALE", 0, "plate scale at source (arcsec/pixel)", pltScale); 1089 1748 1090 1749 // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_XX", 0, "second moment in x", extPars->Mxx); … … 1103 1762 1104 1763 // XXX these should be major and minor, not 'x' and 'y' 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1764 if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) { 1765 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (major axis), length for trail", PAR[PM_PAR_LENGTH]); 1766 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (minor axis), sigma for trail", PAR[PM_PAR_SIGMA]); // this is not fitted 1767 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", PAR[PM_PAR_THETA]); 1768 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)", dPAR[PM_PAR_LENGTH]); 1769 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)", NAN); // this is not fitted, so error is NAN 1770 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_THETA]); 1771 } else { 1772 if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY])) { 1773 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (SXX, isnan)", PAR[PM_PAR_SXX]); 1774 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (SYY, isnan)", PAR[PM_PAR_SYY]); 1775 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT angle (SXY, isnan)", PAR[PM_PAR_SXY]); 1776 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR", 0, "EXT width err (SXX, isnan)", dPAR[PM_PAR_SXX]); 1777 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR", 0, "EXT width err (SYY, isnan)", dPAR[PM_PAR_SYY]); 1778 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]); 1779 } else { 1780 psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->class->useReff); 1781 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (major axis), length for trail", axes.major); 1782 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (minor axis), sigma for trail", axes.minor); 1783 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); 1784 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)", dPAR[PM_PAR_SXX]); 1785 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)", dPAR[PM_PAR_SYY]); 1786 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_SXY]); 1787 } 1788 } 1130 1789 1131 1790 // write out the other generic parameters … … 1140 1799 1141 1800 snprintf (name, 64, "EXT_PAR_%02d", k); 1142 1801 1143 1802 if (k < model->params->n) { 1144 1803 psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]); … … 1148 1807 } 1149 1808 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 }1160 1809 // optionally, write out the covariance matrix values 1810 // XXX do I need to pad this to match the biggest covar matrix? 1811 if (false && model->covar) { 1812 for (int iy = 0; iy < model->covar->numCols; iy++) { 1813 for (int ix = iy; ix < model->covar->numCols; ix++) { 1814 snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix); 1815 psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]); 1816 1817 } 1818 } 1819 } 1161 1820 psArrayAdd (table, 100, row); 1162 1821 psFree (row); … … 1280 1939 PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07"); 1281 1940 // XXX add an error: 1282 1283 } 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1941 // dPAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07_"); 1942 } 1943 1944 // NOTE: we no longer write out the covariance matrix 1945 if (false) { 1946 // read the covariance matrix 1947 int nparams = model->params->n; 1948 psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32); 1949 for (int y = 0; y < nparams; y++) { 1950 for (int x = 0; x < nparams; x++) { 1951 char name[64]; 1952 snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x); 1953 covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name); 1954 } 1955 } 1956 model->covar = covar; 1957 } 1958 1959 // we are only saving the values stored in dPAR[SXX,etc] 1301 1960 dPAR[PM_PAR_SXX] = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ_ERR"); 1302 1961 dPAR[PM_PAR_SYY] = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN_ERR"); 1303 1962 dPAR[PM_PAR_SXY] = psMetadataLookupF32(&status, row, "EXT_THETA_ERR"); 1304 1963 1305 1964 // other parameters that we need to read 1306 1965 PAR[PM_PAR_SKY] = psMetadataLookupF32(&status, row, "SKY_EXT"); 1307 1966 … … 1341 2000 // perform full non-linear fits / extended source analysis? 1342 2001 if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) { 1343 1344 2002 psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n"); 2003 return true; 1345 2004 } 1346 2005 … … 1385 2044 for (int i = 0; i < sources->n; i++) { 1386 2045 1387 2046 // this is the source associated with this image 1388 2047 pmSource *thisSource = sources->data[i]; 1389 2048 1390 1391 2049 // this is the "real" version of this source 2050 pmSource *source = thisSource->parent ? thisSource->parent : thisSource; 1392 2051 1393 2052 // skip sources without radial aper measurements (or insufficient) 1394 2053 if (source->radialAper == NULL) continue; 1395 2054 1396 2055 // psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set"); 1397 2056 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 2057 for (int entry = 0; entry < source->radialAper->n; entry++) { 2058 2059 // choose the convolved EXT model, if available, otherwise the simple one 2060 pmSourceRadialApertures *radialAper = source->radialAper->data[entry]; 2061 assert (radialAper); 2062 2063 if (pmSourcePositionUseMoments(source)) { 2064 xPos = source->moments->Mx; 2065 yPos = source->moments->My; 2066 } else { 2067 xPos = source->peak->xf; 2068 yPos = source->peak->yf; 2069 } 2070 2071 row = psMetadataAlloc (); 2072 2073 // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) 2074 // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation 2075 // This set of psMetadataAdd Entries marks the "----" "Start of the XRAD segment" 2076 psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); 2077 psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER", 0, "Center of aperture measurements", xPos); 2078 psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER", 0, "Center of aperture measurements", yPos); 2079 if (fwhmValues) { 2080 psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM", 0, "FWHM of matched PSF", fwhmValues->data.F32[entry]); 2081 } else { 2082 psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM", 0, "image is not FWHM-matched", NAN); 2083 } 2084 2085 // XXX if we have raw radial apertures, write them out here 2086 psVector *radFlux = psVectorAlloc(radMax->n, PS_TYPE_F32); 2087 psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32); 2088 psVector *radFill = psVectorAlloc(radMax->n, PS_TYPE_F32); 2089 psVector *radFluxStdev = psVectorAlloc(radMax->n, PS_TYPE_F32); 2090 psVectorInit (radFlux, NAN); 2091 psVectorInit (radFluxErr, NAN); 2092 psVectorInit (radFill, NAN); 2093 if (!radialAper->flux) goto write_annuli; 2094 if (!radialAper->fill) goto write_annuli; 2095 psAssert (radialAper->flux->n <= radFlux->n, "inconsistent vector lengths"); 2096 psAssert (radialAper->fill->n <= radFlux->n, "inconsistent vector lengths"); 2097 2098 // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux 2099 for (int j = 0; j < radialAper->flux->n; j++) { 2100 radFlux->data.F32[j] = radialAper->flux->data.F32[j]; 2101 radFluxErr->data.F32[j] = radialAper->fluxErr->data.F32[j]; 2102 radFluxStdev->data.F32[j] = radialAper->fluxStdev->data.F32[j]; 2103 radFill->data.F32[j] = radialAper->fill->data.F32[j]; 2104 } 2105 2106 write_annuli: 2107 psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX", PS_META_REPLACE, "flux within annuli", radFlux); 2108 psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_META_REPLACE, "flux error in annuli", radFluxErr); 2109 psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_META_REPLACE, "flux standard deviation", radFluxStdev); 2110 psMetadataAddVector (row, PS_LIST_TAIL, "APER_FILL", PS_META_REPLACE, "fill factor of annuli", radFill); 2111 psFree (radFlux); 2112 psFree (radFluxErr); 2113 psFree (radFluxStdev); 2114 psFree (radFill); 2115 2116 psArrayAdd (table, 100, row); 2117 psFree (row); 2118 } 1460 2119 } 1461 2120 … … 1581 2240 // perform full non-linear fits / extended source analysis? 1582 2241 if (!psMetadataLookupBool (&status, recipe, "GALAXY_SHAPES")) { 1583 1584 2242 psLogMsg ("psphot", PS_LOG_INFO, "galaxy shapes were not measured, skipping\n"); 2243 return true; 1585 2244 } 1586 2245 … … 1609 2268 pmSource *thisSource = sources->data[i]; 1610 2269 1611 1612 1613 1614 2270 // this is the "real" version of this source 2271 pmSource *source = thisSource->parent ? thisSource->parent : thisSource; 2272 2273 // if we did not fit the galaxy model, modelFits will be NULL 1615 2274 if (source->modelFits == NULL) continue; 1616 2275 1617 2276 // if we did not fit the galaxy model, galaxyFits will also be NULL 1618 2277 if (source->galaxyFits == NULL) continue; 1619 2278
Note:
See TracChangeset
for help on using the changeset viewer.
