Index: branches/eam_branches/ipp-20110710/psModules/src/objects/Makefile.am =================================================================== --- branches/eam_branches/ipp-20110710/psModules/src/objects/Makefile.am (revision 31926) +++ branches/eam_branches/ipp-20110710/psModules/src/objects/Makefile.am (revision 31927) @@ -133,11 +133,11 @@ 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 + mksource.pl pmSourceIO_CMF.c.in PS1_V1 pmSourceIO_CMF_PS1_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 + mksource.pl pmSourceIO_CMF.c.in PS1_V2 pmSourceIO_CMF_PS1_V2.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 + mksource.pl pmSourceIO_CMF.c.in PS1_V3 pmSourceIO_CMF_PS1_V3.c # EXTRA_DIST = pmErrorCodes.h.in pmErrorCodes.dat pmErrorCodes.c.in Index: branches/eam_branches/ipp-20110710/psModules/src/objects/mksource.pl =================================================================== --- branches/eam_branches/ipp-20110710/psModules/src/objects/mksource.pl (revision 31926) +++ branches/eam_branches/ipp-20110710/psModules/src/objects/mksource.pl (revision 31927) @@ -14,7 +14,10 @@ # see if we can add in PS1_DV* and PS1_SV* as well... -@cmfmodes = ("PS1_V1", 1, +%cmfmodes = ("PS1_V1", 1, "PS1_V2", 2, "PS1_V3", 3); + +print "1: $cmfmodes{1}\n"; +print "PS1_V1: $cmfmodes{'PS1_V1'}\n"; open (FILE, "$template") || die "failed to open template $template\n"; @@ -50,15 +53,19 @@ if ($gtMode) { + # print "gtMode : $line\n"; $thisLevel = $cmfmodes{$gtMode}; $myLevel = $cmfmodes{$cmfmode}; - if ($thisLevel <= $myLevel) { next; } - $line =~ s|\@<\S*\@\s*||; + print "gtMode : $gtMode vs $cmfmode, $thisLevel, $myLevel\n"; + if ($myLevel <= $thisLevel) { next; } + $line =~ s|\@>\S*\@\s*||; } if ($ltMode) { + # print "ltMode : $line\n"; $thisLevel = $cmfmodes{$ltMode}; $myLevel = $cmfmodes{$cmfmode}; - if ($thisLevel >= $myLevel) { next; } - $line =~ s|\@>\S*\@\s*||; + print "ltMode : $ltMode vs $cmfmode, $thisLevel, $myLevel\n"; + if ($myLevel >= $thisLevel) { next; } + $line =~ s|\@<\S*\@\s*||; } Index: branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceIO_CMF.c.in (revision 31926) +++ branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceIO_CMF.c.in (revision 31927) @@ -222,5 +222,5 @@ // read in a readout from the fits file -psArray *pmSourcesRead_CMF_PS1_V3 (psFits *fits, psMetadata *header) +psArray *pmSourcesRead_CMF_@CMFMODE@ (psFits *fits, psMetadata *header) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -281,6 +281,11 @@ // 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->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); @ALL@ source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); + @=PS1_V3@ source->apMagRaw = psMetadataLookupF32 (&status, row, "AP_MAG_RAW"); + + // XXX use these to determine PAR[PM_PAR_I0] if they exist? + @=PS1_V3@ source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); + @=PS1_V3@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); // XXX this scaling is incorrect: does not include the 2 \pi AREA factor @@ -353,5 +358,5 @@ } -bool pmSourcesWrite_CMF_PS1_V3_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) +bool pmSourcesWrite_CMF_@CMFMODE@_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) { bool status; @@ -541,5 +546,5 @@ // 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) +bool pmSourcesWrite_CMF_@CMFMODE@_XFIT (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) { @@ -673,5 +678,5 @@ } -bool pmSourcesWrite_CMF_PS1_V3_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) +bool pmSourcesWrite_CMF_@CMFMODE@_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) { return true; Index: branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c =================================================================== --- branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 31926) +++ (revision ) @@ -1,605 +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 - -bool pmSourcesWrite_CMF_PS1_V1 (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 (); - 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, "RA_PSF", PS_DATA_F32, "PSF RA coordinate (degrees)", outputs.ra); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F32, "PSF DEC coordinate (degrees)", outputs.dec); - 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, "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, "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, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - - // 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; - } - - psMetadata *header = psMetadataCopy(NULL, tableHeader); - pmSourceMasksHeader(header); - - if (table->n == 0) { - if (!psFitsWriteBlank(fits, header, extname)) { - psError(psErrorCodeLast(), false, "Unable to write empty 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_V1 (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->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); - - 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); - - sources->data[i] = source; - psFree(row); - } - - return sources; -} - -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V1_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) { - if (!psFitsWriteBlank (fits, outhead, extname)) { - psError(psErrorCodeLast(), false, "Unable to write empty sources."); - 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_V1_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."); - 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_V1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - return true; -} Index: branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c =================================================================== --- branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 31926) +++ (revision ) @@ -1,669 +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 - -bool pmSourcesWrite_CMF_PS1_V2 (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 (); - 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, "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, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - - // 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_V2 (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); - - 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); - - sources->data[i] = source; - psFree(row); - } - - return sources; -} - -bool pmSourcesWrite_CMF_PS1_V2_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, 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"); - // 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_V2_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_V2_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - return true; -} Index: branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c =================================================================== --- branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 31926) +++ (revision ) @@ -1,680 +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 - -bool pmSourcesWrite_CMF_PS1_V3 (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 (); - 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, "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); - - // NOTE: RA & DEC (both double) need to be on an 8-byte boundary... - 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, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); - psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); - - 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); - - // 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); - 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, "PADDING2", PS_DATA_S32, "more padding", 0); - - // 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_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; - - 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]; - - source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); - 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 - 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); - - 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]; - - // we no longer sort by S/N, only flux - // if (isfinite (source->psfMagErr) && (source->psfMagErr > 0.0)) { - // source->peak->SN = 1.0 / source->psfMagErr; - // } else { - // source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N - // } - - 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_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: branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceOutputs.c =================================================================== --- branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceOutputs.c (revision 31926) +++ branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceOutputs.c (revision 31927) @@ -181,2 +181,62 @@ return true; } + +bool pmSourceLocalAstrometry (psSphere *ptSky, float *posAngle, float *pltScale, pmChip *chip, float xPos, float yPos) { + + pmFPA *fpa = chip->parent; + + if (!chip->toFPA) goto escape; + if (!fpa->toTPA) goto escape; + if (!fpa->toSky) goto escape; + + // generate RA,DEC + psPlane ptCH, ptFP, ptTP_o, ptTP_x, ptTP_y; + + // calculate the astrometry for the coordinate of interest + ptCH.x = xPos; + ptCH.y = yPos; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP_o, fpa->toTPA, &ptFP); + psDeproject (ptSky, &ptTP_o, fpa->toSky); + + // calculate the astrometry for the coordinate + 1pix in X + ptCH.x = xPos + 1.0; + ptCH.y = yPos; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP_x, fpa->toTPA, &ptFP); + + // calculate the astrometry for the coordinate + 1pix in Y + ptCH.x = xPos; + ptCH.y = yPos + 1.0; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP_y, fpa->toTPA, &ptFP); + + // the resulting Tangent Plane coordinates are in TP pixels; convert to local Tangent Plane + // degrees + + float dTPx_dCHx = fpa->toSky->Xs * (ptTP_x.x - ptTP_o.x); + float dTPy_dCHx = fpa->toSky->Ys * (ptTP_x.y - ptTP_o.y); + + float dTPx_dCHy = fpa->toSky->Xs * (ptTP_y.x - ptTP_o.x); + float dTPy_dCHy = fpa->toSky->Ys * (ptTP_y.y - ptTP_o.y); + + float pltScale_x = hypot(dTPx_dCHx, dTPy_dCHx); + float pltScale_y = hypot(dTPx_dCHy, dTPy_dCHy); + *pltScale = 0.5*(pltScale_x + pltScale_y); + + float posAngle_x = atan2 (+dTPy_dCHx, +dTPx_dCHx); + float posAngle_y = atan2 (-dTPy_dCHy, +dTPx_dCHy); + *posAngle = 0.5*(posAngle_x + posAngle_y); + + return true; + +escape: + // no astrometry calibration, give up + ptSky->r = NAN; + ptSky->d = NAN; + *posAngle = NAN; + *pltScale = NAN; + + return false; +} +