IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 17309


Ignore:
Timestamp:
Apr 2, 2008, 8:01:38 PM (18 years ago)
Author:
eugene
Message:

finish split of ext params and fits, write out all fitted ext models

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r16819 r17309  
    33 *  @author EAM, IfA
    44 *
    5  *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2008-03-05 01:08:08 $
     5 *  @version $Revision: 1.9.2.1 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2008-04-03 06:01:38 $
    77 *
    88 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    256256    psArray *table;
    257257    psMetadata *row;
    258     int i;
    259258    psF32 *PAR, *dPAR;
    260     psEllipseAxes axes;
    261259    psF32 xPos, yPos;
    262260    psF32 xErr, yErr;
    … …  
    274272
    275273    // we write out all sources, regardless of quality.  the source flags tell us the state
    276     for (i = 0; i < sources->n; i++) {
     274    for (int i = 0; i < sources->n; i++) {
    277275        // skip source if it is not a ext sourc
    278276        // XXX we have two places that extended source parameters are measured:
    … …  
    283281        pmSource *source = sources->data[i];
    284282
    285         // choose the convolved EXT model, if available, otherwise the simple one
    286         pmModel *model = source->modelConv;
    287         if (model == NULL) {
    288             model = source->modelEXT;
    289         }
    290         if (model == NULL) continue;
    291 
    292         // XXX do we need to do this?
     283        // skip sources without measurements
    293284        if (source->extpars == NULL) continue;
     285
     286        // we require a PSF model fit (ignore the real crud)
     287        pmModel *model = source->modelPSF;
     288        if (model == NULL) continue;
    294289
    295290        // XXX I need to split the extended models from the extended aperture measurements
    … …  
    301296        yErr = dPAR[PM_PAR_YPOS];
    302297
    303         // XXX for the aperture values, I probably can / should just refer to the psf position and not write any of this
    304         axes = pmPSF_ModelToAxes (PAR, 20.0);
    305 
    306298        row = psMetadataAlloc ();
     299
    307300        // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
    308301        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    … …  
    311304        psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG",        PS_DATA_F32, "Sigma in EXT x coordinate",                  xErr);
    312305        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    313         psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG",     PS_DATA_F32, "EXT fit instrumental magnitude",             source->extMag);
    314         psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
    315 
    316         // XXX these should be major and minor, not 'x' and 'y'
    317         psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    PS_DATA_F32, "EXT width in x coordinate",                  axes.major);
    318         psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    PS_DATA_F32, "EXT width in y coordinate",                  axes.minor);
    319         psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA",        PS_DATA_F32, "EXT orientation angle",                      axes.theta);
    320 
    321         // XXX at the moment, this will be a fixed sized table, but perhaps have multiple output formats depending on what is measured?
    322 
    323         // other values that I need to report:
    324         // R, Mag Petrosian, .....
    325         if (source->extpars) {
    326 
    327             // Petrosian measurements
    328             pmSourcePetrosianValues *petrosian = source->extpars->petrosian;
    329             if (petrosian) {
    330                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       petrosian->mag);
    331                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr);
    332                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          petrosian->rad);
    333                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    petrosian->radErr);
    334             }
    335 
    336             // Kron measurements
    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             }
    344 
    345             // Isophot measurements
    346             // XXX insert header data: isophotal level
    347             pmSourceIsophotalValues *isophot = source->extpars->isophot;
    348             if (isophot) {
    349                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       isophot->mag);
    350                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr);
    351                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          isophot->rad);
    352                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    isophot->radErr);
    353             }
    354 
    355             pmSourceAnnuli *annuli = source->extpars->annuli;
    356             if (annuli) {
    357                 psVector *fluxVal = annuli->flux;
    358                 psVector *fluxErr = annuli->fluxErr;
    359                 psVector *fluxVar = annuli->fluxVar;
    360 
    361                 for (int j = 0; j < fluxVal->n; j++) {
    362                     char name[32];
    363                     sprintf (name, "FLUX_VAL_R_%02d", j);
    364                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]);
    365                     sprintf (name, "FLUX_ERR_R_%02d", j);
    366                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]);
    367                     sprintf (name, "FLUX_VAR_R_%02d", j);
    368                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]);
    369                 }
     306
     307        // 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        }
     315
     316        // 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);
     323        }
     324
     325        // Isophot measurements
     326        // 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]);
    370349            }
    371350        }
    … …  
    397376    psArray *table;
    398377    psMetadata *row;
    399     int i;
    400378    psF32 *PAR, *dPAR;
    401379    psEllipseAxes axes;
    402380    psF32 xPos, yPos;
    403381    psF32 xErr, yErr;
     382    char name[64];
    404383
    405384    // create a header to hold the output data
    … …  
    415394
    416395    // we write out all sources, regardless of quality.  the source flags tell us the state
    417     for (i = 0; i < sources->n; i++) {
    418         // skip source if it is not a ext sourc
    419         // XXX we have two places that extended source parameters are measured:
    420         // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models,
    421         // psphotFitEXT, which does the simple extended source model fit (not psf-convolved)
    422         // should we require both?
     396    for (int i = 0; i < sources->n; i++) {
    423397
    424398        pmSource *source = sources->data[i];
    425399
    426         // XXX need to have an array of model fits; loop over the array and write out all
    427         // which we calculated
    428 
    429         // choose the convolved EXT model, if available, otherwise the simple one
    430         // XXX should not need to choose: write both out
    431         pmModel *model = source->modelConv;
    432         if (model == NULL) {
    433             model = source->modelEXT;
     400        // XXX if no model fits are saved, write out modelEXT?
     401        if (source->modelFits == NULL) continue;
     402
     403        // We have multiple sources : need to flag the one used to subtract the light (the 'best' model)
     404        for (int j = 0; j < source->modelFits->n; j++) {
     405
     406            // choose the convolved EXT model, if available, otherwise the simple one
     407            pmModel *model = source->modelFits->data[j];
     408            assert (model);
     409
     410            PAR = model->params->data.F32;
     411            dPAR = model->dparams->data.F32;
     412            xPos = PAR[PM_PAR_XPOS];
     413            yPos = PAR[PM_PAR_YPOS];
     414            xErr = dPAR[PM_PAR_XPOS];
     415            yErr = dPAR[PM_PAR_YPOS];
     416
     417            axes = pmPSF_ModelToAxes (PAR, 20.0);
     418
     419            row = psMetadataAlloc ();
     420
     421            // 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);
     429
     430            // 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);
     434
     435            // write out the other generic parameters
     436            for (int k = 0; k < model->params->n; k++) {
     437                if (k == PM_PAR_I0) continue;
     438                if (k == PM_PAR_XPOS) continue;
     439                if (k == PM_PAR_YPOS) continue;
     440                if (k == PM_PAR_SXX) continue;
     441                if (k == PM_PAR_SXY) continue;
     442                if (k == PM_PAR_SYY) continue;
     443
     444                snprintf (name, 64, "EXT_PAR_%02d", k);
     445                psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]);
     446            }
     447
     448            // XXX other parameters which may be set.
     449            // XXX flag / value to define the model
     450            // XXX write out the model type, fit status flags
     451
     452            psArrayAdd (table, 100, row);
     453            psFree (row);
    434454        }
    435         if (model == NULL) continue;
    436 
    437         PAR = model->params->data.F32;
    438         dPAR = model->dparams->data.F32;
    439         xPos = PAR[PM_PAR_XPOS];
    440         yPos = PAR[PM_PAR_YPOS];
    441         xErr = dPAR[PM_PAR_XPOS];
    442         yErr = dPAR[PM_PAR_YPOS];
    443 
    444         axes = pmPSF_ModelToAxes (PAR, 20.0);
    445 
    446         row = psMetadataAlloc ();
    447         // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
    448         psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    449         psMetadataAdd (row, PS_LIST_TAIL, "X_EXT",            PS_DATA_F32, "EXT model x coordinate",                     xPos);
    450         psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT",            PS_DATA_F32, "EXT model y coordinate",                     yPos);
    451         psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG",        PS_DATA_F32, "Sigma in EXT x coordinate",                  xErr);
    452         psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    453         psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG",     PS_DATA_F32, "EXT fit instrumental magnitude",             source->extMag);
    454         psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
    455 
    456         // XXX these should be major and minor, not 'x' and 'y'
    457         psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    PS_DATA_F32, "EXT width in x coordinate",                  axes.major);
    458         psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    PS_DATA_F32, "EXT width in y coordinate",                  axes.minor);
    459         psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA",        PS_DATA_F32, "EXT orientation angle",                      axes.theta);
    460 
    461         // XXX other parameters which may be set.
    462         // XXX flag / value to define the model
    463 
    464         psArrayAdd (table, 100, row);
    465         psFree (row);
    466455    }
    467456
Note: See TracChangeset for help on using the changeset viewer.