Index: /branches/eam_branch_20071023/psModules/src/camera/pmFPAfileIO.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/camera/pmFPAfileIO.c (revision 15358) +++ /branches/eam_branch_20071023/psModules/src/camera/pmFPAfileIO.c (revision 15359) @@ -192,4 +192,7 @@ status = pmPSFmodelReadForView (view, file, config); break; + case PM_FPA_FILE_ASTROM: + status = pmAstromReadForView (view, file, config); + break; case PM_FPA_FILE_JPEG: case PM_FPA_FILE_KAPA: @@ -273,4 +276,5 @@ case PM_FPA_FILE_CMF: case PM_FPA_FILE_PSF: + case PM_FPA_FILE_ASTROM: case PM_FPA_FILE_JPEG: case PM_FPA_FILE_KAPA: @@ -336,4 +340,10 @@ // have their data stored on the readout->analysis metadata structure of another // (existing) fpa + if (file->type == PM_FPA_FILE_CMF) { + if (!pmFPAviewCheckDataStatusForSources (view, file)) { + psTrace("psModules.camera", 6, "skip write for %s, no data for this entry", file->name); + return true; + } + } if (file->type == PM_FPA_FILE_PSF) { if (!pmPSFmodelCheckDataStatusForView (view, file)) { @@ -342,6 +352,6 @@ } } - if (file->type == PM_FPA_FILE_CMF) { - if (!pmFPAviewCheckDataStatusForSources (view, file)) { + if (file->type == PM_FPA_FILE_ASTROM) { + if (!pmAstromCheckDataStatusForView (view, file)) { psTrace("psModules.camera", 6, "skip write for %s, no data for this entry", file->name); return true; @@ -449,4 +459,8 @@ case PM_FPA_FILE_PSF: status = pmPSFmodelWriteForView (view, file, config); + break; + + case PM_FPA_FILE_ASTROM: + status = pmAstromWriteForView (view, file, config); break; @@ -508,4 +522,5 @@ case PM_FPA_FILE_CMF: case PM_FPA_FILE_PSF: + case PM_FPA_FILE_ASTROM: psTrace ("psModules.camera", 5, "closing %s (%s) (%d:%d:%d)\n", file->filename, file->name, view->chip, view->cell, view->readout); status = psFitsClose (file->fits); @@ -574,4 +589,5 @@ case PM_FPA_FILE_CMF: case PM_FPA_FILE_PSF: + case PM_FPA_FILE_ASTROM: psTrace ("psModules.camera", 5, "NOT freeing %s (%s) : save for further analysis\n", file->filename, file->name); return true; @@ -712,4 +728,5 @@ case PM_FPA_FILE_CMF: case PM_FPA_FILE_PSF: + case PM_FPA_FILE_ASTROM: psTrace ("psModules.camera", 5, "opening %s (%s) (%d:%d:%d)\n", file->filename, file->name, view->chip, view->cell, view->readout); @@ -861,4 +878,7 @@ status = pmPSFmodelWritePHU (view, file, config); break; + case PM_FPA_FILE_ASTROM: + status = pmAstromWritePHU (view, file, config); + break; case PM_FPA_FILE_SX: case PM_FPA_FILE_RAW: Index: /branches/eam_branch_20071023/psModules/src/objects/pmModel.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/objects/pmModel.c (revision 15358) +++ /branches/eam_branch_20071023/psModules/src/objects/pmModel.c (revision 15359) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-21 00:04:07 $ + * @version $Revision: 1.15.6.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-10-23 20:54:53 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -217,5 +217,5 @@ Ro = psImageInterpolateOptionsAlloc( PS_INTERPOLATE_BILINEAR, - model->residuals->Ro, NULL, NULL, 0, 0.0, 0.0, 1, 0, 0.0); + model->residuals->Ro, NULL, model->mask, 0, 0.0, 0.0, 1, 0, 0.0); Rx = psImageInterpolateOptionsAlloc( PS_INTERPOLATE_BILINEAR, @@ -257,6 +257,6 @@ float oy = yBin*(imageRow + 0.5 - yCenter) + yResidCenter; + psU8 mflux = 0; if (mode & PM_MODEL_OP_RES0) { - psU8 mflux = 0; double Fo = 0.0; psImageInterpolate (&Fo, NULL, &mflux, ox, oy, Ro); @@ -265,6 +265,6 @@ } } - if (mode & PM_MODEL_OP_RES1) { - psU8 mflux = 0; + // skip Rx,Ry if Ro is masked + if (!mflux && (mode & PM_MODEL_OP_RES1)) { double Fx = 0.0; double Fy = 0.0; Index: /branches/eam_branch_20071023/psModules/src/objects/pmPSF_IO.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/objects/pmPSF_IO.c (revision 15358) +++ /branches/eam_branch_20071023/psModules/src/objects/pmPSF_IO.c (revision 15359) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.27 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-09 19:26:25 $ + * @version $Revision: 1.27.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-10-23 20:54:53 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -93,6 +93,5 @@ } -bool pmPSFmodelWriteForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) -{ +bool pmPSFmodelWriteForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) { pmFPA *fpa = file->fpa; Index: /branches/eam_branch_20071023/psModules/src/objects/pmSource.h =================================================================== --- /branches/eam_branch_20071023/psModules/src/objects/pmSource.h (revision 15358) +++ /branches/eam_branch_20071023/psModules/src/objects/pmSource.h (revision 15359) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-27 03:35:29 $ + * @version $Revision: 1.18.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-10-23 20:54:53 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -84,5 +84,46 @@ psRegion region; ///< area on image covered by selected pixels float sky, skyErr; ///< The sky and its error at the center of the object + pmSourceExtendedParameters *extpars; // extended source parameters }; + +typedef struct { + psVector *radius; + psVector *flux; +} pmSourceRadialProfile; + +typedef struct { + psVector *flux; + psVector *fluxVar; // measured variance + psVector *fluxErr; // formal error +} pmSourceAnnuli; + +typedef struct { + float mag; + float magErr; + float rad; + float radErr; +} pmSourceIsophotalValues; + +typedef struct { + float mag; + float magErr; + float rad; + float radErr; +} pmSourcePetrosianValues; + +typedef struct { + float mag; + float magErr; + float rad; + float radErr; +} pmSourceKronValues; + +typedef struct { + pmSourceRadialProfile *profile; + pmSourceAnnuli *annuli; + pmSourceIsophotalValues *isophot; + pmSourcePetrosianValues *petrosian; + pmSourceKronValues *kron; +} pmSourceExtendedParameters; /** pmPSFClump data structure Index: /branches/eam_branch_20071023/psModules/src/objects/pmSourceIO.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/objects/pmSourceIO.c (revision 15358) +++ /branches/eam_branch_20071023/psModules/src/objects/pmSourceIO.c (revision 15359) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.52 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-05 22:46:34 $ + * @version $Revision: 1.52.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-10-23 20:54:53 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -224,10 +224,12 @@ bool status; - char *exttype; - char *dataname; - char *headname; pmHDU *hdu; psMetadata *updates; psMetadata *outhead; + + char *exttype = NULL; + char *dataname = NULL; + char *xsrcname = NULL; + char *headname = NULL; // XXX if sources is NULL, skip the cell or write out empty tables? @@ -312,4 +314,14 @@ } dataname = pmFPAfileNameFromRule (rule, file, view); + + if (XSRC_OUTPUT) { + // EXTNAME for extended source data table + rule = psMetadataLookupStr(&status, menu, "CMF.XSRC"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XSRC in EXTNAME.RULES in camera.config"); + return false; + } + xsrcname = pmFPAfileNameFromRule (rule, file, view); + } } @@ -375,7 +387,6 @@ } if (!strcmp (exttype, "PS1_DEV_1")) { - status = pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname); - } - + status = pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname, xsrcname); + } if (!status) { psError(PS_ERR_IO, false, "writing CMF data to %s with format %s\n", file->filename, exttype); Index: /branches/eam_branch_20071023/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 15358) +++ /branches/eam_branch_20071023/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 15359) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-09 03:10:05 $ + * @version $Revision: 1.4.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-10-23 20:54:53 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -47,5 +47,5 @@ // XXX how do I generate the source tables which I need to send to PSPS? -bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname) +bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, char *xsrcname) { @@ -57,4 +57,9 @@ psF32 xPos, yPos; psF32 xErr, yErr; + + // if we request XSRC output, add the XSRC name to this header + if (xsrcname) { + psMetadataAddStr (tableHeader, PS_LIST_TAIL, "EXTXSRC", PS_META_REPLACE, "name of XSRC table extension", xsrcname); + } // let's write these out in S/N order @@ -144,6 +149,10 @@ return false; } - psFree (table); + + if (xsrcname) { + pmSourcesWriteXSRC_PS1_DEV_1 (file->fits, sources, xsrcname); + } + return true; } @@ -241,2 +250,85 @@ return sources; } + +bool pmSourcesWriteXSRC_PS1_DEV_1 (psFits *fits, psArray *sources, char *extname) +{ + + psArray *table; + psMetadata *row; + int i; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + psF32 xPos, yPos; + psF32 xErr, yErr; + + // create a header to hold the output data + outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + + table = psArrayAllocEmpty (sources->n); + + // we write out all sources, regardless of quality. the source flags tell us the state + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + source->seq = i; + + // skip source if it is not a ext sourc + + // no difference between PSF and non-PSF model + // XXX the PSF output should report the value for the psf, not the ext, model + pmModel *model = source->modelEXT; + if (model == NULL) continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + + axes = pmPSF_ModelToAxes (PAR, 20.0); + + row = psMetadataAlloc (); + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) + psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); + psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "PSF x coordinate", xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "PSF y coordinate", yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); + psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); + psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", PS_MIN (99.0, source->extMag)); + // XXX need to calculate psfMag, psfMagErr, extMag, extMagErr... + psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", PS_MIN (99.0, source->errMag)); + + // XXX these should be major and minor, not 'x' and 'y' + psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_X", PS_DATA_F32, "PSF width in x coordinate", axes.major); + psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_Y", PS_DATA_F32, "PSF width in y coordinate", axes.minor); + psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + + // other values that I need to report: + // R, Mag Petrosian, ..... + + psArrayAdd (table, 100, row); + psFree (row); + } + + if (table->n == 0) { + psFitsWriteBlank (fits, outhead, extname); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable (fits, outhead, table, extname)) { + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree(table); + return false; + } + psFree (table); + + return true; +} Index: /branches/eam_branch_20071023/psastro/src/psastroAnalysis.c =================================================================== --- /branches/eam_branch_20071023/psastro/src/psastroAnalysis.c (revision 15358) +++ /branches/eam_branch_20071023/psastro/src/psastroAnalysis.c (revision 15359) @@ -16,7 +16,4 @@ return true; } - - // XXX test retrun - // return true; // load the reference stars overlapping the data stars Index: /branches/eam_branch_20071023/psastro/src/psastroDefineFiles.c =================================================================== --- /branches/eam_branch_20071023/psastro/src/psastroDefineFiles.c (revision 15358) +++ /branches/eam_branch_20071023/psastro/src/psastroDefineFiles.c (revision 15359) @@ -4,7 +4,6 @@ // these calls bind the I/O handle to the specified fpa - pmFPAfile *output = NULL; - output = pmFPAfileDefineOutput (config, input->fpa, "PSASTRO.OUTPUT"); + pmFPAfile *output = pmFPAfileDefineOutput (config, input->fpa, "PSASTRO.OUTPUT"); if (!output) { psError(PSASTRO_ERR_CONFIG, false, "Failed to build FPA from PSASTRO.INPUT"); @@ -12,4 +11,11 @@ } output->save = true; + + // XXX not sure of this: should this be bind from args? + pmFPAfile *refAstrom = pmFPAfileDefineFromConf (&status, config, "PSASTRO.REF.ASTROM"); + if (!refAstrom) { + psError(PSASTRO_ERR_CONFIG, false, "Failed to load ref astrometry"); + return false; + } # if (0) Index: /branches/eam_branch_20071023/psastro/src/psastroMosaicAstrom.c =================================================================== --- /branches/eam_branch_20071023/psastro/src/psastroMosaicAstrom.c (revision 15358) +++ /branches/eam_branch_20071023/psastro/src/psastroMosaicAstrom.c (revision 15359) @@ -24,4 +24,9 @@ // compare chips with supplied mosaic model // adjust significant outliers to match model + + if (!psastroEnforceChips (fpa, recipe) { + psError(PSASTRO_ERR_UNKNOWN, false, "failed to align problematic chips"); + return false; + } // given the existing per-chip astrometry, determine matches between raw and ref stars Index: /branches/eam_branch_20071023/psphot/src/psphotAnnuli.c =================================================================== --- /branches/eam_branch_20071023/psphot/src/psphotAnnuli.c (revision 15358) +++ /branches/eam_branch_20071023/psphot/src/psphotAnnuli.c (revision 15359) @@ -3,6 +3,55 @@ bool psphotAnnuli (pmSource *source, psMetadata *recipe, psMaskType maskVal) { - psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n"); + assert (source->extpars); + assert (source->extpars->profile); + assert (source->extpars->profile->radius); + assert (source->extpars->profile->flux); + + psVector *radius = source->extpars->profile->radius; + psVector *weight = source->extpars->profile->weight; + psVector *flux = source->extpars->profile->flux; + + // XXX how do I define the radii? we can put a vector in the recipe... + // radialBins defines the bounds or start and stop (we can skip some that way... + psVector *radialBinsLower = psMetadataLookupVector (&status, recipe, RADIAL_ANNULAR_BINS_LOWER); + psVector *radialBinsUpper = psMetadataLookupVector (&status, recipe, RADIAL_ANNULAR_BINS_UPPER); + assert (radialBinsLower->n == radialBinsUpper->n); + + psVector *fluxValues = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32); + psVector *fluxSquare = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32); + psVector *fluxWeight = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32); + psVector *pixelCount = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32); + psVectorInit (fluxValues, 0.0); + psVectorInit (fluxSquare, 0.0); + psVectorInit (fluxWeight, 0.0); + psVectorInit (pixelCount, 0.0); + + // XXX this code assumes the radii are in pixels. convert from arcsec with plate scale + // XXX assume the annulii above are not overlapping? much faster... + // XXX this might be must faster in the reverse order: loop over annulii and use disection to + // skip to the start of the annulus. + for (int i = 0; i < flux->n; i++) { + for (int j = 0; j < radialBinsLower->n; j++) { + if (radius->data.F32[i] < radialBinsLower->data.F32[j]) continue; + if (radius->data.F32[i] > radialBinsUpper->data.F32[j]) continue; + fluxValues->data.F32[j] += flux->data.F32[i]; + fluxSquare->data.F32[j] += PS_SQR(flux->data.F32[i]); + fluxWeight->data.F32[j] += weight->data.F32[i]; + pixelCount->data.F32[j] += 1.0; + } + } + + for (int j = 0; j < radialBinsLower->n; j++) { + fluxValues->data.F32[j] /= pixelCount->data.F32[j]; + fluxSquare->data.F32[j] /= pixelCount->data.F32[j]; + fluxSquare->data.F32[j] -= PS_SQR(fluxValues->data.F32[j]); + } + + source->extpars->annuli = pmSourceAnnuliAlloc (); + source->extpars->annuli->flux = fluxValues; + source->extpars->annuli->fluxErr = fluxWeight; + source->extpars->annuli->fluxVar = fluxSquare; + return true; +} -} Index: /branches/eam_branch_20071023/psphot/src/psphotExtendedSources.c =================================================================== --- /branches/eam_branch_20071023/psphot/src/psphotExtendedSources.c (revision 15358) +++ /branches/eam_branch_20071023/psphot/src/psphotExtendedSources.c (revision 15359) @@ -66,7 +66,18 @@ return false; } else { + // XXX why am I caching the model? pmSourceCacheModel (source, maskVal); // XXX put this in the source model function? psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y); Npsf ++; + } + + // all of these below require the radial profile + // XXX push this as a test and call in each of the functions below? + // XXX this constructs a pmSourceExtendedParameters element + if (doPetrosian || doIsophotal || doAnnuli || doKron) { + if (!psphotRadialProfile (source, recipe, maskVal)) { + psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate radial profile"); + return false; + } } Index: /branches/eam_branch_20071023/psphot/src/psphotIsophotal.c =================================================================== --- /branches/eam_branch_20071023/psphot/src/psphotIsophotal.c (revision 15358) +++ /branches/eam_branch_20071023/psphot/src/psphotIsophotal.c (revision 15359) @@ -3,6 +3,82 @@ bool psphotIsophotal (pmSource *source, psMetadata *recipe, psMaskType maskVal) { - psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n"); + assert (source->extpars); + assert (source->extpars->profile); + assert (source->extpars->profile->radius); + assert (source->extpars->profile->flux); + + psVector *radius = source->extpars->profile->radius; + psVector *flux = source->extpars->profile->flux; + + // flux at which to measure isophotal parameters + // XXX cache this? + float ISOPHOT_FLUX = psMetadataLookupF32 (&status, recipe, "ISOPHOTAL_FLUX"); + + // find the first bin below the flux level and the last above the level + // XXX can this be done faster with disection? + // XXX do I need to worry about crazy outliers? + // XXX should i be smoothing or fitting the curve? + firstBelow = -1; + lastAbove = -1; + for (int i = 0; i < flux->n; i++) { + if (flux->data.F32[i] > ISOPHOT_FLUX) lastAbove = i; + if ((firstBelow < 0) && (flux->data.F32[i] < ISOPHOT_FLUX)) firstBelow = i; + } + + // need to examine pixels in this vicinity + float isophotalFluxFirst = 0; + float isophotalFluxLast = 0; + for (int i = 0; i <= PS_MAX(firstBelow, lastAbove); i++) { + if (i <= firstBelow) { + isophotFluxFirst += flux->data.F32[i]; + } + if (i <= lastAbove) { + isophotFluxLast += flux->data.F32[i]; + } + } + float isophotalFlux = 0.5*(isophotFluxLast + isophotFluxFirst); + float isophotalFluxErr = 0.5*fabs(isophotFluxLast - isophotFluxFirst); + + float isophotalRad = 0.5*(radius->data.F32[firstBelow] + radius->data.F32[lastAbove]); + float isophotalRadErr = 0.5*fabs(radius->data.F32[firstBelow] - radius->data.F32[lastAbove]); + + if (!source->extpars->isophot) { + source->extpars->isophot = pmSourceIsophotalValuesAlloc (); + } + + // these are uncalibrated: instrumental mags and pixel units + source->extpars->isophot->mag = -2.5*log10(isophotalFlux); + source->extpars->isophot->magErr = isophotalFluxErr / isophotalFlux; + + source->extpars->isophot->rad = isophotalRad; + source->extpars->isophot->radErr = isophotalRadErr; + return true; } + +# if (0) + // XXX cache the tmpScalar + psScalar fluxScalar; + fluxScalar.type.type = PS_TYPE_F32; + fluxScalar.data.F32 = ISOPHOT_FLUX; + + // radius and flux must be sorted by radius + + // XXX in general, flux decreases monotonically with radius. however, since exceptions are + // possible we need to do something to smooth or otherwise handle the flux variations + + // get the index of the flux-sorted vector + psVector *fluxIndex = psVectorSortIndex (flux); + + // XXX need to write psVectorIndexBinaryDisect, which operates on a + + // use disection to get in the right vicinity + binS = psVectorIndexBinaryDisect (&result, flux, fluxIndex, fluxScalar, PS_BISECT_FIRST); + binE = psVectorIndexBinaryDisect (&result, flux, fluxIndex, fluxScalar, PS_BISECT_LAST); + + // find the min and max radius bins in this range + + // XXX do something to choose a more accurate radial bin + +# endif Index: /branches/eam_branch_20071023/psphot/src/psphotMakeResiduals.c =================================================================== --- /branches/eam_branch_20071023/psphot/src/psphotMakeResiduals.c (revision 15358) +++ /branches/eam_branch_20071023/psphot/src/psphotMakeResiduals.c (revision 15359) @@ -194,4 +194,10 @@ resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0; //resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev; + + // if value < NRESID_SIGMA*sigma, mask pixel in resid map + if (resid->Ro->data.F32[oy][ox] < NRESID_SIGMA*fluxStats->sampleStdev) { + resid->mask->data.U8[oy][ox] = RESID_SIG_MASK; + } + } else { assert (SPATIAL_ORDER == 1); @@ -227,4 +233,9 @@ resid->Rx->data.F32[oy][ox] = B->data.F64[1]; resid->Ry->data.F32[oy][ox] = B->data.F64[2]; + + dRo = sqrt(A->data.F32[0][0]); + if (resid->Ro->data.F32[oy][ox] < NRESID_SIGMA*dRo) { + resid->mask->data.U8[oy][ox] = RESID_SIG_MASK; + } //resid->weight->data.F32[oy][ox] = XXX; } Index: /branches/eam_branch_20071023/psphot/src/psphotPetrosian.c =================================================================== --- /branches/eam_branch_20071023/psphot/src/psphotPetrosian.c (revision 15358) +++ /branches/eam_branch_20071023/psphot/src/psphotPetrosian.c (revision 15359) @@ -3,5 +3,77 @@ bool psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal) { - psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n"); + assert (source->extpars); + assert (source->extpars->profile); + assert (source->extpars->profile->radius); + assert (source->extpars->profile->flux); + + psVector *radius = source->extpars->profile->radius; + psVector *flux = source->extpars->profile->flux; + + // flux at which to measure isophotal parameters + // XXX cache this? + float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0"); + float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO"); + + // first find flux at R0 + firstBelow = -1; + lastAbove = -1; + for (int i = 0; i < radius->n; i++) { + if (radius->data.F32[i] > PETROSIAN_R0) lastAbove = i; + if ((firstBelow < 0) && (flux->data.F32[i] < PETROSIAN_R0)) firstBelow = i; + } + + // average flux in this range + fluxR0 = 0.0; + fluxRn = 0; + for (int i = PS_MIN(firstBelow, lastAbove); i <= PS_MAX(firstBelow, lastAbove); i++) { + fluxR0 += flux->data.F32[i]; + fluxRn ++; + } + fluxR0 = fluxR0 / (float)(fluxRn); + + // target flux for petrosian radius + fluxRP = fluxR0 * PETROSIAN_RF; + + // find the first bin below the flux level and the last above the level + // XXX can this be done faster with disection? + // XXX do I need to worry about crazy outliers? + // XXX should i be smoothing or fitting the curve? + firstBelow = -1; + lastAbove = -1; + for (int i = 0; i < flux->n; i++) { + if (flux->data.F32[i] > fluxRP) lastAbove = i; + if ((firstBelow < 0) && (flux->data.F32[i] < fluxRP)) firstBelow = i; + } + + // need to examine pixels in this vicinity + float fluxFirst = 0; + float fluxLast = 0; + for (int i = 0; i <= PS_MAX(firstBelow, lastAbove); i++) { + if (i <= firstBelow) { + fluxFirst += flux->data.F32[i]; + } + if (i <= lastAbove) { + fluxLast += flux->data.F32[i]; + } + } + float flux = 0.5*(fluxLast + fluxFirst); + float fluxErr = 0.5*fabs(fluxLast - fluxFirst); + // XXX need to use the weight appropriately here... + + float rad = 0.5*(radius->data.F32[firstBelow] + radius->data.F32[lastAbove]); + float radErr = 0.5*fabs(radius->data.F32[firstBelow] - radius->data.F32[lastAbove]); + + if (!source->extpars->petrosian) { + source->extpars->petrosian = pmSourcePetrosianValuesAlloc (); + } + + // these are uncalibrated: instrumental mags and pixel units + source->extpars->petrosian->mag = -2.5*log10(flux); + source->extpars->petrosian->magErr = fluxErr / flux; + + source->extpars->petrosian->rad = rad; + source->extpars->petrosian->radErr = radErr; + return true;