Index: /branches/eam_branch_20080324/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- /branches/eam_branch_20080324/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 17326) +++ /branches/eam_branch_20080324/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 17327) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.9.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-04-03 06:01:38 $ + * @version $Revision: 1.9.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-05 01:48:48 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -251,7 +251,8 @@ } -bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname) +bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) { + bool status; psArray *table; psMetadata *row; @@ -270,4 +271,14 @@ table = psArrayAllocEmpty (sources->n); + + // which extended source analyses should we perform? + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + + psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + assert (radialBinsLower->n == radialBinsUpper->n); // we write out all sources, regardless of quality. the source flags tell us the state @@ -306,64 +317,102 @@ // Petrosian measurements - pmSourcePetrosianValues *petrosian = source->extpars->petrosian; - if (petrosian) { - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", petrosian->mag); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", petrosian->rad); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", petrosian->radErr); - } + // XXX insert header data: petrosian ref radius, flux ratio + if (doPetrosian) { + pmSourcePetrosianValues *petrosian = source->extpars->petrosian; + if (petrosian) { + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", petrosian->mag); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", petrosian->rad); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", petrosian->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); + } + } // Kron measurements - pmSourceKronValues *kron = source->extpars->kron; - if (kron) { - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", kron->mag); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", kron->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", kron->rad); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", kron->radErr); + if (doKron) { + pmSourceKronValues *kron = source->extpars->kron; + if (kron) { + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", kron->mag); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", kron->magErr); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", kron->rad); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", kron->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", NAN); + } } // Isophot measurements // XXX insert header data: isophotal level - pmSourceIsophotalValues *isophot = source->extpars->isophot; - if (isophot) { - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", isophot->mag); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", isophot->rad); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", isophot->radErr); - } - - pmSourceAnnuli *annuli = source->extpars->annuli; - if (annuli) { - psVector *fluxVal = annuli->flux; - psVector *fluxErr = annuli->fluxErr; - psVector *fluxVar = annuli->fluxVar; - - for (int j = 0; j < fluxVal->n; j++) { - char name[32]; - sprintf (name, "FLUX_VAL_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]); - sprintf (name, "FLUX_ERR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]); - sprintf (name, "FLUX_VAR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]); + if (doIsophotal) { + pmSourceIsophotalValues *isophot = source->extpars->isophot; + if (isophot) { + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", isophot->mag); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", isophot->rad); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", isophot->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", NAN); } } - psArrayAdd (table, 100, row); - psFree (row); + // Flux Annuli + if (doAnnuli) { + pmSourceAnnuli *annuli = source->extpars->annuli; + if (annuli) { + psVector *fluxVal = annuli->flux; + psVector *fluxErr = annuli->fluxErr; + psVector *fluxVar = annuli->fluxVar; + + for (int j = 0; j < fluxVal->n; j++) { + char name[32]; + sprintf (name, "FLUX_VAL_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]); + sprintf (name, "FLUX_ERR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]); + sprintf (name, "FLUX_VAR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]); + } + } else { + for (int j = 0; j < radialBinsLower->n; j++) { + char name[32]; + sprintf (name, "FLUX_VAL_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN); + sprintf (name, "FLUX_ERR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN); + sprintf (name, "FLUX_VAR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN); + } + } + } + + psArrayAdd (table, 100, row); + psFree (row); } if (table->n == 0) { - psFitsWriteBlank (fits, outhead, extname); - psFree (table); - return true; + psFitsWriteBlank (fits, outhead, extname); + psFree (outhead); + psFree (table); + return true; } psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(PS_ERR_IO, false, "writing ext data %s\n", extname); - psFree(table); - return false; - } + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); psFree (table); @@ -391,4 +440,16 @@ sources = psArraySort (sources, pmSourceSortBySN); + // 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++) { + pmSource *source = sources->data[i]; + 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); + } + } + table = psArrayAllocEmpty (sources->n); @@ -420,20 +481,24 @@ // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG", PS_DATA_F32, "EXT fit instrumental magnitude", model->mag); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", model->magErr); + 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); + 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); + + psMetadataAddF32 (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' - psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", PS_DATA_F32, "EXT width in x coordinate", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", PS_DATA_F32, "EXT width in y coordinate", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA", PS_DATA_F32, "EXT orientation angle", axes.theta); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); // write out the other generic parameters - for (int k = 0; k < model->params->n; k++) { + 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; @@ -443,5 +508,10 @@ snprintf (name, 64, "EXT_PAR_%02d", k); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); + + if (k < model->params->n) { + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); + } else { + psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); + } } @@ -456,17 +526,19 @@ if (table->n == 0) { - psFitsWriteBlank (fits, outhead, extname); - psFree (table); - return true; + psFitsWriteBlank (fits, outhead, extname); + psFree (outhead); + psFree (table); + return true; } psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(PS_ERR_IO, false, "writing ext data %s\n", extname); - psFree(table); - return false; - } + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); psFree (table); - return true; }