Index: branches/eam_branches/ipp-20130904/psModules/src/camera/pmFPAfileIO.c =================================================================== --- branches/eam_branches/ipp-20130904/psModules/src/camera/pmFPAfileIO.c (revision 36269) +++ branches/eam_branches/ipp-20130904/psModules/src/camera/pmFPAfileIO.c (revision 36270) @@ -485,4 +485,5 @@ case PM_FPA_FILE_CMP: case PM_FPA_FILE_CMF: + case PM_FPA_FILE_CFF: status = pmFPAviewWriteObjects (view, file, config); break; @@ -515,8 +516,4 @@ case PM_FPA_FILE_WCS: psError(PS_ERR_IO, true, "cannot write type WCS (%s)", file->name); - return false; - - case PM_FPA_FILE_CFF: - psError(PS_ERR_IO, true, "cannot write type CFF (%s)", file->name); return false; Index: branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO.c =================================================================== --- branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO.c (revision 36269) +++ branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO.c (revision 36270) @@ -645,4 +645,17 @@ return false; + case PM_FPA_FILE_CFF: { + hdu = pmFPAviewThisHDU (view, fpa); + pmConfigConformHeader(hdu->header, file->format); + psFitsWriteBlank (file->fits, hdu->header, NULL); + file->header = hdu->header; + file->wrote_phu = true; + if (!pmSourcesWrite_CFF(file->fits, sources, hdu->header, recipe)) { + psError(PS_ERR_UNKNOWN, false, "failed to write CFF"); + return false; + } + break; + } + default: fprintf (stderr, "warning: type mismatch\n"); Index: branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO.h =================================================================== --- branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO.h (revision 36269) +++ branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO.h (revision 36270) @@ -54,4 +54,5 @@ psArray *pmSourcesReadCMP (char *filename, psMetadata *header); psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header); +bool pmSourcesWrite_CFF (psFits *fits, psArray *sources, psMetadata *header, psMetadata *recipe); bool pmSourcesWritePSFs (psArray *sources, char *filename); Index: branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO_CFF.c =================================================================== --- branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO_CFF.c (revision 36269) +++ branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO_CFF.c (revision 36270) @@ -65,5 +65,5 @@ int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); - char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSFMODEL"); if (PSF_NAME != NULL) { modelType = pmModelClassGetType (PSF_NAME); @@ -107,5 +107,6 @@ // XXX: we need to put a lookup table in the cff header to define the correspondence of the // model type values in the cff with our models. (We want to use an interger for efficiency - // but the values depend on the array in pmModelClass.c + // but the value for each model type is set on the organization of the the array in pmModelClass.c + // For now use the input values verbatim and trust the user that this is valid value int galaxyModelType = psMetadataLookupS32(&status, row, "MODEL_TYPE"); float Sindex = psMetadataLookupF32 (&status, row, "INDEX"); // Should this be PAR_07 not sersic index @@ -173,5 +174,7 @@ } - if (fitGalaxy) { + // XXX: should use < 0 as invalid galaxyModelType + + if (fitGalaxy && galaxyModelType > 0) { source->modelFits = psArrayAllocEmpty (1); pmModel *model = pmModelAlloc(galaxyModelType); @@ -192,4 +195,6 @@ xPAR[PM_PAR_7] = 0.5 / Sindex; } + // XXX: set source->modelEXT to this model and flag as extended so that we have a better + // shot at subtracting extended sources? psArrayAdd (source->modelFits, 1, model); @@ -204,2 +209,112 @@ } +bool pmSourcesWrite_CFF(psFits *fits, psArray *sources, psMetadata *header, psMetadata *recipe) { + + char *extname = "SkyChip.cff"; + + psArray *table = psArrayAllocEmpty(sources->n); + + pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSCIC"); + + for (int i = 0; i < sources->n; i++) { + pmSource *thisSource = sources->data[i]; + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + + psF32 xPos, yPos, flux, rMajor, rMinor, theta; + psS32 modelType = 0; + bool fitGalaxy = false; + bool psfStar = false; + psS32 sersicIndex = 0; + if (source->modelFits == NULL) { + pmModel *model = source->modelPSF; + if (model == NULL) continue; + psF32 *PAR = model->params->data.F32; + if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY]) || + !isfinite(source->psfFlux)) { + continue; + } + + xPos = model->params->data.F32[PM_PAR_XPOS]; + yPos = model->params->data.F32[PM_PAR_YPOS]; + flux = source->psfFlux; + psfStar = (source->mode & PM_SOURCE_MODE_PSFSTAR) ? true : false; + rMajor = 0; + rMinor = 0; + theta = 0; + } else { + // Choose the extended model to use for this source. + // For now we use the one set as the modelEXT for this source which had the best fit. + // XXX: make this controllable by recipe + // use best (as implemented here) + // choose specific type + int iModel = -1; + if (source->modelEXT) { + pmModelType ext_model_type = source->modelEXT->type; + for (int j=0; jmodelFits->n; j++) { + pmModel *aModel = source->modelFits->data[j]; + if (aModel->type == ext_model_type) { + iModel = j; + break; + } + } + } + if (iModel == -1) { + // Can this happen? perhaps if model type is qgauss (extended but below s/n for model fits? + // XXX: Should this be an assert? + continue; + } + pmModel *model = source->modelFits->data[iModel]; + psF32 *PAR = model->params->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) { + // XXX: do we need to handle this type + continue; + } else { + if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY]) || + !isfinite(model->mag)) { + // bad model + continue; + } + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type); + rMajor = axes.major; + rMinor = axes.minor; + theta = axes.theta*PS_DEG_RAD; + flux = pow(10.0, -0.4*model->mag); + fitGalaxy = true; + } + modelType = model->type; + if (modelType == sersicModelType) { + sersicIndex = 0.5 / PAR[PM_PAR_7]; + } + } + psMetadata *row = psMetadataAlloc(); + psMetadataAddU32 (row, PS_LIST_TAIL, "ID", 0, "IPP detection identifier", source->seq); + psMetadataAddF32 (row, PS_LIST_TAIL, "X", 0, "x coordinate", xPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y", 0, "y coordinate", yPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "FLUX", 0, "flux", flux); + psMetadataAddF32 (row, PS_LIST_TAIL, "AP_RADIUS", 0, "aperture radius", source->apRadius); + psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_RADIUS", 0, "Kron radius", source->moments->Mrf * 2.5); + psMetadataAddF32 (row, PS_LIST_TAIL, "PETRO_RADIUS", 0, "Petrosian Radius", source->extpars ? source->extpars->petrosianRadius : NAN); + psMetadataAddBool (row, PS_LIST_TAIL, "FIT_GALAXY", 0, "source has xfit", fitGalaxy); + psMetadataAddBool (row, PS_LIST_TAIL, "PSF_STAR", 0, "source was psf star", psfStar); + psMetadataAddF32 (row, PS_LIST_TAIL, "R_MAJOR", 0, "radius of major axis", rMajor); + psMetadataAddF32 (row, PS_LIST_TAIL, "R_MINOR", 0, "radius of minor axis", rMinor); + psMetadataAddF32 (row, PS_LIST_TAIL, "THETA", 0, "theta", theta); + psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "model type", modelType); + psMetadataAddS32 (row, PS_LIST_TAIL, "INDEX", 0, "sersic index", sersicIndex); + + psArrayAdd(table, 100, row); + } + + if (!psFitsWriteTable(fits, NULL, table, extname)) { + psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); + psFree(table); + psFree(header); + return false; + } + psFree(table); + // psFree(header); + + return true; +}