Index: /branches/czw_branch/20100427/psModules/src/camera/pmFPAMosaic.c =================================================================== --- /branches/czw_branch/20100427/psModules/src/camera/pmFPAMosaic.c (revision 28029) +++ /branches/czw_branch/20100427/psModules/src/camera/pmFPAMosaic.c (revision 28030) @@ -595,6 +595,6 @@ if (cell->readouts->n > 1) { - psWarning("Cell contains more than one readout (%ld) --- mosaicking only the first.\n", - cell->readouts->n); + psWarning("Skipping video cell for mosaic.\n"); + return true; } Index: /branches/czw_branch/20100427/psModules/src/objects/Makefile.am =================================================================== --- /branches/czw_branch/20100427/psModules/src/objects/Makefile.am (revision 28029) +++ /branches/czw_branch/20100427/psModules/src/objects/Makefile.am (revision 28030) @@ -43,4 +43,5 @@ pmSourceIO_CMF_PS1_V1.c \ pmSourceIO_CMF_PS1_V2.c \ + pmSourceIO_CMF_PS1_SV1.c \ pmSourceIO_CMF_PS1_DV1.c \ pmSourceIO_MatchedRefs.c \ Index: /branches/czw_branch/20100427/psModules/src/objects/pmPhotObj.c =================================================================== --- /branches/czw_branch/20100427/psModules/src/objects/pmPhotObj.c (revision 28029) +++ /branches/czw_branch/20100427/psModules/src/objects/pmPhotObj.c (revision 28030) @@ -58,8 +58,26 @@ if (!object->sources) { object->sources = psArrayAllocEmpty(1); - object->x = source->peak->xf; - object->y = source->peak->yf; + object->x = source->peak->xf; + object->y = source->peak->yf; + object->SN = source->peak->SN; + } else { + object->SN = PS_MAX(object->SN, source->peak->SN); } psArrayAdd (object->sources, 1, source); return true; } + +// sort by SN (descending) +int pmPhotObjSortBySN (const void **a, const void **b) +{ + pmPhotObj *objA = *(pmPhotObj **)a; + pmPhotObj *objB = *(pmPhotObj **)b; + + psF32 fA = objA->SN; + psF32 fB = objB->SN; + + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (-1); + if (diff < FLT_EPSILON) return (+1); + return (0); +} Index: /branches/czw_branch/20100427/psModules/src/objects/pmPhotObj.h =================================================================== --- /branches/czw_branch/20100427/psModules/src/objects/pmPhotObj.h (revision 28029) +++ /branches/czw_branch/20100427/psModules/src/objects/pmPhotObj.h (revision 28030) @@ -39,8 +39,12 @@ float x; float y; + float SN; // max of peak->SN for all matched sources } pmPhotObj; +pmPhotObj *pmPhotObjAlloc(void); + bool pmPhotObjAddSource(pmPhotObj *object, pmSource *source); -pmPhotObj *pmPhotObjAlloc(void); + +int pmPhotObjSortBySN (const void **a, const void **b); /// @} Index: /branches/czw_branch/20100427/psModules/src/objects/pmSource.c =================================================================== --- /branches/czw_branch/20100427/psModules/src/objects/pmSource.c (revision 28029) +++ /branches/czw_branch/20100427/psModules/src/objects/pmSource.c (revision 28030) @@ -54,4 +54,5 @@ psFree(tmp->moments); psFree(tmp->diffStats); + psFree(tmp->radial); psFree(tmp->blends); psTrace("psModules.objects", 10, "---- end ----\n"); @@ -116,4 +117,5 @@ source->extpars = NULL; source->diffStats = NULL; + source->radial = NULL; source->region = psRegionSet(NAN, NAN, NAN, NAN); Index: /branches/czw_branch/20100427/psModules/src/objects/pmSource.h =================================================================== --- /branches/czw_branch/20100427/psModules/src/objects/pmSource.h (revision 28029) +++ /branches/czw_branch/20100427/psModules/src/objects/pmSource.h (revision 28030) @@ -91,4 +91,5 @@ pmSourceExtendedPars *extpars; ///< extended source parameters pmSourceDiffStats *diffStats; ///< extra parameters for difference detections + pmSourceRadialApertures *radial; ///< radial flux in circular apertures int imageID; }; Index: /branches/czw_branch/20100427/psModules/src/objects/pmSourceExtendedPars.c =================================================================== --- /branches/czw_branch/20100427/psModules/src/objects/pmSourceExtendedPars.c (revision 28029) +++ /branches/czw_branch/20100427/psModules/src/objects/pmSourceExtendedPars.c (revision 28030) @@ -65,4 +65,30 @@ } +// pmSourceRadialApertures carries the raw radial flux information, including angular bins +static void pmSourceRadialAperturesFree(pmSourceRadialApertures *radial) +{ + if (!radial) return; + psFree(radial->flux); + psFree(radial->fluxErr); + psFree(radial->fill); +} + +pmSourceRadialApertures *pmSourceRadialAperturesAlloc() +{ + pmSourceRadialApertures *radial = (pmSourceRadialApertures *)psAlloc(sizeof(pmSourceRadialApertures)); + psMemSetDeallocator(radial, (psFreeFunc) pmSourceRadialAperturesFree); + + radial->flux = NULL; + radial->fluxErr = NULL; + radial->fill = NULL; + return radial; +} + +bool psMemCheckSourceRadialApertures(psPtr ptr) +{ + PS_ASSERT_PTR(ptr, false); + return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceRadialAperturesFree); +} + // pmSourceEllipticalFlux carries the elliptical renormalized radial flux info static void pmSourceEllipticalFluxFree(pmSourceEllipticalFlux *flux) Index: /branches/czw_branch/20100427/psModules/src/objects/pmSourceExtendedPars.h =================================================================== --- /branches/czw_branch/20100427/psModules/src/objects/pmSourceExtendedPars.h (revision 28029) +++ /branches/czw_branch/20100427/psModules/src/objects/pmSourceExtendedPars.h (revision 28030) @@ -20,4 +20,10 @@ psVector *isophotalRadii; // isophotal radius for the above angles } pmSourceRadialFlux; + +typedef struct { + psVector *flux; // fluxes measured at above radii + psVector *fluxErr; // formal error on the fluxes (sqrt\sum(variance)) + psVector *fill; // angles corresponding to above radial profiles +} pmSourceRadialApertures; typedef struct { @@ -62,4 +68,7 @@ bool psMemCheckSourceRadialFlux(psPtr ptr); +pmSourceRadialApertures *pmSourceRadialAperturesAlloc(); +bool psMemCheckSourceRadialApertures(psPtr ptr); + pmSourceEllipticalFlux *pmSourceEllipticalFluxAlloc(); bool psMemCheckSourceEllipticalFlux(psPtr ptr); Index: /branches/czw_branch/20100427/psModules/src/objects/pmSourceIO.c =================================================================== --- /branches/czw_branch/20100427/psModules/src/objects/pmSourceIO.c (revision 28029) +++ /branches/czw_branch/20100427/psModules/src/objects/pmSourceIO.c (revision 28030) @@ -539,4 +539,7 @@ status &= pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname); } + if (!strcmp (exttype, "PS1_SV1")) { + status &= pmSourcesWrite_CMF_PS1_SV1 (file->fits, readout, sources, file->header, outhead, dataname, recipe); + } if (!strcmp (exttype, "PS1_DV1")) { status &= pmSourcesWrite_CMF_PS1_DV1 (file->fits, readout, sources, file->header, outhead, dataname); @@ -556,4 +559,7 @@ status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); } + if (!strcmp (exttype, "PS1_SV1")) { + status &= pmSourcesWrite_CMF_PS1_SV1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); + } if (!strcmp (exttype, "PS1_DV1")) { status &= pmSourcesWrite_CMF_PS1_DV1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); @@ -573,4 +579,7 @@ status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, readout, sources, xfitname); } + if (!strcmp (exttype, "PS1_SV1")) { + status &= pmSourcesWrite_CMF_PS1_SV1_XFIT (file->fits, readout, sources, xfitname); + } if (!strcmp (exttype, "PS1_DV1")) { status &= pmSourcesWrite_CMF_PS1_DV1_XFIT (file->fits, readout, sources, xfitname); @@ -981,5 +990,5 @@ // we need to find the corresponding table EXTNAME. // first check the header - char *extdata = psMetadataLookupStr (NULL, hdu->header, "EXTDATA"); + char *extdata = psMetadataLookupStr (&status, hdu->header, "EXTDATA"); if (extdata) { // if EXTDATA is defined in the header, use that value for 'dataname' @@ -1023,4 +1032,7 @@ if (!strcmp (exttype, "PS1_V2")) { sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header); + } + if (!strcmp (exttype, "PS1_SV1")) { + sources = pmSourcesRead_CMF_PS1_SV1 (file->fits, hdu->header); } if (!strcmp (exttype, "PS1_DV1")) { Index: /branches/czw_branch/20100427/psModules/src/objects/pmSourceIO.h =================================================================== --- /branches/czw_branch/20100427/psModules/src/objects/pmSourceIO.h (revision 28029) +++ /branches/czw_branch/20100427/psModules/src/objects/pmSourceIO.h (revision 28030) @@ -43,4 +43,8 @@ bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname); +bool pmSourcesWrite_CMF_PS1_SV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_SV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_SV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname); + bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); @@ -57,4 +61,5 @@ psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header); psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header); +psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header); psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header); Index: /branches/czw_branch/20100427/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c =================================================================== --- /branches/czw_branch/20100427/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 28029) +++ /branches/czw_branch/20100427/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 28030) @@ -175,6 +175,6 @@ 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", errMag); - 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, "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", apRadius); Index: /branches/czw_branch/20100427/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c =================================================================== --- /branches/czw_branch/20100427/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 28030) +++ /branches/czw_branch/20100427/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 28030) @@ -0,0 +1,745 @@ +/** @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 "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmResiduals.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmModelClass.h" +#include "pmSourceIO.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 in 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; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + psF32 xPos, yPos; + psF32 xErr, yErr; + psF32 errMag, chisq, apRadius; + psS32 nPix, nDOF; + + pmChip *chip = readout->parent->parent; + pmFPA *fpa = chip->parent; + + bool status1 = false; + bool status2 = false; + float magOffset = NAN; + 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); + } + float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); + + // 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, pmSourceSortBySN); + } else { + sources = psArraySort (sources, pmSourceSortBySeq); + } + } + + table = psArrayAllocEmpty (sources->n); + + // we use this just to define the output vectors (which must be present for all objects) + bool status = false; + 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!"); + + // 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; + } + + // no difference between PSF and non-PSF model + pmModel *model = source->modelPSF; + + if (model != NULL) { + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + } else { + // in linear-fit mode, there is no error on the centroid + xErr = source->peak->dx; + yErr = source->peak->dy; + } + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { + axes = pmPSF_ModelToAxes (PAR, 20.0); + } else { + axes.major = NAN; + axes.minor = NAN; + axes.theta = NAN; + } + chisq = model->chisq; + nDOF = model->nDOF; + nPix = model->nPix; + apRadius = source->apRadius; + errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; + } else { + xPos = source->peak->xf; + yPos = source->peak->yf; + xErr = source->peak->dx; + yErr = source->peak->dy; + axes.major = NAN; + axes.minor = NAN; + axes.theta = NAN; + chisq = NAN; + nDOF = 0; + nPix = 0; + apRadius = NAN; + errMag = NAN; + } + + float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; + float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; + psS16 nImageOverlap = 1; + + psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; + float posAngle = 0.0; + float pltScale = 0.0; + pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + + 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", xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", 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", yErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + 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", errMag); + 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_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", 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)", ptSky.r*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); + 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", 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)", axes.major); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor", source->pixWeight); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", nDOF); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); + + // distinguish moments measure from window vs S/N > XX ?? + float mxx = source->moments ? source->moments->Mxx : NAN; + float mxy = source->moments ? source->moments->Mxy : NAN; + float myy = source->moments ? source->moments->Myy : NAN; + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", myy); + + psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); + + psVector *radFlux = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVector *radFill = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVectorInit (radFlux, NAN); + psVectorInit (radFluxErr, NAN); + psVectorInit (radFill, NAN); + if (!source->radial) goto empty_annuli; + if (!source->radial->flux) goto empty_annuli; + if (!source->radial->fill) goto empty_annuli; + psAssert (source->radial->flux->n <= radFlux->n, "inconsistent vector lengths"); + psAssert (source->radial->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 < source->radial->flux->n; j++) { + radFlux->data.F32[j] = source->radial->flux->data.F32[j]; + radFluxErr->data.F32[j] = source->radial->fluxErr->data.F32[j]; + radFill->data.F32[j] = source->radial->fill->data.F32[j]; + } + + empty_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_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); + psFree (radFlux); + psFree (radFluxErr); + psFree (radFill); + + // 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", nImageOverlap); + 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->errMag = 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->errMag : NAN; + + pmPSF_AxesToModel (PAR, axes); + + 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->flux = peakFlux; + source->peak->dx = dPAR[PM_PAR_XPOS]; + source->peak->dy = dPAR[PM_PAR_YPOS]; + source->peak->SN = sqrt(source->peak->flux); // XXX a proxy: various functions sort by peak S/N + + source->pixWeight = 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->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); + + 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; + + // 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, pmSourceSortBySN); + + 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"); + // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + + 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"); + + // int nRadialBins = 0; + // if (doAnnuli) { + // // get the max count of radial bins + // for (int i = 0; i < sources->n; i++) { + // pmSource *source = sources->data[i]; + // if (!source->extpars) continue; + // if (!source->extpars->radProfile ) continue; + // if (!source->extpars->radProfile->binSB) continue; + // nRadialBins = PS_MAX(nRadialBins, source->extpars->radProfile->binSB->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); + + 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); + } 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); + } + } + +# if (0) + // 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); + } + } +# endif + + // 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; +} + +// XXX this layout is still the same as PS1_DEV_1 +bool pmSourcesWrite_CMF_PS1_SV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, 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, pmSourceSortBySN); + + // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams + int nParamMax = 0; + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (source->modelFits == NULL) continue; + for (int j = 0; j < source->modelFits->n; j++) { + pmModel *model = source->modelFits->data[j]; + assert (model); + nParamMax = PS_MAX (nParamMax, model->params->n); + } + } + + table = psArrayAllocEmpty (sources->n); + + // 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); + + 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); + + 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) { + 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; +}