IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 21528


Ignore:
Timestamp:
Feb 17, 2009, 4:44:19 PM (17 years ago)
Author:
Paul Price
Message:

Catching case of !source->moments.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c

    r21516 r21528  
    33 *  @author EAM, IfA
    44 *
    5  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2009-02-16 22:30:50 $
     5 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2009-02-18 02:44:19 $
    77 *
    88 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4545// followed by a zero-size matrix, followed by the table data
    4646
    47 bool pmSourcesWrite_CMF_PS1_V1 (psFits *fits, pmReadout *readout, psArray *sources, 
    48                                 psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
     47bool pmSourcesWrite_CMF_PS1_V1 (psFits *fits, pmReadout *readout, psArray *sources,
     48                                psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
    4949{
    5050    PS_ASSERT_PTR_NON_NULL(fits, false);
    … …  
    7272    float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
    7373    if (status1 && status2 && (exptime > 0.0)) {
    74         magOffset = zeropt + 2.5*log10(exptime);
     74        magOffset = zeropt + 2.5*log10(exptime);
    7575    }
    7676
    … …  
    7979    if (sources->n > 0) {
    8080        pmSource *source = (pmSource *) sources->data[0];
    81         if (source->seq == -1) {
    82           // let's write these out in S/N order
    83           sources = psArraySort (sources, pmSourceSortBySN);
    84         } else {
    85           sources = psArraySort (sources, pmSourceSortBySeq);
    86         }
     81        if (source->seq == -1) {
     82          // let's write these out in S/N order
     83          sources = psArraySort (sources, pmSourceSortBySN);
     84        } else {
     85          sources = psArraySort (sources, pmSourceSortBySeq);
     86        }
    8787    }
    8888
    … …  
    9494        pmSource *source = (pmSource *) sources->data[i];
    9595
    96         // If source->seq is -1, source was generated in this analysis.  If source->seq is
    97         // not -1, source was read from elsewhere: in the latter case, preserve the source
    98         // ID.  source.seq is used instead of source.id since the latter is a const
    99         // generated on Alloc, and would thus be wrong for read in sources.
    100         if (source->seq == -1) {
    101           source->seq = i;
    102         }
     96        // If source->seq is -1, source was generated in this analysis.  If source->seq is
     97        // not -1, source was read from elsewhere: in the latter case, preserve the source
     98        // ID.  source.seq is used instead of source.id since the latter is a const
     99        // generated on Alloc, and would thus be wrong for read in sources.
     100        if (source->seq == -1) {
     101          source->seq = i;
     102        }
    103103
    104104        // no difference between PSF and non-PSF model
    … …  
    110110            xPos = PAR[PM_PAR_XPOS];
    111111            yPos = PAR[PM_PAR_YPOS];
    112             if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
    113               xErr = dPAR[PM_PAR_XPOS];
    114               yErr = dPAR[PM_PAR_YPOS];
    115             } else {
    116               // in linear-fit mode, there is no error on the centroid
    117               xErr = source->peak->dx;
    118               yErr = source->peak->dy;
    119             }         
    120             if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
    121                 axes = pmPSF_ModelToAxes (PAR, 20.0);
    122             } else {
    123                 axes.major = NAN;
    124                 axes.minor = NAN;
    125                 axes.theta = NAN;
    126             }
    127             chisq = model->chisq;
    128             nDOF = model->nDOF;
    129             nPix = model->nPix;
    130             apRadius = model->radiusFit; // XXX should we really use the fitRadius for aperture Radius?
    131             errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
     112            if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
     113              xErr = dPAR[PM_PAR_XPOS];
     114              yErr = dPAR[PM_PAR_YPOS];
     115            } else {
     116              // in linear-fit mode, there is no error on the centroid
     117              xErr = source->peak->dx;
     118              yErr = source->peak->dy;
     119            }
     120            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
     121                axes = pmPSF_ModelToAxes (PAR, 20.0);
     122            } else {
     123                axes.major = NAN;
     124                axes.minor = NAN;
     125                axes.theta = NAN;
     126            }
     127            chisq = model->chisq;
     128            nDOF = model->nDOF;
     129            nPix = model->nPix;
     130            apRadius = model->radiusFit; // XXX should we really use the fitRadius for aperture Radius?
     131            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    132132        } else {
    133133            xPos = source->peak->xf;
    … …  
    138138            axes.minor = NAN;
    139139            axes.theta = NAN;
    140             chisq = NAN;
    141             nDOF = 0;
    142             nPix = 0;
    143             apRadius = NAN;
    144             errMag = NAN;
    145         }
    146 
    147         float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
     140            chisq = NAN;
     141            nDOF = 0;
     142            nPix = 0;
     143            apRadius = NAN;
     144            errMag = NAN;
     145        }
     146
     147        float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
    148148        float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
    149149        psS16 nImageOverlap = 1;
    150150
    151         psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
    152         float posAngle = 0.0;
    153         float pltScale = 0.0;
    154         pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
     151        psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
     152        float posAngle = 0.0;
     153        float pltScale = 0.0;
     154        pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
    155155
    156156        row = psMetadataAlloc ();
    … …  
    185185        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
    186186
    187         // distinguish moments measure from window vs S/N > XX ??
    188         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      source->moments->Mxx);
    189         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      source->moments->Mxy);
    190         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      source->moments->Myy);
     187        // distinguish moments measure from window vs S/N > XX ??
     188        float mxx = source->moments ? source->moments->Mxx : NAN;
     189        float mxy = source->moments ? source->moments->Mxy : NAN;
     190        float myy = source->moments ? source->moments->Myy : NAN;
     191        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      mxx);
     192        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      mxy);
     193        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      myy);
    191194
    192195        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
    … …  
    199202        psFree (row);
    200203
    201         // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not
    202         // subtracted
    203 
    204         // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image
    205         // edge; 3) any pixels in the 3x3 peak region are masked;
     204        // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not
     205        // subtracted
     206
     207        // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image
     208        // edge; 3) any pixels in the 3x3 peak region are masked;
    206209    }
    207210
    … …  
    281284        source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
    282285
    283         // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
    284         PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
    285         dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
     286        // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
     287        PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
     288        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
    286289
    287290        pmPSF_AxesToModel (PAR, axes);
    … …  
    290293        float peakFlux    = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN;
    291294
    292         // recreate the peak to match (xPos, yPos) +/- (xErr, yErr)
     295        // recreate the peak to match (xPos, yPos) +/- (xErr, yErr)
    293296        source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
    294297        source->peak->flux = peakFlux;
    … …  
    300303        source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
    301304
    302         // note that some older versions used PSF_PROBABILITY: this was not well defined.
    303         model->chisq      = psMetadataLookupF32 (&status, row, "PSF_CHISQ");
    304         model->nDOF       = psMetadataLookupS32 (&status, row, "PSF_NDOF");
    305         model->nPix       = psMetadataLookupS32 (&status, row, "PSF_NPIX");
     305        // note that some older versions used PSF_PROBABILITY: this was not well defined.
     306        model->chisq      = psMetadataLookupF32 (&status, row, "PSF_CHISQ");
     307        model->nDOF       = psMetadataLookupS32 (&status, row, "PSF_NDOF");
     308        model->nPix       = psMetadataLookupS32 (&status, row, "PSF_NPIX");
    306309        model->radiusFit  = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS");
    307310
    308         source->moments = pmMomentsAlloc ();
     311        source->moments = pmMomentsAlloc ();
    309312        source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX");
    310313        source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY");
    … …  
    321324}
    322325
    323 // XXX this layout is still the same as PS1_DEV_1 
     326// XXX this layout is still the same as PS1_DEV_1
    324327bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe)
    325328{
    … …  
    355358    // we write out all sources, regardless of quality.  the source flags tell us the state
    356359    for (int i = 0; i < sources->n; i++) {
    357         // skip source if it is not a ext sourc
    358         // XXX we have two places that extended source parameters are measured:
    359         // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models,
    360         // psphotFitEXT, which does the simple extended source model fit (not psf-convolved)
    361         // should we require both?
    362 
    363         pmSource *source = sources->data[i];
    364 
    365         // skip sources without measurements
    366         if (source->extpars == NULL) continue;
    367 
    368         // we require a PSF model fit (ignore the real crud)
    369         pmModel *model = source->modelPSF;
    370         if (model == NULL) continue;
    371 
    372         // XXX I need to split the extended models from the extended aperture measurements
    373         PAR = model->params->data.F32;
    374         dPAR = model->dparams->data.F32;
    375         xPos = PAR[PM_PAR_XPOS];
    376         yPos = PAR[PM_PAR_YPOS];
    377         xErr = dPAR[PM_PAR_XPOS];
    378         yErr = dPAR[PM_PAR_YPOS];
     360        // skip source if it is not a ext sourc
     361        // XXX we have two places that extended source parameters are measured:
     362        // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models,
     363        // psphotFitEXT, which does the simple extended source model fit (not psf-convolved)
     364        // should we require both?
     365
     366        pmSource *source = sources->data[i];
     367
     368        // skip sources without measurements
     369        if (source->extpars == NULL) continue;
     370
     371        // we require a PSF model fit (ignore the real crud)
     372        pmModel *model = source->modelPSF;
     373        if (model == NULL) continue;
     374
     375        // XXX I need to split the extended models from the extended aperture measurements
     376        PAR = model->params->data.F32;
     377        dPAR = model->dparams->data.F32;
     378        xPos = PAR[PM_PAR_XPOS];
     379        yPos = PAR[PM_PAR_YPOS];
     380        xErr = dPAR[PM_PAR_XPOS];
     381        yErr = dPAR[PM_PAR_YPOS];
    379382
    380383        row = psMetadataAlloc ();
    … …  
    387390        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    388391
    389         // Petrosian measurements
    390         // XXX insert header data: petrosian ref radius, flux ratio
    391         if (doPetrosian) {
    392             pmSourcePetrosianValues *petrosian = source->extpars->petrosian;
    393             if (petrosian) {
    394                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       petrosian->mag);
    395                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr);
    396                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          petrosian->rad);
    397                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    petrosian->radErr);
    398             } else {
    399                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
    400                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", NAN);
    401                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          NAN);
    402                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    NAN);
    403             }
    404         }
    405 
    406         // Kron measurements
    407         if (doKron) {
    408             pmSourceKronValues *kron = source->extpars->kron;
    409             if (kron) {
    410                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       kron->mag);
    411                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", kron->magErr);
    412                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          kron->rad);
    413                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    kron->radErr);
    414             } else {
    415                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       NAN);
    416                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", NAN);
    417                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          NAN);
    418                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    NAN);
    419             }
    420         }
    421 
    422         // Isophot measurements
    423         // XXX insert header data: isophotal level
    424         if (doIsophotal) {
    425             pmSourceIsophotalValues *isophot = source->extpars->isophot;
    426             if (isophot) {
    427                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       isophot->mag);
    428                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr);
    429                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          isophot->rad);
    430                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    isophot->radErr);
    431             } else {
    432                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       NAN);
    433                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", NAN);
    434                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          NAN);
    435                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    NAN);
    436             }
    437         }
    438 
    439         // Flux Annuli
    440         if (doAnnuli) {
    441             pmSourceAnnuli *annuli = source->extpars->annuli;
    442             if (annuli) {
    443                 psVector *fluxVal = annuli->flux;
    444                 psVector *fluxErr = annuli->fluxErr;
    445                 psVector *fluxVar = annuli->fluxVar;
    446 
    447                 for (int j = 0; j < fluxVal->n; j++) {
    448                     char name[32];
    449                     sprintf (name, "FLUX_VAL_R_%02d", j);
    450                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]);
    451                     sprintf (name, "FLUX_ERR_R_%02d", j);
    452                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]);
    453                     sprintf (name, "FLUX_VAR_R_%02d", j);
    454                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]);
    455                 }
    456             } else {
    457                 for (int j = 0; j < radialBinsLower->n; j++) {
    458                     char name[32];
    459                     sprintf (name, "FLUX_VAL_R_%02d", j);
    460                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN);
    461                     sprintf (name, "FLUX_ERR_R_%02d", j);
    462                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN);
    463                     sprintf (name, "FLUX_VAR_R_%02d", j);
    464                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN);
    465                 }
    466             }
    467         }
    468 
    469         psArrayAdd (table, 100, row);
    470         psFree (row);
     392        // Petrosian measurements
     393        // XXX insert header data: petrosian ref radius, flux ratio
     394        if (doPetrosian) {
     395            pmSourcePetrosianValues *petrosian = source->extpars->petrosian;
     396            if (petrosian) {
     397                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       petrosian->mag);
     398                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr);
     399                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          petrosian->rad);
     400                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    petrosian->radErr);
     401            } else {
     402                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
     403                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", NAN);
     404                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          NAN);
     405                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    NAN);
     406            }
     407        }
     408
     409        // Kron measurements
     410        if (doKron) {
     411            pmSourceKronValues *kron = source->extpars->kron;
     412            if (kron) {
     413                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       kron->mag);
     414                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", kron->magErr);
     415                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          kron->rad);
     416                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    kron->radErr);
     417            } else {
     418                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       NAN);
     419                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", NAN);
     420                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          NAN);
     421                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    NAN);
     422            }
     423        }
     424
     425        // Isophot measurements
     426        // XXX insert header data: isophotal level
     427        if (doIsophotal) {
     428            pmSourceIsophotalValues *isophot = source->extpars->isophot;
     429            if (isophot) {
     430                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       isophot->mag);
     431                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr);
     432                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          isophot->rad);
     433                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    isophot->radErr);
     434            } else {
     435                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       NAN);
     436                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", NAN);
     437                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          NAN);
     438                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    NAN);
     439            }
     440        }
     441
     442        // Flux Annuli
     443        if (doAnnuli) {
     444            pmSourceAnnuli *annuli = source->extpars->annuli;
     445            if (annuli) {
     446                psVector *fluxVal = annuli->flux;
     447                psVector *fluxErr = annuli->fluxErr;
     448                psVector *fluxVar = annuli->fluxVar;
     449
     450                for (int j = 0; j < fluxVal->n; j++) {
     451                    char name[32];
     452                    sprintf (name, "FLUX_VAL_R_%02d", j);
     453                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]);
     454                    sprintf (name, "FLUX_ERR_R_%02d", j);
     455                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]);
     456                    sprintf (name, "FLUX_VAR_R_%02d", j);
     457                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]);
     458                }
     459            } else {
     460                for (int j = 0; j < radialBinsLower->n; j++) {
     461                    char name[32];
     462                    sprintf (name, "FLUX_VAL_R_%02d", j);
     463                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN);
     464                    sprintf (name, "FLUX_ERR_R_%02d", j);
     465                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN);
     466                    sprintf (name, "FLUX_VAR_R_%02d", j);
     467                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN);
     468                }
     469            }
     470        }
     471
     472        psArrayAdd (table, 100, row);
     473        psFree (row);
    471474    }
    472475
    473476    if (table->n == 0) {
    474         psFitsWriteBlank (fits, outhead, extname);
    475         psFree (outhead);
    476         psFree (table);
    477         return true;
     477        psFitsWriteBlank (fits, outhead, extname);
     478        psFree (outhead);
     479        psFree (table);
     480        return true;
    478481    }
    479482
    480483    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
    481484    if (!psFitsWriteTable (fits, outhead, table, extname)) {
    482         psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
    483         psFree (outhead);
    484         psFree(table);
    485         return false;
     485        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
     486        psFree (outhead);
     487        psFree(table);
     488        return false;
    486489    }
    487490    psFree (outhead);
    … …  
    491494}
    492495
    493 // XXX this layout is still the same as PS1_DEV_1 
     496// XXX this layout is still the same as PS1_DEV_1
    494497bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, psArray *sources, char *extname)
    495498{
    … …  
    515518    int nParamMax = 0;
    516519    for (int i = 0; i < sources->n; i++) {
    517         pmSource *source = sources->data[i];
    518         if (source->modelFits == NULL) continue;
    519         for (int j = 0; j < source->modelFits->n; j++) {
    520             pmModel *model = source->modelFits->data[j];
    521             assert (model);
    522             nParamMax = PS_MAX (nParamMax, model->params->n);
    523         }
     520        pmSource *source = sources->data[i];
     521        if (source->modelFits == NULL) continue;
     522        for (int j = 0; j < source->modelFits->n; j++) {
     523            pmModel *model = source->modelFits->data[j];
     524            assert (model);
     525            nParamMax = PS_MAX (nParamMax, model->params->n);
     526        }
    524527    }
    525528
    … …  
    529532    for (int i = 0; i < sources->n; i++) {
    530533
    531         pmSource *source = sources->data[i];
    532 
    533         // XXX if no model fits are saved, write out modelEXT?
    534         if (source->modelFits == NULL) continue;
    535 
    536         // We have multiple sources : need to flag the one used to subtract the light (the 'best' model)
    537         for (int j = 0; j < source->modelFits->n; j++) {
    538 
    539             // choose the convolved EXT model, if available, otherwise the simple one
    540             pmModel *model = source->modelFits->data[j];
    541             assert (model);
    542 
    543             PAR = model->params->data.F32;
    544             dPAR = model->dparams->data.F32;
    545             xPos = PAR[PM_PAR_XPOS];
    546             yPos = PAR[PM_PAR_YPOS];
    547             xErr = dPAR[PM_PAR_XPOS];
    548             yErr = dPAR[PM_PAR_YPOS];
    549 
    550             axes = pmPSF_ModelToAxes (PAR, 20.0);
    551 
    552             row = psMetadataAlloc ();
    553 
    554             // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
    555             psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
    556             psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT",            0, "EXT model x coordinate",                     xPos);
    557             psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT",            0, "EXT model y coordinate",                     yPos);
    558             psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG",        0, "Sigma in EXT x coordinate",                  xErr);
    559             psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG",        0, "Sigma in EXT y coordinate",                  yErr);
    560             psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG",     0, "EXT fit instrumental magnitude",             model->mag);
    561             psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude",        model->magErr);
    562 
    563             psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS",          0, "number of model parameters",                 model->params->n);
    564             psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "name of model",                              pmModelClassGetName (model->type));
    565 
    566             // XXX these should be major and minor, not 'x' and 'y'
    567             psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width in x coordinate",                  axes.major);
    568             psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width in y coordinate",                  axes.minor);
    569             psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                      axes.theta);
    570 
    571             // write out the other generic parameters
    572             for (int k = 0; k < nParamMax; k++) {
    573                 if (k == PM_PAR_I0) continue;
    574                 if (k == PM_PAR_SKY) continue;
    575                 if (k == PM_PAR_XPOS) continue;
    576                 if (k == PM_PAR_YPOS) continue;
    577                 if (k == PM_PAR_SXX) continue;
    578                 if (k == PM_PAR_SXY) continue;
    579                 if (k == PM_PAR_SYY) continue;
    580 
    581                 snprintf (name, 64, "EXT_PAR_%02d", k);
    582 
    583                 if (k < model->params->n) {
    584                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]);
    585                 } else {
    586                     psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN);
    587                 }
    588             }
    589 
    590             // XXX other parameters which may be set.
    591             // XXX flag / value to define the model
    592             // XXX write out the model type, fit status flags
    593 
    594             psArrayAdd (table, 100, row);
    595             psFree (row);
    596         }
     534        pmSource *source = sources->data[i];
     535
     536        // XXX if no model fits are saved, write out modelEXT?
     537        if (source->modelFits == NULL) continue;
     538
     539        // We have multiple sources : need to flag the one used to subtract the light (the 'best' model)
     540        for (int j = 0; j < source->modelFits->n; j++) {
     541
     542            // choose the convolved EXT model, if available, otherwise the simple one
     543            pmModel *model = source->modelFits->data[j];
     544            assert (model);
     545
     546            PAR = model->params->data.F32;
     547            dPAR = model->dparams->data.F32;
     548            xPos = PAR[PM_PAR_XPOS];
     549            yPos = PAR[PM_PAR_YPOS];
     550            xErr = dPAR[PM_PAR_XPOS];
     551            yErr = dPAR[PM_PAR_YPOS];
     552
     553            axes = pmPSF_ModelToAxes (PAR, 20.0);
     554
     555            row = psMetadataAlloc ();
     556
     557            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
     558            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
     559            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT",            0, "EXT model x coordinate",                     xPos);
     560            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT",            0, "EXT model y coordinate",                     yPos);
     561            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG",        0, "Sigma in EXT x coordinate",                  xErr);
     562            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG",        0, "Sigma in EXT y coordinate",                  yErr);
     563            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG",     0, "EXT fit instrumental magnitude",             model->mag);
     564            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude",        model->magErr);
     565
     566            psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS",          0, "number of model parameters",                 model->params->n);
     567            psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "name of model",                              pmModelClassGetName (model->type));
     568
     569            // XXX these should be major and minor, not 'x' and 'y'
     570            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width in x coordinate",                  axes.major);
     571            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width in y coordinate",                  axes.minor);
     572            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                      axes.theta);
     573
     574            // write out the other generic parameters
     575            for (int k = 0; k < nParamMax; k++) {
     576                if (k == PM_PAR_I0) continue;
     577                if (k == PM_PAR_SKY) continue;
     578                if (k == PM_PAR_XPOS) continue;
     579                if (k == PM_PAR_YPOS) continue;
     580                if (k == PM_PAR_SXX) continue;
     581                if (k == PM_PAR_SXY) continue;
     582                if (k == PM_PAR_SYY) continue;
     583
     584                snprintf (name, 64, "EXT_PAR_%02d", k);
     585
     586                if (k < model->params->n) {
     587                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]);
     588                } else {
     589                    psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN);
     590                }
     591            }
     592
     593            // XXX other parameters which may be set.
     594            // XXX flag / value to define the model
     595            // XXX write out the model type, fit status flags
     596
     597            psArrayAdd (table, 100, row);
     598            psFree (row);
     599        }
    597600    }
    598601
    599602    if (table->n == 0) {
    600         psFitsWriteBlank (fits, outhead, extname);
    601         psFree (outhead);
    602         psFree (table);
    603         return true;
     603        psFitsWriteBlank (fits, outhead, extname);
     604        psFree (outhead);
     605        psFree (table);
     606        return true;
    604607    }
    605608
    606609    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
    607610    if (!psFitsWriteTable (fits, outhead, table, extname)) {
    608         psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
    609         psFree (outhead);
    610         psFree(table);
    611         return false;
     611        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
     612        psFree (outhead);
     613        psFree(table);
     614        return false;
    612615    }
    613616    psFree (outhead);
    … …  
    645648    psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
    646649    psPlaneTransformApply (&ptTP_y, fpa->toTPA, &ptFP);
    647    
     650
    648651    // the resulting Tangent Plane coordinates are in TP pixels; convert to local Tangent Plane
    649652    // degrees
    … …  
    671674    *posAngle = NAN;
    672675    *pltScale = NAN;
    673    
     676
    674677    return false;
    675678}
Note: See TracChangeset for help on using the changeset viewer.