IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 30704


Ignore:
Timestamp:
Feb 19, 2011, 10:32:56 AM (15 years ago)
Author:
eugene
Message:

allow radial apertures for non-psf-matched images; optionally save the covar matrix for extended pars; save the petrosian fill factor

File:
1 edited

Legend:

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

    r30621 r30704  
    556556        }
    557557        psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO",       PS_DATA_F32, "Axial Ratio of radial profile",              AxialRatio);
    558         psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA",        PS_DATA_F32, "Angle of radial profile ellipse",                  AxialTheta);
     558        psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA",        PS_DATA_F32, "Angle of radial profile ellipse",            AxialTheta);
    559559
    560560        // Petrosian measurements
     
    567567                float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point
    568568                float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point
    569                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude", mag);
    570                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", magErr);
    571                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius);
    572                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr);
     569                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",           PS_DATA_F32, "Petrosian Magnitude", mag);
     570                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",       PS_DATA_F32, "Petrosian Magnitude Error", magErr);
     571                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",        PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius);
     572                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR",    PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr);
    573573                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50",     PS_DATA_F32, "Petrosian R50 (pix)", extpars->petrosianR50);
    574574                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)", extpars->petrosianR50Err);
    575575                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90);
    576576                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err);
     577                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL",          PS_DATA_F32, "Petrosian Fill Factor", extpars->petrosianFill);
    577578            } else {
    578                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
    579                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", NAN);
    580                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          NAN);
    581                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    NAN);
     579                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",           PS_DATA_F32, "Petrosian Magnitude",       NAN);
     580                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",       PS_DATA_F32, "Petrosian Magnitude Error", NAN);
     581                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",        PS_DATA_F32, "Petrosian Radius",          NAN);
     582                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR",    PS_DATA_F32, "Petrosian Radius Error",    NAN);
    582583                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50",     PS_DATA_F32, "Petrosian R50 (pix)", NAN);
    583584                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)",NAN);
    584585                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", NAN);
    585586                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN);
     587                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL",          PS_DATA_F32, "Petrosian Fill Factor", NAN);
    586588            }
    587589        }
     
    747749
    748750                if (k < model->params->n) {
    749                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]);
     751                    psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]);
    750752                } else {
    751                     psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN);
     753                    psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", NAN);
    752754                }
    753755            }
    754756
    755             // XXX other parameters which may be set.
    756             // XXX flag / value to define the model
    757             // XXX write out the model type, fit status flags
    758 
     757            // optionally, write out the covariance matrix values
     758            // XXX do I need to pad this to match the biggest covar matrix?
     759            if (model->covar) {
     760                for (int iy = 0; iy < model->covar->numCols; iy++) {
     761                    for (int ix = iy; ix < model->covar->numCols; ix++) {
     762                        snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix);
     763                        psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]);
     764
     765                    }
     766                }                   
     767            }
    759768            psArrayAdd (table, 100, row);
    760769            psFree (row);
     
    790799bool pmSourcesWrite_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
    791800{
    792 
     801    bool status = false;
    793802    psArray *table;
    794803    psMetadata *row;
     
    796805    char keyword1[80], keyword2[80];
    797806
     807    // perform full non-linear fits / extended source analysis?
     808    if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) {
     809        psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n");
     810        return true;
     811    }
     812
    798813    // create a header to hold the output data
    799814    psMetadata *outhead = psMetadataAlloc ();
     
    803818
    804819    // we use this just to define the output vectors (which must be present for all objects)
    805     bool status = false;
    806820    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    807821    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
     
    818832    }
    819833
     834    // the FWHM values are available if we measured a psf-matched convolved set
    820835    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
    821     if (!fwhmValues) {
    822         psError (PM_ERR_CONFIG, true, "convolved or measured FWHM is not defined for this readout");
    823         return false;
    824     }
    825836
    826837    // let's write these out in S/N order
     
    836847        // skip sources without radial aper measurements (or insufficient)
    837848        if (source->radialAper == NULL) continue;
    838         psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set");
    839 
    840         for (int entry = 0; entry < fwhmValues->n; entry++) {
     849        // psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set");
     850
     851        for (int entry = 0; entry < source->radialAper->n; entry++) {
    841852
    842853            // choose the convolved EXT model, if available, otherwise the simple one
     
    863874            psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER",           0, "Center of aperture measurements",            xPos);
    864875            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER",           0, "Center of aperture measurements",            yPos);
    865             psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "FWHM of matched PSF",                        fwhmValues->data.F32[entry]);
     876            if (fwhmValues) {
     877                psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "FWHM of matched PSF",                    fwhmValues->data.F32[entry]);
     878            } else {
     879                psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "image is not FWHM-matched",              NAN);
     880            }
    866881
    867882            // XXX if we have raw radial apertures, write them out here
    868             psVector *radFlux    = psVectorAlloc(radMax->n, PS_TYPE_F32);
    869             psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32);
    870             psVector *radFill    = psVectorAlloc(radMax->n, PS_TYPE_F32);
     883            psVector *radFlux      = psVectorAlloc(radMax->n, PS_TYPE_F32);
     884            psVector *radFluxErr   = psVectorAlloc(radMax->n, PS_TYPE_F32);
     885            psVector *radFill      = psVectorAlloc(radMax->n, PS_TYPE_F32);
     886            psVector *radFluxStdev = psVectorAlloc(radMax->n, PS_TYPE_F32);
    871887            psVectorInit (radFlux,    NAN);
    872888            psVectorInit (radFluxErr, NAN);
     
    879895            // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
    880896            for (int j = 0; j < radialAper->flux->n; j++) {
    881                 radFlux->data.F32[j]    = radialAper->flux->data.F32[j];
    882                 radFluxErr->data.F32[j] = radialAper->fluxErr->data.F32[j];
    883                 radFill->data.F32[j]    = radialAper->fill->data.F32[j];
     897                radFlux->data.F32[j]      = radialAper->flux->data.F32[j];
     898                radFluxErr->data.F32[j]   = radialAper->fluxErr->data.F32[j];
     899                radFluxStdev->data.F32[j] = radialAper->fluxStdev->data.F32[j];
     900                radFill->data.F32[j]      = radialAper->fill->data.F32[j];
    884901            }
    885902
    886903        write_annuli:
    887             psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX",     PS_DATA_VECTOR, "flux within annuli",    radFlux);
    888             psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli",  radFluxErr);
    889             psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL",     PS_DATA_VECTOR, "fill factor of annuli", radFill);
     904            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX",       PS_DATA_VECTOR, "flux within annuli",       radFlux);
     905            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR",   PS_DATA_VECTOR, "flux error in annuli",     radFluxErr);
     906            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_DATA_VECTOR, "flux standard deviation",  radFluxStdev);
     907            psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL",       PS_DATA_VECTOR, "fill factor of annuli",    radFill);
    890908            psFree (radFlux);
    891909            psFree (radFluxErr);
     910            psFree (radFluxStdev);
    892911            psFree (radFill);
    893912
Note: See TracChangeset for help on using the changeset viewer.