IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 17327


Ignore:
Timestamp:
Apr 4, 2008, 3:48:48 PM (18 years ago)
Author:
eugene
Message:

fix up the output tables of XSRC / XFIT

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branch_20080324/psModules/src/objects/pmSourceIO_PS1_DEV_1.c

    r17309 r17327  
    33 *  @author EAM, IfA
    44 *
    5  *  @version $Revision: 1.9.2.1 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2008-04-03 06:01:38 $
     5 *  @version $Revision: 1.9.2.2 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2008-04-05 01:48:48 $
    77 *
    88 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    251251}
    252252
    253 bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname)
     253bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe)
    254254{
    255255
     256    bool status;
    256257    psArray *table;
    257258    psMetadata *row;
    … …  
    270271
    271272    table = psArrayAllocEmpty (sources->n);
     273
     274    // which extended source analyses should we perform?
     275    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     276    bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
     277    bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
     278    bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
     279
     280    psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
     281    psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
     282    assert (radialBinsLower->n == radialBinsUpper->n);
    272283
    273284    // we write out all sources, regardless of quality.  the source flags tell us the state
    … …  
    306317
    307318        // Petrosian measurements
    308         pmSourcePetrosianValues *petrosian = source->extpars->petrosian;
    309         if (petrosian) {
    310             psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       petrosian->mag);
    311             psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr);
    312             psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          petrosian->rad);
    313             psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    petrosian->radErr);
    314         }
     319        // XXX insert header data: petrosian ref radius, flux ratio
     320        if (doPetrosian) {
     321            pmSourcePetrosianValues *petrosian = source->extpars->petrosian;
     322            if (petrosian) {
     323                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       petrosian->mag);
     324                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr);
     325                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          petrosian->rad);
     326                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    petrosian->radErr);
     327            } else {
     328                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
     329                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", NAN);
     330                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          NAN);
     331                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    NAN);
     332            }
     333        }
    315334
    316335        // Kron measurements
    317         pmSourceKronValues *kron = source->extpars->kron;
    318         if (kron) {
    319             psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       kron->mag);
    320             psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", kron->magErr);
    321             psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          kron->rad);
    322             psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    kron->radErr);
     336        if (doKron) {
     337            pmSourceKronValues *kron = source->extpars->kron;
     338            if (kron) {
     339                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       kron->mag);
     340                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", kron->magErr);
     341                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          kron->rad);
     342                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    kron->radErr);
     343            } else {
     344                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       NAN);
     345                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", NAN);
     346                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          NAN);
     347                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    NAN);
     348            }
    323349        }
    324350
    325351        // Isophot measurements
    326352        // XXX insert header data: isophotal level
    327         pmSourceIsophotalValues *isophot = source->extpars->isophot;
    328         if (isophot) {
    329             psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       isophot->mag);
    330             psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr);
    331             psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          isophot->rad);
    332             psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    isophot->radErr);
    333         }
    334 
    335         pmSourceAnnuli *annuli = source->extpars->annuli;
    336         if (annuli) {
    337             psVector *fluxVal = annuli->flux;
    338             psVector *fluxErr = annuli->fluxErr;
    339             psVector *fluxVar = annuli->fluxVar;
    340 
    341             for (int j = 0; j < fluxVal->n; j++) {
    342                 char name[32];
    343                 sprintf (name, "FLUX_VAL_R_%02d", j);
    344                 psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]);
    345                 sprintf (name, "FLUX_ERR_R_%02d", j);
    346                 psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]);
    347                 sprintf (name, "FLUX_VAR_R_%02d", j);
    348                 psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]);
     353        if (doIsophotal) {
     354            pmSourceIsophotalValues *isophot = source->extpars->isophot;
     355            if (isophot) {
     356                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       isophot->mag);
     357                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr);
     358                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          isophot->rad);
     359                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    isophot->radErr);
     360            } else {
     361                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       NAN);
     362                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", NAN);
     363                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          NAN);
     364                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    NAN);
    349365            }
    350366        }
    351367
    352         psArrayAdd (table, 100, row);
    353         psFree (row);
     368        // Flux Annuli
     369        if (doAnnuli) {
     370            pmSourceAnnuli *annuli = source->extpars->annuli;
     371            if (annuli) {
     372                psVector *fluxVal = annuli->flux;
     373                psVector *fluxErr = annuli->fluxErr;
     374                psVector *fluxVar = annuli->fluxVar;
     375
     376                for (int j = 0; j < fluxVal->n; j++) {
     377                    char name[32];
     378                    sprintf (name, "FLUX_VAL_R_%02d", j);
     379                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]);
     380                    sprintf (name, "FLUX_ERR_R_%02d", j);
     381                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]);
     382                    sprintf (name, "FLUX_VAR_R_%02d", j);
     383                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]);
     384                }
     385            } else {
     386                for (int j = 0; j < radialBinsLower->n; j++) {
     387                    char name[32];
     388                    sprintf (name, "FLUX_VAL_R_%02d", j);
     389                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN);
     390                    sprintf (name, "FLUX_ERR_R_%02d", j);
     391                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN);
     392                    sprintf (name, "FLUX_VAR_R_%02d", j);
     393                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN);
     394                }
     395            }
     396        }
     397
     398        psArrayAdd (table, 100, row);
     399        psFree (row);
    354400    }
    355401
    356402    if (table->n == 0) {
    357         psFitsWriteBlank (fits, outhead, extname);
    358         psFree (table);
    359         return true;
     403        psFitsWriteBlank (fits, outhead, extname);
     404        psFree (outhead);
     405        psFree (table);
     406        return true;
    360407    }
    361408
    362409    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
    363410    if (!psFitsWriteTable (fits, outhead, table, extname)) {
    364         psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
    365         psFree(table);
    366         return false;
    367     }
     411        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
     412        psFree (outhead);
     413        psFree(table);
     414        return false;
     415    }
     416    psFree (outhead);
    368417    psFree (table);
    369418
    … …  
    391440    sources = psArraySort (sources, pmSourceSortBySN);
    392441
     442    // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams
     443    int nParamMax = 0;
     444    for (int i = 0; i < sources->n; i++) {
     445        pmSource *source = sources->data[i];
     446        if (source->modelFits == NULL) continue;
     447        for (int j = 0; j < source->modelFits->n; j++) {
     448            pmModel *model = source->modelFits->data[j];
     449            assert (model);
     450            nParamMax = PS_MAX (nParamMax, model->params->n);
     451        }
     452    }
     453
    393454    table = psArrayAllocEmpty (sources->n);
    394455
    … …  
    420481
    421482            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
    422             psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    423             psMetadataAdd (row, PS_LIST_TAIL, "X_EXT",            PS_DATA_F32, "EXT model x coordinate",                     xPos);
    424             psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT",            PS_DATA_F32, "EXT model y coordinate",                     yPos);
    425             psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG",        PS_DATA_F32, "Sigma in EXT x coordinate",                  xErr);
    426             psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    427             psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG",     PS_DATA_F32, "EXT fit instrumental magnitude",             model->mag);
    428             psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        model->magErr);
     483            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
     484            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT",            0, "EXT model x coordinate",                     xPos);
     485            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT",            0, "EXT model y coordinate",                     yPos);
     486            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG",        0, "Sigma in EXT x coordinate",                  xErr);
     487            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG",        0, "Sigma in EXT y coordinate",                  yErr);
     488            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG",     0, "EXT fit instrumental magnitude",             model->mag);
     489            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude",        model->magErr);
     490
     491            psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS",          0, "number of model parameters",                 model->params->n);
     492            psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "name of model",                              pmModelClassGetName (model->type));
    429493
    430494            // XXX these should be major and minor, not 'x' and 'y'
    431             psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    PS_DATA_F32, "EXT width in x coordinate",                  axes.major);
    432             psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    PS_DATA_F32, "EXT width in y coordinate",                  axes.minor);
    433             psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA",        PS_DATA_F32, "EXT orientation angle",                      axes.theta);
     495            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width in x coordinate",                  axes.major);
     496            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width in y coordinate",                  axes.minor);
     497            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                      axes.theta);
    434498
    435499            // write out the other generic parameters
    436             for (int k = 0; k < model->params->n; k++) {
     500            for (int k = 0; k < nParamMax; k++) {
    437501                if (k == PM_PAR_I0) continue;
     502                if (k == PM_PAR_SKY) continue;
    438503                if (k == PM_PAR_XPOS) continue;
    439504                if (k == PM_PAR_YPOS) continue;
    … …  
    443508
    444509                snprintf (name, 64, "EXT_PAR_%02d", k);
    445                 psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]);
     510
     511                if (k < model->params->n) {
     512                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]);
     513                } else {
     514                    psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN);
     515                }
    446516            }
    447517
    … …  
    456526
    457527    if (table->n == 0) {
    458         psFitsWriteBlank (fits, outhead, extname);
    459         psFree (table);
    460         return true;
     528        psFitsWriteBlank (fits, outhead, extname);
     529        psFree (outhead);
     530        psFree (table);
     531        return true;
    461532    }
    462533
    463534    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
    464535    if (!psFitsWriteTable (fits, outhead, table, extname)) {
    465         psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
    466         psFree(table);
    467         return false;
    468     }
     536        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
     537        psFree (outhead);
     538        psFree(table);
     539        return false;
     540    }
     541    psFree (outhead);
    469542    psFree (table);
    470 
    471543    return true;
    472544}
Note: See TracChangeset for help on using the changeset viewer.