Index: branches/pap/psModules/src/astrom/pmAstrometryWCS.c =================================================================== --- branches/pap/psModules/src/astrom/pmAstrometryWCS.c (revision 27708) +++ branches/pap/psModules/src/astrom/pmAstrometryWCS.c (revision 28003) @@ -972,4 +972,8 @@ // FPA -> TPA : identidy + // NOTE: streaksremove's linearizeTransform function passes pointers to skeleton outFPA and outChip structs + // Only outFPA->toSky is valid. No other memebers in these structs should be read. The resulting output transforms + // are copied to inFPA and inChip by the caller. + int nSamples = 10; // 10 samples in each dimension Index: branches/pap/psModules/src/camera/pmFPAConstruct.c =================================================================== --- branches/pap/psModules/src/camera/pmFPAConstruct.c (revision 27708) +++ branches/pap/psModules/src/camera/pmFPAConstruct.c (revision 28003) @@ -100,27 +100,4 @@ return num; -} - -// Get the name of a PHU chip or cell from the header -static psString phuNameFromHeader(const char *name, // The name to lookup: "CELL.NAME" or "CHIP.NAME" - const psMetadata *fileInfo, // FILE within the camera format description - const psMetadata *header // Primary header - ) -{ - assert(name && strlen(name) > 0); - assert(fileInfo); - assert(header); - - bool mdok = true; // Result of MD lookup - psString keyword = psMetadataLookupStr(&mdok, fileInfo, name); - if (!mdok || strlen(keyword) == 0) { - return false; - } - psMetadataItem *resultItem = psMetadataLookup(header, keyword); - if (!resultItem) { - psError(PS_ERR_IO, true, "Unable to find %s in primary header to identify %s.\n", keyword, name); - return NULL; - } - return psMetadataItemParseString(resultItem); } @@ -1102,5 +1079,4 @@ // It returns a view corresponding to the PHU static pmFPAview *addSource(pmFPA *fpa, // The FPA - const char *fpaObs, // The desired FPA observation name const pmFPAview *phuView, // The view corresponding to the PHU, or NULL const psMetadata *header, // The PHU header, or NULL @@ -1114,19 +1090,4 @@ bool mdok; // Status of MD lookup - - // If FPA.OBS is already defined, new name must match it; otherwise, warn the user that something - // potentially bad is happening. - if (fpaObs && install) { - const char *currentName = psMetadataLookupStr(&mdok, fpa->concepts, "FPA.OBS"); // Current name - if (mdok && currentName && strlen(currentName) > 0 && strcmp(currentName, fpaObs) != 0) { - psWarning("FPA.OBS for new source (%s) doesn't match FPA.OBS for current fpa (%s).", - fpaObs, currentName); - } - psMetadataAddStr(fpa->concepts, PS_LIST_HEAD, "FPA.OBS", PS_META_REPLACE, "Observation identifier", - fpaObs); - } else if (!psMetadataLookup(fpa->concepts, "FPA.OBS")) { - // Make sure there is an FPA.OBS - psMetadataAddStr(fpa->concepts, PS_LIST_HEAD, "FPA.OBS", 0, "Observation identifier", "UNKNOWN"); - } psMetadata *fileInfo = psMetadataLookupMetadata(&mdok, format, "FILE"); // The file information @@ -1350,5 +1311,5 @@ } -bool pmFPAAddSourceFromFormat(pmFPA *fpa, const char *fpaObs, const psMetadata *format) +bool pmFPAAddSourceFromFormat(pmFPA *fpa, const psMetadata *format) { PS_ASSERT_PTR_NON_NULL(fpa, false); @@ -1359,5 +1320,5 @@ pmFPAview *view = pmFPAviewAlloc(0);// View for current level if (phuLevel == PM_FPA_LEVEL_FPA) { - if (!pmFPAAddSourceFromView(fpa, fpaObs, view, format)) { + if (!pmFPAAddSourceFromView(fpa, view, format)) { psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); psFree(view); @@ -1368,5 +1329,5 @@ while ((chip = pmFPAviewNextChip(view, fpa, 1))) { if (phuLevel == PM_FPA_LEVEL_CHIP) { - if (!pmFPAAddSourceFromView(fpa, fpaObs, view, format)) { + if (!pmFPAAddSourceFromView(fpa, view, format)) { psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); psFree(view); @@ -1377,5 +1338,5 @@ while ((cell = pmFPAviewNextCell(view, fpa, 1))) { if (phuLevel == PM_FPA_LEVEL_CELL) { - if (!pmFPAAddSourceFromView(fpa, fpaObs, view, format)) { + if (!pmFPAAddSourceFromView(fpa, view, format)) { psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); psFree(view); @@ -1392,5 +1353,5 @@ } -bool pmFPAAddSourceFromView(pmFPA *fpa, const char *fpaObs, const pmFPAview *phuView, +bool pmFPAAddSourceFromView(pmFPA *fpa, const pmFPAview *phuView, const psMetadata *format) { @@ -1399,5 +1360,5 @@ PS_ASSERT_PTR_NON_NULL(format, false); - pmFPAview *view = addSource(fpa, fpaObs, phuView, NULL, format, true); + pmFPAview *view = addSource(fpa, phuView, NULL, format, true); bool status = (view != NULL); psFree(view); @@ -1418,12 +1379,5 @@ } - // Check the name of the FPA - psString fpaObs = phuNameFromHeader("FPA.OBS", fileInfo, phu); // New observation name for the FPA - if (!fpaObs || strlen(fpaObs) == 0) { - psWarning("Unable to determine FPA.OBS: check for FPA.OBS in FILE in camera format"); - } - - pmFPAview *view = addSource(fpa, fpaObs, NULL, phu, format, true); // View of PHU, to return - psFree(fpaObs); + pmFPAview *view = addSource(fpa, NULL, phu, format, true); // View of PHU, to return return view; @@ -1437,5 +1391,5 @@ PS_ASSERT_PTR_NON_NULL(format, NULL); - return addSource(fpa, NULL, NULL, phu, format, false); + return addSource(fpa, NULL, phu, format, false); } Index: branches/pap/psModules/src/camera/pmFPAConstruct.h =================================================================== --- branches/pap/psModules/src/camera/pmFPAConstruct.h (revision 27708) +++ branches/pap/psModules/src/camera/pmFPAConstruct.h (revision 28003) @@ -29,5 +29,4 @@ /// This is suitable for generating an output FPA given the desired format. bool pmFPAAddSourceFromFormat(pmFPA *fpa, ///< The FPA - const char *fpaObs, ///< FPA.NAME for the source const psMetadata *format ///< Format of file ); @@ -38,5 +37,4 @@ /// configuration is required in order to describe how the FPA is laid out in terms of disk files. bool pmFPAAddSourceFromView(pmFPA *fpa, ///< The FPA - const char *fpaObs, ///< FPA.NAME for the source const pmFPAview *phuView, ///< The view, corresponding to the PHU const psMetadata *format ///< Format of file Index: branches/pap/psModules/src/camera/pmFPAWrite.c =================================================================== --- branches/pap/psModules/src/camera/pmFPAWrite.c (revision 27708) +++ branches/pap/psModules/src/camera/pmFPAWrite.c (revision 28003) @@ -165,21 +165,21 @@ // generate the HDU, but only copies the structure. if (!blank && !hdu->blankPHU && !*imageArray) { - if (!pmHDUGenerateForCell(cell)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate HDU for cell --- likely programming error."); - return false; - } - if (!*imageArray) { - if (type == FPA_WRITE_TYPE_IMAGE) { - psError(PS_ERR_UNKNOWN, false, "Expected to write an image, but it is missing...programming error?."); - return false; - } - if (type == FPA_WRITE_TYPE_MASK) { - psWarning("No mask image for this cell; skipping"); - } - if (type == FPA_WRITE_TYPE_VARIANCE) { - psWarning("No variance image for this cell; skipping"); - } - return true; - } + if (!pmHDUGenerateForCell(cell)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate HDU for cell --- likely programming error."); + return false; + } + if (!*imageArray) { + if (type == FPA_WRITE_TYPE_IMAGE) { + psError(PS_ERR_UNKNOWN, false, "Expected to write an image, but it is missing...programming error?."); + return false; + } + if (type == FPA_WRITE_TYPE_MASK) { + psWarning("No mask image for this cell; skipping"); + } + if (type == FPA_WRITE_TYPE_VARIANCE) { + psWarning("No variance image for this cell; skipping"); + } + return true; + } } @@ -307,5 +307,5 @@ if (writeBlank || writeImage) { - if (!pmConceptsWriteFPA(fpa, true, config)) { + if (!pmConceptsWriteFPA(fpa, true, config)) { psError(PS_ERR_IO, false, "Unable to write concepts for FPA.\n"); return false; @@ -341,5 +341,5 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// Update the FPA.OBS, CHIP.NAME and CELL.NAME in the FITS header, if required +// Update the CHIP.NAME and CELL.NAME in the FITS header, if required bool pmFPAUpdateNames(pmFPA *fpa, pmChip *chip, pmCell *cell, psS64 imageId, psS64 sourceId) { @@ -375,12 +375,4 @@ psError(PS_ERR_UNKNOWN, true, "Unable to find FILE information in camera format.\n"); return false; - } - if (fpa) { - const char *fpaObsHdr = psMetadataLookupStr(&mdok, fileData, "FPA.OBS"); - if (mdok && fpaObsHdr && strlen(fpaObsHdr) > 0) { - const char *fpaObs = psMetadataLookupStr(NULL, fpa->concepts, "FPA.OBS"); - psMetadataAddStr(hduHigh->header, PS_LIST_TAIL, fpaObsHdr, PS_META_REPLACE, - "Observation identifier", fpaObs); - } } Index: branches/pap/psModules/src/camera/pmFPAfileFitsIO.c =================================================================== --- branches/pap/psModules/src/camera/pmFPAfileFitsIO.c (revision 27708) +++ branches/pap/psModules/src/camera/pmFPAfileFitsIO.c (revision 28003) @@ -75,13 +75,6 @@ } - pmFPA *nameSource = file->src; // Source of FPA.OBS - if (!nameSource) { - nameSource = file->fpa; - } - bool mdok; // Status of MD lookup - const char *fpaObs = psMetadataLookupStr(&mdok, nameSource->concepts, "FPA.OBS"); // Observation id - pmFPA *copy = pmFPAConstruct(file->camera, file->cameraName); // FPA to return - if (!pmFPAAddSourceFromView(copy, fpaObs, phuView, file->format)) { + if (!pmFPAAddSourceFromView(copy, phuView, file->format)) { psError(PS_ERR_UNKNOWN, false, "Unable to insert HDU into FPA for writing.\n"); psFree(copy); Index: branches/pap/psModules/src/camera/pmFPAfileIO.c =================================================================== --- branches/pap/psModules/src/camera/pmFPAfileIO.c (revision 27708) +++ branches/pap/psModules/src/camera/pmFPAfileIO.c (revision 28003) @@ -286,12 +286,5 @@ } - pmFPA *nameSource = file->src; // Source of FPA.OBS - if (!nameSource) { - nameSource = file->fpa; - } - bool mdok; // Status of MD lookup - const char *fpaObs = psMetadataLookupStr(&mdok, nameSource->concepts, "FPA.OBS"); // Obs. id - - pmFPAAddSourceFromView(file->fpa, fpaObs, view, format); + pmFPAAddSourceFromView(file->fpa, view, format); psTrace ("psModules.camera", 5, "created fpa data elements for %s (%s) (%d:%d:%d)\n", file->name, file->name, view->chip, view->cell, view->readout); Index: branches/pap/psModules/src/detrend/pmShutterCorrection.c =================================================================== --- branches/pap/psModules/src/detrend/pmShutterCorrection.c (revision 27708) +++ branches/pap/psModules/src/detrend/pmShutterCorrection.c (revision 28003) @@ -1015,11 +1015,15 @@ pmShutterCorrection *corr = pmShutterCorrectionFullFit(newtimes, newcounts, newerrors, guess); // The actual correction - psTrace("psModules.detrend", 5, "Shutter correction fit: scale: %f, offset: %f, offref: %f\n", corr->scale, corr->offset, corr->offref); - - if (corr && isfinite(corr->offref) && corr->valid) { - psTrace("psModules.detrend", 5, "Sample reference value: %f\n", corr->offref); - meanRef += corr->offref; - numGood++; - } + + if (corr) { + psTrace("psModules.detrend", 5, "Shutter correction fit: scale: %f, offset: %f, offref: %f\n", corr->scale, corr->offset, corr->offref); + if (isfinite(corr->offref) && corr->valid) { + psTrace("psModules.detrend", 5, "Sample reference value: %f\n", corr->offref); + meanRef += corr->offref; + numGood++; + } + } else { + psTrace("psModules.detrend", 5, "failed Shutter correction fit\n"); + } psFree(corr); Index: branches/pap/psModules/src/extras/Makefile.am =================================================================== --- branches/pap/psModules/src/extras/Makefile.am (revision 27708) +++ branches/pap/psModules/src/extras/Makefile.am (revision 28003) @@ -17,4 +17,5 @@ pmKapaPlots.h \ pmVisual.h \ + ippDiffMode.h \ ippStages.h Index: branches/pap/psModules/src/extras/ippDiffMode.h =================================================================== --- branches/pap/psModules/src/extras/ippDiffMode.h (revision 28003) +++ branches/pap/psModules/src/extras/ippDiffMode.h (revision 28003) @@ -0,0 +1,20 @@ +/* @file ippDiffMode.h + * @brief some macro defintions for the stages of the pipeline + * @author Bill Sweeney, IfA + * + * Copyright 2010 Institute for Astronomy, University of Hawaii + */ + +#ifndef IPP_DIFF_MODE_H +#define IPP_DIFF_MODE_H + +typedef enum { + IPP_DIFF_MODE_UNDEFINED = 0, + IPP_DIFF_MODE_WARP_WARP = 1, + IPP_DIFF_MODE_WARP_STACK = 2, + IPP_DIFF_MODE_STACK_WARP = 3, + IPP_DIFF_MODE_STACK_STACK = 4 +} ippDiffMode; + + +#endif Index: branches/pap/psModules/src/extras/ippStages.h =================================================================== --- branches/pap/psModules/src/extras/ippStages.h (revision 27708) +++ branches/pap/psModules/src/extras/ippStages.h (revision 28003) @@ -1,3 +1,3 @@ -/* @file psVisual.h +/* @file ippStages.h * @brief some macro defintions for the stages of the pipeline * @author Bill Sweeney, IfA Index: branches/pap/psModules/src/objects/pmPSFtryModel.c =================================================================== --- branches/pap/psModules/src/objects/pmPSFtryModel.c (revision 27708) +++ branches/pap/psModules/src/objects/pmPSFtryModel.c (revision 28003) @@ -107,6 +107,6 @@ // as we loop over orders, we need to refer to the initial selection, but we modify the // option values to match the current guess: save the max values here: - int Nx = options->psfTrendNx; - int Ny = options->psfTrendNy; + int Nx = (options->psfTrendMode == PM_TREND_MAP) ? options->psfTrendNx : options->psfTrendNx + 1; + int Ny = (options->psfTrendMode == PM_TREND_MAP) ? options->psfTrendNy : options->psfTrendNy + 1; for (int i = orderMin; i <= orderMax; i++) { Index: branches/pap/psModules/src/objects/pmSourceExtendedPars.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceExtendedPars.c (revision 27708) +++ branches/pap/psModules/src/objects/pmSourceExtendedPars.c (revision 28003) @@ -36,23 +36,67 @@ #include "pmSourceExtendedPars.h" -// *** pmSourceRadialProfile describes the radial profile of a source in elliptical contours, and -// intermediate data used to measure the profile +// pmSourceRadialFlux carries the raw radial flux information, including angular bins +static void pmSourceRadialFluxFree(pmSourceRadialFlux *flux) +{ + if (!flux) return; + psFree(flux->radii); + psFree(flux->fluxes); + psFree(flux->theta); + psFree(flux->isophotalRadii); +} + +pmSourceRadialFlux *pmSourceRadialFluxAlloc() +{ + pmSourceRadialFlux *flux = (pmSourceRadialFlux *)psAlloc(sizeof(pmSourceRadialFlux)); + psMemSetDeallocator(flux, (psFreeFunc) pmSourceRadialFluxFree); + + flux->radii = NULL; + flux->fluxes = NULL; + flux->theta = NULL; + flux->isophotalRadii = NULL; + + return flux; +} + +bool psMemCheckSourceRadialFlux(psPtr ptr) +{ + PS_ASSERT_PTR(ptr, false); + return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceRadialFluxFree); +} + +// pmSourceEllipticalFlux carries the elliptical renormalized radial flux info +static void pmSourceEllipticalFluxFree(pmSourceEllipticalFlux *flux) +{ + if (!flux) return; + psFree(flux->radiusElliptical); + psFree(flux->fluxElliptical); +} + +pmSourceEllipticalFlux *pmSourceEllipticalFluxAlloc() +{ + pmSourceEllipticalFlux *flux = (pmSourceEllipticalFlux *)psAlloc(sizeof(pmSourceEllipticalFlux)); + psMemSetDeallocator(flux, (psFreeFunc) pmSourceEllipticalFluxFree); + + flux->radiusElliptical = NULL; + flux->fluxElliptical = NULL; + + return flux; +} + +bool psMemCheckSourceEllipticalFlux(psPtr ptr) +{ + PS_ASSERT_PTR(ptr, false); + return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceEllipticalFluxFree); +} + +// pmSourceRadialProfile defines flux information in radial bins static void pmSourceRadialProfileFree(pmSourceRadialProfile *profile) { - if (!profile) { - return; - } - psFree(profile->radii); - psFree(profile->fluxes); - psFree(profile->theta); - psFree(profile->isophotalRadii); - - psFree(profile->radiusElliptical); - psFree(profile->fluxElliptical); - + if (!profile) return; psFree(profile->binSB); psFree(profile->binSBstdev); psFree(profile->binSBerror); - + psFree(profile->binSum); + psFree(profile->binFill); psFree(profile->radialBins); psFree(profile->area); @@ -63,20 +107,12 @@ pmSourceRadialProfile *profile = (pmSourceRadialProfile *)psAlloc(sizeof(pmSourceRadialProfile)); psMemSetDeallocator(profile, (psFreeFunc) pmSourceRadialProfileFree); - - profile->radii = NULL; - profile->fluxes = NULL; - profile->theta = NULL; - profile->isophotalRadii = NULL; - - profile->radiusElliptical = NULL; - profile->fluxElliptical = NULL; profile->binSB = NULL; profile->binSBstdev = NULL; profile->binSBerror = NULL; - + profile->binSum = NULL; + profile->binFill = NULL; profile->radialBins = NULL; profile->area = NULL; - return profile; } @@ -88,6 +124,6 @@ } - -// *** pmSourceRadialProfileFreeVectors frees the intermediate data values +# if (0) +// pmSourceRadialProfileFreeVectors frees the intermediate data values bool pmSourceRadialProfileFreeVectors(pmSourceRadialProfile *profile) { @@ -124,4 +160,5 @@ return true; } +# endif // *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors @@ -150,7 +187,8 @@ if (!pars) return; - psFree(pars->profile); - psFree(pars->petrosian_50); - psFree(pars->petrosian_80); + psFree(pars->radFlux); + psFree(pars->ellipticalFlux); + psFree(pars->radProfile); + psFree(pars->petProfile); return; } @@ -160,8 +198,17 @@ psMemSetDeallocator(pars, (psFreeFunc) pmSourceExtendedParsFree); - pars->profile = NULL; - pars->petrosian_50 = NULL; - pars->petrosian_80 = NULL; - + pars->radFlux = NULL; + pars->ellipticalFlux = NULL; + pars->radProfile = NULL; + pars->petProfile = NULL; + + pars->petrosianFlux = NAN; + pars->petrosianFluxErr = NAN; + pars->petrosianRadius = NAN; + pars->petrosianRadiusErr = NAN; + pars->petrosianR90 = NAN; + pars->petrosianR90Err = NAN; + pars->petrosianR50 = NAN; + pars->petrosianR50Err = NAN; return pars; } Index: branches/pap/psModules/src/objects/pmSourceExtendedPars.h =================================================================== --- branches/pap/psModules/src/objects/pmSourceExtendedPars.h (revision 27708) +++ branches/pap/psModules/src/objects/pmSourceExtendedPars.h (revision 28003) @@ -19,30 +19,49 @@ psVector *theta; // angles corresponding to above radial profiles psVector *isophotalRadii; // isophotal radius for the above angles +} pmSourceRadialFlux; +typedef struct { psVector *radiusElliptical; // normalized radial coordinates for all relevant pixels psVector *fluxElliptical; // flux for the above radial coordinates +} pmSourceEllipticalFlux; +typedef struct { psVector *binSB; // mean surface brightness within radial bins psVector *binSBstdev; // scatter of mean surface brightness within radial bins psVector *binSBerror; // formal error on mean surface brightness within radial bins - - psVector *radialBins; // radii corresponding to above binnedBlux + psVector *binSum; // sum of flux within radial bins + psVector *binFill; // fraction of area actually lit + psVector *radialBins; // radii corresponding to above binnedFlux psVector *area; // differential area of the non-overlapping radial bins - - psEllipseAxes axes; // shape of elliptical contour } pmSourceRadialProfile; typedef struct { - float flux; - float fluxErr; - float radius; - float radiusErr; + float flux; + float fluxErr; + float radius; + float radiusErr; } pmSourceExtendedFlux; typedef struct { - pmSourceRadialProfile *profile; - pmSourceExtendedFlux *petrosian_50; - pmSourceExtendedFlux *petrosian_80; + pmSourceRadialFlux *radFlux; // raw radial flux information + pmSourceEllipticalFlux *ellipticalFlux; // flux for elliptically-renormalized radii + pmSourceRadialProfile *radProfile; // surface brightness profile in specified fixed bins + pmSourceRadialProfile *petProfile; // surface brightness profile in petrosian bins + psEllipseAxes axes; // shape of elliptical contour + float petrosianFlux; + float petrosianFluxErr; + float petrosianRadius; + float petrosianRadiusErr; + float petrosianR90; + float petrosianR90Err; + float petrosianR50; + float petrosianR50Err; } pmSourceExtendedPars; + +pmSourceRadialFlux *pmSourceRadialFluxAlloc(); +bool psMemCheckSourceRadialFlux(psPtr ptr); + +pmSourceEllipticalFlux *pmSourceEllipticalFluxAlloc(); +bool psMemCheckSourceEllipticalFlux(psPtr ptr); // *** pmSourceRadialProfile describes the radial profile of a source in elliptical contours, and @@ -50,7 +69,4 @@ pmSourceRadialProfile *pmSourceRadialProfileAlloc(); bool psMemCheckSourceRadialProfile(psPtr ptr); - -// *** pmSourceRadialProfileFreeVectors frees the intermediate data values -bool pmSourceRadialProfileFreeVectors(pmSourceRadialProfile *profile); // *** pmSourceExtendedPars describes the possible collection of extended flux measurements for a source @@ -65,4 +81,6 @@ bool pmSourceRadialProfileSortPair(psVector *index, psVector *extra); + + /// @} # endif /* PM_SOURCE_H */ Index: branches/pap/psModules/src/objects/pmSourceFitSet.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceFitSet.c (revision 27708) +++ branches/pap/psModules/src/objects/pmSourceFitSet.c (revision 28003) @@ -62,5 +62,5 @@ // before pmSourceFitSetInit is called for (int i = 0; i < fitSets->n; i++) { - psAssert (fitSets->data[i] == NULL, "failure to init or clear fitSets?"); + psAssert(fitSets->data[i] == NULL, "failure to init or clear fitSets?"); } return true; @@ -68,5 +68,6 @@ void pmSourceFitSetDone (void) { - psFree (fitSets); + psFree(fitSets); + fitSets = NULL; } @@ -113,38 +114,31 @@ } -pmSourceFitSetData *pmSourceFitSetDataSet (psArray *modelSet) { - - pmSourceFitSetData *thisSet = NULL; - - psAssert (fitSets, "pmSourceFitSetInit not called"); +pmSourceFitSetData *pmSourceFitSetDataSet (psArray *modelSet) +{ + psAssert(fitSets, "pmSourceFitSetInit not called"); // find the fitSet used by this thread pthread_t id = pthread_self(); - // is our ID is already on the stack, abort. - // we do not need to lock to do this.... + // If our ID is already on the stack, abort. + // We do need to lock on this because someone might pull one of the fitSets out from under us + pthread_mutex_lock(&fitSetInitMutex); for (int i = 0; i < fitSets->n; i++) { - thisSet = fitSets->data[i]; + pmSourceFitSetData *thisSet = fitSets->data[i]; if (!thisSet) continue; if (thisSet->thread == id) { - break; - } - thisSet = NULL; - } - psAssert (thisSet == NULL, "pmSourceFitSetDataSet() called but previous entry not cleared"); - - thisSet = pmSourceFitSetDataAlloc(modelSet); - - pthread_mutex_lock (&fitSetInitMutex); - - // find an entry that is NULL and place it there + psAbort("pmSourceFitSetDataSet() called but previous entry not cleared"); + } + } + + // Find an open slot for (int i = 0; i < fitSets->n; i++) { if (fitSets->data[i]) continue; - fitSets->data[i] = thisSet; - pthread_mutex_unlock (&fitSetInitMutex); + pmSourceFitSetData *thisSet = fitSets->data[i] = pmSourceFitSetDataAlloc(modelSet); + pthread_mutex_unlock(&fitSetInitMutex); return thisSet; } - pthread_mutex_unlock (&fitSetInitMutex); - psAbort ("no empty slot for new pmSourceFitSetData"); + pthread_mutex_unlock(&fitSetInitMutex); + psAbort("no empty slot for new pmSourceFitSetData"); return NULL; } @@ -152,47 +146,41 @@ pmSourceFitSetData *pmSourceFitSetDataGet (void) { - psAssert (fitSets, "pmSourceFitSetInit not called"); - - // find the fitSet used by this thread + psAssert(fitSets, "pmSourceFitSetInit not called"); + + // Find the fitSet used by this thread + // We do need to lock on this because someone might pull one of the fitSets out from under us pthread_t id = pthread_self(); - - // can we find our fitSet? we do not need to lock to do this.... - pmSourceFitSetData *thisSet = NULL; + pthread_mutex_lock(&fitSetInitMutex); for (int i = 0; i < fitSets->n; i++) { - thisSet = fitSets->data[i]; + pmSourceFitSetData *thisSet = fitSets->data[i]; if (!thisSet) continue; if (thisSet->thread == id) { - break; - } - thisSet = NULL; - } - psAssert (thisSet != NULL, "pmSourceFitSetDataGet() called, but no entry found"); - - return thisSet; -} - -void pmSourceFitSetDataClear (void) { - - int i; - + pthread_mutex_unlock(&fitSetInitMutex); + return thisSet; + } + } + psAbort("pmSourceFitSetDataGet() called, but no entry found"); +} + +void pmSourceFitSetDataClear (void) +{ psAssert (fitSets, "pmSourceFitSetInit not called"); - // find the fitSet used by this thread + // Find the fitSet used by this thread + // We do need to lock on this because someone might pull one of the fitSets out from under us pthread_t id = pthread_self(); - - // can we find our fitSet? we do not need to lock to do this.... - pmSourceFitSetData *thisSet = NULL; - for (i = 0; i < fitSets->n; i++) { - thisSet = fitSets->data[i]; + pthread_mutex_lock(&fitSetInitMutex); + for (int i = 0; i < fitSets->n; i++) { + pmSourceFitSetData *thisSet = fitSets->data[i]; if (!thisSet) continue; if (thisSet->thread == id) { - break; - } - thisSet = NULL; - } - psAssert (thisSet != NULL, "pmSourceFitSetDataClear() called, but no entry found"); - - psFree (thisSet); - fitSets->data[i] = NULL; + psFree(thisSet); + fitSets->data[i] = NULL; + pthread_mutex_unlock(&fitSetInitMutex); + return; + } + } + psAbort("pmSourceFitSetDataClear() called, but no entry found"); + return; } @@ -481,5 +469,5 @@ // Convert i/j to image space: - // 0.5 PIX: the coordinate values must be in pixel coords, not index + // 0.5 PIX: the coordinate values must be in pixel coords, not index coord->data.F32[0] = (psF32) (j + 0.5 + source->pixels->col0); coord->data.F32[1] = (psF32) (i + 0.5 + source->pixels->row0); Index: branches/pap/psModules/src/objects/pmSourceIO.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceIO.c (revision 27708) +++ branches/pap/psModules/src/objects/pmSourceIO.c (revision 28003) @@ -551,11 +551,11 @@ } if (!strcmp (exttype, "PS1_V1")) { - status &= pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe); + status &= pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); } if (!strcmp (exttype, "PS1_V2")) { - status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe); + status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); } if (!strcmp (exttype, "PS1_DV1")) { - status &= pmSourcesWrite_CMF_PS1_DV1_XSRC (file->fits, sources, xsrcname, recipe); + status &= pmSourcesWrite_CMF_PS1_DV1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); } } @@ -568,11 +568,11 @@ } if (!strcmp (exttype, "PS1_V1")) { - status &= pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname); + status &= pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, readout, sources, xfitname); } if (!strcmp (exttype, "PS1_V2")) { - status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname); + status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, readout, sources, xfitname); } if (!strcmp (exttype, "PS1_DV1")) { - status &= pmSourcesWrite_CMF_PS1_DV1_XFIT (file->fits, sources, xfitname); + status &= pmSourcesWrite_CMF_PS1_DV1_XFIT (file->fits, readout, sources, xfitname); } } Index: branches/pap/psModules/src/objects/pmSourceIO.h =================================================================== --- branches/pap/psModules/src/objects/pmSourceIO.h (revision 27708) +++ branches/pap/psModules/src/objects/pmSourceIO.h (revision 28003) @@ -36,14 +36,14 @@ bool pmSourcesWrite_CMF_PS1_V1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); -bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe); -bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, psArray *sources, char *extname); +bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname); bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); -bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe); -bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname); +bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname); bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); -bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe); -bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, psArray *sources, char *extname); +bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname); bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config); Index: branches/pap/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 27708) +++ branches/pap/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 28003) @@ -364,5 +364,5 @@ // XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) { @@ -536,5 +536,5 @@ // XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, psArray *sources, char *extname) +bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname) { Index: branches/pap/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 27708) +++ branches/pap/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 28003) @@ -339,5 +339,5 @@ // XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) { @@ -516,5 +516,5 @@ // XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, psArray *sources, char *extname) +bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname) { Index: branches/pap/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c =================================================================== --- branches/pap/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 27708) +++ branches/pap/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 28003) @@ -45,6 +45,5 @@ // followed by a zero-size matrix, followed by the table data -bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, - psMetadata *imageHeader, psMetadata *tableHeader, char *extname) +bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -54,5 +53,4 @@ psArray *table; psMetadata *row; - int i; psF32 *PAR, *dPAR; psEllipseAxes axes; @@ -83,8 +81,8 @@ pmSource *source = (pmSource *) sources->data[0]; if (source->seq == -1) { - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); } else { - sources = psArraySort (sources, pmSourceSortBySeq); + sources = psArraySort (sources, pmSourceSortBySeq); } } @@ -94,5 +92,5 @@ // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. - for (i = 0; i < sources->n; i++) { + for (int i = 0; i < sources->n; i++) { pmSource *source = (pmSource *) sources->data[i]; @@ -102,5 +100,5 @@ // generated on Alloc, and would thus be wrong for read in sources. if (source->seq == -1) { - source->seq = i; + source->seq = i; } @@ -114,10 +112,10 @@ yPos = PAR[PM_PAR_YPOS]; if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; } else { - // in linear-fit mode, there is no error on the centroid - xErr = source->peak->dx; - yErr = source->peak->dy; + // in linear-fit mode, there is no error on the centroid + xErr = source->peak->dx; + yErr = source->peak->dy; } if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { @@ -212,4 +210,5 @@ } + // XXX why do we make a copy here to be supplemented with the masks? why not do this in the calling function? psMetadata *header = psMetadataCopy(NULL, tableHeader); pmSourceMasksHeader(header); @@ -344,6 +343,5 @@ } -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) { @@ -354,4 +352,5 @@ psF32 xPos, yPos; psF32 xErr, yErr; + int nRow = -1; // create a header to hold the output data @@ -361,4 +360,20 @@ psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); + pmChip *chip = readout->parent->parent; + pmFPA *fpa = chip->parent; + + // zero point corrections + bool status1 = false; + bool status2 = false; + float magOffset = 0.0; + float exptime = psMetadataLookupF32(&status1, fpa->concepts, "FPA.EXPOSURE"); + float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); + if (!isfinite(zeropt)) { + zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); + } + if (status1 && status2 && (exptime > 0.0)) { + magOffset = zeropt + 2.5*log10(exptime); + } + // let's write these out in S/N order sources = psArraySort (sources, pmSourceSortBySN); @@ -367,12 +382,24 @@ // which extended source analyses should we perform? - // bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); - // bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); - psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); - psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); - assert (radialBinsLower->n == radialBinsUpper->n); + psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + psAssert (radMin->n == radMax->n, "inconsistent annular bins"); + + // int nRadialBins = 0; + // if (doAnnuli) { + // // get the max count of radial bins + // for (int i = 0; i < sources->n; i++) { + // pmSource *source = sources->data[i]; + // if (!source->extpars) continue; + // if (!source->extpars->radProfile ) continue; + // if (!source->extpars->radProfile->binSB) continue; + // nRadialBins = PS_MAX(nRadialBins, source->extpars->radProfile->binSB->n); + // } + // } // we write out all sources, regardless of quality. the source flags tell us the state @@ -410,14 +437,31 @@ psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); -# if (0) + float AxialRatio = NAN; + float AxialTheta = NAN; + pmSourceExtendedPars *extpars = source->extpars; + if (extpars) { + AxialRatio = extpars->axes.minor / extpars->axes.major; + AxialTheta = extpars->axes.theta; + } + psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO", PS_DATA_F32, "Axial Ratio of radial profile", AxialRatio); + psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA", PS_DATA_F32, "Angle of radial profile ellipse", AxialTheta); + // Petrosian measurements // XXX insert header data: petrosian ref radius, flux ratio + // XXX check flags to see if Pet was measured if (doPetrosian) { - pmSourcePetrosianValues *petrosian = source->extpars->petrosian; - if (petrosian) { - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", petrosian->mag); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", petrosian->rad); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", petrosian->radErr); + pmSourceExtendedPars *extpars = source->extpars; + if (extpars) { + // XXX note that this mag is either calibrated or instrumental depending on existence of zero point + float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point + float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", mag); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", magErr); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", extpars->petrosianR50); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)", extpars->petrosianR50Err); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err); } else { psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); @@ -425,7 +469,12 @@ psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)",NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN); } } +# if (0) // Kron measurements if (doKron) { @@ -460,67 +509,74 @@ } } - - // Flux Annuli +# endif + + // Flux Annuli (if we have extended source measurements, we have these. only optionally save them) if (doAnnuli) { - pmSourceAnnuli *annuli = source->extpars->annuli; - if (annuli) { - psVector *fluxVal = annuli->flux; - psVector *fluxErr = annuli->fluxErr; - psVector *fluxVar = annuli->fluxVar; - - for (int j = 0; j < fluxVal->n; j++) { - char name[32]; - sprintf (name, "FLUX_VAL_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]); - sprintf (name, "FLUX_ERR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]); - sprintf (name, "FLUX_VAR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]); - } - } else { - for (int j = 0; j < radialBinsLower->n; j++) { - char name[32]; - sprintf (name, "FLUX_VAL_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN); - sprintf (name, "FLUX_ERR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN); - sprintf (name, "FLUX_VAR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN); - } - } - } - -# endif - psArrayAdd (table, 100, row); - psFree (row); - } - + psVector *radSB = psVectorAlloc(radMin->n, PS_TYPE_F32); + psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32); + psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32); + psVectorInit (radSB, NAN); + psVectorInit (radFlux, NAN); + psVectorInit (radFill, NAN); + if (!source->extpars) goto empty_annuli; + if (!source->extpars->radProfile) goto empty_annuli; + if (!source->extpars->radProfile->binSB) goto empty_annuli; + psAssert (source->extpars->radProfile->binSum, "programming error"); + psAssert (source->extpars->radProfile->binFill, "programming error"); + psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths"); + psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths"); + psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths"); + + // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux + for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) { + radSB->data.F32[j] = source->extpars->radProfile->binSB->data.F32[j]; + radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j]; + radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j]; + } + + empty_annuli: + psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB); + psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); + psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); + psFree (radSB); + psFree (radFlux); + psFree (radFill); + } + if (nRow < 0) { + nRow = row->list->n; + } else { + psAssert (nRow == row->list->n, "inconsistent row lengths"); + } + psArrayAdd (table, 100, row); + psFree (row); + } + if (table->n == 0) { - if (!psFitsWriteBlank (fits, outhead, extname)) { - psError(psErrorCodeLast(), false, "Unable to write empty sources file."); - psFree(outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; - } - + if (!psFitsWriteBlank (fits, outhead, extname)) { + psError(psErrorCodeLast(), false, "Unable to write empty sources file."); + psFree(outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; + } + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); - psFree (outhead); - psFree(table); - return false; + psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; } psFree (outhead); psFree (table); - + return true; } // XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname) +bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname) { Index: branches/pap/psModules/src/psmodules.h =================================================================== --- branches/pap/psModules/src/psmodules.h (revision 27708) +++ branches/pap/psModules/src/psmodules.h (revision 28003) @@ -10,4 +10,6 @@ #include #include +#include +#include // XXX the following headers define constructs needed by the elements below