Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/Makefile.am (revision 30762) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/Makefile.am (revision 30763) @@ -32,4 +32,5 @@ pmSourceFitSet.c \ pmSourcePhotometry.c \ + pmSourceOutputs.c \ pmSourceIO.c \ pmSourceIO_RAW.c \ @@ -102,4 +103,5 @@ pmSourceFitSet.h \ pmSourcePhotometry.h \ + pmSourceOutputs.h \ pmSourceIO.h \ pmSourcePlots.h \ Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmPCMdata.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmPCMdata.c (revision 30762) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmPCMdata.c (revision 30763) @@ -341,5 +341,5 @@ return false; } - pcm->nDOF = pcm->nPix - nParams - 1; + pcm->nDOF = pcm->nPix - nParams; // has the source pixel window changed? Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitModel.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitModel.c (revision 30762) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitModel.c (revision 30763) @@ -40,4 +40,5 @@ #include "pmSourceDiffStats.h" #include "pmSource.h" +#include "pmSourcePhotometry.h" #include "pmSourceFitModel.h" @@ -238,9 +239,15 @@ // save the resulting chisq, nDOF, nIter - model->chisq = myMin->value; + if (options->poissonErrors) { + model->chisq = myMin->value; + model->nPix = y->n; + model->nDOF = y->n - nParams; + model->chisqNorm = model->chisq / model->nDOF; + } else { + pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, options->covarFactor, nParams); + } model->nIter = myMin->iter; - model->nPix = y->n; - model->nDOF = y->n - nParams; - model->chisqNorm = model->chisq / model->nDOF; + + // set the model success or failure status model->flags |= PM_MODEL_STATUS_FITTED; if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE; Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitModel.h =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitModel.h (revision 30762) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitModel.h (revision 30763) @@ -31,4 +31,5 @@ float maxChisqDOF; ///< convergence criterion float weight; ///< use this weight for constant-weight fits + float covarFactor; ///< covariance factor for calculating the chisq bool poissonErrors; ///< use poisson errors for fits? bool saveCovariance; Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitPCM.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitPCM.c (revision 30762) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitPCM.c (revision 30763) @@ -38,4 +38,5 @@ #include "pmSourceDiffStats.h" #include "pmSource.h" +#include "pmSourcePhotometry.h" #include "pmSourceFitModel.h" #include "pmPCMdata.h" @@ -85,9 +86,15 @@ // save the resulting chisq, nDOF, nIter - pcm->modelConv->chisq = myMin->value; + if (fitOptions->poissonErrors) { + pcm->modelConv->chisq = myMin->value; + pcm->modelConv->nPix = pcm->nPix; + pcm->modelConv->nDOF = pcm->nDOF; + pcm->modelConv->chisqNorm = pcm->modelConv->chisq / pcm->modelConv->nDOF; + } else { + pmSourceChisq (pcm->modelConv, source->pixels, source->maskObj, source->variance, maskVal, fitOptions->covarFactor, pcm->nPix - pcm->nDOF - 1); + } pcm->modelConv->nIter = myMin->iter; - pcm->modelConv->nPix = pcm->nPix; - pcm->modelConv->nDOF = pcm->nDOF; - pcm->modelConv->chisqNorm = pcm->modelConv->chisq / pcm->modelConv->nDOF; + + // set the model success or failure status pcm->modelConv->flags |= PM_MODEL_STATUS_FITTED; if (!fitStatus) pcm->modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE; Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitSet.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitSet.c (revision 30762) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitSet.c (revision 30763) @@ -39,4 +39,5 @@ #include "pmSourceDiffStats.h" #include "pmSource.h" +#include "pmSourcePhotometry.h" #include "pmSourceFitModel.h" @@ -299,5 +300,5 @@ const psVector *dparam, const psVector *param, const psImage *covar, pmSource *source, psMinimization *myMin, int nPix, - bool fitStatus, bool saveCovariance) + bool fitStatus, pmSourceFitOptions *options, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(set, false); @@ -322,5 +323,5 @@ psTrace ("psModules.objects", 4, "%f +/- %f", param->data.F32[n], dparam->data.F32[n]); } - if (saveCovariance) { + if (options->saveCovariance) { // we only save the covar matrix for this object with itself (ignore cross terms between objects) model->covar = psImageAlloc(model->params->n, model->params->n, PS_TYPE_F32); @@ -336,7 +337,13 @@ // save the resulting chisq, nDOF, nIter // these are not unique for any one source - model->chisq = myMin->value; - model->nIter = myMin->iter; - model->nDOF = nPix - model->params->n; + if (options->poissonErrors) { + model->chisq = myMin->value; + model->nPix = nPix; + model->nDOF = nPix - model->params->n; + model->chisqNorm = model->chisq / model->nDOF; + } else { + pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, options->covarFactor, model->params->n); + } + model->nIter = myMin->iter; // set the model success or failure status @@ -592,5 +599,5 @@ } - pmSourceFitSetValues (thisSet, dparams, params, covar, source, myMin, y->n, fitStatus, options->saveCovariance); + pmSourceFitSetValues (thisSet, dparams, params, covar, source, myMin, y->n, fitStatus, options, maskVal); psTrace ("psModules.objects", 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nPix: %ld\n", onPic, fitStatus, myMin->iter, myMin->value, y->n); Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitSet.h =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitSet.h (revision 30762) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitSet.h (revision 30763) @@ -44,5 +44,5 @@ const psVector *dparam, const psVector *param, const psImage *covar, pmSource *source, psMinimization *myMin, int nPix, - bool fitStatus, bool saveCovariance); + bool fitStatus, pmSourceFitOptions *options, psImageMaskType maskVal); psF32 pmSourceFitSetFunction(psVector *deriv, const psVector *param, const psVector *x); Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 30762) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 30763) @@ -49,4 +49,5 @@ #include "pmSourceIO.h" +#include "pmSourceOutputs.h" // panstarrs-style FITS table output (header + table in 1st extension) @@ -62,29 +63,9 @@ PS_ASSERT_PTR_NON_NULL(extname, false); - int i; psArray *table; psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - psF32 errMag, chisq, apRadius; - psS32 nPix, nDOF; pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - bool status1 = false; - bool status2 = false; - float magOffset = NAN; - float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); + // if the sequence is defined, write these in seq order; otherwise @@ -100,9 +81,16 @@ } + short nImageOverlap; + float magOffset; + float zeroptErr; + float fwhmMajor; + float fwhmMinor; + pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); + table = psArrayAllocEmpty (sources->n); // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. - for (i = 0; i < sources->n; i++) { + for (int i = 0; i < sources->n; i++) { pmSource *source = (pmSource *) sources->data[i]; @@ -115,55 +103,10 @@ } - // no difference between PSF and non-PSF model - pmModel *model = source->modelPSF; - - if (model != NULL) { - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - } else { - // in linear-fit mode, there is no error on the centroid - xErr = source->peak->dx; - yErr = source->peak->dy; - } - if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0); - } else { - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - } - chisq = model->chisq; - nDOF = model->nDOF; - nPix = model->nPix; - apRadius = source->apRadius; - errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; - } else { - xPos = source->peak->xf; - yPos = source->peak->yf; - xErr = source->peak->dx; - yErr = source->peak->dy; - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - chisq = NAN; - nDOF = 0; - nPix = 0; - apRadius = NAN; - errMag = NAN; - } - - float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; - float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; - psS16 nImageOverlap = 1; - - psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; - float posAngle = 0.0; - float pltScale = 0.0; - pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + // set the 'best' values for various output fields: + pmSourceOutputs outputs; + pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); + + pmSourceOutputsMoments moments; + pmSourceOutputsSetMoments (&moments, source); pmSourceDiffStats diffStats; @@ -175,48 +118,44 @@ row = psMetadataAlloc (); psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "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", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfFlux); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfFluxErr); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + psMetadataAdd (row, PS_LIST_TAIL, "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)", ptSky.r*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "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", chisq); + 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)", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_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", nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); - - // distinguish moments measure from window vs S/N > XX ?? - float mxx = source->moments ? source->moments->Mxx : NAN; - float mxy = source->moments ? source->moments->Mxy : NAN; - float myy = source->moments ? source->moments->Myy : NAN; - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", myy); - - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); + + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c (revision 30762) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c (revision 30763) @@ -49,4 +49,5 @@ #include "pmSourceIO.h" +#include "pmSourceOutputs.h" // panstarrs-style FITS table output (header + table in 1st extension) @@ -62,45 +63,8 @@ PS_ASSERT_PTR_NON_NULL(extname, false); - int i; psArray *table; psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - psF32 chisq, apRadius; - psS32 nPix, nDOF; pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - bool status1 = false; - bool status2 = false; - float magOffset = NAN; - float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); - - // we need a measure of the image quality (FWHM) for this image, in order to get the positional errors - float fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MAJ"); - if (!status1) { - fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW1"); - if (!status1) { - fwhmMajor = 5.0; // XXX just a guess! - } - } - float fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MIN"); - if (!status1) { - fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW2"); - if (!status1) { - fwhmMinor = 5.0; // XXX just a guess! - } - } // if the sequence is defined, write these in seq order; otherwise @@ -118,7 +82,14 @@ 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 (i = 0; i < sources->n; i++) { + for (int i = 0; i < sources->n; i++) { pmSource *source = (pmSource *) sources->data[i]; @@ -131,65 +102,10 @@ } - // no difference between PSF and non-PSF model - pmModel *model = source->modelPSF; - - if (model != NULL) { - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - } else { - // in linear-fit mode, there is no error on the centroid - xErr = fwhmMajor * source->errMag / 2.35; - yErr = fwhmMinor * source->errMag / 2.35; - } - if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0); - } else { - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - } - chisq = model->chisq; - nDOF = model->nDOF; - nPix = model->nPix; - apRadius = source->apRadius; - } else { - bool useMoments = true; - useMoments = (useMoments && source->moments); // can't if there are no moments - useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured - useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed... - - if (useMoments) { - xPos = source->moments->Mx; - yPos = source->moments->My; - xErr = fwhmMajor * source->errMag / 2.35; - yErr = fwhmMinor * source->errMag / 2.35; - } else { - xPos = source->peak->xf; - yPos = source->peak->yf; - xErr = source->peak->dx; - yErr = source->peak->dy; - } - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - chisq = NAN; - nDOF = 0; - nPix = 0; - apRadius = NAN; - } - - float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; - float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; - psS16 nImageOverlap = 1; - - psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; - float posAngle = 0.0; - float pltScale = 0.0; - pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + // set the 'best' values for various output fields: + pmSourceOutputs outputs; + pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); + + pmSourceOutputsMoments moments; + pmSourceOutputsSetMoments (&moments, source); pmSourceDiffStats diffStats; @@ -201,10 +117,10 @@ row = psMetadataAlloc (); psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "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->errMag); @@ -214,57 +130,44 @@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in real aperture", source->apMagRaw); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", source->apFlux); psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", source->apFluxErr); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + 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)", ptSky.r*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "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", chisq); + 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)", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_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", nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); - - // distinguish moments measure from window vs S/N > XX ?? - float mxx = source->moments ? source->moments->Mxx : NAN; - float mxy = source->moments ? source->moments->Mxy : NAN; - float myy = source->moments ? source->moments->Myy : NAN; - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", myy); - - float Mrf = source->moments ? source->moments->Mrf : NAN; - float Mrh = source->moments ? source->moments->Mrh : NAN; - float Krf = source->moments ? source->moments->KronFlux : NAN; - float dKrf = source->moments ? source->moments->KronFluxErr : NAN; - - float Kinner = source->moments ? source->moments->KronFinner : NAN; - float Kouter = source->moments ? source->moments->KronFouter : NAN; - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", Mrf); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", Mrh); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", Krf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", dKrf); - - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 1.0 R1)", Kinner); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 4.0 R1)", Kouter); - - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 1.0 R1)", moments.Kinner); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 4.0 R1)", moments.Kouter); + + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P", PS_DATA_F32, "distance to positive match source", diffStats.Rp); Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 30762) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 30763) @@ -50,4 +50,5 @@ #include "pmSourceIO.h" +#include "pmSourceOutputs.h" // panstarrs-style FITS table output (header + table in 1st extension) @@ -66,26 +67,6 @@ psArray *table; psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - psF32 errMag, chisq, apRadius; - psS32 nPix, nDOF; pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - bool status1 = false; - bool status2 = false; - float magOffset = NAN; - float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); // if the sequence is defined, write these in seq order; otherwise @@ -103,21 +84,10 @@ table = psArrayAllocEmpty (sources->n); -# if (0) - // we use this just to define the output vectors (which must be present for all objects) - bool status = false; - psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); - psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); - psAssert (radMax, "this must have been defined and tested earlier!"); - psAssert (radMax->n, "this must have been defined and tested earlier!"); - psAssert (radMin->n == radMax->n, "inconsistent annular bins"); - - // write the radial profile apertures to header - for (int i = 0; i < radMax->n; i++) { - sprintf (keyword1, "RMIN_%02d", i); - sprintf (keyword2, "RMAX_%02d", i); - psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]); - psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]); - } -# endif + 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 @@ -134,158 +104,64 @@ } - // no difference between PSF and non-PSF model - pmModel *model = source->modelPSF; - - if (model != NULL) { - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - } else { - // in linear-fit mode, there is no error on the centroid - xErr = source->peak->dx; - yErr = source->peak->dy; - } - if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0); - } else { - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - } - chisq = model->chisq; - nDOF = model->nDOF; - nPix = model->nPix; - apRadius = source->apRadius; - errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; - } else { - xPos = source->peak->xf; - yPos = source->peak->yf; - xErr = source->peak->dx; - yErr = source->peak->dy; - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - chisq = NAN; - nDOF = 0; - nPix = 0; - apRadius = NAN; - errMag = NAN; - } - - float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; - float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; - psS16 nImageOverlap = 1; - - psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; - float posAngle = 0.0; - float pltScale = 0.0; - pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + // 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", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "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", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux", source->psfFluxErr); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_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", apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + psMetadataAdd (row, PS_LIST_TAIL, "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)", ptSky.r*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "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", chisq); + 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)", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_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", nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); - - // distinguish moments measure from window vs S/N > XX ?? - float Mxx = source->moments ? source->moments->Mxx : NAN; - float Mxy = source->moments ? source->moments->Mxy : NAN; - float Myy = source->moments ? source->moments->Myy : NAN; - - float Mrf = source->moments ? source->moments->Mrf : NAN; - float Mrh = source->moments ? source->moments->Mrh : NAN; - float Krf = source->moments ? source->moments->KronFlux : NAN; - float dKrf = source->moments ? source->moments->KronFluxErr : NAN; - - float Kinner = source->moments ? source->moments->KronFinner : NAN; - float Kouter = source->moments ? source->moments->KronFouter : NAN; - - float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN; - float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN; - float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN; - float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN; - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", Mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", Mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", Myy); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", M_c3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", M_s3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", M_c4); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", M_s4); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", Mrf); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", Mrh); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", Krf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", dKrf); - - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", Kinner); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", Kouter); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kinner); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kouter); psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); - -# if (0) - // XXX if we have raw radial apertures, write them out here - psVector *radFlux = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVector *radFill = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVectorInit (radFlux, NAN); - psVectorInit (radFluxErr, NAN); - psVectorInit (radFill, NAN); - if (!source->radial) goto empty_annuli; - if (!source->radial->flux) goto empty_annuli; - if (!source->radial->fill) goto empty_annuli; - psAssert (source->radial->flux->n <= radFlux->n, "inconsistent vector lengths"); - psAssert (source->radial->fill->n <= radFlux->n, "inconsistent vector lengths"); - - // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux - for (int j = 0; j < source->radial->flux->n; j++) { - radFlux->data.F32[j] = source->radial->flux->data.F32[j]; - radFluxErr->data.F32[j] = source->radial->fluxErr->data.F32[j]; - radFill->data.F32[j] = source->radial->fill->data.F32[j]; - } - - empty_annuli: - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli", radFluxErr); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); - psFree (radFlux); - psFree (radFluxErr); - psFree (radFill); -# endif // XXX not sure how to get this : need to load Nimages with weight? Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 30762) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 30763) @@ -49,4 +49,5 @@ #include "pmSourceIO.h" +#include "pmSourceOutputs.h" // panstarrs-style FITS table output (header + table in 1st extension) @@ -62,27 +63,6 @@ psArray *table; psMetadata *row; - int i; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - psF32 errMag, chisq, apRadius; - psS32 nPix, nDOF; pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - bool status1 = false; - bool status2 = false; - float magOffset = NAN; - float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); // if the sequence is defined, write these in seq order; otherwise @@ -100,7 +80,14 @@ 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 (i = 0; i < sources->n; i++) { + for (int i = 0; i < sources->n; i++) { pmSource *source = (pmSource *) sources->data[i]; @@ -113,94 +100,45 @@ } - // no difference between PSF and non-PSF model - pmModel *model = source->modelPSF; - - if (model != NULL) { - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - } else { - // in linear-fit mode, there is no error on the centroid - xErr = source->peak->dx; - yErr = source->peak->dy; - } - if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0); - } else { - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - } - chisq = model->chisq; - nDOF = model->nDOF; - nPix = model->nPix; - apRadius = source->apRadius; - errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; - } else { - xPos = source->peak->xf; - yPos = source->peak->yf; - xErr = source->peak->dx; - yErr = source->peak->dy; - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - chisq = NAN; - nDOF = 0; - nPix = 0; - apRadius = NAN; - errMag = NAN; - } - - float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; - float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; - psS16 nImageOverlap = 1; - - psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; - float posAngle = 0.0; - float pltScale = 0.0; - pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + // 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", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F32, "PSF RA coordinate (degrees)", ptSky.r*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F32, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "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", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + psMetadataAdd (row, PS_LIST_TAIL, "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", chisq); + 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)", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_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", nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); - - // distinguish moments measure from window vs S/N > XX ?? - float mxx = source->moments ? source->moments->Mxx : NAN; - float mxy = source->moments ? source->moments->Mxy : NAN; - float myy = source->moments ? source->moments->Myy : NAN; - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", myy); + psMetadataAdd (row, PS_LIST_TAIL, "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); Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 30762) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 30763) @@ -49,4 +49,5 @@ #include "pmSourceIO.h" +#include "pmSourceOutputs.h" // panstarrs-style FITS table output (header + table in 1st extension) @@ -62,26 +63,6 @@ psArray *table; psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - psF32 errMag, chisq, apRadius; - psS32 nPix, nDOF; pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - bool status1 = false; - bool status2 = false; - float magOffset = NAN; - float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); // if the sequence is defined, write these in seq order; otherwise @@ -99,4 +80,11 @@ 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. @@ -112,94 +100,45 @@ } - // no difference between PSF and non-PSF model - pmModel *model = source->modelPSF; - - if (model != NULL) { - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - } else { - // in linear-fit mode, there is no error on the centroid - xErr = source->peak->dx; - yErr = source->peak->dy; - } - if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0); - } else { - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - } - chisq = model->chisq; - nDOF = model->nDOF; - nPix = model->nPix; - apRadius = source->apRadius; - errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; - } else { - xPos = source->peak->xf; - yPos = source->peak->yf; - xErr = source->peak->dx; - yErr = source->peak->dy; - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - chisq = NAN; - nDOF = 0; - nPix = 0; - apRadius = NAN; - errMag = NAN; - } - - float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; - float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; - psS16 nImageOverlap = 1; - - psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; - float posAngle = 0.0; - float pltScale = 0.0; - pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + // 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", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "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", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + psMetadataAdd (row, PS_LIST_TAIL, "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)", ptSky.r*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "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", chisq); + 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)", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_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", nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); - - // distinguish moments measure from window vs S/N > XX ?? - float mxx = source->moments ? source->moments->Mxx : NAN; - float mxy = source->moments ? source->moments->Mxy : NAN; - float myy = source->moments ? source->moments->Myy : NAN; - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", myy); + psMetadataAdd (row, PS_LIST_TAIL, "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); Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 30762) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 30763) @@ -49,4 +49,5 @@ #include "pmSourceIO.h" +#include "pmSourceOutputs.h" // panstarrs-style FITS table output (header + table in 1st extension) @@ -62,42 +63,6 @@ psArray *table; psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - psF32 chisq, apRadius; - psS32 nPix, nDOF; pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - bool status1 = false; - bool status2 = false; - float magOffset = NAN; - float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); - - // we need a measure of the image quality (FWHM) for this image, in order to get the positional errors - float fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MAJ"); - if (!status1) { - fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW1"); - if (!status1) { - fwhmMajor = 5.0; // XXX just a guess! - } - } - float fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MIN"); - if (!status1) { - fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW2"); - if (!status1) { - fwhmMinor = 5.0; // XXX just a guess! - } - } // if the sequence is defined, write these in seq order; otherwise @@ -115,4 +80,11 @@ 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. @@ -128,73 +100,19 @@ } - // no difference between PSF and non-PSF model - pmModel *model = source->modelPSF; - - if (model != NULL) { - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - } else { - xErr = fwhmMajor * source->errMag / 2.35; - yErr = fwhmMinor * source->errMag / 2.35; - } - if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0); - } else { - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - } - chisq = model->chisq; - nDOF = model->nDOF; - nPix = model->nPix; - apRadius = source->apRadius; - } else { - bool useMoments = true; - useMoments = (useMoments && source->moments); // can't if there are no moments - useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured - useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed... - - if (useMoments) { - xPos = source->moments->Mx; - yPos = source->moments->My; - xErr = fwhmMajor * source->errMag / 2.35; - yErr = fwhmMinor * source->errMag / 2.35; - } else { - xPos = source->peak->xf; - yPos = source->peak->yf; - xErr = source->peak->dx; - yErr = source->peak->dy; - } - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - chisq = NAN; - nDOF = 0; - nPix = 0; - apRadius = NAN; - } - - float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; - float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; - psS16 nImageOverlap = 1; - - psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; - float posAngle = 0.0; - float pltScale = 0.0; - pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + // 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", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "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->errMag); @@ -203,66 +121,46 @@ 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", apRadius); - - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + 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)", ptSky.r*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); - - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); + 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", chisq); + 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)", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_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", nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); - - // distinguish moments measure from window vs S/N > XX ?? - float Mxx = source->moments ? source->moments->Mxx : NAN; - float Mxy = source->moments ? source->moments->Mxy : NAN; - float Myy = source->moments ? source->moments->Myy : NAN; - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", Mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", Mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", Myy); - - float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN; - float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN; - float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN; - float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN; - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", M_c3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", M_s3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", M_c4); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", M_s4); - - float Mrf = source->moments ? source->moments->Mrf : NAN; - float Mrh = source->moments ? source->moments->Mrh : NAN; - float Krf = source->moments ? source->moments->KronFlux : NAN; - float dKrf = source->moments ? source->moments->KronFluxErr : NAN; - - float Kinner = source->moments ? source->moments->KronFinner : NAN; - float Kouter = source->moments ? source->moments->KronFouter : NAN; - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", Mrf); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", Mrh); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", Krf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", dKrf); - - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", Kinner); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", Kouter); - - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kinner); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kouter); + psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); + psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); psMetadataAdd (row, PS_LIST_TAIL, "PADDING2", PS_DATA_S32, "more padding", 0); Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceOutputs.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceOutputs.c (revision 30763) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceOutputs.c (revision 30763) @@ -0,0 +1,174 @@ +#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" + +bool pmSourceOutputsCommonValues (short *nImageOverlap, float *magOffset, float *zeroptErr, float *fwhmMajor, float *fwhmMinor, pmReadout *readout, psMetadata *header) { + + pmFPA *fpa = readout->parent->parent->parent; + + bool status1 = false; + bool status2 = false; + float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); + float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); + if (!isfinite(zeropt)) { + zeropt = psMetadataLookupF32 (&status2, header, "ZPT_OBS"); + } + if (status1 && status2 && (exptime > 0.0)) { + *magOffset = zeropt + 2.5*log10(exptime); + } + *zeroptErr = psMetadataLookupF32 (&status2, header, "ZPT_ERR"); + + // we need a measure of the image quality (FWHM) for this image, in order to get the positional errors + *fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MAJ"); + if (!status1) { + *fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW1"); + if (!status1) { + *fwhmMajor = 5.0; // XXX just a guess! + } + } + *fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MIN"); + if (!status1) { + *fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW2"); + if (!status1) { + *fwhmMinor = 5.0; // XXX just a guess! + } + } + + *nImageOverlap = 1; + return true; + +} + +bool pmSourceOutputsSetValues (pmSourceOutputs *outputs, pmSource *source, pmChip *chip, float fwhmMajor, float fwhmMinor, float magOffset) { + + psF32 *PAR, *dPAR; + psEllipseAxes axes; + + // no difference between PSF and non-PSF model + pmModel *model = source->modelPSF; + + if (model != NULL) { + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + outputs->xPos = PAR[PM_PAR_XPOS]; + outputs->yPos = PAR[PM_PAR_YPOS]; + if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { + outputs->xErr = dPAR[PM_PAR_XPOS]; + outputs->yErr = dPAR[PM_PAR_YPOS]; + } else { + outputs->xErr = fwhmMajor * source->errMag / 2.35; + outputs->yErr = fwhmMinor * source->errMag / 2.35; + } + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { + axes = pmPSF_ModelToAxes (PAR, 20.0); + outputs->psfMajor = axes.major; + outputs->psfMinor = axes.minor; + outputs->psfTheta = axes.theta*PS_DEG_RAD; + } else { + outputs->psfMajor = NAN; + outputs->psfMinor = NAN; + outputs->psfTheta = NAN; + } + outputs->chisq = model->chisq; + outputs->nDOF = model->nDOF; + outputs->nPix = model->nPix; + outputs->apRadius = source->apRadius; + } else { + bool useMoments = true; + useMoments = (useMoments && source->moments); // can't if there are no moments + useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured + useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed... + + if (useMoments) { + outputs->xPos = source->moments->Mx; + outputs->yPos = source->moments->My; + outputs->xErr = fwhmMajor * source->errMag / 2.35; + outputs->yErr = fwhmMinor * source->errMag / 2.35; + } else { + outputs->xPos = source->peak->xf; + outputs->yPos = source->peak->yf; + outputs->xErr = source->peak->dx; + outputs->yErr = source->peak->dy; + } + outputs->psfMajor = NAN; + outputs->psfMinor = NAN; + outputs->psfTheta = NAN; + outputs->chisq = NAN; + outputs->nDOF = 0; + outputs->nPix = 0; + outputs->apRadius = NAN; + } + + outputs->calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; + outputs->peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; + + psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; + float posAngle = 0.0; + float pltScale = 0.0; + pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, outputs->xPos, outputs->yPos); + + outputs->posAngle = posAngle*PS_DEG_RAD; + outputs->pltScale = pltScale*PS_DEG_RAD*3600.0; + + outputs->ra = ptSky.r*PS_DEG_RAD; + outputs->dec = ptSky.d*PS_DEG_RAD; + + return true; +} + +bool pmSourceOutputsSetMoments (pmSourceOutputsMoments *moments, pmSource *source) { + + // distinguish moments measure from window vs S/N > XX ?? + moments->Mxx = source->moments ? source->moments->Mxx : NAN; + moments->Mxy = source->moments ? source->moments->Mxy : NAN; + moments->Myy = source->moments ? source->moments->Myy : NAN; + moments->M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN; + moments->M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN; + moments->M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN; + moments->M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN; + moments->Mrf = source->moments ? source->moments->Mrf : NAN; + moments->Mrh = source->moments ? source->moments->Mrh : NAN; + moments->Krf = source->moments ? source->moments->KronFlux : NAN; + moments->dKrf = source->moments ? source->moments->KronFluxErr : NAN; + moments->Kinner = source->moments ? source->moments->KronFinner : NAN; + moments->Kouter = source->moments ? source->moments->KronFouter : NAN; + + return true; +} Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceOutputs.h =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceOutputs.h (revision 30763) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceOutputs.h (revision 30763) @@ -0,0 +1,60 @@ +/* @file pmSourceOutputs.h + * @brief functions to perform common I/O conversions + * + * @author EAM, IfA + * + * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-16 22:30:50 $ + * Copyright 2011 Institute for Astronomy, University of Hawaii + * + */ + +# ifndef PM_SOURCE_OUTPUTS_H +# define PM_SOURCE_OUTPUTS_H + +/// @addtogroup Objects Object Detection / Analysis Functions +/// @{ + +typedef struct { + float xPos; + float yPos; + float xErr; + float yErr; + double ra; + double dec; + float psfMajor; + float psfMinor; + float psfTheta; + float chisq; + float apRadius; + int nPix; + int nDOF; + float calMag; + float peakMag; + float posAngle; + float pltScale; +} pmSourceOutputs; + +typedef struct { + float Mxx; + float Mxy; + float Myy; + float M_c3; + float M_s3; + float M_c4; + float M_s4; + float Mrf; + float Mrh; + float Krf; + float dKrf; + float Kinner; + float Kouter; +} pmSourceOutputsMoments; + +bool pmSourceOutputsCommonValues (short *nImageOverlap, float *magOffset, float *zeroptErr, float *fwhmMajor, float *fwhmMinor, pmReadout *readout, psMetadata *header); + +bool pmSourceOutputsSetValues (pmSourceOutputs *outputs, pmSource *source, pmChip *chip, float fwhmMajor, float fwhmMinor, float magOffset); + +bool pmSourceOutputsSetMoments (pmSourceOutputsMoments *moments, pmSource *source); + +# endif Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.c (revision 30762) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.c (revision 30763) @@ -422,4 +422,8 @@ *pixWeightNotPoor = notPoorSum / modelSum; + if (false && isfinite(*pixWeightNotBad) && isfinite(*pixWeightNotPoor)) { + psAssert (*pixWeightNotBad <= *pixWeightNotPoor, "error: all bad pixels should also be poor"); + } + return (true); } @@ -674,5 +678,5 @@ // determine chisq, etc for linear normalization-only fit -bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal, const float covarFactor) +bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal, const float covarFactor, int nParams) { PS_ASSERT_PTR_NON_NULL(model, false); @@ -689,11 +693,14 @@ if (variance->data.F32[j][i] <= 0) continue; - dC += PS_SQR (image->data.F32[j][i]) / (covarFactor * variance->data.F32[j][i]); + // dC += PS_SQR (image->data.F32[j][i]) / (covarFactor * variance->data.F32[j][i]); + dC += PS_SQR (image->data.F32[j][i]) / variance->data.F32[j][i]; Npix ++; } } + model->nPix = Npix; - model->nDOF = Npix - 1; + model->nDOF = Npix - nParams - 1; model->chisq = dC; + model->chisqNorm = dC / model->nDOF; return (true); Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.h =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.h (revision 30762) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.h (revision 30763) @@ -69,5 +69,5 @@ bool pmSourcePixelWeight (float *pixWeightNotBad, float *pixWeightNotPoor, pmModel *model, psImage *mask, psImageMaskType maskVal, psImageMaskType markVal); -bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal, const float covarFactor); +bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal, const float covarFactor, int nParams); bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);