IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 18, 2010, 12:49:05 PM (16 years ago)
Author:
eugene
Message:

merging changes from trunk into branches/pap

Location:
branches/pap
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/pap

  • branches/pap/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c

    r27177 r28003  
    4545// followed by a zero-size matrix, followed by the table data
    4646
    47 bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources,
    48                                 psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
     47bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
    4948{
    5049    PS_ASSERT_PTR_NON_NULL(fits, false);
    … …  
    5453    psArray *table;
    5554    psMetadata *row;
    56     int i;
    5755    psF32 *PAR, *dPAR;
    5856    psEllipseAxes axes;
    … …  
    8381        pmSource *source = (pmSource *) sources->data[0];
    8482        if (source->seq == -1) {
    85           // let's write these out in S/N order
    86           sources = psArraySort (sources, pmSourceSortBySN);
     83            // let's write these out in S/N order
     84            sources = psArraySort (sources, pmSourceSortBySN);
    8785        } else {
    88           sources = psArraySort (sources, pmSourceSortBySeq);
     86            sources = psArraySort (sources, pmSourceSortBySeq);
    8987        }
    9088    }
    … …  
    9492    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    9593    // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
    96     for (i = 0; i < sources->n; i++) {
     94    for (int i = 0; i < sources->n; i++) {
    9795        pmSource *source = (pmSource *) sources->data[i];
    9896
    … …  
    102100        // generated on Alloc, and would thus be wrong for read in sources.
    103101        if (source->seq == -1) {
    104           source->seq = i;
     102            source->seq = i;
    105103        }
    106104
    … …  
    114112            yPos = PAR[PM_PAR_YPOS];
    115113            if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
    116               xErr = dPAR[PM_PAR_XPOS];
    117               yErr = dPAR[PM_PAR_YPOS];
     114                xErr = dPAR[PM_PAR_XPOS];
     115                yErr = dPAR[PM_PAR_YPOS];
    118116            } else {
    119               // in linear-fit mode, there is no error on the centroid
    120               xErr = source->peak->dx;
    121               yErr = source->peak->dy;
     117                // in linear-fit mode, there is no error on the centroid
     118                xErr = source->peak->dx;
     119                yErr = source->peak->dy;
    122120            }
    123121            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
    … …  
    212210    }
    213211
     212    // XXX why do we make a copy here to be supplemented with the masks?  why not do this in the calling function?
    214213    psMetadata *header = psMetadataCopy(NULL, tableHeader);
    215214    pmSourceMasksHeader(header);
    … …  
    344343}
    345344
    346 // XXX this layout is still the same as PS1_DEV_1
    347 bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe)
     345bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
    348346{
    349347
    … …  
    354352    psF32 xPos, yPos;
    355353    psF32 xErr, yErr;
     354    int nRow = -1;
    356355
    357356    // create a header to hold the output data
    … …  
    361360    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname);
    362361
     362    pmChip *chip = readout->parent->parent;
     363    pmFPA  *fpa  = chip->parent;
     364
     365    // zero point corrections
     366    bool status1 = false;
     367    bool status2 = false;
     368    float magOffset = 0.0;
     369    float exptime   = psMetadataLookupF32(&status1, fpa->concepts, "FPA.EXPOSURE");
     370    float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
     371    if (!isfinite(zeropt)) {
     372        zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
     373    }
     374    if (status1 && status2 && (exptime > 0.0)) {
     375        magOffset = zeropt + 2.5*log10(exptime);
     376    }
     377
    363378    // let's write these out in S/N order
    364379    sources = psArraySort (sources, pmSourceSortBySN);
    … …  
    367382
    368383    // which extended source analyses should we perform?
    369     // bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     384    bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
     385    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    370386    // bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    371     // bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    372387    // bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
    373388
    374     psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    375     psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
    376     assert (radialBinsLower->n == radialBinsUpper->n);
     389    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
     390    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
     391    psAssert (radMin->n == radMax->n, "inconsistent annular bins");
     392
     393    // int nRadialBins = 0;
     394    // if (doAnnuli) {
     395    //  // get the max count of radial bins
     396    //  for (int i = 0; i < sources->n; i++) {
     397    //      pmSource *source = sources->data[i];
     398    //      if (!source->extpars) continue;
     399    //      if (!source->extpars->radProfile ) continue;
     400    //         if (!source->extpars->radProfile->binSB) continue;
     401    //      nRadialBins = PS_MAX(nRadialBins, source->extpars->radProfile->binSB->n);
     402    //  }
     403    // }
    377404
    378405    // we write out all sources, regardless of quality.  the source flags tell us the state
    … …  
    410437        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    411438
    412 # if (0)
     439        float AxialRatio = NAN;
     440        float AxialTheta = NAN;
     441        pmSourceExtendedPars *extpars = source->extpars;
     442        if (extpars) {
     443            AxialRatio = extpars->axes.minor / extpars->axes.major;
     444            AxialTheta = extpars->axes.theta;
     445        }
     446        psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO",       PS_DATA_F32, "Axial Ratio of radial profile",              AxialRatio);
     447        psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA",        PS_DATA_F32, "Angle of radial profile ellipse",                  AxialTheta);
     448
    413449        // Petrosian measurements
    414450        // XXX insert header data: petrosian ref radius, flux ratio
     451        // XXX check flags to see if Pet was measured
    415452        if (doPetrosian) {
    416             pmSourcePetrosianValues *petrosian = source->extpars->petrosian;
    417             if (petrosian) {
    418                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       petrosian->mag);
    419                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr);
    420                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          petrosian->rad);
    421                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    petrosian->radErr);
     453            pmSourceExtendedPars *extpars = source->extpars;
     454            if (extpars) {
     455                // XXX note that this mag is either calibrated or instrumental depending on existence of zero point
     456                float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point
     457                float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point
     458                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude", mag);
     459                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", magErr);
     460                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius);
     461                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr);
     462                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50",     PS_DATA_F32, "Petrosian R50 (pix)", extpars->petrosianR50);
     463                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)", extpars->petrosianR50Err);
     464                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90);
     465                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err);
    422466            } else {
    423467                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
    … …  
    425469                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          NAN);
    426470                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    NAN);
     471                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50",     PS_DATA_F32, "Petrosian R50 (pix)", NAN);
     472                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)",NAN);
     473                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", NAN);
     474                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN);
    427475            }
    428476        }
    429477
     478# if (0)
    430479        // Kron measurements
    431480        if (doKron) {
    … …  
    460509            }
    461510        }
    462 
    463         // Flux Annuli
     511# endif
     512
     513        // Flux Annuli (if we have extended source measurements, we have these.  only optionally save them)
    464514        if (doAnnuli) {
    465             pmSourceAnnuli *annuli = source->extpars->annuli;
    466             if (annuli) {
    467                 psVector *fluxVal = annuli->flux;
    468                 psVector *fluxErr = annuli->fluxErr;
    469                 psVector *fluxVar = annuli->fluxVar;
    470 
    471                 for (int j = 0; j < fluxVal->n; j++) {
    472                     char name[32];
    473                     sprintf (name, "FLUX_VAL_R_%02d", j);
    474                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]);
    475                     sprintf (name, "FLUX_ERR_R_%02d", j);
    476                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]);
    477                     sprintf (name, "FLUX_VAR_R_%02d", j);
    478                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]);
    479                 }
    480             } else {
    481                 for (int j = 0; j < radialBinsLower->n; j++) {
    482                     char name[32];
    483                     sprintf (name, "FLUX_VAL_R_%02d", j);
    484                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN);
    485                     sprintf (name, "FLUX_ERR_R_%02d", j);
    486                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN);
    487                     sprintf (name, "FLUX_VAR_R_%02d", j);
    488                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN);
    489                 }
    490             }
    491         }
    492 
    493 # endif
    494         psArrayAdd (table, 100, row);
    495         psFree (row);
    496     }
    497 
     515            psVector *radSB   = psVectorAlloc(radMin->n, PS_TYPE_F32);
     516            psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32);
     517            psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32);
     518            psVectorInit (radSB, NAN);
     519            psVectorInit (radFlux, NAN);
     520            psVectorInit (radFill, NAN);
     521            if (!source->extpars) goto empty_annuli;
     522            if (!source->extpars->radProfile) goto empty_annuli;
     523            if (!source->extpars->radProfile->binSB) goto empty_annuli;
     524            psAssert (source->extpars->radProfile->binSum, "programming error");
     525            psAssert (source->extpars->radProfile->binFill, "programming error");
     526            psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths");
     527            psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths");
     528            psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths");
     529
     530            // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
     531            for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) {
     532                radSB->data.F32[j]   = source->extpars->radProfile->binSB->data.F32[j];
     533                radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j];
     534                radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j];
     535            }
     536
     537        empty_annuli:
     538            psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB);
     539            psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux);
     540            psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill);
     541            psFree (radSB);
     542            psFree (radFlux);
     543            psFree (radFill);
     544        }
     545        if (nRow < 0) {
     546            nRow = row->list->n;
     547        } else {
     548            psAssert (nRow == row->list->n, "inconsistent row lengths");
     549        }
     550        psArrayAdd (table, 100, row);
     551        psFree (row);
     552    }
     553   
    498554    if (table->n == 0) {
    499         if (!psFitsWriteBlank (fits, outhead, extname)) {
    500             psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
    501             psFree(outhead);
    502             psFree(table);
    503             return false;
    504         }
    505         psFree (outhead);
    506         psFree (table);
    507         return true;
    508     }
    509 
     555        if (!psFitsWriteBlank (fits, outhead, extname)) {
     556            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
     557            psFree(outhead);
     558            psFree(table);
     559            return false;
     560        }
     561        psFree (outhead);
     562        psFree (table);
     563        return true;
     564    }
     565   
    510566    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
    511567    if (!psFitsWriteTable (fits, outhead, table, extname)) {
    512         psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
    513         psFree (outhead);
    514         psFree(table);
    515         return false;
     568        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
     569        psFree (outhead);
     570    psFree(table);
     571    return false;
    516572    }
    517573    psFree (outhead);
    518574    psFree (table);
    519 
     575   
    520576    return true;
    521577}
    522578
    523579// XXX this layout is still the same as PS1_DEV_1
    524 bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname)
     580bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname)
    525581{
    526582
Note: See TracChangeset for help on using the changeset viewer.