Index: anches/eam_branches/ipp-20120627/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c =================================================================== --- /branches/eam_branches/ipp-20120627/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 34185) +++ (revision ) @@ -1,618 +1,0 @@ -/** @file pmSourceIO.c - * - * @author EAM, IfA - * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-18 02:44:19 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include - -#include "pmConfig.h" -#include "pmDetrendDB.h" - -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmFPALevel.h" -#include "pmFPAview.h" -#include "pmFPAfile.h" - -#include "pmTrend2D.h" -#include "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmSpan.h" -#include "pmFootprintSpans.h" -#include "pmFootprint.h" -#include "pmPeaks.h" -#include "pmMoments.h" -#include "pmModelFuncs.h" -#include "pmModel.h" -#include "pmModelUtils.h" -#include "pmModelClass.h" -#include "pmSourceMasks.h" -#include "pmSourceExtendedPars.h" -#include "pmSourceDiffStats.h" -#include "pmSource.h" -#include "pmSourceFitModel.h" -#include "pmPSF.h" -#include "pmPSFtry.h" - -#include "pmSourceIO.h" -#include "pmSourceOutputs.h" - -// panstarrs-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 version has elements intended for difference images & forced photometry: -// diffStats entries (good for dipoles); flux + flux error (for insignificant detections) -bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(sources, false); - PS_ASSERT_PTR_NON_NULL(extname, false); - - psArray *table; - psMetadata *row; - - pmChip *chip = readout->parent->parent; - - - // if the sequence is defined, write these in seq order; otherwise - // write them in S/N order: - if (sources->n > 0) { - pmSource *source = (pmSource *) sources->data[0]; - if (source->seq == -1) { - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - } else { - sources = psArraySort (sources, pmSourceSortBySeq); - } - } - - float magOffset; - float zeroptErr; - float fwhmMajor; - float fwhmMinor; - pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); - - table = psArrayAllocEmpty (sources->n); - - // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state - // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. - for (int i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - - // If source->seq is -1, source was generated in this analysis. If source->seq is - // not -1, source was read from elsewhere: in the latter case, preserve the source - // ID. source.seq is used instead of source.id since the latter is a const - // generated on Alloc, and would thus be wrong for read in sources. - if (source->seq == -1) { - source->seq = i; - } - - // set the 'best' values for various output fields: - pmSourceOutputs outputs; - pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); - - pmSourceOutputsMoments moments; - pmSourceOutputsSetMoments (&moments, source); - - pmSourceDiffStats diffStats; - pmSourceDiffStatsInit(&diffStats); - if (source->diffStats) { - diffStats = *source->diffStats; - } - - row = psMetadataAlloc (); - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfFlux); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfFluxErr); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); - 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); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); - psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor", source->pixWeightNotBad); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); - - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); - - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - - psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", source->nFrames); - psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); - - psArrayAdd (table, 100, row); - psFree (row); - - // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not - // subtracted - - // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image - // edge; 3) any pixels in the 3x3 peak region are masked; - } - - psMetadata *header = psMetadataCopy(NULL, tableHeader); - pmSourceMasksHeader(header); - - if (table->n == 0) { - psFitsWriteBlank(fits, header, extname); - psFree(table); - psFree(header); - 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); - psFree(header); - return false; - } - psFree(table); - psFree(header); - - return true; -} - -// read in a readout from the fits file -psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(header, false); - - bool status; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - - // define PSF model type - // XXX need to carry the extra model parameters - int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); - - char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); - if (PSF_NAME != NULL) { - modelType = pmModelClassGetType (PSF_NAME); - } - assert (modelType > -1); - - // We get the size of the table, and allocate the array of sources first because the table - // is large and ephemeral --- when the table gets blown away, whatever is allocated after - // the table is read blocks the free. In fact, it's better to read the table row by row. - long numSources = psFitsTableSize(fits); // Number of sources in table - psArray *sources = psArrayAlloc(numSources); // Array of sources, to return - - // convert the table to the pmSource entriesa - for (int i = 0; i < numSources; i++) { - psMetadata *row = psFitsReadTableRow(fits, i); // Table row - - pmSource *source = pmSourceAlloc (); - pmModel *model = pmModelAlloc (modelType); - source->modelPSF = model; - source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags - - // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - source->seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); - 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_MAJOR"); - axes.minor = psMetadataLookupF32 (&status, row, "PSF_MINOR"); - axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); - - PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); - dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); - source->sky = PAR[PM_PAR_SKY]; - source->skyErr = dPAR[PM_PAR_SKY]; - - // XXX use these to determine PAR[PM_PAR_I0]? - source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); - source->psfFluxErr = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); - - source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); - source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); - source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); - - // XXX this scaling is incorrect: does not include the 2 \pi AREA factor - PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; - dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - - pmPSF_AxesToModel (PAR, axes, modelType); - - float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); - float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; - - // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) - source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->rawFlux = peakFlux; - source->peak->smoothFlux = peakFlux; - source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position - source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position - source->peak->dx = dPAR[PM_PAR_XPOS]; - source->peak->dy = dPAR[PM_PAR_YPOS]; - - source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); - source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); - source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); - source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); - - // note that some older versions used PSF_PROBABILITY: this was not well defined. - model->chisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); - model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); - model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); - - source->moments = pmMomentsAlloc (); - source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf - source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf - - source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); - source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); - source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); - - source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); - assert (status); - - int nPos = psMetadataLookupS32 (&status, row, "DIFF_NPOS"); - if (nPos) { - source->diffStats = pmSourceDiffStatsAlloc(); - source->diffStats->nGood = nPos; - source->diffStats->fRatio = psMetadataLookupF32 (&status, row, "DIFF_FRATIO"); - source->diffStats->nRatioBad = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_BAD"); - source->diffStats->nRatioMask = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_MASK"); - source->diffStats->nRatioAll = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_ALL"); - } - - sources->data[i] = source; - psFree(row); - } - - return sources; -} - -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - - bool status; - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psF32 xPos, yPos; - psF32 xErr, yErr; - - // 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, "xsrc table extension", extname); - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - 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 - for (int i = 0; i < sources->n; i++) { - // skip source if it is not a ext sourc - // XXX we have two places that extended source parameters are measured: - // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models, - // psphotFitEXT, which does the simple extended source model fit (not psf-convolved) - // should we require both? - - pmSource *source = sources->data[i]; - - // skip sources without measurements - if (source->extpars == NULL) continue; - - // we require a PSF model fit (ignore the real crud) - pmModel *model = source->modelPSF; - if (model == NULL) continue; - - // XXX I need to split the extended models from the extended aperture measurements - 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]; - - row = psMetadataAlloc (); - - // 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); - -# if (0) - // Petrosian measurements - // 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 - 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 - 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); - } - } - - // 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); - } - } - } - -# endif - psArrayAdd (table, 100, row); - psFree (row); - } - - if (table->n == 0) { - 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 (outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - - return true; -} - -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_DV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) -{ - - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - 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, "xsrc 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++) { - 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); - - // we write out all sources, regardless of quality. the source flags tell us the state - for (int i = 0; i < sources->n; i++) { - - pmSource *source = sources->data[i]; - - // 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); - - // 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]; - - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - 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' - 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 < 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) { - 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); - } - } - - // XXX other parameters which may be set. - // XXX flag / value to define the model - // XXX write out the model type, fit status flags - - psArrayAdd (table, 100, row); - psFree (row); - } - } - - if (table->n == 0) { - 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 (outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; -} - -bool pmSourcesWrite_CMF_PS1_DV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - return true; -} Index: anches/eam_branches/ipp-20120627/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c =================================================================== --- /branches/eam_branches/ipp-20120627/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c (revision 34185) +++ (revision ) @@ -1,644 +1,0 @@ -/** @file pmSourceIO.c - * - * @author EAM, IfA - * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-18 02:44:19 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include - -#include "pmConfig.h" -#include "pmDetrendDB.h" - -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmFPALevel.h" -#include "pmFPAview.h" -#include "pmFPAfile.h" - -#include "pmTrend2D.h" -#include "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmSpan.h" -#include "pmFootprintSpans.h" -#include "pmFootprint.h" -#include "pmPeaks.h" -#include "pmMoments.h" -#include "pmModelFuncs.h" -#include "pmModel.h" -#include "pmModelUtils.h" -#include "pmModelClass.h" -#include "pmSourceMasks.h" -#include "pmSourceExtendedPars.h" -#include "pmSourceDiffStats.h" -#include "pmSource.h" -#include "pmSourceFitModel.h" -#include "pmPSF.h" -#include "pmPSFtry.h" - -#include "pmSourceIO.h" -#include "pmSourceOutputs.h" - -// panstarrs-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 version has elements intended for difference images & forced photometry: -// diffStats entries (good for dipoles); flux + flux error (for insignificant detections) -bool pmSourcesWrite_CMF_PS1_DV2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(sources, false); - PS_ASSERT_PTR_NON_NULL(extname, false); - - psArray *table; - psMetadata *row; - - pmChip *chip = readout->parent->parent; - - // if the sequence is defined, write these in seq order; otherwise - // write them in S/N order: - if (sources->n > 0) { - pmSource *source = (pmSource *) sources->data[0]; - if (source->seq == -1) { - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - } else { - sources = psArraySort (sources, pmSourceSortBySeq); - } - } - - table = psArrayAllocEmpty (sources->n); - - float magOffset; - float zeroptErr; - float fwhmMajor; - float fwhmMinor; - pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); - - // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state - // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. - for (int i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - - // If source->seq is -1, source was generated in this analysis. If source->seq is - // not -1, source was read from elsewhere: in the latter case, preserve the source - // ID. source.seq is used instead of source.id since the latter is a const - // generated on Alloc, and would thus be wrong for read in sources. - if (source->seq == -1) { - source->seq = i; - } - - // set the 'best' values for various output fields: - pmSourceOutputs outputs; - pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); - - pmSourceOutputsMoments moments; - pmSourceOutputsSetMoments (&moments, source); - - pmSourceDiffStats diffStats; - pmSourceDiffStatsInit(&diffStats); - if (source->diffStats) { - diffStats = *source->diffStats; - } - - row = psMetadataAlloc (); - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfFlux); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfFluxErr); - - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in real aperture", source->apMagRaw); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", source->apFlux); - psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", source->apFluxErr); - - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); - 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); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); - psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", source->pixWeightNotPoor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 1.0 R1)", moments.Kinner); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 4.0 R1)", moments.Kouter); - - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); - - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P", PS_DATA_F32, "distance to positive match source", diffStats.Rp); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P", PS_DATA_F32, "signal-to-noise of pos match src", diffStats.SNp); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M", PS_DATA_F32, "distance to negative match source", diffStats.Rm); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M", PS_DATA_F32, "signal-to-noise of neg match src", diffStats.SNm); - - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags (group 2)", source->mode2); - - // XXX not sure how to get this : need to load Nimages with weight? - psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", source->nFrames); - psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); - - psArrayAdd (table, 100, row); - psFree (row); - - // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not - // subtracted - - // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image - // edge; 3) any pixels in the 3x3 peak region are masked; - } - - psMetadata *header = psMetadataCopy(NULL, tableHeader); - pmSourceMasksHeader(header); - - if (table->n == 0) { - psFitsWriteBlank(fits, header, extname); - psFree(table); - psFree(header); - 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); - psFree(header); - return false; - } - psFree(table); - psFree(header); - - return true; -} - -// read in a readout from the fits file -psArray *pmSourcesRead_CMF_PS1_DV2 (psFits *fits, psMetadata *header) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(header, false); - - bool status; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - - // define PSF model type - // XXX need to carry the extra model parameters - int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); - - char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); - if (PSF_NAME != NULL) { - modelType = pmModelClassGetType (PSF_NAME); - } - assert (modelType > -1); - - // We get the size of the table, and allocate the array of sources first because the table - // is large and ephemeral --- when the table gets blown away, whatever is allocated after - // the table is read blocks the free. In fact, it's better to read the table row by row. - long numSources = psFitsTableSize(fits); // Number of sources in table - psArray *sources = psArrayAlloc(numSources); // Array of sources, to return - - // convert the table to the pmSource entriesa - for (int i = 0; i < numSources; i++) { - psMetadata *row = psFitsReadTableRow(fits, i); // Table row - - pmSource *source = pmSourceAlloc (); - pmModel *model = pmModelAlloc (modelType); - source->modelPSF = model; - source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags - - // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - source->seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); - 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_MAJOR"); - axes.minor = psMetadataLookupF32 (&status, row, "PSF_MINOR"); - axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); - - PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); - dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); - source->sky = PAR[PM_PAR_SKY]; - source->skyErr = dPAR[PM_PAR_SKY]; - - // XXX use these to determine PAR[PM_PAR_I0]? - source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); - source->psfFluxErr = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); - - source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); - source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); - source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); - - // XXX this scaling is incorrect: does not include the 2 \pi AREA factor - PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; - dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - - pmPSF_AxesToModel (PAR, axes, modelType); - - float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); - float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; - - // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) - source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->rawFlux = peakFlux; - source->peak->smoothFlux = peakFlux; - source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position - source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position - source->peak->dx = dPAR[PM_PAR_XPOS]; - source->peak->dy = dPAR[PM_PAR_YPOS]; - - source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); - source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT"); - source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); - source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); - source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); - - // note that some older versions used PSF_PROBABILITY: this was not well defined. - model->chisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); - model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); - model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); - - source->moments = pmMomentsAlloc (); - source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf - source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf - - source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); - source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); - source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); - - source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); - source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2"); - assert (status); - - int nPos = psMetadataLookupS32 (&status, row, "DIFF_NPOS"); - if (nPos) { - source->diffStats = pmSourceDiffStatsAlloc(); - source->diffStats->nGood = nPos; - source->diffStats->fRatio = psMetadataLookupF32 (&status, row, "DIFF_FRATIO"); - source->diffStats->nRatioBad = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_BAD"); - source->diffStats->nRatioMask = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_MASK"); - source->diffStats->nRatioAll = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_ALL"); - - source->diffStats->Rp = psMetadataLookupF32 (&status, row, "DIFF_R_P"); - source->diffStats->SNp = psMetadataLookupF32 (&status, row, "DIFF_SN_P"); - source->diffStats->Rm = psMetadataLookupF32 (&status, row, "DIFF_R_M"); - source->diffStats->SNm = psMetadataLookupF32 (&status, row, "DIFF_SN_M"); - } - - sources->data[i] = source; - psFree(row); - } - - return sources; -} - -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_DV2_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - - bool status; - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psF32 xPos, yPos; - psF32 xErr, yErr; - - // 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, "xsrc table extension", extname); - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - 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 - for (int i = 0; i < sources->n; i++) { - // skip source if it is not a ext sourc - // XXX we have two places that extended source parameters are measured: - // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models, - // psphotFitEXT, which does the simple extended source model fit (not psf-convolved) - // should we require both? - - pmSource *source = sources->data[i]; - - // skip sources without measurements - if (source->extpars == NULL) continue; - - // we require a PSF model fit (ignore the real crud) - pmModel *model = source->modelPSF; - if (model == NULL) continue; - - // XXX I need to split the extended models from the extended aperture measurements - 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]; - - row = psMetadataAlloc (); - - // 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); - -# if (0) - // Petrosian measurements - // 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 - 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 - 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); - } - } - - // 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); - } - } - } - -# endif - psArrayAdd (table, 100, row); - psFree (row); - } - - if (table->n == 0) { - 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 (outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - - return true; -} - -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_DV2_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) -{ - - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - 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, "xsrc 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++) { - 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); - - // we write out all sources, regardless of quality. the source flags tell us the state - for (int i = 0; i < sources->n; i++) { - - pmSource *source = sources->data[i]; - - // 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); - - // 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]; - - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - 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' - 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 < 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) { - 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); - } - } - - // XXX other parameters which may be set. - // XXX flag / value to define the model - // XXX write out the model type, fit status flags - - psArrayAdd (table, 100, row); - psFree (row); - } - } - - if (table->n == 0) { - 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 (outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; -} - -bool pmSourcesWrite_CMF_PS1_DV2_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - return true; -} Index: anches/eam_branches/ipp-20120627/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c =================================================================== --- /branches/eam_branches/ipp-20120627/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 34185) +++ (revision ) @@ -1,1122 +1,0 @@ -/** @file pmSourceIO.c - * - * @author EAM, IfA - * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-18 02:44:19 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include - -#include "pmConfig.h" -#include "pmErrorCodes.h" -#include "pmDetrendDB.h" - -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmFPALevel.h" -#include "pmFPAview.h" -#include "pmFPAfile.h" - -#include "pmTrend2D.h" -#include "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmSpan.h" -#include "pmFootprintSpans.h" -#include "pmFootprint.h" -#include "pmPeaks.h" -#include "pmMoments.h" -#include "pmModelFuncs.h" -#include "pmModel.h" -#include "pmModelUtils.h" -#include "pmModelClass.h" -#include "pmSourceMasks.h" -#include "pmSourceExtendedPars.h" -#include "pmSourceDiffStats.h" -#include "pmSource.h" -#include "pmSourceFitModel.h" -#include "pmPSF.h" -#include "pmPSFtry.h" - -#include "pmSourceIO.h" -#include "pmSourceOutputs.h" - -// panstarrs-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 - -// NOTE: this output function is intended for psphotStack analysis: it includes per-psf radial fluxes -// XXX currently, the 'read' function is NOT consistent with the 'write' function (does not read radial fluxes) - -bool pmSourcesWrite_CMF_PS1_SV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(sources, false); - PS_ASSERT_PTR_NON_NULL(extname, false); - - psArray *table; - psMetadata *row; - - pmChip *chip = readout->parent->parent; - - // if the sequence is defined, write these in seq order; otherwise write them in S/N order. - // Careful: if we are working with child sources, then we need to sort by the parent info, - // not our info - if (sources->n > 0) { - pmSource *source = sources->data[0]; - if (source->seq == -1) { - sources = psArraySort (sources, pmSourceSortByFlux); - } else { - sources = psArraySort (sources, pmSourceSortBySeq); - } - } - - table = psArrayAllocEmpty (sources->n); - - float magOffset; - float zeroptErr; - float fwhmMajor; - float fwhmMinor; - pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); - - // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state - // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. - 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->seq is -1, source is unique and generated in this analysis. If - // source->seq is not -1, source was read from elsewhere or tied to other source (eg - // from another image): in the latter case, preserve the source ID. source.seq is used - // instead of source.id since the latter is a const generated on Alloc, and would thus - // be wrong for read in sources. - if (source->seq == -1) { - source->seq = i; - } - - // set the 'best' values for various output fields: - pmSourceOutputs outputs; - pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); - - pmSourceOutputsMoments moments; - pmSourceOutputsSetMoments (&moments, source); - - row = psMetadataAlloc (); - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux", source->psfFluxErr); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in reported aperture", source->apMagRaw); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); - 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); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); - psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", source->pixWeightNotPoor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kinner); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kouter); - - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); - - psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", source->nFrames); - psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); - - psArrayAdd (table, 100, row); - psFree (row); - - // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not - // subtracted - - // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image - // edge; 3) any pixels in the 3x3 peak region are masked; - } - - // XXX why do we make a copy here to be supplemented with the masks? why not do this in the calling function? - psMetadata *header = psMetadataCopy(NULL, tableHeader); - pmSourceMasksHeader(header); - - if (table->n == 0) { - if (!psFitsWriteBlank(fits, header, extname)) { - psError(psErrorCodeLast(), false, "Unable to write blank sources file."); - psFree(table); - psFree(header); - return false; - } - psFree(table); - psFree(header); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable(fits, header, table, extname)) { - psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); - psFree(table); - psFree(header); - return false; - } - psFree(table); - psFree(header); - - return true; -} - -// read in a readout from the fits file -psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(header, false); - - bool status; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - - // define PSF model type - // XXX need to carry the extra model parameters - int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); - - char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); - if (PSF_NAME != NULL) { - modelType = pmModelClassGetType (PSF_NAME); - } - assert (modelType > -1); - - // We get the size of the table, and allocate the array of sources first because the table - // is large and ephemeral --- when the table gets blown away, whatever is allocated after - // the table is read blocks the free. In fact, it's better to read the table row by row. - long numSources = psFitsTableSize(fits); // Number of sources in table - psArray *sources = psArrayAlloc(numSources); // Array of sources, to return - - // convert the table to the pmSource entriesa - for (int i = 0; i < numSources; i++) { - psMetadata *row = psFitsReadTableRow(fits, i); // Table row - if (!row) { - psError(psErrorCodeLast(), false, "Unable to read row %d of sources", i); - psFree(sources); - return NULL; - } - - pmSource *source = pmSourceAlloc (); - pmModel *model = pmModelAlloc (modelType); - source->modelPSF = model; - source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags - - // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - source->seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); - 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_MAJOR"); - axes.minor = psMetadataLookupF32 (&status, row, "PSF_MINOR"); - axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); - - PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); - dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); - source->sky = PAR[PM_PAR_SKY]; - source->skyErr = dPAR[PM_PAR_SKY]; - - // XXX use these to determine PAR[PM_PAR_I0]? - source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); - source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); - source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); - - // XXX this scaling is incorrect: does not include the 2 \pi AREA factor - PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; - dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - - pmPSF_AxesToModel (PAR, axes, modelType); - - float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); - float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; - - // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) - source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->rawFlux = peakFlux; - source->peak->smoothFlux = peakFlux; - source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position - source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position - source->peak->dx = dPAR[PM_PAR_XPOS]; - source->peak->dy = dPAR[PM_PAR_YPOS]; - - source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); - source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT"); - source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); - source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); - source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); - - // note that some older versions used PSF_PROBABILITY: this was not well defined. - model->chisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); - model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); - model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); - - source->moments = pmMomentsAlloc (); - source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf - source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf - - source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); - source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); - source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); - - source->moments->Mrf = psMetadataLookupF32 (&status, row, "MOMENTS_R1"); - source->moments->Mrh = psMetadataLookupF32 (&status, row, "MOMENTS_RH"); - source->moments->KronFlux = psMetadataLookupF32 (&status, row, "KRON_FLUX"); - source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR"); - - source->moments->KronFinner = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER"); - source->moments->KronFouter = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER"); - - // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data, - // we are storing enough information so the output will be consistent with the input - source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C"); - source->moments->Mxxy = 0.0; - source->moments->Mxyy = 0.0; - source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S"); - - source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C"); - source->moments->Mxxxy = 0.0; - source->moments->Mxxyy = 0.0; - source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S"); - source->moments->Myyyy = 0.0; - - source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); - source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2"); - assert (status); - - sources->data[i] = source; - psFree(row); - } - - return sources; -} - -bool pmSourcesWrite_CMF_PS1_SV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - - bool status; - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psF32 xPos, yPos; - psF32 xErr, yErr; - int nRow = -1; - char keyword1[80], keyword2[80]; - - // 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, "xsrc table extension", extname); - - pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - // zero point corrections - bool status1 = false; - bool status2 = false; - float magOffset = 0.0; - float exptime = psMetadataLookupF32(&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - table = psArrayAllocEmpty (sources->n); - - // which extended source analyses should we perform? - bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - - psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); - psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); - psAssert (radMin->n == radMax->n, "inconsistent annular bins"); - - // write the radial profile apertures to header - for (int i = 0; i < radMax->n; i++) { - sprintf (keyword1, "RMIN_%02d", i); - sprintf (keyword2, "RMAX_%02d", i); - psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]); - psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]); - } - - // we write out all sources, regardless of quality. the source flags tell us the state - 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; - - // skip sources without measurements - if (source->extpars == NULL) continue; - - // we require a PSF model fit (ignore the real crud) - pmModel *model = source->modelPSF; - if (model == NULL) continue; - - // XXX I need to split the extended models from the extended aperture measurements - 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]; - - row = psMetadataAlloc (); - - // 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); - - float AxialRatio = NAN; - float AxialTheta = NAN; - pmSourceExtendedPars *extpars = source->extpars; - if (extpars) { - AxialRatio = extpars->axes.minor / extpars->axes.major; - AxialTheta = extpars->axes.theta; - } - psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO", PS_DATA_F32, "Axial Ratio of radial profile", AxialRatio); - psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA", PS_DATA_F32, "Angle of radial profile ellipse", AxialTheta); - - // Petrosian measurements - // XXX insert header data: petrosian ref radius, flux ratio - // XXX check flags to see if Pet was measured - if (doPetrosian) { - pmSourceExtendedPars *extpars = source->extpars; - if (extpars) { - // XXX note that this mag is either calibrated or instrumental depending on existence of zero point - float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point - float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", mag); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", magErr); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", extpars->petrosianR50); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)", extpars->petrosianR50Err); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL", PS_DATA_F32, "Petrosian Fill Factor", extpars->petrosianFill); - } 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); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)",NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL", PS_DATA_F32, "Petrosian Fill Factor", NAN); - } - } - - // Flux Annuli (if we have extended source measurements, we have these. only optionally save them) - if (doAnnuli) { - psVector *radSB = psVectorAlloc(radMin->n, PS_TYPE_F32); - psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32); - psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32); - psVectorInit (radSB, NAN); - psVectorInit (radFlux, NAN); - psVectorInit (radFill, NAN); - if (!source->extpars) goto empty_annuli; - if (!source->extpars->radProfile) goto empty_annuli; - if (!source->extpars->radProfile->binSB) goto empty_annuli; - psAssert (source->extpars->radProfile->binSum, "programming error"); - psAssert (source->extpars->radProfile->binFill, "programming error"); - psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths"); - psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths"); - psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths"); - - // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux - for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) { - radSB->data.F32[j] = source->extpars->radProfile->binSB->data.F32[j]; - radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j]; - radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j]; - } - - empty_annuli: - psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB); - psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); - psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); - psFree (radSB); - psFree (radFlux); - psFree (radFill); - } - if (nRow < 0) { - nRow = row->list->n; - } else { - psAssert (nRow == row->list->n, "inconsistent row lengths"); - } - 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; -} - -bool pmSourcesRead_CMF_PS1_SV1_XSRC(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(sources, false); - - bool status; - long numSources = psFitsTableSize(fits); // Number of sources in table - if (numSources == 0) { - psError(psErrorCodeLast(), false, "XSRC Table contains no entries\n"); - return false; - } - - // petrosian mags are not saved, we need to calculate fluxes. For this we need exptime and zero point - float zeropt = psMetadataLookupF32(&status, hduHeader, "FPA.ZP"); - float exptime = psMetadataLookupF32(&status, hduHeader, "EXPTIME"); - float magOffset = zeropt + 2.5*log10(exptime); - - for (long i = 0; i < numSources; i++) { - psMetadata *row = psFitsReadTableRow(fits, i); // Table row - if (!row) { - psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); - psFree(row); - return false; - } - // Find the source with this sequence number. - // XXX: I am assuming that sources is sorted in order of seq - long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); - pmSource *source = NULL; -#ifndef ASSUME_SORTED - long j = seq < sources->n ? seq : sources->n - 1; - for (; j >= 0; j--) { - source = sources->data[j]; - if (source->seq == seq) { - break; - } - } -#else - long j = sourceIndex[seq]; - psAssert(j >= 0 && j < sources->n, "invalid sourceIndex"); - source = sources->data[j]; -#endif - if (!source) { - psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); - psFree(row); - return false; - } - - if (!source->extpars) { - source->extpars = pmSourceExtendedParsAlloc (); - } - pmSourceExtendedPars *extpars = source->extpars; - - // Assume that X_EXT Y_EXT and sigmas match the psf src so skip - - // We don't have enough information to calculate the major and minor axis. Set major to 1. Should we scale this by - // psf size or something? - extpars->axes.major = 1.0; - extpars->axes.minor = extpars->axes.major * psMetadataLookupF32(&status, row, "F25_ARATIO"); - extpars->axes.theta = psMetadataLookupF32(&status, row, "F25_THETA"); - - float mag = psMetadataLookupF32(&status, row, "PETRO_MAG"); - float magErr = psMetadataLookupF32(&status, row, "PETRO_MAG_ERR"); - if (isfinite(mag)) { - extpars->petrosianFlux = pow(10., (magOffset - mag) / 2.5); - if (isfinite(magErr)) { - extpars->petrosianFluxErr = extpars->petrosianFlux / magErr; - } - } - - extpars->petrosianRadius = psMetadataLookupF32(&status, row, "PETRO_RADIUS"); - extpars->petrosianRadiusErr= psMetadataLookupF32(&status, row, "PETRO_RADIUS_ERR"); - extpars->petrosianR50 = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50"); - extpars->petrosianR50Err = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50_ERR"); - extpars->petrosianR90 = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90"); - extpars->petrosianR90Err = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90_ERR"); - extpars->petrosianFill = psMetadataLookupF32(&status, row, "PETRO_FILL"); - - psVector *radSB = psMetadataLookupVector(&status, row, "PROF_SB"); - psVector *radFlux = psMetadataLookupVector(&status, row, "PROF_FLUX"); - psVector *radFill = psMetadataLookupVector(&status, row, "PROF_FILL"); - - if (radSB && radSB->n > 0) { - extpars->radProfile = pmSourceRadialProfileAlloc(); - extpars->radProfile->binSB = psMemIncrRefCounter(radSB); - extpars->radProfile->binSum = psMemIncrRefCounter(radFlux); - extpars->radProfile->binFill = psMemIncrRefCounter(radFill); - } - - psFree(row); - } - - return true; -} - -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_SV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) -{ - - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - 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, "xsrc 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); - } - } - - 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); - - // 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]; - - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - 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' - 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]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", PAR[PM_PAR_THETA]); - } else { - 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); - } - - // 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; -} - -bool pmSourcesRead_CMF_PS1_SV1_XFIT(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(sources, false); - - bool status; - long numSources = psFitsTableSize(fits); // Number of sources in table - if (numSources == 0) { - psError(psErrorCodeLast(), false, "XFIT Table contains no entries\n"); - return false; - } - - for (long i = 0; i < numSources; i++) { - psMetadata *row = psFitsReadTableRow(fits, i); // Table row - if (!row) { - psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); - psFree(row); - return false; - } - // Find the source with this sequence number. - // XXX: I am assuming that sources is sorted in order of seq. - long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); - long j = seq < sources->n ? seq : sources->n - 1; - pmSource *source = NULL; - for (; j >= 0; j--) { - source = sources->data[j]; - if (source->seq == seq) { - break; - } - } - if (!source) { - psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); - psFree(row); - return false; - } - if (!source->modelFits) { - // XXX: where to find the number of models to expect? - source->modelFits = psArrayAllocEmpty(5); - } - psString modelName = psMetadataLookupStr(&status, row, "MODEL_TYPE"); - if (!modelName) { - psError(PS_ERR_UNKNOWN, true, "Failed to find model name for row %ld\n", i); - psFree(row); - return false; - } - pmModelType modelType = pmModelClassGetType(modelName); - if (modelType < 0) { - psError(PS_ERR_UNKNOWN, true, "Failed to find model type for %s\n", modelName); - psFree(row); - return false; - } - pmModel *model = pmModelAlloc(modelType); - - psF32 *PAR = model->params->data.F32; - psF32 *dPAR = model->dparams->data.F32; - - PAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT"); - PAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT"); - dPAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT_SIG"); - dPAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT_SIG"); - - model->mag = psMetadataLookupF32(&status, row, "EXT_INST_MAG"); - model->magErr = psMetadataLookupF32(&status, row, "EXT_INST_MAG_SIG"); - - psEllipseAxes axes; - axes.major = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ"); - axes.minor = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN"); - axes.theta = psMetadataLookupF32(&status, row, "EXT_THETA"); - if (!pmPSF_AxesToModel(PAR, axes, modelType)) { - // Do we need to fail here or can this happen? - psError(PS_ERR_UNKNOWN, false, "Failed to convert psf axes to model"); - psFree(model); - psFree(row); - return false; - } - // XXX: clean this up - if (model->params->n > 7) { - PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07"); - } - // read the covariance matrix - int nparams = model->params->n; - psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32); - for (int y = 0; y < nparams; y++) { - for (int x = 0; x < nparams; x++) { - char name[64]; - snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x); - covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name); - } - } - model->covar = covar; - - psArrayAdd(source->modelFits, 1, model); - psFree(model); - - psFree(row); - } - - return true; -} - -// **** write out the radial flux values for the sources for a given matched-PSF image -// **** how do we distinguish the matched-PSF images from the non-matched version -bool pmSourcesWrite_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - bool status = false; - psArray *table; - psMetadata *row; - psF32 xPos, yPos; - char keyword1[80], keyword2[80]; - - // perform full non-linear fits / extended source analysis? - if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) { - psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n"); - return true; - } - - // 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, "radial flux table extension", extname); - - // we use this just to define the output vectors (which must be present for all objects) - psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); - psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); - psAssert (radMax, "this must have been defined and tested earlier!"); - psAssert (radMax->n, "this must have been defined and tested earlier!"); - psAssert (radMin->n == radMax->n, "inconsistent annular bins"); - - // write the radial profile apertures to header - for (int i = 0; i < radMax->n; i++) { - sprintf (keyword1, "RMIN_%02d", i); - sprintf (keyword2, "RMAX_%02d", i); - psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]); - psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]); - } - - // the FWHM values are available if we measured a psf-matched convolved set - psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES"); - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - 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++) { - - // 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; - - // skip sources without radial aper measurements (or insufficient) - if (source->radialAper == NULL) continue; - - // psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set"); - - for (int entry = 0; entry < source->radialAper->n; entry++) { - - // choose the convolved EXT model, if available, otherwise the simple one - pmSourceRadialApertures *radialAper = source->radialAper->data[entry]; - assert (radialAper); - - if (pmSourcePositionUseMoments(source)) { - xPos = source->moments->Mx; - yPos = source->moments->My; - } else { - xPos = source->peak->xf; - yPos = source->peak->yf; - } - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER", 0, "Center of aperture measurements", xPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER", 0, "Center of aperture measurements", yPos); - if (fwhmValues) { - psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM", 0, "FWHM of matched PSF", fwhmValues->data.F32[entry]); - } else { - psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM", 0, "image is not FWHM-matched", NAN); - } - - // XXX if we have raw radial apertures, write them out here - psVector *radFlux = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVector *radFill = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVector *radFluxStdev = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVectorInit (radFlux, NAN); - psVectorInit (radFluxErr, NAN); - psVectorInit (radFill, NAN); - if (!radialAper->flux) goto write_annuli; - if (!radialAper->fill) goto write_annuli; - psAssert (radialAper->flux->n <= radFlux->n, "inconsistent vector lengths"); - psAssert (radialAper->fill->n <= radFlux->n, "inconsistent vector lengths"); - - // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux - for (int j = 0; j < radialAper->flux->n; j++) { - radFlux->data.F32[j] = radialAper->flux->data.F32[j]; - radFluxErr->data.F32[j] = radialAper->fluxErr->data.F32[j]; - radFluxStdev->data.F32[j] = radialAper->fluxStdev->data.F32[j]; - radFill->data.F32[j] = radialAper->fill->data.F32[j]; - } - - write_annuli: - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli", radFluxErr); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_DATA_VECTOR, "flux standard deviation", radFluxStdev); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); - psFree (radFlux); - psFree (radFluxErr); - psFree (radFluxStdev); - psFree (radFill); - - 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; -} - -bool pmSourcesRead_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psMetadata *hduHeader, psArray *sources, long *sourceIndex) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(sources, false); - - bool status; - long numSources = psFitsTableSize(fits); // Number of sources in table - if (numSources == 0) { - psError(psErrorCodeLast(), false, "XRAD Table contains no entries\n"); - return false; - } - - long seq_first = -1; - long seq_last = -1; - psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32); - long max_entries = -1; - long num_entries = -1; - - for (long i = 0; i < numSources; i++) { - psMetadata *row = psFitsReadTableRow(fits, i); // Table row - if (!row) { - psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); - psFree(row); - return false; - } - // Find the source with this sequence number. - // XXX: I am assuming that sources is sorted in order of seq. - long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); - long j = seq < sources->n ? seq : sources->n - 1; - pmSource *source = NULL; - for (; j >= 0; j--) { - source = sources->data[j]; - if (source->seq == seq) { - break; - } - } - if (!source) { - psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); - psFree(row); - return false; - } - if (seq_first == -1) { - seq_first = seq; - } - if (seq == seq_first) { - psF32 value = psMetadataLookupF32(&status, row, "PSF_FWHM"); - psVectorAppend(fwhmValues, value); - } - if (seq == seq_last) { - num_entries++; - } else { - num_entries = 1; - seq_last = seq; - } - if (num_entries > max_entries) { - max_entries = num_entries; - } - - if (!source->radialAper) { - // XXX: where to find the number of models to expect? - source->radialAper = psArrayAllocEmpty(5); - } - pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc(); - - radialAper->flux = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX")); - radialAper->fluxStdev = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_STDEV")); - radialAper->fluxErr = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_ERR")); - radialAper->fill = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FILL")); - - psArrayAdd(source->radialAper, 1, radialAper); - - psFree(radialAper); - psFree(row); - } - - // check for consistency between the length of fwhmValues and the maximum number of entries for each row - if (fwhmValues->n != max_entries) { - psError(PS_ERR_PROGRAMMING, true, "number of PSF_FWHM values found %ld does not match expected number: %ld\n", - fwhmValues->n, max_entries); - psAssert(0, "fixme"); - } - - if (!readout->analysis) { - readout->analysis = psMetadataAlloc(); - } - - psMetadataAddVector(readout->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues); - psFree(fwhmValues); - - return true; -}