IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 30763


Ignore:
Timestamp:
Feb 28, 2011, 2:44:39 PM (15 years ago)
Author:
eugene
Message:

unified function to generate output parameters for output cmf files; add chisqNorm calculation to pmSourceChisq, pass in nParams; re-calculate the chisq if we use constant weights; make sure we correctly set nDOF = nPoints - nParams

Location:
branches/eam_branches/ipp-20110213/psModules/src/objects
Files:
2 added
15 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20110213/psModules/src/objects/Makefile.am

    r29004 r30763  
    3232        pmSourceFitSet.c \
    3333        pmSourcePhotometry.c \
     34        pmSourceOutputs.c \
    3435        pmSourceIO.c \
    3536        pmSourceIO_RAW.c \
    … …  
    102103        pmSourceFitSet.h \
    103104        pmSourcePhotometry.h \
     105        pmSourceOutputs.h \
    104106        pmSourceIO.h \
    105107        pmSourcePlots.h \

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmPCMdata.c

    r29546 r30763  
    341341        return false;
    342342    }
    343     pcm->nDOF = pcm->nPix - nParams - 1;
     343    pcm->nDOF = pcm->nPix - nParams;
    344344
    345345    // has the source pixel window changed?

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitModel.c

    r30705 r30763  
    4040#include "pmSourceDiffStats.h"
    4141#include "pmSource.h"
     42#include "pmSourcePhotometry.h"
    4243#include "pmSourceFitModel.h"
    4344
    … …  
    238239
    239240    // save the resulting chisq, nDOF, nIter
    240     model->chisq = myMin->value;
     241    if (options->poissonErrors) {
     242        model->chisq = myMin->value;
     243        model->nPix  = y->n;
     244        model->nDOF  = y->n - nParams;
     245        model->chisqNorm = model->chisq / model->nDOF;
     246    } else {
     247        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, options->covarFactor, nParams);
     248    }
    241249    model->nIter = myMin->iter;
    242     model->nPix  = y->n;
    243     model->nDOF  = y->n - nParams;
    244     model->chisqNorm = model->chisq / model->nDOF;
     250
     251    // set the model success or failure status
    245252    model->flags |= PM_MODEL_STATUS_FITTED;
    246253    if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE;

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitModel.h

    r30705 r30763  
    3131    float maxChisqDOF;                  ///< convergence criterion
    3232    float weight;                       ///< use this weight for constant-weight fits
     33    float covarFactor;                  ///< covariance factor for calculating the chisq
    3334    bool poissonErrors;                 ///< use poisson errors for fits?
    3435    bool saveCovariance;

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitPCM.c

    r30719 r30763  
    3838#include "pmSourceDiffStats.h"
    3939#include "pmSource.h"
     40#include "pmSourcePhotometry.h"
    4041#include "pmSourceFitModel.h"
    4142#include "pmPCMdata.h"
    … …  
    8586
    8687    // save the resulting chisq, nDOF, nIter
    87     pcm->modelConv->chisq = myMin->value;
     88    if (fitOptions->poissonErrors) {
     89        pcm->modelConv->chisq = myMin->value;
     90        pcm->modelConv->nPix = pcm->nPix;
     91        pcm->modelConv->nDOF = pcm->nDOF;
     92        pcm->modelConv->chisqNorm = pcm->modelConv->chisq / pcm->modelConv->nDOF;
     93    } else {
     94        pmSourceChisq (pcm->modelConv, source->pixels, source->maskObj, source->variance, maskVal, fitOptions->covarFactor, pcm->nPix - pcm->nDOF - 1);
     95    }
    8896    pcm->modelConv->nIter = myMin->iter;
    89     pcm->modelConv->nPix = pcm->nPix;
    90     pcm->modelConv->nDOF = pcm->nDOF;
    91     pcm->modelConv->chisqNorm = pcm->modelConv->chisq / pcm->modelConv->nDOF;
     97
     98    // set the model success or failure status
    9299    pcm->modelConv->flags |= PM_MODEL_STATUS_FITTED;
    93100    if (!fitStatus) pcm->modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE;

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitSet.c

    r30705 r30763  
    3939#include "pmSourceDiffStats.h"
    4040#include "pmSource.h"
     41#include "pmSourcePhotometry.h"
    4142
    4243#include "pmSourceFitModel.h"
    … …  
    299300                           const psVector *dparam, const psVector *param, const psImage *covar,
    300301                           pmSource *source, psMinimization *myMin, int nPix,
    301                            bool fitStatus, bool saveCovariance)
     302                           bool fitStatus, pmSourceFitOptions *options, psImageMaskType maskVal)
    302303{
    303304    PS_ASSERT_PTR_NON_NULL(set, false);
    … …  
    322323            psTrace ("psModules.objects", 4, "%f +/- %f", param->data.F32[n], dparam->data.F32[n]);
    323324        }
    324         if (saveCovariance) {
     325        if (options->saveCovariance) {
    325326            // we only save the covar matrix for this object with itself (ignore cross terms between objects)
    326327            model->covar = psImageAlloc(model->params->n, model->params->n, PS_TYPE_F32);
    … …  
    336337        // save the resulting chisq, nDOF, nIter
    337338        // these are not unique for any one source
    338         model->chisq = myMin->value;
    339         model->nIter = myMin->iter;
    340         model->nDOF  = nPix - model->params->n;
     339        if (options->poissonErrors) {
     340            model->chisq = myMin->value;
     341            model->nPix  = nPix;
     342            model->nDOF  = nPix - model->params->n;
     343            model->chisqNorm = model->chisq / model->nDOF;
     344        } else {
     345            pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, options->covarFactor, model->params->n);
     346        }
     347        model->nIter = myMin->iter;
    341348
    342349        // set the model success or failure status
    … …  
    592599    }
    593600
    594     pmSourceFitSetValues (thisSet, dparams, params, covar, source, myMin, y->n, fitStatus, options->saveCovariance);
     601    pmSourceFitSetValues (thisSet, dparams, params, covar, source, myMin, y->n, fitStatus, options, maskVal);
    595602    psTrace ("psModules.objects", 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nPix: %ld\n", onPic, fitStatus, myMin->iter, myMin->value, y->n);
    596603

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitSet.h

    r30705 r30763  
    4444                           const psVector *dparam, const psVector *param, const psImage *covar,
    4545                           pmSource *source, psMinimization *myMin, int nPix,
    46                            bool fitStatus, bool saveCovariance);
     46                           bool fitStatus, pmSourceFitOptions *options, psImageMaskType maskVal);
    4747
    4848psF32 pmSourceFitSetFunction(psVector *deriv, const psVector *param, const psVector *x);

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c

    r30621 r30763  
    4949
    5050#include "pmSourceIO.h"
     51#include "pmSourceOutputs.h"
    5152
    5253// panstarrs-style FITS table output (header + table in 1st extension)
    … …  
    6263    PS_ASSERT_PTR_NON_NULL(extname, false);
    6364
    64     int i;
    6565    psArray *table;
    6666    psMetadata *row;
    67     psF32 *PAR, *dPAR;
    68     psEllipseAxes axes;
    69     psF32 xPos, yPos;
    70     psF32 xErr, yErr;
    71     psF32 errMag, chisq, apRadius;
    72     psS32 nPix, nDOF;
    7367
    7468    pmChip *chip = readout->parent->parent;
    75     pmFPA  *fpa  = chip->parent;
    76 
    77     bool status1 = false;
    78     bool status2 = false;
    79     float magOffset = NAN;
    80     float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
    81     float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
    82     if (!isfinite(zeropt)) {
    83         zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
    84     }
    85     if (status1 && status2 && (exptime > 0.0)) {
    86         magOffset = zeropt + 2.5*log10(exptime);
    87     }
    88     float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
     69
    8970
    9071    // if the sequence is defined, write these in seq order; otherwise
    … …  
    10081    }
    10182
     83    short nImageOverlap;
     84    float magOffset;
     85    float zeroptErr;
     86    float fwhmMajor;
     87    float fwhmMinor;
     88    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
     89
    10290    table = psArrayAllocEmpty (sources->n);
    10391
    10492    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    10593    // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
    106     for (i = 0; i < sources->n; i++) {
     94    for (int i = 0; i < sources->n; i++) {
    10795        pmSource *source = (pmSource *) sources->data[i];
    10896
    … …  
    115103        }
    116104
    117         // no difference between PSF and non-PSF model
    118         pmModel *model = source->modelPSF;
    119 
    120         if (model != NULL) {
    121             PAR = model->params->data.F32;
    122             dPAR = model->dparams->data.F32;
    123             xPos = PAR[PM_PAR_XPOS];
    124             yPos = PAR[PM_PAR_YPOS];
    125             if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
    126               xErr = dPAR[PM_PAR_XPOS];
    127               yErr = dPAR[PM_PAR_YPOS];
    128             } else {
    129               // in linear-fit mode, there is no error on the centroid
    130               xErr = source->peak->dx;
    131               yErr = source->peak->dy;
    132             }
    133             if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
    134                 axes = pmPSF_ModelToAxes (PAR, 20.0);
    135             } else {
    136                 axes.major = NAN;
    137                 axes.minor = NAN;
    138                 axes.theta = NAN;
    139             }
    140             chisq = model->chisq;
    141             nDOF = model->nDOF;
    142             nPix = model->nPix;
    143             apRadius = source->apRadius;
    144             errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    145         } else {
    146             xPos = source->peak->xf;
    147             yPos = source->peak->yf;
    148             xErr = source->peak->dx;
    149             yErr = source->peak->dy;
    150             axes.major = NAN;
    151             axes.minor = NAN;
    152             axes.theta = NAN;
    153             chisq = NAN;
    154             nDOF = 0;
    155             nPix = 0;
    156             apRadius = NAN;
    157             errMag = NAN;
    158         }
    159 
    160         float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
    161         float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
    162         psS16 nImageOverlap = 1;
    163 
    164         psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
    165         float posAngle = 0.0;
    166         float pltScale = 0.0;
    167         pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
     105        // set the 'best' values for various output fields:
     106        pmSourceOutputs outputs;
     107        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
     108
     109        pmSourceOutputsMoments moments;
     110        pmSourceOutputsSetMoments (&moments, source);
    168111
    169112        pmSourceDiffStats diffStats;
    … …  
    175118        row = psMetadataAlloc ();
    176119        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    177         psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
    178         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
    179         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
    180         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
    181         psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
    182         psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
     120        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
     121        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
     122        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
     123        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
     124        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
     125        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
    183126        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
    184         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        errMag);
     127        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
    185128        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfFlux);
    186129        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfFluxErr);
    187130        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
    188         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
    189         psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
    190         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
     131        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
     132        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
     133        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
    191134        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
    192         psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
    193         psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
     135        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
     136        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
    194137        psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
    195138        psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
    196139
    197         psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
     140        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
    198141        psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
    199142        psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
    200143
    201         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
    202         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
    203         psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
     144        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
     145        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
     146        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
    204147        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeightNotBad);
    205         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
    206         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
    207 
    208         // distinguish moments measure from window vs S/N > XX ??
    209         float mxx = source->moments ? source->moments->Mxx : NAN;
    210         float mxy = source->moments ? source->moments->Mxy : NAN;
    211         float myy = source->moments ? source->moments->Myy : NAN;
    212         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      mxx);
    213         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      mxy);
    214         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      myy);
    215 
    216         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS",        PS_DATA_S32, "nPos (n pix > 3 sigma)",                    diffStats.nGood);
    217         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO",      PS_DATA_F32, "fPos / (fPos + fNeg)",                      diffStats.fRatio);
    218         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD",  PS_DATA_F32, "nPos / (nPos + nNeg)",                      diffStats.nRatioBad);
    219         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)",                     diffStats.nRatioMask);
    220         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL",  PS_DATA_F32, "nPos / (nGood + nMask + nBad)",             diffStats.nRatioAll);
     148        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
     149        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
     150
     151        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
     152        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
     153        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
     154
     155        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS",        PS_DATA_S32, "nPos (n pix > 3 sigma)",                     diffStats.nGood);
     156        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO",      PS_DATA_F32, "fPos / (fPos + fNeg)",                       diffStats.fRatio);
     157        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD",  PS_DATA_F32, "nPos / (nPos + nNeg)",                       diffStats.nRatioBad);
     158        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)",                      diffStats.nRatioMask);
     159        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL",  PS_DATA_F32, "nPos / (nGood + nMask + nBad)",              diffStats.nRatioAll);
    221160
    222161        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c

    r30621 r30763  
    4949
    5050#include "pmSourceIO.h"
     51#include "pmSourceOutputs.h"
    5152
    5253// panstarrs-style FITS table output (header + table in 1st extension)
    … …  
    6263    PS_ASSERT_PTR_NON_NULL(extname, false);
    6364
    64     int i;
    6565    psArray *table;
    6666    psMetadata *row;
    67     psF32 *PAR, *dPAR;
    68     psEllipseAxes axes;
    69     psF32 xPos, yPos;
    70     psF32 xErr, yErr;
    71     psF32 chisq, apRadius;
    72     psS32 nPix, nDOF;
    7367
    7468    pmChip *chip = readout->parent->parent;
    75     pmFPA  *fpa  = chip->parent;
    76 
    77     bool status1 = false;
    78     bool status2 = false;
    79     float magOffset = NAN;
    80     float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
    81     float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
    82     if (!isfinite(zeropt)) {
    83         zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
    84     }
    85     if (status1 && status2 && (exptime > 0.0)) {
    86         magOffset = zeropt + 2.5*log10(exptime);
    87     }
    88     float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
    89 
    90     // we need a measure of the image quality (FWHM) for this image, in order to get the positional errors
    91     float fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MAJ");
    92     if (!status1) {
    93         fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW1");
    94         if (!status1) {
    95             fwhmMajor = 5.0; // XXX just a guess!
    96         }
    97     }
    98     float fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MIN");
    99     if (!status1) {
    100         fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW2");
    101         if (!status1) {
    102             fwhmMinor = 5.0; // XXX just a guess!
    103         }
    104     }
    10569
    10670    // if the sequence is defined, write these in seq order; otherwise
    … …  
    11882    table = psArrayAllocEmpty (sources->n);
    11983
     84    short nImageOverlap;
     85    float magOffset;
     86    float zeroptErr;
     87    float fwhmMajor;
     88    float fwhmMinor;
     89    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
     90
    12091    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    12192    // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
    122     for (i = 0; i < sources->n; i++) {
     93    for (int i = 0; i < sources->n; i++) {
    12394        pmSource *source = (pmSource *) sources->data[i];
    12495
    … …  
    131102        }
    132103
    133         // no difference between PSF and non-PSF model
    134         pmModel *model = source->modelPSF;
    135 
    136         if (model != NULL) {
    137             PAR = model->params->data.F32;
    138             dPAR = model->dparams->data.F32;
    139             xPos = PAR[PM_PAR_XPOS];
    140             yPos = PAR[PM_PAR_YPOS];
    141             if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
    142               xErr = dPAR[PM_PAR_XPOS];
    143               yErr = dPAR[PM_PAR_YPOS];
    144             } else {
    145               // in linear-fit mode, there is no error on the centroid
    146                 xErr = fwhmMajor * source->errMag / 2.35;
    147                 yErr = fwhmMinor * source->errMag / 2.35;
    148             }
    149             if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
    150                 axes = pmPSF_ModelToAxes (PAR, 20.0);
    151             } else {
    152                 axes.major = NAN;
    153                 axes.minor = NAN;
    154                 axes.theta = NAN;
    155             }
    156             chisq = model->chisq;
    157             nDOF = model->nDOF;
    158             nPix = model->nPix;
    159             apRadius = source->apRadius;
    160         } else {
    161             bool useMoments = true;
    162             useMoments = (useMoments && source->moments);          // can't if there are no moments
    163             useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured
    164             useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed...
    165 
    166             if (useMoments) {
    167                 xPos = source->moments->Mx;
    168                 yPos = source->moments->My;
    169                 xErr = fwhmMajor * source->errMag / 2.35;
    170                 yErr = fwhmMinor * source->errMag / 2.35;
    171             } else {
    172                 xPos = source->peak->xf;
    173                 yPos = source->peak->yf;
    174                 xErr = source->peak->dx;
    175                 yErr = source->peak->dy;
    176             }
    177             axes.major = NAN;
    178             axes.minor = NAN;
    179             axes.theta = NAN;
    180             chisq = NAN;
    181             nDOF = 0;
    182             nPix = 0;
    183             apRadius = NAN;
    184         }
    185 
    186         float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
    187         float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
    188         psS16 nImageOverlap = 1;
    189 
    190         psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
    191         float posAngle = 0.0;
    192         float pltScale = 0.0;
    193         pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
     104        // set the 'best' values for various output fields:
     105        pmSourceOutputs outputs;
     106        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
     107
     108        pmSourceOutputsMoments moments;
     109        pmSourceOutputsSetMoments (&moments, source);
    194110
    195111        pmSourceDiffStats diffStats;
    … …  
    201117        row = psMetadataAlloc ();
    202118        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    203         psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
    204         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
    205         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
    206         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
    207         psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
    208         psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
     119        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
     120        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
     121        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
     122        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
     123        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
     124        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
    209125        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
    210126        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
    … …  
    214130        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
    215131        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in real aperture",                 source->apMagRaw);
    216         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
     132        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
    217133        psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",     source->apFlux);
    218134        psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",                        source->apFluxErr);
    219135
    220         psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
    221         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
     136        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
     137        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
    222138        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
    223         psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
    224         psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
     139        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
     140        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
    225141        psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
    226142        psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
    227143
    228         psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
     144        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
    229145        psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
    230146        psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
    231147
    232         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
    233         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
    234         psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
     148        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
     149        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
     150        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
    235151        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
    236152        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
    237         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
    238         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
    239 
    240         // distinguish moments measure from window vs S/N > XX ??
    241         float mxx = source->moments ? source->moments->Mxx : NAN;
    242         float mxy = source->moments ? source->moments->Mxy : NAN;
    243         float myy = source->moments ? source->moments->Myy : NAN;
    244         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      mxx);
    245         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      mxy);
    246         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      myy);
    247 
    248         float Mrf  = source->moments ? source->moments->Mrf : NAN;
    249         float Mrh  = source->moments ? source->moments->Mrh : NAN;
    250         float Krf  = source->moments ? source->moments->KronFlux : NAN;
    251         float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
    252 
    253         float Kinner = source->moments ? source->moments->KronFinner : NAN;
    254         float Kouter = source->moments ? source->moments->KronFouter : NAN;
    255 
    256         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                       Mrf);
    257         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                        Mrh);
    258         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Krf);
    259         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                          dKrf);
    260 
    261         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 1.0 R1)",                     Kinner);
    262         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 4.0 R1)",                     Kouter);
    263 
    264         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS",        PS_DATA_S32, "nPos (n pix > 3 sigma)",                    diffStats.nGood);
    265         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO",      PS_DATA_F32, "fPos / (fPos + fNeg)",                      diffStats.fRatio);
    266         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD",  PS_DATA_F32, "nPos / (nPos + nNeg)",                      diffStats.nRatioBad);
    267         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)",                     diffStats.nRatioMask);
    268         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL",  PS_DATA_F32, "nPos / (nGood + nMask + nBad)",             diffStats.nRatioAll);
     153        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
     154        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
     155
     156        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
     157        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
     158        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
     159
     160        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
     161        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                         moments.Mrh);
     162        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Krf);
     163        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                            moments.dKrf);
     164        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 1.0 R1)",                      moments.Kinner);
     165        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 4.0 R1)",                      moments.Kouter);
     166
     167        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS",        PS_DATA_S32, "nPos (n pix > 3 sigma)",                     diffStats.nGood);
     168        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO",      PS_DATA_F32, "fPos / (fPos + fNeg)",                       diffStats.fRatio);
     169        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD",  PS_DATA_F32, "nPos / (nPos + nNeg)",                       diffStats.nRatioBad);
     170        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)",                      diffStats.nRatioMask);
     171        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL",  PS_DATA_F32, "nPos / (nGood + nMask + nBad)",              diffStats.nRatioAll);
    269172
    270173        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P",         PS_DATA_F32, "distance to positive match source",          diffStats.Rp);

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c

    r30704 r30763  
    5050
    5151#include "pmSourceIO.h"
     52#include "pmSourceOutputs.h"
    5253
    5354// panstarrs-style FITS table output (header + table in 1st extension)
    … …  
    6667    psArray *table;
    6768    psMetadata *row;
    68     psF32 *PAR, *dPAR;
    69     psEllipseAxes axes;
    70     psF32 xPos, yPos;
    71     psF32 xErr, yErr;
    72     psF32 errMag, chisq, apRadius;
    73     psS32 nPix, nDOF;
    7469
    7570    pmChip *chip = readout->parent->parent;
    76     pmFPA  *fpa  = chip->parent;
    77 
    78     bool status1 = false;
    79     bool status2 = false;
    80     float magOffset = NAN;
    81     float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
    82     float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
    83     if (!isfinite(zeropt)) {
    84         zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
    85     }
    86     if (status1 && status2 && (exptime > 0.0)) {
    87         magOffset = zeropt + 2.5*log10(exptime);
    88     }
    89     float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
    9071
    9172    // if the sequence is defined, write these in seq order; otherwise
    … …  
    10384    table = psArrayAllocEmpty (sources->n);
    10485
    105 # if (0)
    106     // we use this just to define the output vectors (which must be present for all objects)
    107     bool status = false;
    108     psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    109     psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
    110     psAssert (radMax, "this must have been defined and tested earlier!");
    111     psAssert (radMax->n, "this must have been defined and tested earlier!");
    112     psAssert (radMin->n == radMax->n, "inconsistent annular bins");
    113 
    114     // write the radial profile apertures to header
    115     for (int i = 0; i < radMax->n; i++) {
    116       sprintf (keyword1, "RMIN_%02d", i);
    117       sprintf (keyword2, "RMAX_%02d", i);
    118       psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
    119       psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
    120     }
    121 # endif
     86    short nImageOverlap;
     87    float magOffset;
     88    float zeroptErr;
     89    float fwhmMajor;
     90    float fwhmMinor;
     91    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
    12292
    12393    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    … …  
    134104        }
    135105
    136         // no difference between PSF and non-PSF model
    137         pmModel *model = source->modelPSF;
    138 
    139         if (model != NULL) {
    140             PAR = model->params->data.F32;
    141             dPAR = model->dparams->data.F32;
    142             xPos = PAR[PM_PAR_XPOS];
    143             yPos = PAR[PM_PAR_YPOS];
    144             if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
    145                 xErr = dPAR[PM_PAR_XPOS];
    146                 yErr = dPAR[PM_PAR_YPOS];
    147             } else {
    148                 // in linear-fit mode, there is no error on the centroid
    149                 xErr = source->peak->dx;
    150                 yErr = source->peak->dy;
    151             }
    152             if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
    153                 axes = pmPSF_ModelToAxes (PAR, 20.0);
    154             } else {
    155                 axes.major = NAN;
    156                 axes.minor = NAN;
    157                 axes.theta = NAN;
    158             }
    159             chisq = model->chisq;
    160             nDOF = model->nDOF;
    161             nPix = model->nPix;
    162             apRadius = source->apRadius;
    163             errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    164         } else {
    165             xPos = source->peak->xf;
    166             yPos = source->peak->yf;
    167             xErr = source->peak->dx;
    168             yErr = source->peak->dy;
    169             axes.major = NAN;
    170             axes.minor = NAN;
    171             axes.theta = NAN;
    172             chisq = NAN;
    173             nDOF = 0;
    174             nPix = 0;
    175             apRadius = NAN;
    176             errMag = NAN;
    177         }
    178 
    179         float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
    180         float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
    181         psS16 nImageOverlap = 1;
    182 
    183         psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
    184         float posAngle = 0.0;
    185         float pltScale = 0.0;
    186         pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
     106        // set the 'best' values for various output fields:
     107        pmSourceOutputs outputs;
     108        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
     109
     110        pmSourceOutputsMoments moments;
     111        pmSourceOutputsSetMoments (&moments, source);
    187112
    188113        row = psMetadataAlloc ();
    189114        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    190         psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
    191         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
    192         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
    193         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
    194         psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
    195         psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
     115        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
     116        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
     117        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr);
     118        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr);
     119        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
     120        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
    196121        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
    197         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        errMag);
     122        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
    198123        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental flux (counts)",         source->psfFlux);
    199124        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux",             source->psfFluxErr);
    200125        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
    201126        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
    202         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
    203         psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
    204         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
     127        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
     128        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
     129        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
    205130        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
    206         psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
    207         psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
     131        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
     132        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
    208133        psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
    209134        psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
    210135
    211         psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
     136        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
    212137        psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
    213138        psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
    214139
    215         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
    216         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
    217         psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
     140        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
     141        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
     142        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
    218143        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
    219144        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
    220         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
    221         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
    222 
    223         // distinguish moments measure from window vs S/N > XX ??
    224         float Mxx = source->moments ? source->moments->Mxx : NAN;
    225         float Mxy = source->moments ? source->moments->Mxy : NAN;
    226         float Myy = source->moments ? source->moments->Myy : NAN;
    227 
    228         float Mrf  = source->moments ? source->moments->Mrf : NAN;
    229         float Mrh  = source->moments ? source->moments->Mrh : NAN;
    230         float Krf  = source->moments ? source->moments->KronFlux : NAN;
    231         float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
    232 
    233         float Kinner = source->moments ? source->moments->KronFinner : NAN;
    234         float Kouter = source->moments ? source->moments->KronFouter : NAN;
    235 
    236         float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN;
    237         float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN;
    238         float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN;
    239         float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN;
    240 
    241         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      Mxx);
    242         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      Mxy);
    243         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      Myy);
    244 
    245         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                    M_c3);
    246         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                    M_s3);
    247         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                   M_c4);
    248         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                   M_s4);
    249 
    250         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                       Mrf);
    251         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                        Mrh);
    252         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Krf);
    253         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                          dKrf);
    254 
    255         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kinner);
    256         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kouter);
     145        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
     146        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
     147
     148        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
     149        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
     150        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
     151
     152        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
     153        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
     154        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
     155        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
     156
     157        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
     158        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                         moments.Mrh);
     159        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Krf);
     160        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                            moments.dKrf);
     161        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kinner);
     162        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kouter);
    257163
    258164        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                     source->mode);
    259165        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                     source->mode2);
    260 
    261 # if (0)
    262         // XXX if we have raw radial apertures, write them out here
    263         psVector *radFlux    = psVectorAlloc(radMax->n, PS_TYPE_F32);
    264         psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32);
    265         psVector *radFill    = psVectorAlloc(radMax->n, PS_TYPE_F32);
    266         psVectorInit (radFlux,    NAN);
    267         psVectorInit (radFluxErr, NAN);
    268         psVectorInit (radFill,    NAN);
    269         if (!source->radial) goto empty_annuli;
    270         if (!source->radial->flux) goto empty_annuli;
    271         if (!source->radial->fill) goto empty_annuli;
    272         psAssert (source->radial->flux->n <= radFlux->n, "inconsistent vector lengths");
    273         psAssert (source->radial->fill->n <= radFlux->n, "inconsistent vector lengths");
    274 
    275         // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
    276         for (int j = 0; j < source->radial->flux->n; j++) {
    277             radFlux->data.F32[j]    = source->radial->flux->data.F32[j];
    278             radFluxErr->data.F32[j] = source->radial->fluxErr->data.F32[j];
    279             radFill->data.F32[j]    = source->radial->fill->data.F32[j];
    280         }
    281 
    282     empty_annuli:
    283         psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX",     PS_DATA_VECTOR, "flux within annuli",    radFlux);
    284         psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli",  radFluxErr);
    285         psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL",     PS_DATA_VECTOR, "fill factor of annuli", radFill);
    286         psFree (radFlux);
    287         psFree (radFluxErr);
    288         psFree (radFill);
    289 # endif
    290166
    291167        // XXX not sure how to get this : need to load Nimages with weight?

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c

    r30621 r30763  
    4949
    5050#include "pmSourceIO.h"
     51#include "pmSourceOutputs.h"
    5152
    5253// panstarrs-style FITS table output (header + table in 1st extension)
    … …  
    6263    psArray *table;
    6364    psMetadata *row;
    64     int i;
    65     psF32 *PAR, *dPAR;
    66     psEllipseAxes axes;
    67     psF32 xPos, yPos;
    68     psF32 xErr, yErr;
    69     psF32 errMag, chisq, apRadius;
    70     psS32 nPix, nDOF;
    7165
    7266    pmChip *chip = readout->parent->parent;
    73     pmFPA  *fpa  = chip->parent;
    74 
    75     bool status1 = false;
    76     bool status2 = false;
    77     float magOffset = NAN;
    78     float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
    79     float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
    80     if (!isfinite(zeropt)) {
    81         zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
    82     }
    83     if (status1 && status2 && (exptime > 0.0)) {
    84         magOffset = zeropt + 2.5*log10(exptime);
    85     }
    86     float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
    8767
    8868    // if the sequence is defined, write these in seq order; otherwise
    … …  
    10080    table = psArrayAllocEmpty (sources->n);
    10181
     82    short nImageOverlap;
     83    float magOffset;
     84    float zeroptErr;
     85    float fwhmMajor;
     86    float fwhmMinor;
     87    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
     88
    10289    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    10390    // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
    104     for (i = 0; i < sources->n; i++) {
     91    for (int i = 0; i < sources->n; i++) {
    10592        pmSource *source = (pmSource *) sources->data[i];
    10693
    … …  
    113100        }
    114101
    115         // no difference between PSF and non-PSF model
    116         pmModel *model = source->modelPSF;
    117 
    118         if (model != NULL) {
    119             PAR = model->params->data.F32;
    120             dPAR = model->dparams->data.F32;
    121             xPos = PAR[PM_PAR_XPOS];
    122             yPos = PAR[PM_PAR_YPOS];
    123             if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
    124               xErr = dPAR[PM_PAR_XPOS];
    125               yErr = dPAR[PM_PAR_YPOS];
    126             } else {
    127               // in linear-fit mode, there is no error on the centroid
    128               xErr = source->peak->dx;
    129               yErr = source->peak->dy;
    130             }
    131             if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
    132                 axes = pmPSF_ModelToAxes (PAR, 20.0);
    133             } else {
    134                 axes.major = NAN;
    135                 axes.minor = NAN;
    136                 axes.theta = NAN;
    137             }
    138             chisq = model->chisq;
    139             nDOF = model->nDOF;
    140             nPix = model->nPix;
    141             apRadius = source->apRadius;
    142             errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    143         } else {
    144             xPos = source->peak->xf;
    145             yPos = source->peak->yf;
    146             xErr = source->peak->dx;
    147             yErr = source->peak->dy;
    148             axes.major = NAN;
    149             axes.minor = NAN;
    150             axes.theta = NAN;
    151             chisq = NAN;
    152             nDOF = 0;
    153             nPix = 0;
    154             apRadius = NAN;
    155             errMag = NAN;
    156         }
    157 
    158         float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
    159         float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
    160         psS16 nImageOverlap = 1;
    161 
    162         psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
    163         float posAngle = 0.0;
    164         float pltScale = 0.0;
    165         pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
     102        // set the 'best' values for various output fields:
     103        pmSourceOutputs outputs;
     104        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
     105
     106        pmSourceOutputsMoments moments;
     107        pmSourceOutputsSetMoments (&moments, source);
    166108
    167109        row = psMetadataAlloc ();
    168110        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    169         psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
    170         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
    171         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
    172         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
    173         psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F32, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
    174         psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F32, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
    175         psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
    176         psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
     111        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
     112        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
     113        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
     114        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
     115        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F32, "PSF RA coordinate (degrees)",                outputs.ra);
     116        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F32, "PSF DEC coordinate (degrees)",               outputs.dec);
     117        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
     118        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
    177119        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
    178         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        errMag);
     120        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
    179121        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
    180         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
    181         psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
    182         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
     122        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
     123        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
     124        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
    183125        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
    184126        psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
    185127        psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
    186128
    187         psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
     129        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
    188130        psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
    189131        psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
    190132
    191         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
    192         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
    193         psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
     133        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
     134        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
     135        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
    194136        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeightNotBad);
    195         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
    196         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
    197 
    198         // distinguish moments measure from window vs S/N > XX ??
    199         float mxx = source->moments ? source->moments->Mxx : NAN;
    200         float mxy = source->moments ? source->moments->Mxy : NAN;
    201         float myy = source->moments ? source->moments->Myy : NAN;
    202         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      mxx);
    203         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      mxy);
    204         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      myy);
     137        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
     138        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
     139
     140        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
     141        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
     142        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
    205143
    206144        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c

    r30621 r30763  
    4949
    5050#include "pmSourceIO.h"
     51#include "pmSourceOutputs.h"
    5152
    5253// panstarrs-style FITS table output (header + table in 1st extension)
    … …  
    6263    psArray *table;
    6364    psMetadata *row;
    64     psF32 *PAR, *dPAR;
    65     psEllipseAxes axes;
    66     psF32 xPos, yPos;
    67     psF32 xErr, yErr;
    68     psF32 errMag, chisq, apRadius;
    69     psS32 nPix, nDOF;
    7065
    7166    pmChip *chip = readout->parent->parent;
    72     pmFPA  *fpa  = chip->parent;
    73 
    74     bool status1 = false;
    75     bool status2 = false;
    76     float magOffset = NAN;
    77     float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
    78     float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
    79     if (!isfinite(zeropt)) {
    80         zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
    81     }
    82     if (status1 && status2 && (exptime > 0.0)) {
    83         magOffset = zeropt + 2.5*log10(exptime);
    84     }
    85     float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
    8667
    8768    // if the sequence is defined, write these in seq order; otherwise
    … …  
    9980    table = psArrayAllocEmpty (sources->n);
    10081
     82    short nImageOverlap;
     83    float magOffset;
     84    float zeroptErr;
     85    float fwhmMajor;
     86    float fwhmMinor;
     87    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
     88
    10189    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    10290    // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
    … …  
    112100        }
    113101
    114         // no difference between PSF and non-PSF model
    115         pmModel *model = source->modelPSF;
    116 
    117         if (model != NULL) {
    118             PAR = model->params->data.F32;
    119             dPAR = model->dparams->data.F32;
    120             xPos = PAR[PM_PAR_XPOS];
    121             yPos = PAR[PM_PAR_YPOS];
    122             if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
    123                 xErr = dPAR[PM_PAR_XPOS];
    124                 yErr = dPAR[PM_PAR_YPOS];
    125             } else {
    126                 // in linear-fit mode, there is no error on the centroid
    127                 xErr = source->peak->dx;
    128                 yErr = source->peak->dy;
    129             }
    130             if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
    131                 axes = pmPSF_ModelToAxes (PAR, 20.0);
    132             } else {
    133                 axes.major = NAN;
    134                 axes.minor = NAN;
    135                 axes.theta = NAN;
    136             }
    137             chisq = model->chisq;
    138             nDOF = model->nDOF;
    139             nPix = model->nPix;
    140             apRadius = source->apRadius;
    141             errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    142         } else {
    143             xPos = source->peak->xf;
    144             yPos = source->peak->yf;
    145             xErr = source->peak->dx;
    146             yErr = source->peak->dy;
    147             axes.major = NAN;
    148             axes.minor = NAN;
    149             axes.theta = NAN;
    150             chisq = NAN;
    151             nDOF = 0;
    152             nPix = 0;
    153             apRadius = NAN;
    154             errMag = NAN;
    155         }
    156 
    157         float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
    158         float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
    159         psS16 nImageOverlap = 1;
    160 
    161         psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
    162         float posAngle = 0.0;
    163         float pltScale = 0.0;
    164         pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
     102        // set the 'best' values for various output fields:
     103        pmSourceOutputs outputs;
     104        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
     105
     106        pmSourceOutputsMoments moments;
     107        pmSourceOutputsSetMoments (&moments, source);
    165108
    166109        row = psMetadataAlloc ();
    167110        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    168         psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
    169         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
    170         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
    171         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
    172         psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
    173         psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
     111        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
     112        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
     113        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
     114        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
     115        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
     116        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
    174117        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
    175         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        errMag);
     118        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
    176119        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
    177         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
    178         psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
    179         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
     120        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
     121        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
     122        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
    180123        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
    181         psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
    182         psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
     124        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
     125        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
    183126        psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
    184127        psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
    185128
    186         psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
     129        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
    187130        psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
    188131        psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
    189132
    190         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
    191         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
    192         psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
     133        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
     134        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
     135        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
    193136        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeightNotBad);
    194         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
    195         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
    196 
    197         // distinguish moments measure from window vs S/N > XX ??
    198         float mxx = source->moments ? source->moments->Mxx : NAN;
    199         float mxy = source->moments ? source->moments->Mxy : NAN;
    200         float myy = source->moments ? source->moments->Myy : NAN;
    201         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      mxx);
    202         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      mxy);
    203         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      myy);
     137        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
     138        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
     139
     140        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
     141        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
     142        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
    204143
    205144        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c

    r30621 r30763  
    4949
    5050#include "pmSourceIO.h"
     51#include "pmSourceOutputs.h"
    5152
    5253// panstarrs-style FITS table output (header + table in 1st extension)
    … …  
    6263    psArray *table;
    6364    psMetadata *row;
    64     psF32 *PAR, *dPAR;
    65     psEllipseAxes axes;
    66     psF32 xPos, yPos;
    67     psF32 xErr, yErr;
    68     psF32 chisq, apRadius;
    69     psS32 nPix, nDOF;
    7065
    7166    pmChip *chip = readout->parent->parent;
    72     pmFPA  *fpa  = chip->parent;
    73 
    74     bool status1 = false;
    75     bool status2 = false;
    76     float magOffset = NAN;
    77     float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
    78     float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
    79     if (!isfinite(zeropt)) {
    80         zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
    81     }
    82     if (status1 && status2 && (exptime > 0.0)) {
    83         magOffset = zeropt + 2.5*log10(exptime);
    84     }
    85     float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
    86 
    87     // we need a measure of the image quality (FWHM) for this image, in order to get the positional errors
    88     float fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MAJ");
    89     if (!status1) {
    90         fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW1");
    91         if (!status1) {
    92             fwhmMajor = 5.0; // XXX just a guess!
    93         }
    94     }
    95     float fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MIN");
    96     if (!status1) {
    97         fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW2");
    98         if (!status1) {
    99             fwhmMinor = 5.0; // XXX just a guess!
    100         }
    101     }
    10267
    10368    // if the sequence is defined, write these in seq order; otherwise
    … …  
    11580    table = psArrayAllocEmpty (sources->n);
    11681
     82    short nImageOverlap;
     83    float magOffset;
     84    float zeroptErr;
     85    float fwhmMajor;
     86    float fwhmMinor;
     87    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
     88
    11789    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    11890    // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
    … …  
    128100        }
    129101
    130         // no difference between PSF and non-PSF model
    131         pmModel *model = source->modelPSF;
    132 
    133         if (model != NULL) {
    134             PAR = model->params->data.F32;
    135             dPAR = model->dparams->data.F32;
    136             xPos = PAR[PM_PAR_XPOS];
    137             yPos = PAR[PM_PAR_YPOS];
    138             if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
    139                 xErr = dPAR[PM_PAR_XPOS];
    140                 yErr = dPAR[PM_PAR_YPOS];
    141             } else {
    142                 xErr = fwhmMajor * source->errMag / 2.35;
    143                 yErr = fwhmMinor * source->errMag / 2.35;
    144             }
    145             if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
    146                 axes = pmPSF_ModelToAxes (PAR, 20.0);
    147             } else {
    148                 axes.major = NAN;
    149                 axes.minor = NAN;
    150                 axes.theta = NAN;
    151             }
    152             chisq = model->chisq;
    153             nDOF = model->nDOF;
    154             nPix = model->nPix;
    155             apRadius = source->apRadius;
    156         } else {
    157             bool useMoments = true;
    158             useMoments = (useMoments && source->moments);          // can't if there are no moments
    159             useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured
    160             useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed...
    161 
    162             if (useMoments) {
    163                 xPos = source->moments->Mx;
    164                 yPos = source->moments->My;
    165                 xErr = fwhmMajor * source->errMag / 2.35;
    166                 yErr = fwhmMinor * source->errMag / 2.35;
    167             } else {
    168                 xPos = source->peak->xf;
    169                 yPos = source->peak->yf;
    170                 xErr = source->peak->dx;
    171                 yErr = source->peak->dy;
    172             }
    173             axes.major = NAN;
    174             axes.minor = NAN;
    175             axes.theta = NAN;
    176             chisq = NAN;
    177             nDOF = 0;
    178             nPix = 0;
    179             apRadius = NAN;
    180         }
    181 
    182         float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
    183         float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
    184         psS16 nImageOverlap = 1;
    185 
    186         psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
    187         float posAngle = 0.0;
    188         float pltScale = 0.0;
    189         pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
     102        // set the 'best' values for various output fields:
     103        pmSourceOutputs outputs;
     104        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
     105
     106        pmSourceOutputsMoments moments;
     107        pmSourceOutputsSetMoments (&moments, source);
    190108
    191109        row = psMetadataAlloc ();
    192110        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    193         psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
    194         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
    195         psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  xErr);
    196         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr);
    197         psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
    198         psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
     111        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
     112        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
     113        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr);
     114        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr);
     115        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
     116        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
    199117        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
    200118        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
    … …  
    203121        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
    204122        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
    205         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
    206 
    207         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
     123        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
     124
     125        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
    208126        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
    209127       
    210128        // NOTE: RA & DEC (both double) need to be on an 8-byte boundary...
    211         psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
    212         psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
    213 
    214         psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
     129        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
     130        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
     131
     132        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
    215133        psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
    216134        psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
    217135
    218         psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
     136        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
    219137        psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
    220138        psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
    221139
    222         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
    223         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
    224         psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
     140        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
     141        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
     142        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
    225143        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
    226144        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
    227         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
    228         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
    229 
    230         // distinguish moments measure from window vs S/N > XX ??
    231         float Mxx = source->moments ? source->moments->Mxx : NAN;
    232         float Mxy = source->moments ? source->moments->Mxy : NAN;
    233         float Myy = source->moments ? source->moments->Myy : NAN;
    234 
    235         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      Mxx);
    236         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      Mxy);
    237         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      Myy);
    238 
    239         float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN;
    240         float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN;
    241         float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN;
    242         float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN;
    243 
    244         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                    M_c3);
    245         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                    M_s3);
    246         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                   M_c4);
    247         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                   M_s4);
    248 
    249         float Mrf  = source->moments ? source->moments->Mrf : NAN;
    250         float Mrh  = source->moments ? source->moments->Mrh : NAN;
    251         float Krf  = source->moments ? source->moments->KronFlux : NAN;
    252         float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
    253 
    254         float Kinner = source->moments ? source->moments->KronFinner : NAN;
    255         float Kouter = source->moments ? source->moments->KronFouter : NAN;
    256 
    257         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                       Mrf);
    258         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                        Mrh);
    259         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Krf);
    260         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                          dKrf);
    261 
    262         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kinner);
    263         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kouter);
    264 
    265         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                     source->mode);
    266         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                     source->mode2);
     145        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
     146        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
     147
     148        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
     149        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
     150        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
     151
     152        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
     153        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
     154        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
     155        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
     156
     157        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
     158        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                         moments.Mrh);
     159        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Krf);
     160        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                            moments.dKrf);
     161        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kinner);
     162        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kouter);
     163        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
     164        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                      source->mode2);
    267165        psMetadataAdd (row, PS_LIST_TAIL, "PADDING2",         PS_DATA_S32, "more padding", 0);
    268166

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.c

    r29935 r30763  
    422422    *pixWeightNotPoor = notPoorSum / modelSum;
    423423
     424    if (false && isfinite(*pixWeightNotBad) && isfinite(*pixWeightNotPoor)) {
     425        psAssert (*pixWeightNotBad <= *pixWeightNotPoor, "error: all bad pixels should also be poor");
     426    }
     427
    424428    return (true);
    425429}
    … …  
    674678
    675679// determine chisq, etc for linear normalization-only fit
    676 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal, const float covarFactor)
     680bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal, const float covarFactor, int nParams)
    677681{
    678682    PS_ASSERT_PTR_NON_NULL(model, false);
    … …  
    689693            if (variance->data.F32[j][i] <= 0)
    690694                continue;
    691             dC += PS_SQR (image->data.F32[j][i]) / (covarFactor * variance->data.F32[j][i]);
     695            // dC += PS_SQR (image->data.F32[j][i]) / (covarFactor * variance->data.F32[j][i]);
     696            dC += PS_SQR (image->data.F32[j][i]) / variance->data.F32[j][i];
    692697            Npix ++;
    693698        }
    694699    }
     700
    695701    model->nPix = Npix;
    696     model->nDOF = Npix - 1;
     702    model->nDOF = Npix - nParams - 1;
    697703    model->chisq = dC;
     704    model->chisqNorm = dC / model->nDOF;
    698705
    699706    return (true);

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.h

    r29546 r30763  
    6969bool pmSourcePixelWeight (float *pixWeightNotBad, float *pixWeightNotPoor, pmModel *model, psImage *mask, psImageMaskType maskVal, psImageMaskType markVal);
    7070
    71 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal, const float covarFactor);
     71bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal, const float covarFactor, int nParams);
    7272
    7373bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
Note: See TracChangeset for help on using the changeset viewer.