Index: trunk/psphot/src/Makefile.am =================================================================== --- trunk/psphot/src/Makefile.am (revision 36375) +++ trunk/psphot/src/Makefile.am (revision 36441) @@ -25,5 +25,5 @@ libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -bin_PROGRAMS = psphot psphotForced psphotFullForce psphotMinimal psphotMakePSF psphotStack psphotModelTest psmakecff +bin_PROGRAMS = psphot psphotForced psphotFullForce psphotFullForceSummary psphotMinimal psphotMakePSF psphotStack psphotModelTest psmakecff # bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy @@ -39,4 +39,8 @@ psphotFullForce_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) psphotFullForce_LDADD = libpsphot.la + +psphotFullForceSummary_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +psphotFullForceSummary_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) +psphotFullForceSummary_LDADD = libpsphot.la psphotMinimal_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) @@ -99,4 +103,9 @@ psphotMosaicChip.c \ psphotCleanup.c + +# combine full force results from several inputs +psphotFullForceSummary_SOURCES = \ + psphotFullForceSummary.c \ + psphotFullForceSummaryReadout.c # forced photometry of specified positions given a specified psf Index: trunk/psphot/src/psmakecff.c =================================================================== --- trunk/psphot/src/psmakecff.c (revision 36375) +++ trunk/psphot/src/psmakecff.c (revision 36441) @@ -2,4 +2,7 @@ # include # endif + +// psmakecff : A program to make read a cmf file and write a cff file +// The real work is done in psModules. #include @@ -8,4 +11,5 @@ #include "psphot.h" +// For simplicilty, this program's (simple) functions are all contained in this file. static pmConfig* psmakecffArguments(int, char**); static bool psmakecffParseCamera(pmConfig *); @@ -25,5 +29,4 @@ // psphotVersionPrint(); - // load input data (config and images (signal, noise, mask) if (!psmakecffParseCamera (config)) { psErrorStackPrint(stderr, "Error setting up the camera\n"); @@ -31,5 +34,4 @@ } - // call psphot for each readout if (!psmakecffImageLoop (config)) { psErrorStackPrint(stderr, "Error in the psphot image loop\n"); @@ -37,8 +39,4 @@ } -// psLogMsg ("psphot", PS_LOG_WARN, "complete psphot run: %f sec\n", psTimerMark ("complete")); - -// psErrorCode exit_status = psphotGetExitStatus(); -// psphotCleanup (config); exit (0); } @@ -73,4 +71,5 @@ return config; } + static bool psmakecffParseCamera(pmConfig *config) { bool status = false; @@ -127,28 +126,24 @@ psAssert (recipe, "missing recipe?"); -// psImageMaskType maskTest = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // for psphot, we force data to be read at the chip level while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { - psLogMsg ("psmakecmf", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + psLogMsg ("psmakecff", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); if (! chip->process || ! chip->file_exists) { continue; } -#ifdef notdef - pmFPAfileActivate (config->files, false, NULL); - pmFPAfileActivate (config->files, true, "PSPHOT.LOAD"); -#endif if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psmakecff."); // there is now only a single chip (multiple readouts?). loop over it and process while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { - psLogMsg ("psmakecmf", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + psLogMsg ("psmakecff", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); // process each of the readouts while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { - psLogMsg ("psmakecmf", 6, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + psLogMsg ("psmakecff", 6, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); if (! readout->data_exists) { continue; } pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell); if (hdu && hdu != lastHDU) { + // XXX: probably should do this // psphotVersionHeaderFull(hdu->header); lastHDU = hdu; @@ -157,12 +152,9 @@ } - // Defer output and closing of files until we've (possibly) done the NFrames analysis below -// pmFPAfileActivate (config->files, false, NULL); - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for Chip in psmakecff."); } if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for FPA in psphot."); - // fail if we failed to handle an error + // fail if we encountered an unhandled error if (psErrorCodeLast() != PS_ERR_NONE) psAbort ("failed to handle an error!"); Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 36375) +++ trunk/psphot/src/psphot.h (revision 36441) @@ -361,4 +361,6 @@ bool psphotModelTestReadout(pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotFullForceSummaryReadout (pmConfig * config, const pmFPAview *view); + int psphotFileruleCount(const pmConfig *config, const char *filerule); bool psphotFileruleCountSet(const pmConfig *config, const char *filerule, int num); Index: trunk/psphot/src/psphotChooseAnalysisOptions.c =================================================================== --- trunk/psphot/src/psphotChooseAnalysisOptions.c (revision 36375) +++ trunk/psphot/src/psphotChooseAnalysisOptions.c (revision 36441) @@ -60,5 +60,5 @@ } -bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos); +static bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos); // this function use an internal flag to mark sources which have already been measured @@ -121,9 +121,21 @@ psSphereRot *toGal = NULL; float GAL_LIMIT = 0.0; + float GAL_LIMIT_BULGE = 0.0; + float GAL_LIMIT_SIGMA2 = 0.0; bool useGAL_LIMIT = psMetadataLookupBool(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.USE"); assert (status); if (useGAL_LIMIT) { + toGal = psSphereRotICRSToGalactic(); GAL_LIMIT = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT"); - toGal = psSphereRotICRSToGalactic(); + assert (status); + GAL_LIMIT = DEG_TO_RAD(GAL_LIMIT); + GAL_LIMIT_BULGE = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.BULGE"); + assert (status); + GAL_LIMIT_BULGE = DEG_TO_RAD(GAL_LIMIT_BULGE); + float GAL_LIMIT_SIGMA = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.BULGE.SIGMA"); + assert (status); + assert(GAL_LIMIT_SIGMA > 0); + GAL_LIMIT_SIGMA = DEG_TO_RAD(GAL_LIMIT_SIGMA); + GAL_LIMIT_SIGMA2 = GAL_LIMIT_SIGMA * GAL_LIMIT_SIGMA; } @@ -135,4 +147,5 @@ } + float petroFluxLim = NAN; float extFitFluxLim = NAN; @@ -140,4 +153,6 @@ float extFitMagLimDefault = NAN; float extFitMagLim = NAN; + float petroMagLimDefault = NAN; + float petroMagLim = NAN; // match to the given filter @@ -156,17 +171,26 @@ // find a matching filter or default to 'any' if (!strcasecmp (thisFilter, "any")) { - psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT"); - psAssert(extFitMagLimStr, "missing MAG.LIMIT"); + psString petroMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.PETRO"); + psAssert(petroMagLimStr, "missing MAG.LIMIT.PETRO"); + petroMagLimDefault = atof (petroMagLimStr); + + psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.EXTFIT"); + psAssert(extFitMagLimStr, "missing MAG.LIMIT.EXTFIT"); extFitMagLimDefault = atof (extFitMagLimStr); } - // find a matching filter or default to 'any' if (!strcasecmp (thisFilter, filterID)) { - psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT"); - psAssert(extFitMagLimStr, "missing MAG.LIMIT"); + psString petroMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.PETRO"); + psAssert(petroMagLimStr, "missing MAG.LIMIT.PETRO"); + petroMagLim = atof (petroMagLimStr); + + psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.EXTFIT"); + psAssert(extFitMagLimStr, "missing MAG.LIMIT.EXTFIT"); extFitMagLim = atof (extFitMagLimStr); break; } } + psFree(iter); + if (!isfinite (petroMagLim)) petroMagLim = petroMagLimDefault; if (!isfinite (extFitMagLim)) extFitMagLim = extFitMagLimDefault; @@ -185,4 +209,5 @@ psAssert (status, "missing FPA.ZP?"); + petroFluxLim = exptime * pow (10.0, 0.4*(zeropt - petroMagLim)); extFitFluxLim = exptime * pow (10.0, 0.4*(zeropt - extFitMagLim)); } @@ -225,18 +250,27 @@ psSphere ptGal, ptSky; GetGalacticCoords (&ptGal, &ptSky, toGal, chip, source->peak->xf, source->peak->yf); - if (fabs(ptGal.d) < GAL_LIMIT) continue; + float b = ptGal.r; + float limit = GAL_LIMIT + GAL_LIMIT_BULGE * exp(-0.5*(b*b/GAL_LIMIT_SIGMA2)); + if (fabs(ptGal.d) < limit) continue; // include an exception for low density skycells below the limit? } // for petro and extFit, we will either use the mag limits or the S/N + if (doPetrosian) { + if (isfinite(petroFluxLim)) { + if (source->moments->KronFlux > petroFluxLim) { + source->tmpFlags |= PM_SOURCE_TMPF_PETRO; + } + } else if (source->moments->KronFlux > SN_LIM * source->moments->KronFluxErr) { + source->tmpFlags |= PM_SOURCE_TMPF_PETRO; + } + } if (isfinite(extFitFluxLim)) { if (source->moments->KronFlux > extFitFluxLim) { source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT; - if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO; } } else { if (source->moments->KronFlux > SN_LIM * source->moments->KronFluxErr) { source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT; - if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO; } } @@ -244,4 +278,6 @@ psLogMsg ("psphot.options", PS_LOG_WARN, "choose analysis options for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.options")); + + psFree(toGal); return true; @@ -287,9 +323,21 @@ psSphereRot *toGal = NULL; float GAL_LIMIT = 0.0; + float GAL_LIMIT_BULGE = 0.0; + float GAL_LIMIT_SIGMA2 = 0.0; bool useGAL_LIMIT = psMetadataLookupBool(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.USE"); assert (status); if (useGAL_LIMIT) { + toGal = psSphereRotICRSToGalactic(); GAL_LIMIT = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT"); - toGal = psSphereRotICRSToGalactic(); + assert (status); + GAL_LIMIT = DEG_TO_RAD(GAL_LIMIT); + GAL_LIMIT_BULGE = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.BULGE"); + assert (status); + GAL_LIMIT_BULGE = DEG_TO_RAD(GAL_LIMIT_BULGE); + float GAL_LIMIT_SIGMA = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.BULGE.SIGMA"); + assert (status); + assert(GAL_LIMIT_SIGMA > 0); + GAL_LIMIT_SIGMA = DEG_TO_RAD(GAL_LIMIT_SIGMA); + GAL_LIMIT_SIGMA2 = GAL_LIMIT_SIGMA * GAL_LIMIT_SIGMA; } @@ -365,9 +413,13 @@ // find extFitFluxLim->data.F32[i] for i == image number psVector *extFitFluxLim = NULL; + psVector *petroFluxLim = NULL; if (magLimits) { extFitFluxLim = psVectorAlloc (inputFilters->n, PS_TYPE_F32); + petroFluxLim = psVectorAlloc (inputFilters->n, PS_TYPE_F32); psVectorInit (extFitFluxLim, NAN); + psVectorInit (petroFluxLim, NAN); float extFitMagLimDefault = NAN; + float petroMagLimDefault = NAN; // match mag limits (flux limits) to the filters @@ -383,7 +435,11 @@ // save the default value to assign to unset filters if (!strcasecmp (thisFilter, "any")) { - psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT"); - psAssert(magString, "missing MAG.LIMIT"); + psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.EXTFIT"); + psAssert(magString, "missing MAG.LIMIT.EXTFIT"); extFitMagLimDefault = atof (magString); + + magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.PETRO"); + psAssert(magString, "missing MAG.LIMIT.PETRO"); + petroMagLimDefault = atof (magString); continue; } @@ -393,7 +449,12 @@ if (i == chisqNum) continue; if (!strcasecmp (thisFilter, inputFilters->data[i])) { - psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT"); - psAssert(magString, "missing MAG.LIMIT"); + psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.PETRO"); + psAssert(magString, "missing MAG.LIMIT.PETRO"); float magvalue = atof (magString); + petroFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - magvalue)); + + magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.EXTFIT"); + psAssert(magString, "missing MAG.LIMIT.EXTFIT"); + magvalue = atof (magString); extFitFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - magvalue)); break; @@ -401,9 +462,14 @@ } } + psFree(iter); for (int i = 0; i < num; i++) { if (i == chisqNum) continue; - if (isfinite(extFitFluxLim->data.F32[i])) continue; - extFitFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - extFitMagLimDefault)); + if (!isfinite(petroFluxLim->data.F32[i])) { + petroFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - petroMagLimDefault)); + } + if (!isfinite(extFitFluxLim->data.F32[i])) { + extFitFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - extFitMagLimDefault)); + } } } @@ -420,5 +486,6 @@ bool doObjectRadial = false; bool doObjectExtFit = false; - for (int j = 0; !doObjectExtFit && !doObjectRadial && (j < object->sources->n); j++) { + bool doObjectPetrosian = false; + for (int j = 0; !doObjectExtFit && !doObjectRadial && !doObjectPetrosian && (j < object->sources->n); j++) { pmSource *source = object->sources->data[j]; @@ -451,4 +518,5 @@ if (extFitAll) { doObjectExtFit = true; + doObjectPetrosian = doPetrosian; continue; } @@ -465,9 +533,24 @@ psSphere ptGal, ptSky; GetGalacticCoords (&ptGal, &ptSky, toGal, chips->data[imageID], source->peak->xf, source->peak->yf); - if (fabs(ptGal.d) < GAL_LIMIT) continue; - // include an exception for low density skycells below the limit? - } - - float fluxLim = extFitFluxLim ? extFitFluxLim->data.F32[imageID] : NAN; + float b = ptGal.r; + float limit = GAL_LIMIT + GAL_LIMIT_BULGE * exp(-0.5*(b*b/GAL_LIMIT_SIGMA2)); + if (fabs(ptGal.d) < limit) continue; + } + + float fluxLim = NAN; + // for petro and extFit, we will either use the mag limits or the S/N + if (doPetrosian) { + fluxLim = petroFluxLim ? petroFluxLim->data.F32[imageID] : NAN; + if (isfinite(fluxLim)) { + if (source->moments->KronFlux > extFitFluxLim->data.F32[imageID]) { + doObjectPetrosian = true; + } + } else { + if (source->moments->KronFlux > SN_LIM * source->moments->KronFluxErr) { + doObjectPetrosian = true; + } + } + } + fluxLim = extFitFluxLim ? extFitFluxLim->data.F32[imageID] : NAN; // for petro and extFit, we will either use the mag limits or the S/N if (isfinite(fluxLim)) { @@ -503,9 +586,10 @@ if (source->modelPSF == NULL) continue; - // Do the fits if the recipe requests we do extended source fits to everything if (doObjectExtFit) { source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT; - if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO; - } + } + if (doObjectPetrosian) { + source->tmpFlags |= PM_SOURCE_TMPF_PETRO; + } // Do the fits if the recipe requests we do extended source fits to everything @@ -522,8 +606,10 @@ psLogMsg ("psphot.options", PS_LOG_WARN, "choose analysis options for %ld objects: %f sec\n", objects->n, psTimerMark ("psphot.options")); + psFree(toGal); + return true; } -bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos) { +static bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos) { pmFPA *fpa = chip->parent; @@ -543,4 +629,9 @@ psSphereRotApply (ptGal, toGal, ptSky); + // psSphereRotApply insures that 0 < r < 2PI. We want -PI < b <= PI + if (ptGal->r > M_PI) { + ptGal->r -= 2.0 * M_PI; + } + return true; Index: trunk/psphot/src/psphotFullForceSummary.c =================================================================== --- trunk/psphot/src/psphotFullForceSummary.c (revision 36441) +++ trunk/psphot/src/psphotFullForceSummary.c (revision 36441) @@ -0,0 +1,266 @@ +# ifdef HAVE_CONFIG_H +# include +# endif + +#include +#include +#include +#include "psphot.h" + +// For simplicilty, this program's (simple) functions are all contained in this file. +static pmConfig* psphotFullForceSummaryArguments(int, char**); +static bool psphotFullForceSummaryParseCamera(pmConfig *); +static bool psphotFullForceSummaryImageLoop(pmConfig*); + +int main (int argc, char **argv) { + + psMemInit(); + psTimerStart ("complete"); + pmErrorRegister(); // register psModule's error codes/messages + psphotInit(); + + // load command-line arguments, options, and system config data + pmConfig *config = psphotFullForceSummaryArguments (argc, argv); + assert(config); + +// psphotVersionPrint(); + + if (!psphotFullForceSummaryParseCamera (config)) { + psErrorStackPrint(stderr, "Error setting up the camera\n"); + exit (1); + } + + if (!psphotFullForceSummaryImageLoop (config)) { + psErrorStackPrint(stderr, "Error in the psphot image loop\n"); + exit(1); + } + + // XXX:check for memory leaks + exit (0); +} + +// all functions which return to this level must raise one of the top-level error codes if they +// exit with an error. these error codes are used to specify the program exit status + +void usage() { + fprintf(stderr, "usage: psphotFullForceSummary -inputs \n"); + exit (1); +} + +static pmConfig* psphotFullForceSummaryArguments(int argc, char **argv) { + + pmConfig *config = pmConfigRead(&argc, argv, PSPHOT_RECIPE); + int N; + + if (config == NULL) { + psErrorStackPrint(stderr, "Can't read site configuration"); + exit(PS_EXIT_CONFIG_ERROR); + } + if ((N = psArgumentGet(argc, argv, "-input"))) { + if (argc <= N+1) { + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); + usage(); + } + psArgumentRemove(N, &argc, argv); + + unsigned int numBad = 0; // Number of bad lines + psMetadata *inputs = psMetadataConfigRead(NULL, &numBad, argv[N], false); // Input file info + if (!inputs || numBad > 0) { + psErrorStackPrint(stderr, "Unable to cleanly read MDC file with inputs."); + exit(PS_EXIT_CONFIG_ERROR); + } + psMetadataAddMetadata(config->arguments, PS_LIST_TAIL, "INPUTS", 0, "Metadata with input details", inputs); + psFree(inputs); + + psArgumentRemove(N, &argc, argv); + } else { + psErrorStackPrint(stderr, "-input must be supplied."); + usage(); + } + + if (argc < 2) { + psErrorStackPrint(stderr, "Output is required."); + usage(); + } + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", PS_DATA_STRING, "", argv[1]); + + return config; +} + +static bool psphotFullForceSummaryParseCamera(pmConfig *config) { + bool status = false; + + psMetadata *inputs = psMetadataLookupMetadata(&status, config->arguments, "INPUTS"); // The inputs info + if (!inputs) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find inputs."); + return false; + } + + int nInputs = inputs->list->n; + + for (int i = 0; i < nInputs; i++) { + psMetadataItem *item = psMetadataGet(inputs, i); + + if (item->type != PS_DATA_STRING) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Component %s of the input metadata is not of type STRING", item->name); + return false; + + } + psString sourcesFilename = item->data.str; + + psArray *dummy = psArrayAlloc(1); // dummy array of filenames + dummy->data[0] = psStringCopy(sourcesFilename); + + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "FILENAMES", PS_META_REPLACE, + "Filenames for file rule definition", dummy); + psFree(dummy); + + bool found = false; + pmFPAfile *file = pmFPAfileDefineFromArgs(&found, config, "PSPHOT.INPUT.CMF", "FILENAMES"); + if (!file || !found) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file %s from %s", "PSPHOT.INPUT.CMF", sourcesFilename); + return false; + } + if (file->type != PM_FPA_FILE_CMF) { + psError(PS_ERR_IO, true, "%s is not of type %s", sourcesFilename, pmFPAfileStringFromType(PM_FPA_FILE_CMF)); + return false; + } + } + psMetadataRemoveKey(config->arguments, "FILENAMES"); + psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.CMF.NUM", PS_META_REPLACE, "number of inputs", + nInputs); + + + pmFPA *outputFPA = pmFPAConstruct(config->camera, config->cameraName); + if (!outputFPA) { + psError(psErrorCodeLast(), false, "Unable to construct an FPA from camera configuration."); + return false; + } + pmFPAfile *output = pmFPAfileDefineOutput(config, outputFPA, "PSPHOT.FULLFORCE.OUTPUT"); + psFree(outputFPA); // Drop reference + if (!output) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PSPHOT.FULLFORCE.OUTPUT")); + return false; + } + if (output->type != PM_FPA_FILE_CMF) { + psError(PSPHOT_ERR_CONFIG, true, "PSPHOT.FULLFORCE.OUTPUT is not of type CMF"); + return false; + } + output->save = true; + + return true; +} + +# define ESCAPE(MESSAGE) { \ + psError(PSPHOT_ERR_DATA, false, MESSAGE); \ + psFree (view); \ + return false; \ + } + +bool psphotFullForceSummaryImageLoop (pmConfig *config) { + + bool status; + pmChip *chip; + pmCell *cell; + pmReadout *readout; + + pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT.CMF"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "Can't find input data!"); + return false; + } + + pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSPHOT.FULLFORCE.OUTPUT"); + if (!output) { + psError(PSPHOT_ERR_PROG, false, "Can't find output data!"); + return false; + } + + pmFPAview *view = pmFPAviewAlloc (0); + if (!pmFPAAddSourceFromView(output->fpa, view, output->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to insert HDU into FPA for writing.\n"); + psFree(view); + return NULL; + } + + // files associated with the science image + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot."); + + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { +// psLogMsg ("psphotFullForceSummary", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (! chip->process || ! chip->file_exists) { continue; } + + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psphotFullForceSummary."); + + // We read the WCS from the first input + { + pmHDU *hduLow = pmHDUGetLowest(input->fpa, chip, NULL); + if (hduLow && !pmAstromReadWCS(input->fpa, chip, hduLow->header, 1.0)) { + psWarning("Unable to read WCS astrometry from header."); + psErrorClear(); + pmHDU *hduHigh = pmHDUGetHighest(input->fpa, chip, NULL); + if (hduHigh && hduHigh != hduLow && + !pmAstromReadWCS(input->fpa, chip, hduHigh->header, 1.0)) { + psWarning("Unable to read WCS astrometry from primary header."); + psErrorClear(); + } + } + } + + // Copy the transformations from the input + output->fpa->fromTPA = psMemIncrRefCounter(input->fpa->fromTPA); + output->fpa->toTPA = psMemIncrRefCounter(input->fpa->toTPA); + output->fpa->toSky = psMemIncrRefCounter(input->fpa->toSky); + pmChip *outputChip = pmFPAviewThisChip(view, output->fpa); + outputChip->toFPA = psMemIncrRefCounter(chip->toFPA); + outputChip->fromFPA = psMemIncrRefCounter(chip->fromFPA); + if (output->fpa->hdu->header == NULL) { + output->fpa->hdu->header = psMetadataAlloc(); + } + // XXX: how come psphot and psphotStack don't have to do this? + if (!pmAstromWriteWCS(output->fpa->hdu->header, output->fpa, outputChip, 0.001)) { + ESCAPE("failure to copy WCS to header"); + } + + // there is now only a single chip (multiple readouts?). loop over it and process + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + // psLogMsg ("psphotFullForceSummary", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { +// psLogMsg ("psphotFullForceSummary", 6, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (! readout->data_exists) { continue; } + + if (!psphotFullForceSummaryReadout(config, view)) { + ESCAPE ("failure in psphotFullForceSummaryReadout"); + } + } + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for Chip in psphotFullForceSummary."); + } + + + // If these keywords are not set we get a warning message on output. Since we are not copying the input header + // and do not have an image pmFPA stuff doesn't have values for these + // XXX: What other keywords and concepts should be set (or copied)? + int numCols = psMetadataLookupS32(&status, input->fpa->hdu->header, "IMNAXIS1"); + int numRows = psMetadataLookupS32(&status, input->fpa->hdu->header, "IMNAXIS2"); + psMetadataAddS32(output->fpa->hdu->header, PS_LIST_TAIL, "IMNAXIS1", PS_META_REPLACE, "", numCols); + psMetadataAddS32(output->fpa->hdu->header, PS_LIST_TAIL, "IMNAXIS2", PS_META_REPLACE, "", numRows); + + // XXX: Also add psphot version information + + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for FPA in psphot."); + + // fail if we encountered an unhandled error + if (psErrorCodeLast() != PS_ERR_NONE) psAbort ("failed to handle an error!"); + + psFree (view); + + return true; +} Index: trunk/psphot/src/psphotFullForceSummaryReadout.c =================================================================== --- trunk/psphot/src/psphotFullForceSummaryReadout.c (revision 36441) +++ trunk/psphot/src/psphotFullForceSummaryReadout.c (revision 36441) @@ -0,0 +1,393 @@ +#include "psphotInternal.h" + + +typedef struct { + int numTrials; + psF32 fRmajorMin; + psF32 fRmajorMax; + psF32 fRmajorDel; + psF32 fRminorMin; + psF32 fRminorMax; + psF32 fRminorDel; + psVector *rMajor; + psVector *rMinor; + psArray *zeroPt; + psArray *exptime; +} galaxyShapeOptions; + +static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeOptions *options); +static pmPhotObj *findObjectForSource(psArray **pObjects, pmSource *source); + +static bool setOptions(galaxyShapeOptions *options, pmReadout *readout, psMetadata *recipe, bool saveVectors); +static bool checkOptions(galaxyShapeOptions *options, pmReadout *readout, psMetadata *recipe); + + +bool psphotFullForceSummaryReadout (pmConfig *config, const pmFPAview *view) { + + bool status; + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psphotFileruleCount(config, "PSPHOT.INPUT.CMF"); + + pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSPHOT.FULLFORCE.OUTPUT"); + pmCell *outputCell = pmFPAviewThisCell(view, output->fpa); + pmReadout *outputReadout = pmFPAviewThisReadout(view, output->fpa); + if (!outputReadout) { + outputReadout = pmReadoutAlloc(outputCell); + } + + // Get the exposure parameters for the output from recipe and set them on the output + psF32 outputZeroPoint = psMetadataLookupF32(&status, recipe, "PSPHOT.FULLFORCE.ZERO_PT"); + psF32 outputExptime = psMetadataLookupF32(&status, recipe, "PSPHOT.FULLFORCE.EXPTIME"); + + psMetadataAddF32(output->fpa->concepts, PS_LIST_TAIL, "FPA.ZP", PS_META_REPLACE, "Magnitude zero point", + outputZeroPoint); + psMetadataAddF32(outputCell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE, "Exposure time (sec)", + outputExptime); + // don't think this one matters +// psMetadataAddF32(output->fpa->concepts, PS_LIST_TAIL, "FPA.EXPOSURE", PS_META_REPLACE, "Exposure time (sec)", + // outputExptime); + + // Create objects from the various input's sources + // loop over the available readouts +// psArray *readouts = psArrayAlloc(num); + psVector *fluxScaleFactor = psVectorAlloc(num, PS_TYPE_F32); + galaxyShapeOptions options; + psArray *objects = NULL; + for (int index = 0; index < num; index++) { + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT.CMF", index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + if (index == 0) { + // Get the galaxy shape recipe values, from the analysis if present + // or from the recipe if not + if (!setOptions(&options, readout, recipe, true)) { + psError (PS_ERR_UNKNOWN, false, "problem determining galaxy shape options."); + return false; + } + } else { + // Make sure that this input was created with the same galaxy shapes recipe + if (!checkOptions(&options, readout, recipe)) { + psError (PS_ERR_UNKNOWN, false, "galaxy shape options do not match for input %d", index); + return false; + } + } + + // look up zero point + psF32 zero_point = psMetadataLookupF32(&status, readout->parent->parent->parent->concepts, "FPA.ZP"); + psF32 exptime = psMetadataLookupF32(&status, readout->parent->parent->parent->concepts, "FPA.EXPOSURE"); + + psF32 scaleFactor = pow(10, 0.4 * (outputZeroPoint - zero_point)) * outputExptime / exptime; + fluxScaleFactor->data.F32[index] = scaleFactor; + +// readouts->data[index] = readout; + + // find detections + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + sources = psArraySort (sources, pmSourceSortBySeq); + + if (objects == NULL) { + pmSource *lastSource = sources->data[sources->n - 1]; + psAssert(lastSource, "last source is null!"); + objects = psArrayAlloc(lastSource->seq + 1); + } + + for (int i=0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + source->imageID = index; + findObjectForSource(&objects, source); + } + } + + pmDetections *outputDetections = pmDetectionsAlloc(); + if (!psMetadataAddPtr (outputReadout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", + PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", outputDetections)) { + psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout"); + return false; + } + + psArray *outputSources = psArrayAllocEmpty (objects->n); + + // Loop over objects and compute the summaries + long nObjects = 0; + long nSources = 0; + for (int i=0 ; in; i++) { + pmPhotObj *obj = objects->data[i]; + if (!obj) continue; + + ++nObjects; + pmSource *source = psphotFullForceSummarizeObject(config, obj, fluxScaleFactor, &options); + if (source) { + psArrayAdd (outputSources, 100, source); + ++nSources; + } + } + + psLogMsg("psphot", PS_LOG_INFO, "constructed %ld output sources, from %ld objects.", nSources, nObjects); + + psFree(fluxScaleFactor); + + if (nSources) { + // We have data + outputDetections->allSources = outputSources; + outputReadout->data_exists = true; + outputReadout->parent->data_exists = true; + outputReadout->parent->parent->data_exists = true; + } else { + // XXX: tooo set a quality or fault code + return false; + } + + return true; +} + +static pmPhotObj *findObjectForSource(psArray **pObjects, pmSource *source) { + int seq = source->seq; + + psArray *objects = *pObjects; + if (seq >= objects->n) { + // We need to expand the object array. Kind of suprising. + objects = *pObjects = psArrayRealloc(objects, seq+1); + } + + // Look up object for this seq + pmPhotObj *obj = objects->data[seq]; + if (!obj) { + // not found allocate one + obj = pmPhotObjAlloc(); + objects->data[seq] = obj; + } + pmPhotObjAddSource(obj, source); + + return obj; +} + +static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeOptions *options) { + + pmSource *outSrc = NULL; + + psVector *sumWeightedFlux = NULL; + psVector *sumInvSig2 = NULL; + psVector *numerator = NULL; +// psVector *tmp = NULL; + psF32 totalNPix = 0; + long vectorLength = 0; + + for (int i=0; i < obj->sources->n; i++) { + pmSource *source = obj->sources->data[i]; + + // XXX: get cut from recipe + if (source->pixWeightNotPoor < .9) continue; + + // For now just start the output source as a copy of the first input source that makes cuts + if (!outSrc) { + outSrc = pmSourceCopy(source); + // This copies + // the peak + // the moments which are mostly nan except for Mrf Mx, My, and some of the kron parameters + // + // type, mode, flags + // magnitudes + outSrc->imageID = 0; + outSrc->seq = source->seq; + + if (source->modelPSF) { + outSrc->modelPSF = psMemIncrRefCounter(source->modelPSF); + } + if (source->extpars) { + outSrc->extpars = psMemIncrRefCounter(source->extpars); + } + if (source->modelEXT) { + outSrc->modelEXT = psMemIncrRefCounter(source->modelEXT); + } + if (source->modelFits) { + outSrc->modelFits = psMemIncrRefCounter(source->modelFits); + } + } + + if (source->galaxyFits && isfinite(source->galaxyFits->nPix) && source->galaxyFits->chisq->n) { + if (numerator == NULL) { + vectorLength = source->galaxyFits->chisq->n; + // tmp = psVectorAlloc(vectorLength, PS_TYPE_F32); + sumWeightedFlux = psVectorAlloc(vectorLength, PS_TYPE_F32); + psVectorInit(sumWeightedFlux, 0.0); + sumInvSig2 = psVectorAlloc(vectorLength, PS_TYPE_F32); + psVectorInit(sumInvSig2, 0.0); + numerator = psVectorAlloc(vectorLength, PS_TYPE_F32); + psVectorInit(numerator, 0.0); + } + + psAssert(vectorLength == options->numTrials, "length of chisq vector %ld does not match options %d", + vectorLength, options->numTrials); + psAssert(source->galaxyFits->chisq->n == vectorLength, "length of chisq vectors do not match %ld %ld", + source->galaxyFits->chisq->n, vectorLength); + + + psF32 scaleFactor = fluxScaleFactor->data.F32[source->imageID]; + + for (int k = 0; k < vectorLength; k++) { + + psF32 chisq = source->galaxyFits->chisq->data.F32[k]; + psF32 flux = source->galaxyFits->Flux->data.F32[k] * scaleFactor; + psF32 dFlux = source->galaxyFits->dFlux->data.F32[k] * scaleFactor; + // for chisq + // Numerator = sum( nPix[i] * chisq[i] ) + // denominator will be sum( nPix ) + + +// tmp = (psVector *) psBinaryOp(tmp, + // source->galaxyFits->chisq, "*", psScalarAlloc(source->galaxyFits->nPix, PS_TYPE_F32)); +// numerator = (psVector *) psBinaryOp(numerator, numerator, "+", tmp); + + numerator->data.F32[k] += chisq * source->galaxyFits->nPix; + totalNPix += source->galaxyFits->nPix; + + // sumInvSig2 = sum( 1 / dFlux**2 ) +// tmp = (psVector *) psBinaryOp(tmp, source->galaxyFits->dFlux, "*", source->galaxyFits->dFlux); +// tmp = (psVector *) psBinaryOp(tmp, psScalarAlloc(1.0, PS_TYPE_F32), "/", tmp); +// sumInvSig2 = (psVector *) psBinaryOp(sumInvSig2, sumInvSig2, "+", tmp); +// + psF32 invSig2 = 1.0 / (dFlux * dFlux); + sumInvSig2->data.F32[k] += invSig2; + + // sumWeightedFlux = sum ( Flux / dFlux**2 ) +// tmp = (psVector *) psBinaryOp(tmp, source->galaxyFits->Flux, "*", tmp); +// sumWeightedFlux = (psVector *) psBinaryOp(sumWeightedFlux, sumWeightedFlux, "+", tmp); + sumWeightedFlux->data.F32[k] += flux * invSig2; + } + } + } + + if (outSrc) { + if (vectorLength) { + outSrc->galaxyFits = pmSourceGalaxyFitsAlloc(); + outSrc->galaxyFits->nPix = totalNPix; + outSrc->galaxyFits->Flux = psVectorRecycle(outSrc->galaxyFits->Flux, vectorLength, PS_TYPE_F32); + outSrc->galaxyFits->dFlux = psVectorRecycle(outSrc->galaxyFits->dFlux, vectorLength, PS_TYPE_F32); + outSrc->galaxyFits->chisq = psVectorRecycle(outSrc->galaxyFits->chisq, vectorLength, PS_TYPE_F32); + for (int k = 0; k < vectorLength; k++) { + outSrc->galaxyFits->Flux->data.F32[k] = sumWeightedFlux->data.F32[k] / sumInvSig2->data.F32[k]; + outSrc->galaxyFits->dFlux->data.F32[k] = 1.0 / sumInvSig2->data.F32[k]; + outSrc->galaxyFits->chisq->data.F32[k] = numerator->data.F32[k] / totalNPix; + +#ifdef nodef + outSrc->galaxyFits->Flux = (psVector *) psBinaryOp(outSrc->galaxyFits->Flux, + sumWeightedFlux, "/", sumInvSig2); + + outSrc->galaxyFits->dFlux = (psVector *) psBinaryOp(outSrc->galaxyFits->dFlux, + psScalarAlloc(1.0, PS_TYPE_F32), "/", sumInvSig2); + + outSrc->galaxyFits->chisq = (psVector *) psBinaryOp(outSrc->galaxyFits->chisq, + numerator, "/", psScalarAlloc(totalNPix, PS_TYPE_F32)); +#endif + } + + psFree(numerator); + psFree(sumInvSig2); + psFree(sumWeightedFlux); +// psFree(tmp); + + int min_j = -1; + psF32 minChisq = NAN; + psVector *chisq = outSrc->galaxyFits->chisq; + for (int j=0; j < chisq->n; j++) { + psF32 thischisq = chisq->data.F32[j]; + if (isfinite(thischisq) && (!isfinite(minChisq) || thischisq < minChisq)) { + min_j = j; + minChisq = thischisq; + } + } + if (min_j >= 0 && isfinite(minChisq) && outSrc->modelEXT) { + // copy the best fit params to the model + psEllipseAxes axes = pmPSF_ModelToAxes(outSrc->modelEXT->params->data.F32, outSrc->modelEXT->type); + axes.major *= options->rMajor->data.F32[min_j]; + axes.minor *= options->rMinor->data.F32[min_j]; + pmPSF_AxesToModel (outSrc->modelEXT->params->data.F32, axes, outSrc->modelEXT->type); + outSrc->modelEXT->chisq = minChisq; + outSrc->modelEXT->mag = -2.5 * log10(outSrc->galaxyFits->Flux->data.F32[min_j]); + outSrc->modelEXT->magErr = 1.0 / + (outSrc->galaxyFits->Flux->data.F32[min_j]/outSrc->galaxyFits->dFlux->data.F32[min_j]); // 1 / SN + } + } + } + + return outSrc; +} + +#define GETVAL(member, key) \ + opt->member = psMetadataLookupF32(&status, md, key); \ + if (!status) { \ + psError (PSPHOT_ERR_CONFIG, true, "failed to looup value for %s in %s", key, \ + useAnalysis ? "readout->analysis" : "recipe"); \ + } + +static bool setOptions(galaxyShapeOptions *opt, pmReadout *readout, psMetadata *recipe, bool makeVectors) { + bool status; + bool useAnalysis; + + psMetadataLookupF32(&useAnalysis, readout->analysis, "GALAXY_SHAPES_FR_MAJOR_MIN"); + psMetadata *md = useAnalysis ? readout->analysis : recipe; + + GETVAL(fRmajorMin, "GALAXY_SHAPES_FR_MAJOR_MIN"); + GETVAL(fRmajorMax, "GALAXY_SHAPES_FR_MAJOR_MAX"); + GETVAL(fRmajorDel, "GALAXY_SHAPES_FR_MAJOR_DEL"); + GETVAL(fRminorMin, "GALAXY_SHAPES_FR_MINOR_MIN"); + GETVAL(fRminorMax, "GALAXY_SHAPES_FR_MINOR_MAX"); + GETVAL(fRminorDel, "GALAXY_SHAPES_FR_MINOR_DEL"); + + opt->numTrials = ceil((opt->fRmajorMax - opt->fRmajorMin + 0.5*opt->fRmajorDel) / opt->fRmajorDel) * + ceil((opt->fRminorMax - opt->fRminorMin + 0.5*opt->fRminorDel) / opt->fRminorDel) ; + + if (makeVectors) { + opt->rMinor = psVectorAlloc(opt->numTrials, PS_TYPE_F32); + opt->rMajor = psVectorAlloc(opt->numTrials, PS_TYPE_F32); + int i = 0; + for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel; + fRmajor += opt->fRmajorDel) { + for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel; + fRminor += opt->fRminorDel) { + opt->rMinor->data.F32[i] = fRminor; + opt->rMajor->data.F32[i] = fRmajor; + i++; + } + } + psAssert(i == opt->numTrials, "Something's wrong with my loop got %d entries expected %d", i, opt->numTrials); + } else { + opt->rMinor = NULL; + opt->rMajor = NULL; + } + + return true; +} + +#define CHECKVAL(left, right, val, message) \ + if (left->val != right.val) { \ + psError (PSPHOT_ERR_CONFIG, true, message); \ + return false; \ + } + +static bool checkOptions(galaxyShapeOptions *options, pmReadout *readout, psMetadata *recipe) { + galaxyShapeOptions thisReadoutsOptions; + + if (!setOptions(&thisReadoutsOptions, readout, recipe, false)) { + psError (PS_ERR_UNKNOWN, false, "problem determining galaxy shape options for readout"); + return false; + } + CHECKVAL(options, thisReadoutsOptions, numTrials, "mismatched number of trials") + CHECKVAL(options, thisReadoutsOptions, fRmajorMin, "mismatched fRmajorMin") + CHECKVAL(options, thisReadoutsOptions, fRmajorMax, "mismatched fRmajorMax") + CHECKVAL(options, thisReadoutsOptions, fRmajorDel, "mismatched fRmajorDel") + CHECKVAL(options, thisReadoutsOptions, fRminorMin, "mismatched fRminorMin") + CHECKVAL(options, thisReadoutsOptions, fRminorMax, "mismatched fRminorMax") + CHECKVAL(options, thisReadoutsOptions, fRminorDel, "mismatched fRminorDel") + + return true; +} Index: trunk/psphot/src/psphotGalaxyShape.c =================================================================== --- trunk/psphot/src/psphotGalaxyShape.c (revision 36375) +++ trunk/psphot/src/psphotGalaxyShape.c (revision 36441) @@ -283,5 +283,13 @@ } +#ifdef SAVE_BEST_MODEL + // Save model with smallest chisq if (isfinite(chisqBest)) { +#else + // Save model with nominal parameters + { + fRmajorBest = 1; + fRminorBest = 1; +#endif // now save the best fitting model as the source's extended model psEllipseAxes testAxes = guessAxes; @@ -293,4 +301,6 @@ psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, false); + // Replace modelEXT with this model + // XXX: only do this if the model is good psFree (source->modelEXT); @@ -299,7 +309,4 @@ source->mode |= PM_SOURCE_MODE_EXTMODEL; source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; - - // adjust the window so the subtraction covers the faint wings - // psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal); // cache the model flux