Index: branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO_CMF.c.in (revision 36157) +++ branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO_CMF.c.in (revision 36158) @@ -1343,2 +1343,213 @@ return true; } + +// XXX this layout is still the same as PS1_DEV_1 +bool pmSourcesWrite_CMF_@CMFMODE@_GAL (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) +{ + + psArray *table; + psMetadata *row; + psF32 *PAR, *dPAR; + psF32 xPos, yPos; + psF32 xErr, yErr; + char name[64]; + + // create a header to hold the output data + psMetadata *outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "galaxy table extension", extname); + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortByFlux); + + // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams + int nParamMax = 0; + for (int i = 0; i < sources->n; i++) { + // this is the source associated with this image + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + + if (source->modelFits == NULL) continue; + for (int j = 0; j < source->modelFits->n; j++) { + pmModel *model = source->modelFits->data[j]; + assert (model); + nParamMax = PS_MAX (nParamMax, model->params->n); + } + } + + @>PS1_DV2@ pmChip *chip = readout->parent->parent; + + table = psArrayAllocEmpty (sources->n); + + // we write out all sources, regardless of quality. the source flags tell us the state + for (int i = 0; i < sources->n; i++) { + + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + + // XXX if no model fits are saved, write out modelEXT? + if (source->modelFits == NULL) continue; + + // We have multiple sources : need to flag the one used to subtract the light (the 'best' model) + for (int j = 0; j < source->modelFits->n; j++) { + + // choose the convolved EXT model, if available, otherwise the simple one + pmModel *model = source->modelFits->data[j]; + assert (model); + + // pmSourceExtFitPars *extPars = source->extFitPars->data[j]; + // assert (extPars); + + // skip models which were not actually fitted + if (model->flags & PM_MODEL_STATUS_BADARGS) continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + + @>PS1_DV2@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; + @>PS1_DV2@ float posAngle = 0.0; + @>PS1_DV2@ float pltScale = 0.0; + @>PS1_DV2@ pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + @>PS1_DV2@ double raPos = ptSky.r*PS_DEG_RAD; + @>PS1_DV2@ double decPos = ptSky.d*PS_DEG_RAD; + @>PS1_DV2@ posAngle *= PS_DEG_RAD; + @>PS1_DV2@ pltScale *= PS_DEG_RAD*3600.0; + + float kronFlux = source->moments ? source->moments->KronFlux : NAN; + float kronMag = isfinite(kronFlux) ? -2.5*log10(kronFlux) : NAN; + + row = psMetadataAlloc (); + + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) + // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation + // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment" + psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT", 0, "EXT model y coordinate", yPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); + @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "RA_EXT", 0, "EXT model ra coordinate", raPos); + @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "DEC_EXT", 0, "EXT model dec coordinate", decPos); + @>PS1_DV2@ float instFlux = isfinite(model->mag) ? pow(10.0, -0.4*model->mag) : NAN; + @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_FLUX", 0, "EXT fit instrumental counts", instFlux); + + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); + + @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN; + @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration", calMag); + @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ", PS_DATA_F32, "EXT Magnitude using supplied calibration", model->chisq); + @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF", PS_DATA_S32, "EXT Magnitude using supplied calibration", model->nDOF); + + @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE", 0, "position angle at source (degrees)", posAngle); + @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "PLTSCALE", 0, "plate scale at source (arcsec/pixel)", pltScale); + + // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_XX", 0, "second moment in x", extPars->Mxx); + // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_XY", 0, "second moment in x,y", extPars->Mxy); + // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_YY", 0, "second moment in y", extPars->Myy); + // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_R1", 0, "first radial moment", extPars->Mrf); + // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_RH", 0, "half radial moment", extPars->Mrh); + + psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_INST_MAG", 0, "PSF fit instrumental magnitude", source->psfMag); + psMetadataAddF32 (row, PS_LIST_TAIL, "AP_MAG", 0, "PSF-sized aperture magnitude", source->apMag); + psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_MAG", 0, "Kron Mag", kronMag); + + @PS1_DV0,PS1_DV1@ psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); + @ALL,!PS1_DV0,!PS1_DV1@ psMetadataAddS32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); + psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "name of model", pmModelClassGetName (model->type)); + + // XXX these should be major and minor, not 'x' and 'y' + if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) { + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (major axis), length for trail", PAR[PM_PAR_LENGTH]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (minor axis), sigma for trail", PAR[PM_PAR_SIGMA]); // this is not fitted + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", PAR[PM_PAR_THETA]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)", dPAR[PM_PAR_LENGTH]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)", NAN); // this is not fitted, so error is NAN + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_THETA]); + } else { + if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY])) { + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (SXX, isnan)", PAR[PM_PAR_SXX]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (SYY, isnan)", PAR[PM_PAR_SYY]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT angle (SXY, isnan)", PAR[PM_PAR_SXY]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR", 0, "EXT width err (SXX, isnan)", dPAR[PM_PAR_SXX]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR", 0, "EXT width err (SYY, isnan)", dPAR[PM_PAR_SYY]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]); + } else { + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (major axis), length for trail", axes.major); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (minor axis), sigma for trail", axes.minor); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)", dPAR[PM_PAR_SXX]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)", dPAR[PM_PAR_SYY]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_SXY]); + } + } + + // write out the other generic parameters + for (int k = 0; k < nParamMax; k++) { + if (k == PM_PAR_I0) continue; + if (k == PM_PAR_SKY) continue; + if (k == PM_PAR_XPOS) continue; + if (k == PM_PAR_YPOS) continue; + if (k == PM_PAR_SXX) continue; + if (k == PM_PAR_SXY) continue; + if (k == PM_PAR_SYY) continue; + + snprintf (name, 64, "EXT_PAR_%02d", k); + + if (k < model->params->n) { + psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]); + } else { + psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", NAN); + } + } + + // optionally, write out the covariance matrix values + // XXX do I need to pad this to match the biggest covar matrix? + if (model->covar) { + for (int iy = 0; iy < model->covar->numCols; iy++) { + for (int ix = iy; ix < model->covar->numCols; ix++) { + snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix); + psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]); + + } + } + } + psArrayAdd (table, 100, row); + psFree (row); + } + } + + if (table->n == 0) { + if (!psFitsWriteBlank (fits, outhead, extname)) { + psError(psErrorCodeLast(), false, "Unable to write empty sources file."); + psFree(outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable (fits, outhead, table, extname)) { + psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; +} +