Index: trunk/psModules/src/objects/Makefile.am =================================================================== --- trunk/psModules/src/objects/Makefile.am (revision 31634) +++ trunk/psModules/src/objects/Makefile.am (revision 31670) @@ -128,2 +128,16 @@ CLEANFILES = *~ + +# pmSourceID_CMF_* functions use a common framework +BUILT_SOURCES = pmSourceIO_CMF_PS1_V1.v1.c pmSourceIO_CMF_PS1_V2.v1.c pmSourceIO_CMF_PS1_V3.v1.c + +pmSourceIO_CMF_PS1_V1.v1.c : pmSourceIO_CMF.c.in mksource.pl + mksource.pl pmSourceIO_CMF.c.in PS1_V1 pmSourceIO_CMF_PS1_V1.v1.c + +pmSourceIO_CMF_PS1_V2.v1.c : pmSourceIO_CMF.c.in mksource.pl + mksource.pl pmSourceIO_CMF.c.in PS1_V2 pmSourceIO_CMF_PS1_V2.v1.c + +pmSourceIO_CMF_PS1_V3.v1.c : pmSourceIO_CMF.c.in mksource.pl + mksource.pl pmSourceIO_CMF.c.in PS1_V2 pmSourceIO_CMF_PS1_V3.v1.c + +# EXTRA_DIST = pmErrorCodes.h.in pmErrorCodes.dat pmErrorCodes.c.in Index: trunk/psModules/src/objects/mksource.pl =================================================================== --- trunk/psModules/src/objects/mksource.pl (revision 31670) +++ trunk/psModules/src/objects/mksource.pl (revision 31670) @@ -0,0 +1,70 @@ +#!/usr/bin/env perl + +# this program takes the pmSourceIO_CMF.in.c template file and generates the .c version based on the given I/O format made + +if (@ARGV != 3) { die "USAGE: mksource (template) (cmfmode) (output)\n"; } + +$template = $ARGV[0]; +$cmfmode = $ARGV[1]; +$output = $ARGV[2]; + +if (! -f $template) { die "missing template file $template\n"; } + +# validate the cmfmode + +# see if we can add in PS1_DV* and PS1_SV* as well... +@cmfmodes = ("PS1_V1", 1, + "PS1_V2", 2, + "PS1_V3", 3); + +open (FILE, "$template") || die "failed to open template $template\n"; +@list = ; +close (FILE); + +open (FILE, ">$output"); + +# operations we can perform: +# @CMFMODE@ : replace with $cmfmode +# @ALL@ : remove and keep the rest of the line +# @=MODE@ : remove and keep if cmfmode == MODE +# @>MODE@ : remove and keep if cmfmode > MODE +# @ MODE + +foreach $line (@list) { + + # replace @CMFMODE@ wherever it appears + $line =~ s|\@CMFMODE\@|$cmfmode|g; + + if ($line =~ m|\@ALL\@|) { + $line =~ s|\@ALL\@\s*||; + } + + ($isMode) = ($line =~ m|\@=(\S*)\@|); + ($gtMode) = ($line =~ m|\@>(\S*)\@|); + ($ltMode) = ($line =~ m|\@<(\S*)\@|); + + if ($isMode) { + if ($isMode ne $cmfmode) { next; } + $line =~ s|\@=\S*\@\s*||; + } + + if ($gtMode) { + $thisLevel = $cmfmodes{$gtMode}; + $myLevel = $cmfmodes{$cmfmode}; + if ($thisLevel <= $myLevel) { next; } + $line =~ s|\@<\S*\@\s*||; + } + + if ($ltMode) { + $thisLevel = $cmfmodes{$ltMode}; + $myLevel = $cmfmodes{$cmfmode}; + if ($thisLevel >= $myLevel) { next; } + $line =~ s|\@>\S*\@\s*||; + } + + print FILE $line; +} + +close (FILE); + +exit 0; Index: trunk/psModules/src/objects/models/pmModel_PS1_V1.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 31634) +++ trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 31670) @@ -296,9 +296,9 @@ psF64 sigma = axes.major; + // we can solve directly if PAR_7 goes to 0.0 if (PAR[PM_PAR_7] == 0.0) { psF32 z = powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA); return ( sigma * sqrt (2.0 * z) ); } - psF64 limit = flux / PAR[PM_PAR_I0]; Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 31634) +++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 31670) @@ -289,5 +289,6 @@ if (PAR[PM_PAR_I0] <= 0) return 1.0; if (flux >= PAR[PM_PAR_I0]) return 1.0; - if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA); + + // if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA); shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; @@ -297,4 +298,10 @@ psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); psF64 sigma = axes.major; + + // we can solve directly if PAR_7 goes to 0.0 + if (PAR[PM_PAR_7] == 0.0) { + psF32 z = powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA); + return ( sigma * sqrt (2.0 * z) ); + } psF64 limit = flux / PAR[PM_PAR_I0]; Index: trunk/psModules/src/objects/pmSource.c =================================================================== --- trunk/psModules/src/objects/pmSource.c (revision 31634) +++ trunk/psModules/src/objects/pmSource.c (revision 31670) @@ -362,5 +362,5 @@ *****************************************************************************/ -pmPSFClump pmSourcePSFClump(psImage **savedImage, psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_GRID_SCALE, psF32 SX_MAX, psF32 SY_MAX, psF32 AR_MAX) +pmPSFClump pmSourcePSFClump(psImage **savedImage, psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_GRID_SCALE, psF32 SX_MAX, psF32 SY_MAX, psF32 SX_MIN, psF32 SY_MIN, psF32 AR_MAX) { psTrace("psModules.objects", 10, "---- begin ----\n"); @@ -422,5 +422,5 @@ // Sx,Sy are limited at 0. a peak at 0,0 is artificial - if (fabs(Mxx) < 0.05 || fabs(Myy < 0.05)) { + if (fabs(Mxx) < SX_MIN || fabs(Myy < SY_MIN)) { psTrace("psModules.objects", 10, "Rejecting source from clump because of low moments (%f,%f)\n", Index: trunk/psModules/src/objects/pmSource.h =================================================================== --- trunk/psModules/src/objects/pmSource.h (revision 31634) +++ trunk/psModules/src/objects/pmSource.h (revision 31670) @@ -219,4 +219,6 @@ psF32 SX_MAX, psF32 SY_MAX, + psF32 SX_MIN, + psF32 SY_MIN, psF32 AR_MAX ); Index: trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO.c (revision 31634) +++ trunk/psModules/src/objects/pmSourceIO.c (revision 31670) @@ -727,9 +727,10 @@ psMetadata *fileData = psMetadataLookupMetadata(NULL, file->format, "FILE"); // File information - const char *fpaNameHdr = psMetadataLookupStr(NULL, fileData, "FPA.OBS"); + const char *fpaNameHdr = psMetadataLookupStr(&status, fileData, "FPA.OBS"); if (fpaNameHdr && strlen(fpaNameHdr) > 0) { - const char *fpaName = psMetadataLookupStr(NULL, fpa->concepts, "FPA.OBS"); - psMetadataAddStr(outhead, PS_LIST_TAIL, fpaNameHdr, PS_META_REPLACE, - "FPA observation identifier", fpaName); + const char *fpaName = psMetadataLookupStr(&status, fpa->concepts, "FPA.OBS"); + if (fpaName) { + psMetadataAddStr(outhead, PS_LIST_TAIL, fpaNameHdr, PS_META_REPLACE, "FPA observation identifier", fpaName); + } } Index: trunk/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 31670) +++ trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 31670) @@ -0,0 +1,678 @@ +/** @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 + +// # define MODE @CMFMODE@ +bool pmSourcesWrite_CMF_@CMFMODE@ (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); + + short nImageOverlap; + float magOffset; + float zeroptErr; + float fwhmMajor; + float fwhmMinor; + pmSourceOutputsCommonValues (&nImageOverlap, &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); + + row = psMetadataAlloc (); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); + + // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision + @=PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F32, "PSF RA coordinate (degrees)", outputs.ra); + @=PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F32, "PSF DEC coordinate (degrees)", outputs.dec); + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); + + @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); + @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux", source->psfFluxErr); + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); + @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in reported aperture", source->apMagRaw); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); + @PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); + @>PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); + + @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); + @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", source->pixWeightNotPoor); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); + + @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); + @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); + @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); + @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); + + @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); + @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); + @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); + @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); + @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kinner); + @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kouter); + + // Do NOT write these : not consistent with the definition of PS1_V3 in Ohana/src/libautocode/dev/cmf-ps1-v3.d + // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_FLUX", PS_DATA_F32, "Kron Flux (in 1.0 R1)", moments.KronCore); + // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_ERROR", PS_DATA_F32, "Kron Error (in 1.0 R1)", moments.KronCoreErr); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); + @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); + @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2", PS_DATA_S32, "more padding", 0); + + // XXX not sure how to get this : need to load Nimages with weight? + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap); + @ALL@ 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_V3 (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; + + @ALL@ source->seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + @ALL@ PAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF"); + @ALL@ PAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF"); + @ALL@ dPAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF_SIG"); + @ALL@ dPAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF_SIG"); + @ALL@ axes.major = psMetadataLookupF32 (&status, row, "PSF_MAJOR"); + @ALL@ axes.minor = psMetadataLookupF32 (&status, row, "PSF_MINOR"); + @ALL@ axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); + + @ALL@ PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); + @ALL@ dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); + @ALL@ source->sky = PAR[PM_PAR_SKY]; + @ALL@ source->skyErr = dPAR[PM_PAR_SKY]; + + // XXX use these to determine PAR[PM_PAR_I0]? + @ALL@ source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); + @ALL@ source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); + @ALL@ source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); + + // XXX this scaling is incorrect: does not include the 2 \pi AREA factor + @ALL@ PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; + @ALL@ dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; + + pmPSF_AxesToModel (PAR, axes); + + @ALL@ float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); + @ALL@ float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; + + // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) + @ALL@ source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); + @ALL@ source->peak->rawFlux = peakFlux; + @ALL@ source->peak->smoothFlux = peakFlux; + @ALL@ source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position + @ALL@ source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position + @ALL@ source->peak->dx = dPAR[PM_PAR_XPOS]; + @ALL@ source->peak->dy = dPAR[PM_PAR_YPOS]; + + @ALL@ source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); + @=PS1_V3@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT"); + @ALL@ source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); + @ALL@ source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); + @ALL@ source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); + + // note that some older versions used PSF_PROBABILITY: this was not well defined. + @ALL@ model->chisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); + @ALL@ model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); + @ALL@ model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); + + @ALL@ source->moments = pmMomentsAlloc (); + @ALL@ source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf + @ALL@ source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf + + @ALL@ source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); + @ALL@ source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); + @ALL@ source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); + + @=PS1_V3@ source->moments->Mrf = psMetadataLookupF32 (&status, row, "MOMENTS_R1"); + @=PS1_V3@ source->moments->Mrh = psMetadataLookupF32 (&status, row, "MOMENTS_RH"); + @=PS1_V3@ source->moments->KronFlux = psMetadataLookupF32 (&status, row, "KRON_FLUX"); + @=PS1_V3@ source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR"); + + @=PS1_V3@ source->moments->KronFinner = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER"); + @=PS1_V3@ 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 + @=PS1_V3@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C"); + @=PS1_V3@ source->moments->Mxxy = 0.0; + @=PS1_V3@ source->moments->Mxyy = 0.0; + @=PS1_V3@ source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S"); + + @=PS1_V3@ source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C"); + @=PS1_V3@ source->moments->Mxxxy = 0.0; + @=PS1_V3@ source->moments->Mxxyy = 0.0; + @=PS1_V3@ source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S"); + @=PS1_V3@ source->moments->Myyyy = 0.0; + + @ALL@ source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); + @=PS1_V3@ source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2"); + assert (status); + + sources->data[i] = source; + psFree(row); + } + + return sources; +} + +bool pmSourcesWrite_CMF_PS1_V3_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++) { + 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); + } + } + + // 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_V3_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); + + 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; +} + +bool pmSourcesWrite_CMF_PS1_V3_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) +{ + return true; +} Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 31634) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 31670) @@ -113,4 +113,5 @@ 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); Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 31634) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 31670) @@ -270,8 +270,12 @@ 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->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); + source->apMagRaw = psMetadataLookupF32 (&status, row, "AP_MAG_RAW"); + + // XXX use these to determine PAR[PM_PAR_I0] if they exist? + source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); + source->psfFluxErr = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); // XXX this scaling is incorrect: does not include the 2 \pi AREA factor Index: trunk/psModules/src/objects/pmSourceMasks.h =================================================================== --- trunk/psModules/src/objects/pmSourceMasks.h (revision 31634) +++ trunk/psModules/src/objects/pmSourceMasks.h (revision 31670) @@ -53,4 +53,8 @@ PM_SOURCE_MODE2_PASS1_SRC = 0x00000080, ///< source detected in first pass analysis + + PM_SOURCE_MODE2_HAS_BRIGHTER_NEIGHBOR = 0x00000100, ///< peak is not the brightest in its footprint + PM_SOURCE_MODE2_BRIGHT_NEIGHBOR_1 = 0x00000200, ///< flux_n / (r^2 flux_p) > 1 + PM_SOURCE_MODE2_BRIGHT_NEIGHBOR_10 = 0x00000400, ///< flux_n / (r^2 flux_p) > 10 } pmSourceMode2; Index: trunk/psModules/src/objects/pmSourceMoments.c =================================================================== --- trunk/psModules/src/objects/pmSourceMoments.c (revision 31634) +++ trunk/psModules/src/objects/pmSourceMoments.c (revision 31670) @@ -264,7 +264,10 @@ // if Mrf (first radial moment) is very small, we are getting into low-significance - // territory. saturate at minKronRadius. conversely, if Mrf is > radius, we are clearly - // making an error. saturate at radius. - float kronRefRadius = MIN(radius, MAX(minKronRadius, source->moments->Mrf)); + // territory. saturate at minKronRadius. conversely, if Mrf is >> radius for faint + // sources, we are clearly making an error. saturate at radius. + float kronRefRadius = MAX(minKronRadius, source->moments->Mrf); + if (source->moments->SN < 10) { + kronRefRadius = MIN(radius, kronRefRadius); + } float radKinner = 1.0*kronRefRadius; Index: trunk/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- trunk/psModules/src/objects/pmSourcePhotometry.c (revision 31634) +++ trunk/psModules/src/objects/pmSourcePhotometry.c (revision 31670) @@ -134,5 +134,5 @@ // measure PSF model photometry status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, modelPSF); - source->psfFluxErr = source->psfFlux * source->psfMagErr; + source->psfFluxErr = fabs(source->psfFlux * source->psfMagErr); # if (0) @@ -191,4 +191,6 @@ } + pmSourceNeighborFlags (source); + // measure the aperture magnitude, if (SN > AP_MIN_SN) if (!isfinite(SN)) { @@ -264,4 +266,41 @@ */ + +bool pmSourceNeighborFlags (pmSource *source) { + + return false; + + // source must have a peak to have a footprint + if (!source) return false; + if (!source->peak) return false; + if (!source->peak->footprint) return false; + if (!source->peak->footprint->peaks) return false; + if (!source->peak->footprint->peaks->n) return false; + + // find the brightest peak (first peak) + pmPeak *brightPeak = source->peak->footprint->peaks->data[0]; + + // are we the brightest peak? + if (source->peak == brightPeak) return true; + + // if not, raise a flag: + source->mode2 |= PM_SOURCE_MODE2_HAS_BRIGHTER_NEIGHBOR; + + // but, this is a common situation. more interesting is if the ratio flux_n / (r^2 flux_o) is large + + float radius2 = PS_SQR(source->peak->xf - brightPeak->xf) + PS_SQR(source->peak->yf - brightPeak->yf); + + float ratio = brightPeak->rawFlux / (source->peak->rawFlux * radius2); + + if (ratio > 1) { + source->mode2 |= PM_SOURCE_MODE2_BRIGHT_NEIGHBOR_1; + } + + if (ratio > 10) { + source->mode2 |= PM_SOURCE_MODE2_BRIGHT_NEIGHBOR_10; + } + + return true; +} // return source model magnitude Index: trunk/psModules/src/objects/pmSourcePhotometry.h =================================================================== --- trunk/psModules/src/objects/pmSourcePhotometry.h (revision 31634) +++ trunk/psModules/src/objects/pmSourcePhotometry.h (revision 31670) @@ -79,4 +79,6 @@ double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal); +bool pmSourceNeighborFlags (pmSource *source); + // retire these: // double pmSourceCrossProduct(const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum);