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;