Index: /trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO.c (revision 12401) +++ /trunk/psModules/src/objects/pmSourceIO.c (revision 12402) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-03-10 07:03:20 $ + * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-03-11 18:56:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -329,11 +329,31 @@ } - exttype = psStringCopy ("SMPDATA"); // XXX derive this from the recipes + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, "PSPHOT"); + if (!status) { + psError(PS_ERR_UNKNOWN, true, "missing recipe PSPHOT in config data"); + return false; + } + exttype = psMetadataLookupMetadata(&status, recipe, "OUTPUT.FORMAT"); + if (!status) { + psError(PS_ERR_UNKNOWN, true, "missing OUTPUT.FORMAT in PSPHOT recipe"); + return false; + } + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTHEAD", PS_META_REPLACE, "name of image extension w/", headname); psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTTYPE", PS_META_REPLACE, "extension type", exttype); psFree (exttype); - if (!pmSourcesWriteCMF (file->fits, sources, outhead, dataname)) { - psError(PS_ERR_IO, false, "writing CMF data to %s\n", file->filename); + // XXX these are case-sensitive since the EXTYPE is case-sensitive + + status = false; + if (!strcmp (exttype, "SMPDATA")) { + status = pmSourcesWrite_SMPDATA (file->fits, sources, outhead, dataname); + } + if (!strcmp (exttype, "PS1_DEV_0")) { + status = pmSourcesWrite_PS1_DEV_0 (file->fits, sources, outhead, dataname); + } + + if (!status) { + psError(PS_ERR_IO, false, "writing CMF data to %s with format %s\n", file->filename, exttype); psFree (headname); psFree (dataname); @@ -497,4 +517,6 @@ } + XXX need to test for EXTNAME vs SMPDATA, PS1_DEV_0, etc + // we need to find the corresponding table EXTNAME. // first check the header Index: /trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 12402) +++ /trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 12402) @@ -0,0 +1,205 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-03-11 18:56:38 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceIO.h" + +// panstars-style FITS table output (header + table in 1st extension) +// this format consists of a header derived from the image header +// followed by a zero-size matrix, followed by the table data + +// this output format is valid for psphot analysis of an image, and does not include calibrated +// values derived in the DVO database. +// XXX how do I generate the source tables which I need to send to PSPS? + +bool pmSourcesWrite_PS1_DEV_0 (psFits *fits, psArray *sources, psMetadata *header, char *extname) +{ + + psArray *table; + psMetadata *row; + int i, type; + psF32 *PAR, *dPAR; + float lsky; + bool status; + psEllipseShape shape; + psEllipseAxes axes; + + // find config information for output header + float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); + + table = psArrayAllocEmpty (sources->n); + + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + pmModel *model = pmSourceSelectModel (source); + if (model == NULL) + continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + type = pmSourceDophotType (source); + lsky = (source->sky < 1.0) ? 0.0 : log10(source->sky); + + shape.sx = PAR[PM_PAR_SXX]; + shape.sy = PAR[PM_PAR_SYY]; + shape.sxy = PAR[PM_PAR_SXY]; + axes = psEllipseShapeToAxes (shape, 20.0); + + peakMag = (source->peak->flux > 0) ? -2.5*log(source->peak->flux) : NAN; + nImageOverlap = 1; + ID = 0; // XXX need to figure out how to generate this + + row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_S64, "IPP detection identifier index", ID); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", PAR[PM_PAR_XPOS]); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", PAR[PM_PAR_YPOS]); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", dPAR[PM_PAR_XPOS]); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", dPAR[PM_PAR_YPOS]); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", PS_MIN (99.0, source->psfMag)); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", PS_MIN (99.0, source->errMag)); + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", PS_MIN (99.0, peakMag)); + psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); + psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); + // XXX this is called STAR_GALAXY_SEP in the ICD; PSF_PROB is better + // XXX need to set this value in psphotEvalPSF + // XXX can I use the 2d polynomial peak fit to constrain this for the low-sn sources? + psMetadataAdd (row, PS_LIST_TAIL, "PSF_PROBABILITY", PS_DATA_F32, "Probability of PSF-ness", NAN); + // XXX these should be major and minor, not 'x' and 'y' + psMetadataAdd (row, PS_LIST_TAIL, "PSF_WIDTH_X", PS_DATA_F32, "PSF width in x coordinate", axes.major); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_WIDTH_Y", PS_DATA_F32, "PSF width in y coordinate", axes.minor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor", source->pixWeight); + // XXX not sure how to get this : need to load Nimages with weight + psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_S32, "Number of frames overlapping source center", nImageOverlap); + + // XXX these calibrated values are not supplied by psphot analysis + // psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "RA from PSF fit", RA); + // psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "DEC from PSF fit", DEC); + // psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF_SIG", PS_DATA_F32, "Sigma of PSF fit RA", dRA); + // psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF_SIG", PS_DATA_F32, "Sigma of PSF fit DEC", dDEC); + // psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "Calibrated magnitude", calMag); + // psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_sIG", PS_DATA_F32, "Sigma of calibrated magnitude", calMagErr); + + psArrayAdd (table, 100, row); + psFree (row); + } + + if (table->n == 0) { + psFitsWriteBlank (fits, header, extname); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable (fits, header, table, extname)) { + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree(table); + return false; + } + + psFree (table); + return true; +} + +// read in a readout from the fits file +psArray *pmSourcesRead_PS1_DEV_0 (psFits *fits, psMetadata *header) +{ + + bool status; + psF32 *PAR, *dPAR; + psEllipseShape shape; + psEllipseAxes axes; + float lsky; + + // define PSF model type + int modelType = pmModelSetType ("PS_MODEL_GAUSS"); + + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); + if (PSF_NAME != NULL) { + modelType = pmModelSetType (PSF_NAME); + } + + // find config information for output header + float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); + if (!status) + ZERO_POINT = 25.0; + + psArray *table = psFitsReadTable (fits); + // validate a single row of the table (must match SMP) + + psArray *sources = psArrayAlloc (table->n); + + // convert the table to the pmSource entries + // XXX need to chooose PSF vs EXT, based on type? + for (int i = 0; i < table->n; i++) { + pmSource *source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + source->type = PM_SOURCE_TYPE_STAR; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + psMetadata *row = table->data[i]; + + PAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF"); + PAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF"); + dPAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF_SIG"); + dPAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF_SIG"); + axes.major = psMetadataLookupF32 (&status, row, "PSF_WIDTH_X"); + axes.minor = psMetadataLookupF32 (&status, row, "PSF_WIDTH_Y"); + axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); + + PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); + dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); + + // XXX use these to determine PAR[PM_PAR_I0]? + source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); + source->errMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_ERR"); + + shape = psEllipseAxesToShape (axes); + + PAR[PM_PAR_SXX] = shape.sx; + PAR[PM_PAR_SYY] = shape.sy; + PAR[PM_PAR_SXY] = shape.sxy; + + // XXX other values saved but not loaded? + // psMetadataLookupS64 (&status, row, "IPP_IDET"); + // psMetadataLookupF32 (&status, row, "PSF_QF"); + // psMetadataLookupF32 (&status, row, "PSF_PROBABILITY"); + // psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); + // psMetadataLookupF32 (&status, row, "N_FRAMES"); + + sources->data[i] = source; + } + psFree (table); + return (sources); +} Index: /trunk/psModules/src/objects/pmSourceIO_SMPDATA.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 12402) +++ /trunk/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 12402) @@ -0,0 +1,171 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-03-11 18:56:38 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceIO.h" + +// elixir-style FITS table output (header + table in 1st extension) +// this format consists of a header derived from the image header +// followed by a zero-size matrix, followed by the table data +bool pmSourcesWrite_SMPDATA (psFits *fits, psArray *sources, psMetadata *header, char *extname) +{ + + psArray *table; + psMetadata *row; + int i, type; + psF32 *PAR, *dPAR; + float lsky; + bool status; + psEllipseShape shape; + psEllipseAxes axes; + + // find config information for output header + float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); + + table = psArrayAllocEmpty (sources->n); + + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + pmModel *model = pmSourceSelectModel (source); + if (model == NULL) + continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + type = pmSourceDophotType (source); + lsky = (source->sky < 1.0) ? 0.0 : log10(source->sky); + + shape.sx = PAR[PM_PAR_SXX]; + shape.sy = PAR[PM_PAR_SYY]; + shape.sxy = PAR[PM_PAR_SXY]; + axes = psEllipseShapeToAxes (shape, 20.0); + + row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "X_PIX", PS_DATA_F32, "", PAR[PM_PAR_XPOS]); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PIX", PS_DATA_F32, "", PAR[PM_PAR_YPOS]); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_RAW", PS_DATA_F32, "", PS_MIN (99.0, source->psfMag + ZERO_POINT)); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_ERR", PS_DATA_F32, "", PS_MIN (999, 1000*source->errMag)); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_GAL", PS_DATA_F32, "", PS_MIN (99.0, source->extMag + ZERO_POINT)); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_AP", PS_DATA_F32, "", PS_MIN (99.0, source->apMag + ZERO_POINT)); + psMetadataAdd (row, PS_LIST_TAIL, "LOG_SKY", PS_DATA_F32, "", lsky); + psMetadataAdd (row, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32, "", axes.major); + psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32, "", axes.minor); + psMetadataAdd (row, PS_LIST_TAIL, "THETA", PS_DATA_F32, "", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT", PS_DATA_U8, "", type); + psMetadataAdd (row, PS_LIST_TAIL, "WEIGHT", PS_DATA_U8, "", PS_MIN (255, PS_MAX(0, 255*source->pixWeight))); + psMetadataAdd (row, PS_LIST_TAIL, "DUMMY", PS_DATA_U16, "", 0); + + psArrayAdd (table, 100, row); + psFree (row); + } + + if (table->n == 0) { + psFitsWriteBlank (fits, header, extname); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable (fits, header, table, extname)) { + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree(table); + return false; + } + + psFree (table); + return true; +} + +// read in a readout from the fits file +psArray *pmSourcesRead_SMPDATA (psFits *fits, psMetadata *header) +{ + + bool status; + psF32 *PAR, *dPAR; + psEllipseShape shape; + psEllipseAxes axes; + float lsky; + + // define PSF model type + int modelType = pmModelSetType ("PS_MODEL_GAUSS"); + + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); + if (PSF_NAME != NULL) { + modelType = pmModelSetType (PSF_NAME); + } + + // find config information for output header + float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); + if (!status) + ZERO_POINT = 25.0; + + psArray *table = psFitsReadTable (fits); + // validate a single row of the table (must match SMP) + + psArray *sources = psArrayAlloc (table->n); + + // convert the table to the pmSource entries + // XXX need to chooose PSF vs EXT, based on type? + for (int i = 0; i < table->n; i++) { + pmSource *source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + source->type = PM_SOURCE_TYPE_STAR; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + psMetadata *row = table->data[i]; + + lsky = psMetadataLookupF32 (&status, row, "LOG_SKY"); + PAR[PM_PAR_SKY] = pow(lsky, 10.0); + PAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PIX"); + PAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PIX"); + axes.major = psMetadataLookupF32 (&status, row, "FWHM_X"); + axes.minor = psMetadataLookupF32 (&status, row, "FWHM_Y"); + axes.theta = psMetadataLookupF32 (&status, row, "THETA"); + + shape = psEllipseAxesToShape (axes); + + PAR[PM_PAR_SXX] = shape.sx; + PAR[PM_PAR_SYY] = shape.sy; + PAR[PM_PAR_SXY] = shape.sxy; + + source->psfMag = psMetadataLookupF32 (&status, row, "MAG_RAW") - ZERO_POINT; + source->extMag = psMetadataLookupF32 (&status, row, "MAG_GAL") - ZERO_POINT; + source->errMag = psMetadataLookupF32 (&status, row, "MAG_ERR"); + source->apMag = psMetadataLookupF32 (&status, row, "MAG_AP") - ZERO_POINT; + + sources->data[i] = source; + } + psFree (table); + return (sources); +}