Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/Makefile.am =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/Makefile.am (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/Makefile.am (revision 37403) @@ -11,5 +11,6 @@ pmAstrometryRefstars.c \ pmAstrometryWCS.c \ - pmAstrometryVisual.c + pmAstrometryVisual.c \ + pmKHcorrect.c pkginclude_HEADERS = \ @@ -21,5 +22,6 @@ pmAstrometryRefstars.h \ pmAstrometryWCS.h \ - pmAstrometryVisual.h + pmAstrometryVisual.h \ + pmKHcorrect.h CLEANFILES = *~ Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmAstrometryModel.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmAstrometryModel.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmAstrometryModel.c (revision 37403) @@ -39,4 +39,5 @@ #include "pmFPAExtent.h" #include "pmFPAfileFitsIO.h" +#include "pmConcepts.h" #include "pmAstrometryWCS.h" #include "pmAstrometryUtils.h" @@ -452,28 +453,4 @@ } -int pmConceptsChipNumberFromName (pmFPA *fpa, char *name) { - - for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - if (!chip) continue; - char *thisone = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); - if (!thisone) continue; - if (!strcmp (name, thisone)) return (i); - } - return -1; -} - -pmChip *pmConceptsChipFromName (pmFPA *fpa, char *name) { - - for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - if (!chip) continue; - char *thisone = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); - if (!thisone) continue; - if (!strcmp (name, thisone)) return (chip); - } - return NULL; -} - // first layer converts Chip to Focal Plane bool pmAstromModelReadChips (pmFPAfile *file) { Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmAstrometryModel.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmAstrometryModel.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmAstrometryModel.h (revision 37403) @@ -27,7 +27,4 @@ bool pmAstromModelWriteChips (pmFPAfile *file); -int pmConceptsChipNumberFromName (pmFPA *fpa, char *name); -pmChip *pmConceptsChipFromName (pmFPA *fpa, char *name); - bool pmAstromModelReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config); bool pmAstromModelReadFPA (pmFPAfile *file); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmAstrometryObjects.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmAstrometryObjects.c (revision 37403) @@ -571,13 +571,21 @@ psMemSetDeallocator (stats, (psFreeFunc)pmAstromStatsFree); - // stats->center = {0, 0, 0, 0}; - // stats->offset = {0, 0, 0, 0}; - stats->angle = 0.0; - stats->scale = 1.0; - stats->minMetric = 0.0; - stats->minVar = 0.0; - stats->nMatch = 0; - stats->nTest = 0; - stats->nSigma = 0; + stats->center.x = 0; + stats->center.y = 0; + stats->center.xErr = 0; + stats->center.yErr = 0; + + stats->offset.x = 0; + stats->offset.y = 0; + stats->offset.xErr = 0; + stats->offset.yErr = 0; + + stats->angle = 0.0; + stats->scale = 1.0; + stats->minMetric = 0.0; + stats->minVar = 0.0; + stats->nMatch = 0; + stats->nTest = 0; + stats->nSigma = 0; return (stats); @@ -914,6 +922,5 @@ // fprintf (stderr, "sigma: nMatch: %d, nTest: %d, nTen: %d\n", stats->nMatch, stats->nTest, sort->data.U32[sort->n - 10]); - - psFree (sort); + psFree (sort); psFree (listNP); psFree (gridNP); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmAstrometryObjects.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmAstrometryObjects.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmAstrometryObjects.h (revision 37403) @@ -42,4 +42,5 @@ float Color; ///< object color float dMag; ///< error on object magnitude + float SBinst; ///< surface brightness, used for Koppenhoefer correction } pmAstromObj; Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmAstrometryVisual.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmAstrometryVisual.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmAstrometryVisual.c (revision 37403) @@ -946,5 +946,5 @@ KapaSendLabel (kapa2, "X (FP)", KAPA_LABEL_XM); KapaSendLabel (kapa2, "Y (FP)", KAPA_LABEL_YM); - KapaSendLabel (kapa2, "pmAstromGridAngle residuals. Box: Correlation Peak.", KAPA_LABEL_XP); + KapaSendLabel (kapa2, "pmAstromGridAngle red: raw, black: ref.", KAPA_LABEL_XP); // plot the REF data. (also calculate the plot ranges, accumulate the plot vectors) @@ -966,13 +966,18 @@ KapaSetLimits(kapa2, &graphdata); - psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_MIN | PS_STAT_MAX ); psVectorStats (stats, zPlot, NULL, NULL, 0); - float zero = stats->sampleMedian + 3.0; - float range = 6.0; - + float range = stats->max - stats->min; + range = PS_MAX (0.5, PS_MIN (6.0, range)); + float zero = stats->sampleMedian + 0.25*range; + + float maxZ = zPlot->data.F32[0], minZ = zPlot->data.F32[0]; for (int i = 0; i < zPlot->n; i++) { + maxZ = PS_MAX (maxZ, zPlot->data.F32[i]); + minZ = PS_MIN (minZ, zPlot->data.F32[i]); float value = (zero - zPlot->data.F32[i]) / range; zPlot->data.F32[i] = PS_MAX(0.0, PS_MIN(1.0, value)); } + fprintf (stderr, "ref mags: %f to %f (%f median)\n", minZ, maxZ, stats->sampleMedian); // the point size will be scaled from the z vector @@ -998,11 +1003,18 @@ psStatsInit(stats); psVectorStats (stats, zPlot, NULL, NULL, 0); - zero = stats->sampleMedian + 3.0; - range = 6.0; - + range = stats->max - stats->min; + range = PS_MAX (0.5, PS_MIN (6.0, range)); + zero = stats->sampleMedian + 0.25*range; + // zero = stats->sampleMedian + 1.0; + // range = 6.0; + + maxZ = zPlot->data.F32[0], minZ = zPlot->data.F32[0]; for (int i = 0; i < zPlot->n; i++) { + maxZ = PS_MAX (maxZ, zPlot->data.F32[i]); + minZ = PS_MIN (minZ, zPlot->data.F32[i]); float value = (zero - zPlot->data.F32[i]) / range; zPlot->data.F32[i] = PS_MAX(0.0, PS_MIN(1.0, value)); } + fprintf (stderr, "raw mags: %f to %f (%f median)\n", minZ, maxZ, stats->sampleMedian); // the point size will be scaled from the z vector Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmKHcorrect.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmKHcorrect.c (revision 37403) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmKHcorrect.c (revision 37403) @@ -0,0 +1,252 @@ +/** @file pmKHcorrect.c + * @brief Functions to read (and write?) Koppenhoefer correction file + * + * The Koppenhoefer correction is needed for some chips of gpc1 before the camera voltages were adjusted 2011/05/11. + * The correction is a modification of the X (and possibly Y) coordinate of a star which depends on the instrumental + * surface brightness, defined as -2.5 log_10 (DN) + 5.0 log_10 fwhm_maj [XXX be careful about the definition of fwhm_maj] + * + * @ingroup AstroImage + * @author EAM, IfA + * + * Copyright 2014 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +/******************************************************************************/ +/* INCLUDE FILES */ +/******************************************************************************/ +#include +#include +#include +#include +#include +#include // for unlink +#include + +#include "pmConfig.h" +#include "pmDetrendDB.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" +#include "pmFPAExtent.h" +#include "pmFPAfileFitsIO.h" +#include "pmConcepts.h" +#include "pmKHcorrect.h" + +static void KHcorrectDataFree (KHcorrectData *spline) { + + if (!spline) return; + + psFree (spline->xk); + psFree (spline->yk); + psFree (spline->y2); + return; +} + +KHcorrectData *KHcorrectDataAlloc (int Nrow) { + + // allocate xk[Nrow], etc + + KHcorrectData *spline = (KHcorrectData *) psAlloc(sizeof(KHcorrectData)); + psMemSetDeallocator(spline, (psFreeFunc) KHcorrectDataFree); + + spline->N = Nrow; + spline->xk = (float *) psAlloc(Nrow*sizeof(float)); + spline->yk = (float *) psAlloc(Nrow*sizeof(float)); + spline->y2 = (float *) psAlloc(Nrow*sizeof(float)); + + return spline; +} + +// the KH correction is a function of the instrumental surface brightness, SBinst +float KHcorrectApply (KHcorrectData *spline, float X) { + + int N = spline->N; + + // saturate correction at high and low ends + if (X < spline->xk[ 0]) return spline->yk[ 0]; + if (X > spline->xk[N-1]) return spline->yk[N-1]; + + float *xk = spline->xk; + float *yk = spline->yk; + float *y2 = spline->y2; + + /* find correct element in array (x must be sorted) */ + int lo = 0; + int hi = N-1; + while (hi - lo > 1) { + int i = 0.5*(hi+lo); + if (xk[i] > X) { + hi = i; + } else { + lo = i; + } + } + + /* error condition: duplicate abssisca */ + float dx = xk[hi] - xk[lo]; + if (dx == 0.0) { + return (0.0); + } + + /* evaluate spline */ + float a = (xk[hi] - X) / dx; + float b = (X - xk[lo]) / dx; + + float value = a*yk[lo] + b*yk[hi] + ((a*a - 1.0)*a*y2[lo] + (b*b - 1.0)*b*y2[hi])*(dx*dx) / 6.0; + return (value); +} + +/********************* CheckDataStatus functions *****************************/ + +bool pmKHcorrectCheckDataStatusForView (const pmFPAview *view, pmFPAfile *file) { + psError(PS_ERR_IO, false, "Check Data Status not defined"); + return false; +} + +bool pmKHcorrectCheckDataStatusForFPA (const pmFPA *fpa) { + psError(PS_ERR_IO, false, "Check Data Status not defined"); + return false; +} + +bool pmKHcorrectCheckDataStatusForChip (const pmChip *chip) { + psError(PS_ERR_IO, false, "Check Data Status not defined"); + return false; +} + +/********************* Write Data functions *****************************/ + +// NOTE : these are not exposed because I don't think we need them (we do not create KHcorrect +// in psModules based tool, but in DVO. + +bool pmKHcorrectWriteFPA (pmFPAfile *file, const pmFPA *fpa); +bool pmKHcorrectWriteForView (const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + // write the full model in one pass: require the level to be FPA + if (view->chip != -1) { + psError(PS_ERR_IO, false, "Koppenhoefer Correction must be written at the FPA level"); + return false; + } + + pmFPA *fpa = pmFPAfileSuitableFPA(file, view, config, false); // Suitable FPA for writing + + if (!pmKHcorrectWriteFPA(file, fpa)) { + psError(PS_ERR_IO, false, "Failed to write KH Correction for fpa"); + psFree(fpa); + return false; + } + + psFree(fpa); + + return true; +} + +// write out all chip-level KH Correction data for this FPA +bool pmKHcorrectWriteFPA (pmFPAfile *file, const pmFPA *fpa) +{ + psError(PS_ERR_IO, false, "output for KH Correction is not defined"); + return false; +} + +/********************* Read Data functions *****************************/ + +bool pmKHcorrectReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) + { + // read the full model in one pass: require the level to be FPA + if (view->chip != -1) { + psError(PS_ERR_IO, false, "KH Correction must be read at the FPA level"); + return false; + } + + if (!pmKHcorrectReadFPA (file)) { + psError(PS_ERR_IO, false, "Failed to read KH Correction for fpa"); + return false; + } + return true; + } + +// read in all chip-level KH Correction data for this FPA +bool pmKHcorrectReadFPA (pmFPAfile *file) { + + if (!pmKHcorrectReadChips (file)) { + psError(PS_ERR_IO, false, "Failed to read KH Correction for chips"); + return false; + } + + return true; +} + +// read the set of tables, one for each chip +bool pmKHcorrectReadChips (pmFPAfile *file) { + + bool haveData, status; + + // loop over the extensions + // for each extension, use the extname (eg, XY01.DX.T0) to assign to a chip + + // move to the start of the file + haveData = psFitsMoveExtNum (file->fits, 1, false); + if (!haveData) { + psError(PS_ERR_IO, false, "Failed to read even the first extension?"); + return false; + } + + while (haveData) { + + // load the header + psMetadata *header = psFitsReadHeader(NULL, file->fits); // The FITS header + if (!header) psAbort("cannot read model header"); + + // load the full model in one shot + psArray *model = psFitsReadTable (file->fits); + if (!model) psAbort("cannot read model"); + + // determine the chip: + char *extname = psMetadataLookupStr (&status, header, "EXTNAME"); + psLogMsg ("psModules.astrom", 4, "read %ld rows from Koppenhoefer correction file, extname %s\n", model->n, extname); + + // I expect to find a name of the form: chipName.dir.tset (eg, XY01.DX.T0) + // where chipName like 'XY01' + // dir = 'DX' + // tset = 'T0' + psAssert (strlen(extname) == 10, "invalid extension %s", extname); + psAssert (extname[5] == 'D', "invalid extension %s", extname); + psAssert (extname[6] == 'X', "invalid extension %s", extname); + psAssert (extname[8] == 'T', "invalid extension %s", extname); + psAssert (extname[9] == '0', "invalid extension %s", extname); + + char chipName[5]; + strncpy (chipName, extname, 4); + chipName[4] = 0; + + pmChip *chip = pmConceptsChipFromName (file->fpa, chipName); + if (!chip) psAbort ("invalid chip?"); + + KHcorrectData *spline = KHcorrectDataAlloc (model->n); + + // parse the model entries + for (int i = 0; i < model->n; i++) { + psMetadata *row = model->data[i]; + + spline->xk[i] = psMetadataLookupF32(&status, row, "X_KNOT"); + spline->yk[i] = psMetadataLookupF32(&status, row, "Y_KNOT"); + spline->y2[i] = psMetadataLookupF32(&status, row, "DY2_DX"); + } + psMetadataAddUnknown (chip->analysis, PS_LIST_TAIL, "KH.CORRECT", PS_META_REPLACE, "", spline); + psFree (spline); + + psFree (model); + psFree (header); + + // move to the next extension + haveData = psFitsMoveExtNum (file->fits, 1, true); + } + + return true; +} + Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmKHcorrect.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmKHcorrect.h (revision 37403) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/astrom/pmKHcorrect.h (revision 37403) @@ -0,0 +1,35 @@ +/** @file pmKHcorrect.h + * @brief Functions to read (and write?) Koppenhoefer correction file + * + * The Koppenhoefer correction is needed for some chips of gpc1 before the camera voltages were adjusted 2011/05/11. + * The correction is a modification of the X (and possibly Y) coordinate of a star which depends on the instrumental + * surface brightness, defined as -2.5 log_10 (DN) + 5.0 log_10 fwhm_maj [XXX be careful about the definition of fwhm_maj] + * + * @ingroup AstroImage + * @author EAM, IfA + * + * Copyright 2014 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_KH_CORRECT_H +#define PM_KH_CORRECT_H + +/// @addtogroup Astrometry +/// @{ + +typedef struct { + int N; + float *xk; + float *yk; + float *y2; +} KHcorrectData; + +KHcorrectData *KHcorrectDataAlloc (int Nrow); +float KHcorrectApply (KHcorrectData *spline, float X); + +bool pmKHcorrectReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config); +bool pmKHcorrectReadFPA (pmFPAfile *file); +bool pmKHcorrectReadChips (pmFPAfile *file); + +/// @} +#endif Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/camera/pmFPAfile.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/camera/pmFPAfile.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/camera/pmFPAfile.c (revision 37403) @@ -531,7 +531,6 @@ return PM_FPA_FILE_LINEARITY; } - // deprecate this? if (!strcasecmp(type, "ASTROM")) { - return PM_FPA_FILE_ASTROM_MODEL; + return PM_FPA_FILE_ASTROM_MODEL; } if (!strcasecmp(type, "ASTROM.MODEL")) { @@ -540,4 +539,7 @@ if (!strcasecmp(type, "ASTROM.REFSTARS")) { return PM_FPA_FILE_ASTROM_REFSTARS; + } + if (!strcasecmp(type, "KH.CORRECT")) { + return PM_FPA_FILE_KH_CORRECT; } if (!strcasecmp(type, "SUBKERNEL")) { @@ -593,4 +595,6 @@ case PM_FPA_FILE_ASTROM_REFSTARS: return ("ASTROM.REFSTARS"); + case PM_FPA_FILE_KH_CORRECT: + return ("KH.CORRECT"); case PM_FPA_FILE_SUBKERNEL: return ("SUBKERNEL"); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/camera/pmFPAfile.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/camera/pmFPAfile.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/camera/pmFPAfile.h (revision 37403) @@ -48,4 +48,5 @@ PM_FPA_FILE_ASTROM_MODEL, PM_FPA_FILE_ASTROM_REFSTARS, + PM_FPA_FILE_KH_CORRECT, PM_FPA_FILE_SUBKERNEL, PM_FPA_FILE_SRCTEXT, Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/camera/pmFPAfileFitsIO.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/camera/pmFPAfileFitsIO.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/camera/pmFPAfileFitsIO.c (revision 37403) @@ -150,5 +150,7 @@ case PM_FPA_FILE_PSF: case PM_FPA_FILE_ASTROM_MODEL: - case PM_FPA_FILE_ASTROM_REFSTARS: { + case PM_FPA_FILE_ASTROM_REFSTARS: + case PM_FPA_FILE_KH_CORRECT: + { pmHDU *hdu = pmFPAviewThisHDU(view, fpa); if (hdu) { Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/camera/pmFPAfileIO.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/camera/pmFPAfileIO.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/camera/pmFPAfileIO.c (revision 37403) @@ -34,4 +34,5 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmSourceMasks.h" @@ -48,4 +49,5 @@ #include "pmPSF_IO.h" +#include "pmKHcorrect.h" #include "pmAstrometryModel.h" #include "pmAstrometryRefstars.h" @@ -234,4 +236,7 @@ status = pmAstromModelReadForView (view, file, config); break; + case PM_FPA_FILE_KH_CORRECT: + status = pmKHcorrectReadForView (view, file, config); + break; case PM_FPA_FILE_EXPNUM: status = pmExpNumRead(view, file, config); @@ -324,4 +329,5 @@ case PM_FPA_FILE_ASTROM_MODEL: case PM_FPA_FILE_ASTROM_REFSTARS: + case PM_FPA_FILE_KH_CORRECT: case PM_FPA_FILE_JPEG: case PM_FPA_FILE_KAPA: @@ -408,4 +414,8 @@ } } + if (file->type == PM_FPA_FILE_KH_CORRECT) { + psTrace("psModules.camera", 6, "skip write for %s, no write function defined", file->name); + return true; + } // open the file if not yet opened @@ -500,4 +510,8 @@ case PM_FPA_FILE_ASTROM_REFSTARS: status = pmAstromRefstarsWriteForView (view, file, config); + break; + + case PM_FPA_FILE_KH_CORRECT: + psError(PS_ERR_IO, true, "cannot write type KH.CORRECT (%s)", file->name); break; @@ -576,4 +590,5 @@ case PM_FPA_FILE_ASTROM_MODEL: case PM_FPA_FILE_ASTROM_REFSTARS: + case PM_FPA_FILE_KH_CORRECT: case PM_FPA_FILE_LINEARITY: case PM_FPA_FILE_EXPNUM: @@ -653,4 +668,5 @@ case PM_FPA_FILE_ASTROM_MODEL: case PM_FPA_FILE_ASTROM_REFSTARS: + case PM_FPA_FILE_KH_CORRECT: case PM_FPA_FILE_EXPNUM: psTrace ("psModules.camera", 6, "NOT freeing %s (%s) : save for further analysis\n", file->filename, file->name); @@ -815,4 +831,5 @@ case PM_FPA_FILE_ASTROM_MODEL: case PM_FPA_FILE_ASTROM_REFSTARS: + case PM_FPA_FILE_KH_CORRECT: case PM_FPA_FILE_LINEARITY: case PM_FPA_FILE_EXPNUM: @@ -1018,4 +1035,5 @@ case PM_FPA_FILE_EXPNUM: case PM_FPA_FILE_ASTROM_MODEL: + case PM_FPA_FILE_KH_CORRECT: case PM_FPA_FILE_SX: case PM_FPA_FILE_RAW: Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/camera/pmReadoutFake.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/camera/pmReadoutFake.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/camera/pmReadoutFake.c (revision 37403) @@ -20,7 +20,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -55,5 +55,5 @@ psF32 *params = model->params->data.F32; // Model parameters - psEllipseAxes axes = pmPSF_ModelToAxes(params, model->type); // Ellipse axes + psEllipseAxes axes = pmPSF_ModelToAxes(params, model->class->useReff); // Ellipse axes // Curiously, the minor axis can be larger than the major axis, so need to check. if (axes.major >= axes.minor) { @@ -62,5 +62,5 @@ axes.major = axes.minor; } - return pmPSF_AxesToModel(params, axes, model->type); + return pmPSF_AxesToModel(params, axes, model->class->useReff); } @@ -122,5 +122,5 @@ } - flux /= normModel->modelFlux(normModel->params); + flux /= normModel->class->modelFlux(normModel->params); psFree(normModel); } @@ -164,5 +164,5 @@ float fakeRadius = 1.0; // Radius of fake source if (isfinite(minFlux)) { - fakeRadius = PS_MAX(fakeRadius, fakeModel->modelRadius(fakeModel->params, minFlux)); + fakeRadius = PS_MAX(fakeRadius, fakeModel->class->modelRadius(fakeModel->params, minFlux)); } if (radius > 0) { Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/concepts/pmConcepts.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/concepts/pmConcepts.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/concepts/pmConcepts.c (revision 37403) @@ -312,4 +312,7 @@ conceptRegisterEnum("FPA.TIMESYS", "Time system", p_pmConceptParse_TIMESYS, p_pmConceptFormat_TIMESYS, p_pmConceptCopy_TIMESYS, false, PM_FPA_LEVEL_FPA); conceptRegisterTime("FPA.TIME", "Time of exposure", false, PM_FPA_LEVEL_FPA); + + conceptRegisterTime("FPA.SHUTOUTC", "Time of exposure", false, PM_FPA_LEVEL_FPA); + conceptRegisterF32("FPA.TEMP", "Temperature of focal plane", NULL, NULL, NULL, false, PM_FPA_LEVEL_FPA); conceptRegisterF32("FPA.M1X", "Primary Mirror X Position", NULL, NULL, NULL, false, PM_FPA_LEVEL_FPA); @@ -578,2 +581,26 @@ } +int pmConceptsChipNumberFromName (pmFPA *fpa, char *name) { + + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + if (!chip) continue; + char *thisone = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); + if (!thisone) continue; + if (!strcmp (name, thisone)) return (i); + } + return -1; +} + +pmChip *pmConceptsChipFromName (pmFPA *fpa, char *name) { + + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + if (!chip) continue; + char *thisone = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); + if (!thisone) continue; + if (!strcmp (name, thisone)) return (chip); + } + return NULL; +} + Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/concepts/pmConcepts.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/concepts/pmConcepts.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/concepts/pmConcepts.h (revision 37403) @@ -161,4 +161,8 @@ ); +// some utility functions: +int pmConceptsChipNumberFromName (pmFPA *fpa, char *name); +pmChip *pmConceptsChipFromName (pmFPA *fpa, char *name); + /// @} #endif Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/config/pmConfig.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/config/pmConfig.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/config/pmConfig.c (revision 37403) @@ -965,4 +965,42 @@ "Original replaced by -R option", newRule); } + psFree(camerasIter); + } + + // Look for command-line options for files to replace + while ((argNum = psArgumentGet(*argc, argv, "-photcode-rule")) > 0) { + psArgumentRemove(argNum, argc, argv); + if (argNum >= *argc) { + psError(PM_ERR_CONFIG, true, + "-photcode-rule provided without new rule."); + psFree(config); + return NULL; + } + + psString newrule = psStringCopy(argv[argNum]); // The filerule, to be modified + psArgumentRemove(argNum, argc, argv); + + psMetadata *cameras = psMetadataLookupMetadata(NULL, config->system, "CAMERAS"); // List of cameras + if (!cameras) { + psError(PM_ERR_CONFIG, false, "Unable to find CAMERAS in the site configuration.\n"); + return false; + } + + psMetadataIterator *camerasIter = psMetadataIteratorAlloc(cameras, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *cameraItem; // Item from iteration + while ((cameraItem = psMetadataGetAndIncrement(camerasIter))) { + // Silently ignore problems --- they will be caught later, because if the user wants the nominated + // file and it's not available for that camera, then they will know. + + if (cameraItem->type != PS_DATA_METADATA) { + psTrace("psModules.config", 2, + "Entry %s in CAMERAS is not of type METADATA --- ignored.", cameraItem->name); + continue; + } + psMetadata *camera = cameraItem->data.md; // Camera configuration + + psMetadataAddStr (camera, PS_LIST_TAIL, "PHOTCODE.RULE", PS_META_REPLACE, "original replaced by -photcode-rule option", newrule); + } + psFree(newrule); psFree(camerasIter); } Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/detrend/pmDetrendDB.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/detrend/pmDetrendDB.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/detrend/pmDetrendDB.c (revision 37403) @@ -111,4 +111,5 @@ DETREND_STRING_CASE(LINEARITY); DETREND_STRING_CASE(AUXMASK); + DETREND_STRING_CASE(KH_CORRECT); default: return NULL; Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/detrend/pmDetrendDB.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/detrend/pmDetrendDB.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/detrend/pmDetrendDB.h (revision 37403) @@ -40,4 +40,5 @@ PM_DETREND_TYPE_LINEARITY, PM_DETREND_TYPE_AUXMASK, + PM_DETREND_TYPE_KH_CORRECT, } pmDetrendType; Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/extras/pmThreadTools.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/extras/pmThreadTools.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/extras/pmThreadTools.c (revision 37403) @@ -29,7 +29,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/extras/pmVisual.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/extras/pmVisual.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/extras/pmVisual.c (revision 37403) @@ -22,33 +22,5 @@ bool pmSourceVisualClose(void); -// #include "pmHDU.h" -// #include "pmFPA.h" -// #include "pmFPAfile.h" -// #include "pmAstrometryObjects.h" -// #include "pmSubtractionStamps.h" -// #include "pmTrend2D.h" -// #include "pmResiduals.h" -// #include "pmGrowthCurve.h" -// #include "pmSpan.h" -// #include "pmFootprintSpans.h" -// #include "pmFootprint.h" -// #include "pmPeaks.h" -// #include "pmMoments.h" -// #include "pmModelFuncs.h" -// #include "pmModel.h" -// #include "pmSourceMasks.h" -// #include "pmSourceExtendedPars.h" -// #include "pmSourceDiffStats.h" #include "pmSourceSatstar.h" -// #include "pmSourceLensing.h" -// #include "pmSource.h" -// #include "pmSourceFitModel.h" -// #include "pmPSF.h" -// #include "pmPSFtry.h" -// #include "pmFPAExtent.h" -// #include "pmAstrometryVisual.h" -// #include "pmSubtractionVisual.h" -// #include "pmStackVisual.h" -// #include "pmSourceVisual.h" # if (HAVE_KAPA) @@ -161,4 +133,77 @@ } + +// ask the user to continue or not. give up after 2 seconds. +bool pmVisualAskUserOrDump(bool *plotFlag, bool *dumpData) +{ + struct timeval timeout; + fd_set fdSet; + int status; + + if (dumpData) *dumpData = false; + + char key[10]; + if (plotFlag && dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [s]kip the rest of these plots? [d]ump the data? [a]bort all visual plots? (c) "); + } + if (plotFlag && !dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots? (c) "); + } + if (!plotFlag && dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [d]ump the data? [a]bort all visual plots? (c) "); + } + if (!plotFlag && !dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [a]bort all visual plots? (c) "); + } + + /* Wait up to 1.0 second for a response, then continue */ + timeout.tv_sec = 10; + timeout.tv_usec = 0; + + FD_ZERO (&fdSet); + FD_SET (STDIN_FILENO, &fdSet); + + status = select (1, &fdSet, NULL, NULL, &timeout); + if (status <= 0) { + fprintf (stderr, "\n"); + return true; // if no data, give up + } + + while (true) { + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + switch (key[0]) { + case 's': + if (plotFlag) *plotFlag = false; + return true; + case 'd': + if (dumpData) *dumpData = true; + return true; + case 'a': + isVisual = false; + return true; + case 'c': + case '\n': + return true; + default: + break; + } + + if (plotFlag && dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [s]kip the rest of these plots? [d]ump the data? [a]bort all visual plots? (c) "); + } + if (plotFlag && !dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots? (c) "); + } + if (!plotFlag && dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [d]ump the data? [a]bort all visual plots? (c) "); + } + if (!plotFlag && !dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [a]bort all visual plots? (c) "); + } + } + return true; +} bool pmVisualImStats(psImage *image, double *mean, double *stdev, double *min, double *max) { Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/extras/pmVisual.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/extras/pmVisual.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/extras/pmVisual.h (revision 37403) @@ -52,4 +52,12 @@ */ bool pmVisualAskUser(bool *plotFlag); + + +/** Ask the user how to proceed. + * At the user's request, this will disable diagnostic plotting. + * @param plotFlag, set to false if this plot should be disabled in the future + * @param dumpData, set to true if user requests a data dump + */ +bool pmVisualAskUserOrDump(bool *plotFlag, bool *dumpData); @@ -138,4 +146,5 @@ bool pmVisualInitGraph (int kapa, void *section, void *graphdata); bool pmVisualAskUser(bool *plotFlag); +bool pmVisualAskUserOrDump(bool *plotFlag, bool *dumpData); bool pmVisualScaleImage(int kapaFD, psImage *inImage, const char *name, int channel, bool clip); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/imcombine/pmPSFEnvelope.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/imcombine/pmPSFEnvelope.c (revision 37403) @@ -20,7 +20,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -168,9 +168,9 @@ continue; } - model->modelSetLimits(PM_MODEL_LIMITS_MODERATE); + model->class->modelSetLimits(PM_MODEL_LIMITS_MODERATE); bool limits = true; // Model within limits? for (int j = 0; j < model->params->n && limits; j++) { - if (!model->modelLimits(PS_MINIMIZE_PARAM_MIN, j, model->params->data.F32, NULL) || - !model->modelLimits(PS_MINIMIZE_PARAM_MAX, j, model->params->data.F32, NULL)) { + if (!model->class->modelLimits(PS_MINIMIZE_PARAM_MIN, j, model->params->data.F32, NULL) || + !model->class->modelLimits(PS_MINIMIZE_PARAM_MAX, j, model->params->data.F32, NULL)) { limits = false; } @@ -246,5 +246,5 @@ continue; } - float srcRadius = model->modelRadius(model->params, PS_SQR(VARIANCE_VAL)); // Radius for source + float srcRadius = model->class->modelRadius(model->params, PS_SQR(VARIANCE_VAL)); // Radius for source psFree(model); if (srcRadius == 0) { Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/imcombine/pmStack.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/imcombine/pmStack.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/imcombine/pmStack.c (revision 37403) @@ -835,4 +835,5 @@ psImageMaskType goodMask, // Value for good pixels bool safe, // Safe combination? + int nminpix, // Minimum number of input per pixel float invTotalWeight // Inverse of total weight for all inputs ) @@ -854,5 +855,11 @@ CHECKPIX(x, y, "bad vs good : %x %x %x\n", maskValue, badMask, blankMask); - switch (num) { + //MEH -- hackish adding of lower limit for input per pixel + int numN = num; + if (num < nminpix) { + CHECKPIX(x, y, "Nmin (%d) inputs (%d) to combine, pixel %d,%d is manually set bad\n", nminpix, numN, x, y); + numN = 0; + } + switch (numN) { case 0: { // Nothing to combine: it's bad @@ -1518,4 +1525,5 @@ bool useVariance, bool safe, + int nminpix, bool rejection) { @@ -1692,5 +1700,5 @@ psImageMaskType goodMask = 0; // OR of mask bits in all good input pixels combineExtract(&num, &suspect, &badMask, &goodMask, buffer, combinedImage, combinedMask, combinedVariance, input, weights, exps, addVariance, reject, x, y, badMaskBits, suspectMaskBits); - combinePixels(combinedImage, combinedMask, combinedVariance, exp, expnum, expweight, num, buffer, x, y, blankMaskBits, badMask, goodMask, safe, totalExpWeight); + combinePixels(combinedImage, combinedMask, combinedVariance, exp, expnum, expweight, num, buffer, x, y, blankMaskBits, badMask, goodMask, safe, nminpix, totalExpWeight); if (iter > 0) { Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/imcombine/pmStack.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/imcombine/pmStack.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/imcombine/pmStack.h (revision 37403) @@ -61,4 +61,5 @@ bool useVariance, ///< Use variance values for rejection? bool safe, ///< Play safe with small numbers of input pixels (mask if N <= 2)? + int nminpix, ///< Minimum number input per pixel to combine bool rejectInspect ///< Reject pixels instead of marking them for inspection? ); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/imcombine/pmSubtractionStamps.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/imcombine/pmSubtractionStamps.c (revision 37403) @@ -21,4 +21,5 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmSourceMasks.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/Makefile.am =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/Makefile.am (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/Makefile.am (revision 37403) @@ -52,7 +52,10 @@ pmSourceIO_CMF_PS1_SV1.c \ pmSourceIO_CMF_PS1_SV2.c \ + pmSourceIO_CMF_PS1_SV3.c \ + pmSourceIO_CMF_PS1_SV4.c \ pmSourceIO_CMF_PS1_DV1.c \ pmSourceIO_CMF_PS1_DV2.c \ pmSourceIO_CMF_PS1_DV3.c \ + pmSourceIO_CMF_PS1_DV4.c \ pmSourceIO_MatchedRefs.c \ pmSourcePlots.c \ @@ -152,6 +155,9 @@ pmSourceIO_CMF_PS1_DV2.c \ pmSourceIO_CMF_PS1_DV3.c \ +pmSourceIO_CMF_PS1_DV4.c \ pmSourceIO_CMF_PS1_SV1.c \ -pmSourceIO_CMF_PS1_SV2.c +pmSourceIO_CMF_PS1_SV2.c \ +pmSourceIO_CMF_PS1_SV3.c \ +pmSourceIO_CMF_PS1_SV4.c pmSourceIO_CMF_PS1_V1.c : pmSourceIO_CMF.c.in mksource.pl @@ -179,4 +185,7 @@ mksource.pl pmSourceIO_CMF.c.in PS1_DV3 pmSourceIO_CMF_PS1_DV3.c +pmSourceIO_CMF_PS1_DV4.c : pmSourceIO_CMF.c.in mksource.pl + mksource.pl pmSourceIO_CMF.c.in PS1_DV4 pmSourceIO_CMF_PS1_DV4.c + pmSourceIO_CMF_PS1_SV1.c : pmSourceIO_CMF.c.in mksource.pl mksource.pl pmSourceIO_CMF.c.in PS1_SV1 pmSourceIO_CMF_PS1_SV1.c @@ -185,3 +194,9 @@ mksource.pl pmSourceIO_CMF.c.in PS1_SV2 pmSourceIO_CMF_PS1_SV2.c +pmSourceIO_CMF_PS1_SV3.c : pmSourceIO_CMF.c.in mksource.pl + mksource.pl pmSourceIO_CMF.c.in PS1_SV3 pmSourceIO_CMF_PS1_SV3.c + +pmSourceIO_CMF_PS1_SV4.c : pmSourceIO_CMF.c.in mksource.pl + mksource.pl pmSourceIO_CMF.c.in PS1_SV4 pmSourceIO_CMF_PS1_SV4.c + # EXTRA_DIST = pmErrorCodes.h.in pmErrorCodes.dat pmErrorCodes.c.in Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/mksource.pl =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/mksource.pl (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/mksource.pl (revision 37403) @@ -24,7 +24,10 @@ "PS1_DV2", 2, "PS1_DV3", 3, + "PS1_DV4", 4, ); %cmfmodes_sv = ("PS1_SV1", 1, "PS1_SV2", 2, + "PS1_SV3", 3, + "PS1_SV4", 4, ); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/fwhm.sh =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/fwhm.sh (revision 37403) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/fwhm.sh (revision 37403) @@ -0,0 +1,218 @@ + +macro find.fwhm.pgauss + if ($0 != 1) + echo "USAGE: find.fwhm" + break + end + + $dz = 0.01 + create z 0 50 $dz + + # Gaussian taylor expansion + # f = 1 / (1 + z + z^2/2 + z^3/6) + # 1 + z + z^2/2 + z^3/6 = 2 + # f = z + z^2/2 + z^3/6 - 2, find where f == 0.0 + + set f = z + 0.5*z^2 + (1/6.0)*z^3 - 2.0 + set dfdz = 1 + z + 0.5*z^2 + + lim -n 0 z f; clear; box; plot z f + lim -n 1 z dfdz; clear; box; plot z dfdz + + $nZ0 = 0 + $nZ1 = 5 / $dz + + $Zg = 0.5*(z[$nZ0] + z[$nZ1]) + $nZg = int($Zg / $dz) + $dZ = 1.0 + + for i 0 10 + $Fg = f[$nZg] + $dFdz_g = dfdz[$nZg] + + $dZ = $Fg / $dFdz_g + + echo $Zg $Fg $dZ $dFdz_g + + $Zg -= $dZ + $Zg = max ($Zg , 0) + $nZg = int($Zg / $dz) + end + echo $Zg $Fg $dZ $dFdz_g + $Zhm = $Zg + $FWHM = 2*sqrt(2*$Zg) + + echo $Zhm : $FWHM +end + +macro find.fwhm.rgauss + if ($0 != 2) + echo "USAGE: find.fwhm (K)" + break + end + + $K = $1 + + if ($K == 0.0) + $Zhm = (sqrt(5) - 1.0) / 2.0 + $FWHM = 2*sqrt(2*$Zhm) + echo $K : $Zhm : $FWHM + return + end + + $dz = 0.01 + create z 0 50 $dz + + # set f = 1.0 / (1.0 + $K*z + z^1.667) + # f = $K*z + z^1.667 - 1.0, find where f == 0.0 + + set f = z + z^$K - 1.0 + set dfdz = ln(z) * z^$K + 1.0 + + lim -n 0 z f; clear; box; plot z f + lim -n 1 z dfdz; clear; box; plot z dfdz + + $nZ0 = 0 + $nZ1 = 5 / $dz + + $Zg = 0.5*(z[$nZ0] + z[$nZ1]) + $nZg = int($Zg / $dz) + $dZ = 1.0 + + for i 0 10 + $Fg = f[$nZg] + $dFdz_g = dfdz[$nZg] + + $dZ = $Fg / $dFdz_g + + echo $nZg $Zg $Fg $dZ $dFdz_g + + $Zg -= $dZ + $Zg = max ($Zg , 0) + $nZg = int($Zg / $dz) + end + # echo $Zg $Fg $dZ $dFdz_g + $Zhm = $Zg + $FWHM = 2*sqrt(2*$Zg) + + echo $K : $Zhm : $FWHM +end + +macro find.fwhm.ps1v1 + if ($0 != 2) + echo "USAGE: find.fwhm (K)" + break + end + + $K = $1 + + if ($K == 0.0) + $Zhm = 1.0 + $FWHM = 2*sqrt(2) + echo $K : $Zhm : $FWHM + return + end + + $dz = 0.01 + create z 0 50 $dz + + # set f = 1.0 / (1.0 + $K*z + z^1.667) + # f = $K*z + z^1.667 - 1.0, find where f == 0.0 + + set f = $K*z + z^1.667 - 1.0 + set dfdz = $K + z^0.667 + + #lim -n 0 z f; clear; box; plot z f + #lim -n 1 z dfdz; clear; box; plot z dfdz + + $nZ0 = 0 + $nZ1 = 5 / $dz + + $Zg = 0.5*(z[$nZ0] + z[$nZ1]) + $nZg = int($Zg / $dz) + $dZ = 1.0 + + for i 0 10 + $Fg = f[$nZg] + $dFdz_g = dfdz[$nZg] + + $dZ = $Fg / $dFdz_g + + # echo $Zg $Fg $dZ $dFdz_g + + $Zg -= $dZ + $Zg = max ($Zg , 0) + $nZg = int($Zg / $dz) + end + # echo $Zg $Fg $dZ $dFdz_g + $Zhm = $Zg + $FWHM = 2*sqrt(2*$Zg) + + echo $K : $Zhm : $FWHM +end + +# qgauss is like ps1_v1 with z^2.25 +macro find.fwhm.qgauss + if ($0 != 2) + echo "USAGE: find.qgauss (K)" + break + end + + $K = $1 + + if ($K == 0.0) + $Zhm = 1.0 + $FWHM = 2*sqrt(2) + echo $K : $Zhm : $FWHM + return + end + + $dz = 0.01 + create z 0 50 $dz + + # set f = 1.0 / (1.0 + $K*z + z^2.25) + # f = $K*z + z^2.25 - 1.0, find where f == 0.0 + + set f = $K*z + z^2.25 - 1.0 + set dfdz = $K + z^1.25 + + #lim -n 0 z f; clear; box; plot z f + #lim -n 1 z dfdz; clear; box; plot z dfdz + + $nZ0 = 0 + $nZ1 = 5 / $dz + + $Zg = 0.5*(z[$nZ0] + z[$nZ1]) + $nZg = int($Zg / $dz) + $dZ = 1.0 + + for i 0 10 + $Fg = f[$nZg] + $dFdz_g = dfdz[$nZg] + + $dZ = $Fg / $dFdz_g + + # echo $Zg $Fg $dZ $dFdz_g + + $Zg -= $dZ + $Zg = max ($Zg , 0) + $nZg = int($Zg / $dz) + end + # echo $Zg $Fg $dZ $dFdz_g + $Zhm = $Zg + $FWHM = 2*sqrt(2*$Zg) + + echo $K : $Zhm : $FWHM +end + +macro fwhm.trend + + delete fwhm_v k_v + + for k 0 20 1.0 + find.fwhm.qgauss $k + concat $k k_v + concat $FWHM fwhm_v + end +end + Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_DEV.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_DEV.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_DEV.c (revision 37403) @@ -36,7 +36,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -59,4 +59,5 @@ # define PM_MODEL_LIMITS pmModelLimits_DEV # define PM_MODEL_RADIUS pmModelRadius_DEV +# define PM_MODEL_SET_FWHM pmModelSetFWHM_DEV # define PM_MODEL_FROM_PSF pmModelFromPSF_DEV # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_DEV @@ -123,5 +124,5 @@ if (radius <= 1.5) { // Nsub ~ 10*index^2 + 1 - psEllipseAxes axes = pmPSF_ModelToAxes(PAR, pmModelClassGetType ("PS_MODEL_DEV")); + psEllipseAxes axes = pmPSF_ModelToAxes(PAR, true); // DEV uses Reff int Nsub = 2 * ((int)(25 / axes.minor)) + 1; Nsub = PS_MIN (Nsub, 121); @@ -336,4 +337,8 @@ } +psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) { + return (NAN); +} + bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) { @@ -357,5 +362,5 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - bool useReff = pmModelUseReff (modelPSF->type); + bool useReff = modelPSF->class->useReff; if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); @@ -411,5 +416,5 @@ // convert to shape terms (SXX,SYY,SXY) // XXX user-defined value for limit? - bool useReff = pmModelUseReff (model->type); + bool useReff = model->class->useReff; if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_DEV.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_DEV.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_DEV.h (revision 37403) @@ -8,4 +8,5 @@ psF64 pmModelFlux_DEV(const psVector *params); psF64 pmModelRadius_DEV(const psVector *params, psF64 flux); +psF64 pmModelSetFWHM_DEV(const psVector *params, psF64 flux); bool pmModelFromPSF_DEV(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); bool pmModelParamsFromPSF_DEV(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_EXP.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_EXP.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_EXP.c (revision 37403) @@ -33,7 +33,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -56,4 +56,5 @@ # define PM_MODEL_LIMITS pmModelLimits_EXP # define PM_MODEL_RADIUS pmModelRadius_EXP +# define PM_MODEL_SET_FWHM pmModelSetFWHM_EXP # define PM_MODEL_FROM_PSF pmModelFromPSF_EXP # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_EXP @@ -343,4 +344,8 @@ } +psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) { + return (NAN); +} + bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) { @@ -364,5 +369,5 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - bool useReff = pmModelUseReff (modelPSF->type); + bool useReff = modelPSF->class->useReff; if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); @@ -418,5 +423,5 @@ // convert to shape terms (SXX,SYY,SXY) // XXX user-defined value for limit? - bool useReff = pmModelUseReff (model->type); + bool useReff = model->class->useReff; if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_EXP.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_EXP.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_EXP.h (revision 37403) @@ -8,4 +8,5 @@ psF64 pmModelFlux_EXP(const psVector *params); psF64 pmModelRadius_EXP(const psVector *params, psF64 flux); +psF64 pmModelSetFWHM_EXP(const psVector *params, psF64 flux); bool pmModelFromPSF_EXP(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); bool pmModelParamsFromPSF_EXP(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_GAUSS.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_GAUSS.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_GAUSS.c (revision 37403) @@ -33,7 +33,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -54,4 +54,5 @@ # define PM_MODEL_GUESS pmModelGuess_GAUSS # define PM_MODEL_LIMITS pmModelLimits_GAUSS +# define PM_MODEL_SET_FWHM pmModelSetFWHM_GAUSS # define PM_MODEL_RADIUS pmModelRadius_GAUSS # define PM_MODEL_FROM_PSF pmModelFromPSF_GAUSS @@ -257,4 +258,8 @@ } +psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) { + return (2.35482004503*sigma); +} + // construct the PSF model from the FLT model and the psf bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) @@ -279,5 +284,5 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - bool useReff = pmModelUseReff (modelPSF->type); + bool useReff = modelPSF->class->useReff; if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); @@ -331,5 +336,5 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - bool useReff = pmModelUseReff (model->type); + bool useReff = model->class->useReff; if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_GAUSS.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_GAUSS.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_GAUSS.h (revision 37403) @@ -8,4 +8,5 @@ psF64 pmModelFlux_GAUSS(const psVector *params); psF64 pmModelRadius_GAUSS(const psVector *params, psF64 flux); +psF64 pmModelSetFWHM_GAUSS(const psVector *params, psF64 flux); bool pmModelFromPSF_GAUSS(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); bool pmModelParamsFromPSF_GAUSS(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_PGAUSS.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_PGAUSS.c (revision 37403) @@ -33,7 +33,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -55,4 +55,5 @@ # define PM_MODEL_LIMITS pmModelLimits_PGAUSS # define PM_MODEL_RADIUS pmModelRadius_PGAUSS +# define PM_MODEL_SET_FWHM pmModelSetFWHM_PGAUSS # define PM_MODEL_FROM_PSF pmModelFromPSF_PGAUSS # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_PGAUSS @@ -324,4 +325,9 @@ } +// scale factor is constant for PGAUSS, I found it with the fwhm.sh script +psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) { + return (3.0063103*sigma); +} + bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) { @@ -344,5 +350,5 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - bool useReff = pmModelUseReff (modelPSF->type); + bool useReff = modelPSF->class->useReff; if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); @@ -396,5 +402,5 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - bool useReff = pmModelUseReff (model->type); + bool useReff = model->class->useReff; if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_PGAUSS.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_PGAUSS.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_PGAUSS.h (revision 37403) @@ -8,4 +8,5 @@ psF64 pmModelFlux_PGAUSS(const psVector *params); psF64 pmModelRadius_PGAUSS(const psVector *params, psF64 flux); +psF64 pmModelSetFWHM_PGAUSS(const psVector *params, psF64 flux); bool pmModelFromPSF_PGAUSS(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); bool pmModelParamsFromPSF_PGAUSS(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_PS1_V1.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_PS1_V1.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_PS1_V1.c (revision 37403) @@ -35,7 +35,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -57,4 +57,5 @@ # define PM_MODEL_LIMITS pmModelLimits_PS1_V1 # define PM_MODEL_RADIUS pmModelRadius_PS1_V1 +# define PM_MODEL_SET_FWHM pmModelSetFWHM_PS1_V1 # define PM_MODEL_FROM_PSF pmModelFromPSF_PS1_V1 # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_PS1_V1 @@ -336,7 +337,54 @@ } +// I used the script in models/fwhm.sh to generate the trend of FWHM scaling vs the K value +// FWHM = Scale * Sigma (not that PAR[PM_PAR_SXX] = sigma * sqrt(2) +// K : z_hm : FWHM +// 0 : 1.000 : 2.83 +// 1 : 0.597 : 2.19 +// 2 : 0.396 : 1.78 +// 3 : 0.291 : 1.53 +// 4 : 0.232 : 1.36 +// 5 : 0.198 : 1.26 +// 6 : 0.169 : 1.16 +// 7 : 0.142 : 1.07 +// 8 : 0.124 : 0.99 +// 9 : 0.118 : 0.97 +// 10 : 0.106 : 0.92 +// 11 : 0.092 : 0.86 +// 12 : 0.091 : 0.85 +// 13 : 0.080 : 0.80 +// 14 : 0.078 : 0.79 +// 15 : 0.073 : 0.76 +// 16 : 0.063 : 0.71 +// 17 : 0.068 : 0.74 +// 18 : 0.056 : 0.67 +// 19 : 0.058 : 0.68 + +// static float PS1_V1_Core[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0}; +static float PS1_V1_Scale[] = {2.83, 2.19, 1.78, 1.53, 1.36, 1.26, 1.16, 1.07, 0.99, 0.97, 0.92, 0.86, 0.85, 0.80, 0.79, 0.76, 0.71, 0.74, 0.67, 0.68}; + +psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) { + + psF32 *PAR = params->data.F32; + + float core = PAR[PM_PAR_7]; + + if (!isfinite(core)) return (2.0*M_SQRT2*sigma); + + // if PS1_V1_Core is defined as a set of integer steps, so we can simplify: + int binCore = MAX(0, MIN (19, (int)core)); + + float scale = NAN; + if (binCore == 0) { + scale = (core - binCore + 0) * (PS1_V1_Scale[binCore + 1] - PS1_V1_Scale[binCore + 0]) + PS1_V1_Scale[binCore + 0]; + } else { + scale = (core - binCore - 1) * (PS1_V1_Scale[binCore + 0] - PS1_V1_Scale[binCore - 1]) + PS1_V1_Scale[binCore - 1]; + } + + return (scale * sigma); +} + bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) { - psF32 *out = modelPSF->params->data.F32; psF32 *in = modelFLT->params->data.F32; @@ -357,5 +405,5 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - bool useReff = pmModelUseReff (modelPSF->type); + bool useReff = modelPSF->class->useReff; if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); @@ -411,5 +459,5 @@ // convert to shape terms (SXX,SYY,SXY) // XXX user-defined value for limit? - bool useReff = pmModelUseReff (model->type); + bool useReff = model->class->useReff; if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_PS1_V1.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_PS1_V1.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_PS1_V1.h (revision 37403) @@ -8,4 +8,5 @@ psF64 pmModelFlux_PS1_V1(const psVector *params); psF64 pmModelRadius_PS1_V1(const psVector *params, psF64 flux); +psF64 pmModelSetFWHM_PS1_V1(const psVector *params, psF64 flux); bool pmModelFromPSF_PS1_V1(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); bool pmModelParamsFromPSF_PS1_V1(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_QGAUSS.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_QGAUSS.c (revision 37403) @@ -35,7 +35,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -57,4 +57,5 @@ # define PM_MODEL_LIMITS pmModelLimits_QGAUSS # define PM_MODEL_RADIUS pmModelRadius_QGAUSS +# define PM_MODEL_SET_FWHM pmModelSetFWHM_QGAUSS # define PM_MODEL_FROM_PSF pmModelFromPSF_QGAUSS # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_QGAUSS @@ -337,4 +338,52 @@ } +// I used the script in models/fwhm.sh to generate the trend of FWHM scaling vs the K value +// FWHM = Scale * Sigma (not that PAR[PM_PAR_SXX] = sigma * sqrt(2) +// K : z_hm : FWHM +// 0 : 1.000 : 2.83 +// 1 : 0.648 : 2.28 +// 2 : 0.430 : 1.85 +// 3 : 0.310 : 1.58 +// 4 : 0.244 : 1.40 +// 5 : 0.200 : 1.26 +// 6 : 0.165 : 1.15 +// 7 : 0.149 : 1.09 +// 8 : 0.125 : 1.00 +// 9 : 0.116 : 0.96 +// 10 : 0.101 : 0.90 +// 11 : 0.095 : 0.87 +// 12 : 0.083 : 0.82 +// 13 : 0.082 : 0.81 +// 14 : 0.080 : 0.80 +// 15 : 0.074 : 0.77 +// 16 : 0.064 : 0.71 +// 17 : 0.068 : 0.74 +// 18 : 0.057 : 0.67 +// 19 : 0.058 : 0.68 + +// static float QGAUSS_Core[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0}; +static float QGAUSS_Scale[] = {2.83, 2.28, 1.85, 1.58, 1.40, 1.26, 1.15, 1.09, 1.00, 0.96, 0.90, 0.87, 0.82, 0.81, 0.80, 0.77, 0.71, 0.74, 0.67, 0.68}; + +psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) { + + psF32 *PAR = params->data.F32; + + float core = PAR[PM_PAR_7]; + + if (!isfinite(core)) return (2.0*M_SQRT2*sigma); + + // QGAUSS_Core is defined as a set of integer steps, so we can simplify: + int binCore = MAX(0, MIN (19, (int)core)); + + float scale = NAN; + if (binCore == 0) { + scale = (core - binCore + 0) * (QGAUSS_Scale[binCore + 1] - QGAUSS_Scale[binCore + 0]) + QGAUSS_Scale[binCore + 0]; + } else { + scale = (core - binCore - 1) * (QGAUSS_Scale[binCore + 0] - QGAUSS_Scale[binCore - 1]) + QGAUSS_Scale[binCore - 1]; + } + + return (scale * sigma); +} + bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) { @@ -358,5 +407,5 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - bool useReff = pmModelUseReff (modelPSF->type); + bool useReff = modelPSF->class->useReff; if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); @@ -416,5 +465,5 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - bool useReff = pmModelUseReff (model->type); + bool useReff = model->class->useReff; if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_QGAUSS.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_QGAUSS.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_QGAUSS.h (revision 37403) @@ -8,4 +8,5 @@ psF64 pmModelFlux_QGAUSS(const psVector *params); psF64 pmModelRadius_QGAUSS(const psVector *params, psF64 flux); +psF64 pmModelSetFWHM_QGAUSS(const psVector *params, psF64 flux); bool pmModelFromPSF_QGAUSS(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); bool pmModelParamsFromPSF_QGAUSS(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_RGAUSS.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_RGAUSS.c (revision 37403) @@ -34,7 +34,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -56,4 +56,5 @@ # define PM_MODEL_LIMITS pmModelLimits_RGAUSS # define PM_MODEL_RADIUS pmModelRadius_RGAUSS +# define PM_MODEL_SET_FWHM pmModelSetFWHM_RGAUSS # define PM_MODEL_FROM_PSF pmModelFromPSF_RGAUSS # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_RGAUSS @@ -330,4 +331,8 @@ } +psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) { + return (NAN); +} + bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) { @@ -351,5 +356,5 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - bool useReff = pmModelUseReff (modelPSF->type); + bool useReff = modelPSF->class->useReff; if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); @@ -404,5 +409,5 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - bool useReff = pmModelUseReff (model->type); + bool useReff = model->class->useReff; if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_RGAUSS.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_RGAUSS.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_RGAUSS.h (revision 37403) @@ -8,4 +8,5 @@ psF64 pmModelFlux_RGAUSS(const psVector *params); psF64 pmModelRadius_RGAUSS(const psVector *params, psF64 flux); +psF64 pmModelSetFWHM_RGAUSS(const psVector *params, psF64 flux); bool pmModelFromPSF_RGAUSS(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); bool pmModelParamsFromPSF_RGAUSS(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_SERSIC.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_SERSIC.c (revision 37403) @@ -43,7 +43,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -66,4 +66,5 @@ # define PM_MODEL_LIMITS pmModelLimits_SERSIC # define PM_MODEL_RADIUS pmModelRadius_SERSIC +# define PM_MODEL_SET_FWHM pmModelSetFWHM_SERSIC # define PM_MODEL_FROM_PSF pmModelFromPSF_SERSIC # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_SERSIC @@ -124,5 +125,5 @@ if (radius <= 1.5) { // Nsub ~ 10*index^2 + 1 - psEllipseAxes axes = pmPSF_ModelToAxes(PAR, pmModelClassGetType ("PS_MODEL_SERSIC")); + psEllipseAxes axes = pmPSF_ModelToAxes(PAR, true); // SERSIC model uses Reff int Nsub = 2 * ((int)(6.0*Sindex / axes.minor)) + 1; Nsub = PS_MIN (Nsub, 121); @@ -357,4 +358,8 @@ } +psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) { + return (NAN); +} + bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) { @@ -378,5 +383,5 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - bool useReff = pmModelUseReff (modelPSF->type); + bool useReff = modelPSF->class->useReff; if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); @@ -432,5 +437,5 @@ // convert to shape terms (SXX,SYY,SXY) // XXX user-defined value for limit? - bool useReff = pmModelUseReff (model->type); + bool useReff = model->class->useReff; if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_SERSIC.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_SERSIC.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_SERSIC.h (revision 37403) @@ -8,4 +8,5 @@ psF64 pmModelFlux_SERSIC(const psVector *params); psF64 pmModelRadius_SERSIC(const psVector *params, psF64 flux); +psF64 pmModelSetFWHM_SERSIC(const psVector *params, psF64 flux); bool pmModelFromPSF_SERSIC(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); bool pmModelParamsFromPSF_SERSIC(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_TRAIL.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_TRAIL.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_TRAIL.c (revision 37403) @@ -5,4 +5,5 @@ * The meaning of the parameters may thus vary depending on the specifics of the model. * All models which are used as a PSF representations share a few parameters, for which # +#include "pmModelClass.h" * define names are listed in pmModel.h: @@ -33,7 +34,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -55,4 +56,5 @@ # define PM_MODEL_LIMITS pmModelLimits_TRAIL # define PM_MODEL_RADIUS pmModelRadius_TRAIL +# define PM_MODEL_SET_FWHM pmModelSetFWHM_TRAIL # define PM_MODEL_FROM_PSF pmModelFromPSF_TRAIL # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_TRAIL @@ -352,5 +354,5 @@ psF32 *psfPAR = source->modelPSF->params->data.F32; - bool useReff = pmModelUseReff (source->modelPSF->type); + bool useReff = source->modelPSF->class->useReff; psEllipseAxes psfAxes; @@ -413,4 +415,8 @@ // PAR_LENGTH is the unconvolved length. add a bit for safety return (0.5*PAR[PM_PAR_LENGTH] + 2); +} + +psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) { + return (NAN); } Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_TRAIL.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_TRAIL.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_TRAIL.h (revision 37403) @@ -8,4 +8,5 @@ psF64 pmModelFlux_TRAIL(const psVector *params); psF64 pmModelRadius_TRAIL(const psVector *params, psF64 flux); +psF64 pmModelSetFWHM_TRAIL(const psVector *params, psF64 flux); bool pmModelFromPSF_TRAIL(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); bool pmModelParamsFromPSF_TRAIL(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmDetEff.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmDetEff.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmDetEff.c (revision 37403) @@ -17,7 +17,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmGrowthCurve.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmGrowthCurve.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmGrowthCurve.c (revision 37403) @@ -30,7 +30,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmGrowthCurveGenerate.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmGrowthCurveGenerate.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmGrowthCurveGenerate.c (revision 37403) @@ -37,7 +37,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModel.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModel.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModel.c (revision 37403) @@ -33,10 +33,12 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" -#include "pmModelClass.h" static void modelFree(pmModel *tmp) { psTrace("psModules.objects", 10, "---- %s() begin ----\n", __func__); + if (!tmp) return; + psFree(tmp->params); psFree(tmp->dparams); @@ -90,13 +92,15 @@ } - tmp->modelFunc = class->modelFunc; - tmp->modelFlux = class->modelFlux; - tmp->modelRadius = class->modelRadius; - tmp->modelLimits = class->modelLimits; - tmp->modelGuess = class->modelGuess; - tmp->modelFromPSF = class->modelFromPSF; - tmp->modelParamsFromPSF = class->modelParamsFromPSF; - tmp->modelFitStatus = class->modelFitStatus; - tmp->modelSetLimits = class->modelSetLimits; + tmp->class = class; + + // tmp->modelFunc = class->modelFunc; + // tmp->modelFlux = class->modelFlux; + // tmp->modelRadius = class->modelRadius; + // tmp->modelLimits = class->modelLimits; + // tmp->modelGuess = class->modelGuess; + // tmp->modelFromPSF = class->modelFromPSF; + // tmp->modelParamsFromPSF = class->modelParamsFromPSF; + // tmp->modelFitStatus = class->modelFitStatus; + // tmp->modelSetLimits = class->modelSetLimits; psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__); @@ -156,5 +160,5 @@ psF32 tmpF; - tmpF = model->modelFunc (NULL, model->params, x); + tmpF = model->class->modelFunc (NULL, model->params, x); psFree(x); psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__); @@ -176,5 +180,5 @@ psF32 tmpF; - tmpF = model->modelFunc (NULL, model->params, x); + tmpF = model->class->modelFunc (NULL, model->params, x); psFree(x); psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__); @@ -283,5 +287,5 @@ // add in the desired components for this coordinate if (mode & PM_MODEL_OP_FUNC) { - pixelValue += model->modelFunc (NULL, params, x); + pixelValue += model->class->modelFunc (NULL, params, x); } Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModel.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModel.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModel.h (revision 37403) @@ -47,14 +47,17 @@ bool isPCM; ///< is this model fitted with PSF-convolution? + pmModelClass *class; + // functions for this model which depend on the model class - pmModelFunc modelFunc; - pmModelFlux modelFlux; - pmModelRadius modelRadius; - pmModelLimits modelLimits; - pmModelGuessFunc modelGuess; - pmModelFromPSFFunc modelFromPSF; - pmModelParamsFromPSF modelParamsFromPSF; - pmModelFitStatusFunc modelFitStatus; - pmModelSetLimitsFunc modelSetLimits; + + // pmModelFunc modelFunc; + // pmModelFlux modelFlux; + // pmModelRadius modelRadius; + // pmModelLimits modelLimits; + // pmModelGuessFunc modelGuess; + // pmModelFromPSFFunc modelFromPSF; + // pmModelParamsFromPSF modelParamsFromPSF; + // pmModelFitStatusFunc modelFitStatus; + // pmModelSetLimitsFunc modelSetLimits; }; Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModelClass.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModelClass.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModelClass.c (revision 37403) @@ -33,7 +33,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmErrorCodes.h" @@ -54,13 +54,13 @@ static pmModelClass defaultModels[] = { - {"PS_MODEL_GAUSS", 7, (pmModelFunc)pmModelFunc_GAUSS, (pmModelFlux)pmModelFlux_GAUSS, (pmModelRadius)pmModelRadius_GAUSS, (pmModelLimits)pmModelLimits_GAUSS, (pmModelGuessFunc)pmModelGuess_GAUSS, (pmModelFromPSFFunc)pmModelFromPSF_GAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_GAUSS, (pmModelFitStatusFunc)pmModelFitStatus_GAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_GAUSS }, - {"PS_MODEL_PGAUSS", 7, (pmModelFunc)pmModelFunc_PGAUSS, (pmModelFlux)pmModelFlux_PGAUSS, (pmModelRadius)pmModelRadius_PGAUSS, (pmModelLimits)pmModelLimits_PGAUSS, (pmModelGuessFunc)pmModelGuess_PGAUSS, (pmModelFromPSFFunc)pmModelFromPSF_PGAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_PGAUSS, (pmModelFitStatusFunc)pmModelFitStatus_PGAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_PGAUSS }, - {"PS_MODEL_QGAUSS", 8, (pmModelFunc)pmModelFunc_QGAUSS, (pmModelFlux)pmModelFlux_QGAUSS, (pmModelRadius)pmModelRadius_QGAUSS, (pmModelLimits)pmModelLimits_QGAUSS, (pmModelGuessFunc)pmModelGuess_QGAUSS, (pmModelFromPSFFunc)pmModelFromPSF_QGAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_QGAUSS, (pmModelFitStatusFunc)pmModelFitStatus_QGAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_QGAUSS }, - {"PS_MODEL_PS1_V1", 8, (pmModelFunc)pmModelFunc_PS1_V1, (pmModelFlux)pmModelFlux_PS1_V1, (pmModelRadius)pmModelRadius_PS1_V1, (pmModelLimits)pmModelLimits_PS1_V1, (pmModelGuessFunc)pmModelGuess_PS1_V1, (pmModelFromPSFFunc)pmModelFromPSF_PS1_V1, (pmModelParamsFromPSF)pmModelParamsFromPSF_PS1_V1, (pmModelFitStatusFunc)pmModelFitStatus_PS1_V1, (pmModelSetLimitsFunc)pmModelSetLimits_PS1_V1 }, - {"PS_MODEL_RGAUSS", 8, (pmModelFunc)pmModelFunc_RGAUSS, (pmModelFlux)pmModelFlux_RGAUSS, (pmModelRadius)pmModelRadius_RGAUSS, (pmModelLimits)pmModelLimits_RGAUSS, (pmModelGuessFunc)pmModelGuess_RGAUSS, (pmModelFromPSFFunc)pmModelFromPSF_RGAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_RGAUSS, (pmModelFitStatusFunc)pmModelFitStatus_RGAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_RGAUSS }, - {"PS_MODEL_SERSIC", 8, (pmModelFunc)pmModelFunc_SERSIC, (pmModelFlux)pmModelFlux_SERSIC, (pmModelRadius)pmModelRadius_SERSIC, (pmModelLimits)pmModelLimits_SERSIC, (pmModelGuessFunc)pmModelGuess_SERSIC, (pmModelFromPSFFunc)pmModelFromPSF_SERSIC, (pmModelParamsFromPSF)pmModelParamsFromPSF_SERSIC, (pmModelFitStatusFunc)pmModelFitStatus_SERSIC, (pmModelSetLimitsFunc)pmModelSetLimits_SERSIC }, - {"PS_MODEL_EXP", 7, (pmModelFunc)pmModelFunc_EXP, (pmModelFlux)pmModelFlux_EXP, (pmModelRadius)pmModelRadius_EXP, (pmModelLimits)pmModelLimits_EXP, (pmModelGuessFunc)pmModelGuess_EXP, (pmModelFromPSFFunc)pmModelFromPSF_EXP, (pmModelParamsFromPSF)pmModelParamsFromPSF_EXP, (pmModelFitStatusFunc)pmModelFitStatus_EXP, (pmModelSetLimitsFunc)pmModelSetLimits_EXP }, - {"PS_MODEL_DEV", 7, (pmModelFunc)pmModelFunc_DEV, (pmModelFlux)pmModelFlux_DEV, (pmModelRadius)pmModelRadius_DEV, (pmModelLimits)pmModelLimits_DEV, (pmModelGuessFunc)pmModelGuess_DEV, (pmModelFromPSFFunc)pmModelFromPSF_DEV, (pmModelParamsFromPSF)pmModelParamsFromPSF_DEV, (pmModelFitStatusFunc)pmModelFitStatus_DEV, (pmModelSetLimitsFunc)pmModelSetLimits_DEV }, - {"PS_MODEL_TRAIL", 7, (pmModelFunc)pmModelFunc_TRAIL, (pmModelFlux)pmModelFlux_TRAIL, (pmModelRadius)pmModelRadius_TRAIL, (pmModelLimits)pmModelLimits_TRAIL, (pmModelGuessFunc)pmModelGuess_TRAIL, (pmModelFromPSFFunc)pmModelFromPSF_TRAIL, (pmModelParamsFromPSF)pmModelParamsFromPSF_TRAIL, (pmModelFitStatusFunc)pmModelFitStatus_TRAIL, (pmModelSetLimitsFunc)pmModelSetLimits_TRAIL }, + {"PS_MODEL_GAUSS", 7, 0, (pmModelFunc)pmModelFunc_GAUSS, (pmModelFlux)pmModelFlux_GAUSS, (pmModelRadius)pmModelRadius_GAUSS, (pmModelSetFWHM)pmModelSetFWHM_GAUSS, (pmModelLimits)pmModelLimits_GAUSS, (pmModelGuessFunc)pmModelGuess_GAUSS, (pmModelFromPSFFunc)pmModelFromPSF_GAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_GAUSS, (pmModelFitStatusFunc)pmModelFitStatus_GAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_GAUSS }, + {"PS_MODEL_PGAUSS", 7, 0, (pmModelFunc)pmModelFunc_PGAUSS, (pmModelFlux)pmModelFlux_PGAUSS, (pmModelRadius)pmModelRadius_PGAUSS, (pmModelSetFWHM)pmModelSetFWHM_PGAUSS, (pmModelLimits)pmModelLimits_PGAUSS, (pmModelGuessFunc)pmModelGuess_PGAUSS, (pmModelFromPSFFunc)pmModelFromPSF_PGAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_PGAUSS, (pmModelFitStatusFunc)pmModelFitStatus_PGAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_PGAUSS }, + {"PS_MODEL_QGAUSS", 8, 0, (pmModelFunc)pmModelFunc_QGAUSS, (pmModelFlux)pmModelFlux_QGAUSS, (pmModelRadius)pmModelRadius_QGAUSS, (pmModelSetFWHM)pmModelSetFWHM_QGAUSS, (pmModelLimits)pmModelLimits_QGAUSS, (pmModelGuessFunc)pmModelGuess_QGAUSS, (pmModelFromPSFFunc)pmModelFromPSF_QGAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_QGAUSS, (pmModelFitStatusFunc)pmModelFitStatus_QGAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_QGAUSS }, + {"PS_MODEL_PS1_V1", 8, 0, (pmModelFunc)pmModelFunc_PS1_V1, (pmModelFlux)pmModelFlux_PS1_V1, (pmModelRadius)pmModelRadius_PS1_V1, (pmModelSetFWHM)pmModelSetFWHM_PS1_V1, (pmModelLimits)pmModelLimits_PS1_V1, (pmModelGuessFunc)pmModelGuess_PS1_V1, (pmModelFromPSFFunc)pmModelFromPSF_PS1_V1, (pmModelParamsFromPSF)pmModelParamsFromPSF_PS1_V1, (pmModelFitStatusFunc)pmModelFitStatus_PS1_V1, (pmModelSetLimitsFunc)pmModelSetLimits_PS1_V1 }, + {"PS_MODEL_RGAUSS", 8, 0, (pmModelFunc)pmModelFunc_RGAUSS, (pmModelFlux)pmModelFlux_RGAUSS, (pmModelRadius)pmModelRadius_RGAUSS, (pmModelSetFWHM)pmModelSetFWHM_RGAUSS, (pmModelLimits)pmModelLimits_RGAUSS, (pmModelGuessFunc)pmModelGuess_RGAUSS, (pmModelFromPSFFunc)pmModelFromPSF_RGAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_RGAUSS, (pmModelFitStatusFunc)pmModelFitStatus_RGAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_RGAUSS }, + {"PS_MODEL_SERSIC", 8, 1, (pmModelFunc)pmModelFunc_SERSIC, (pmModelFlux)pmModelFlux_SERSIC, (pmModelRadius)pmModelRadius_SERSIC, (pmModelSetFWHM)pmModelSetFWHM_SERSIC, (pmModelLimits)pmModelLimits_SERSIC, (pmModelGuessFunc)pmModelGuess_SERSIC, (pmModelFromPSFFunc)pmModelFromPSF_SERSIC, (pmModelParamsFromPSF)pmModelParamsFromPSF_SERSIC, (pmModelFitStatusFunc)pmModelFitStatus_SERSIC, (pmModelSetLimitsFunc)pmModelSetLimits_SERSIC }, + {"PS_MODEL_EXP", 7, 1, (pmModelFunc)pmModelFunc_EXP, (pmModelFlux)pmModelFlux_EXP, (pmModelRadius)pmModelRadius_EXP, (pmModelSetFWHM)pmModelSetFWHM_EXP, (pmModelLimits)pmModelLimits_EXP, (pmModelGuessFunc)pmModelGuess_EXP, (pmModelFromPSFFunc)pmModelFromPSF_EXP, (pmModelParamsFromPSF)pmModelParamsFromPSF_EXP, (pmModelFitStatusFunc)pmModelFitStatus_EXP, (pmModelSetLimitsFunc)pmModelSetLimits_EXP }, + {"PS_MODEL_DEV", 7, 1, (pmModelFunc)pmModelFunc_DEV, (pmModelFlux)pmModelFlux_DEV, (pmModelRadius)pmModelRadius_DEV, (pmModelSetFWHM)pmModelSetFWHM_DEV, (pmModelLimits)pmModelLimits_DEV, (pmModelGuessFunc)pmModelGuess_DEV, (pmModelFromPSFFunc)pmModelFromPSF_DEV, (pmModelParamsFromPSF)pmModelParamsFromPSF_DEV, (pmModelFitStatusFunc)pmModelFitStatus_DEV, (pmModelSetLimitsFunc)pmModelSetLimits_DEV }, + {"PS_MODEL_TRAIL", 7, 0, (pmModelFunc)pmModelFunc_TRAIL, (pmModelFlux)pmModelFlux_TRAIL, (pmModelRadius)pmModelRadius_TRAIL, (pmModelSetFWHM)pmModelSetFWHM_TRAIL, (pmModelLimits)pmModelLimits_TRAIL, (pmModelGuessFunc)pmModelGuess_TRAIL, (pmModelFromPSFFunc)pmModelFromPSF_TRAIL, (pmModelParamsFromPSF)pmModelParamsFromPSF_TRAIL, (pmModelFitStatusFunc)pmModelFitStatus_TRAIL, (pmModelSetLimitsFunc)pmModelSetLimits_TRAIL }, }; Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModelClass.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModelClass.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModelClass.h (revision 37403) @@ -35,7 +35,9 @@ char *name; int nParams; + bool useReff; pmModelFunc modelFunc; pmModelFlux modelFlux; pmModelRadius modelRadius; + pmModelSetFWHM modelSetFWHM; pmModelLimits modelLimits; pmModelGuessFunc modelGuess; Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModelFuncs.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModelFuncs.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModelFuncs.h (revision 37403) @@ -110,4 +110,7 @@ typedef psF64 (*pmModelRadius)(const psVector *params, double flux); +// This function returns the FWHM given the supplied sigma (major or minor) +typedef psF64 (*pmModelSetFWHM)(const psVector *params, double sigma); + // This function provides the model guess parameters based on the details of // the given source. Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModelUtils.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModelUtils.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModelUtils.c (revision 37403) @@ -32,7 +32,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -46,13 +46,4 @@ #include "pmErrorCodes.h" -// XX static bool useModelVar = false; -// XX -// XX void pmModelSetModelVarOption (bool option) { -// XX useModelVar = option; -// XX } -// XX bool pmModelGetModelVarOption (void) { -// XX return useModelVar; -// XX } - /***************************************************************************** pmModelFromPSF (*modelEXT, *psf): use the model position parameters to @@ -68,5 +59,5 @@ // set model parameters for this source based on PSF information - if (!modelEXT->modelFromPSF (modelPSF, modelEXT, psf)) { + if (!modelEXT->class->modelFromPSF (modelPSF, modelEXT, psf)) { psTrace ("psModules.objects", 3, "Failed to set model params from PSF"); psFree(modelPSF); @@ -89,5 +80,5 @@ // set model parameters for this source based on PSF information - if (!modelPSF->modelParamsFromPSF (modelPSF, psf, Xo, Yo, Io)) { + if (!modelPSF->class->modelParamsFromPSF (modelPSF, psf, Xo, Yo, Io)) { psFree(modelPSF); return NULL; @@ -109,5 +100,5 @@ // determine the normalized flux - float normFlux = model->modelFlux (model->params); + float normFlux = model->class->modelFlux (model->params); assert (isfinite(normFlux)); assert (normFlux > 0); @@ -120,8 +111,8 @@ bool pmModelUseReff (pmModelType type) { - bool useReff = false; - useReff |= (type == pmModelClassGetType ("PS_MODEL_SERSIC")); - useReff |= (type == pmModelClassGetType ("PS_MODEL_DEV")); - useReff |= (type == pmModelClassGetType ("PS_MODEL_EXP")); + + pmModelClass *class = pmModelClassSelect (type); + psAssert (class, "undefined model class?"); + bool useReff = class->useReff; return useReff; } Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPCM_MinimizeChisq.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPCM_MinimizeChisq.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPCM_MinimizeChisq.c (revision 37403) @@ -31,7 +31,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -455,5 +455,5 @@ coord->data.F32[1] = (psF32) (i + 0.5 + source->pixels->row0); - pcm->modelFlux->data.F32[i][j] = pcm->modelConv->modelFunc (deriv, params, coord); + pcm->modelFlux->data.F32[i][j] = pcm->modelConv->class->modelFunc (deriv, params, coord); for (int n = 0; n < params->n; n++) { Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPCMdata.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPCMdata.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPCMdata.c (revision 37403) @@ -31,7 +31,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -293,5 +293,5 @@ psEllipseAxes axes; - bool useReff = pmModelUseReff (modelPSF->type); + bool useReff = modelPSF->class->useReff; psF32 *PAR = modelPSF->params->data.F32; pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff); @@ -302,5 +302,5 @@ // XXX need to do this more carefully if (modelPSF->type == modelType_GAUSS) { - float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]); + float FWHM_MAJOR = 2*modelPSF->class->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]); float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); *sigma = 0.50 * (FWHM_MAJOR + FWHM_MINOR) / 2.35; @@ -330,5 +330,5 @@ psEllipseAxes axes; - bool useReff = pmModelUseReff (modelPSF->type); + bool useReff = modelPSF->class->useReff; psF32 *PAR = modelPSF->params->data.F32; pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff); @@ -339,5 +339,5 @@ // XXX need to do this more carefully if (modelPSF->type == modelType_GAUSS) { - float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]); + float FWHM_MAJOR = 2*modelPSF->class->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]); float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); *sigma = 0.50 * (FWHM_MAJOR + FWHM_MINOR) / 2.35; @@ -393,5 +393,5 @@ psMinConstraint *constraint = psMinConstraintAlloc(); constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_VECTOR_MASK); - constraint->checkLimits = model->modelLimits; + constraint->checkLimits = model->class->modelLimits; int nParams = pmPCMsetParams (constraint, fitOptions->mode); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSF.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSF.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSF.c (revision 37403) @@ -37,7 +37,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -329,5 +329,5 @@ // convert the parameters used in the fitted source model to the psEllipseAxes representation // (major,minor,theta) -psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, pmModelType type) +psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, bool useReff) { psEllipseAxes axes; @@ -338,5 +338,4 @@ PS_ASSERT_PTR_NON_NULL(modelPar, axes); - bool useReff = pmModelUseReff (type); pmModelParamsToAxes (&axes, modelPar[PM_PAR_SXX], modelPar[PM_PAR_SXY], modelPar[PM_PAR_SYY], useReff); return axes; @@ -345,5 +344,5 @@ // convert the psEllipseAxes representation (major,minor,theta) to the parameters used in the // fitted source model -bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes, pmModelType type) +bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes, bool useReff) { PS_ASSERT_PTR_NON_NULL(modelPar, false); @@ -357,5 +356,4 @@ } - bool useReff = pmModelUseReff (type); pmModelAxesToParams (&modelPar[PM_PAR_SXX], &modelPar[PM_PAR_SXY], &modelPar[PM_PAR_SYY], axes, useReff); return true; @@ -420,9 +418,9 @@ // get the model full-width at half-max - float fwhmMajor = 2*model->modelRadius (model->params, 0.5); + float fwhmMajor = 2*model->class->modelRadius (model->params, 0.5); # if (0) psF32 *params = model->params->data.F32; // Model parameters - psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO, model->type); // Ellipse axes + psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO, model->class->useReff); // Ellipse axes // Curiously, the minor axis can be larger than the major axis, so need to check. Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSF.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSF.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSF.h (revision 37403) @@ -106,9 +106,9 @@ pmPSF *pmPSFBuildSimple (char *typeName, float sxx, float syy, float sxy, ...); -bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes, pmModelType type); +bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes, bool useReff); bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis, bool useReff); psEllipsePol pmPSF_ModelToFit (psF32 *modelPar, bool useReff); -psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, pmModelType type); +psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, bool useReff); /// Calculate FWHM value from a PSF Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSF_IO.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSF_IO.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSF_IO.c (revision 37403) @@ -47,7 +47,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -511,11 +511,13 @@ psMetadataAddF32 (header, PS_LIST_TAIL, "SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias); - float PSF_APERTURE = psMetadataLookupF32(&status, roAnalysis, "PSF_APERTURE"); - if (status) { + if (roAnalysis) { + float PSF_APERTURE = psMetadataLookupF32(&status, roAnalysis, "PSF_APERTURE"); + if (status) { psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_APERTURE", PS_DATA_F32, "aperture for psf objects", PSF_APERTURE); - } - float PSF_FIT_RADIUS = psMetadataLookupF32(&status, roAnalysis, "PSF_FIT_RADIUS"); - if (status) { + } + float PSF_FIT_RADIUS = psMetadataLookupF32(&status, roAnalysis, "PSF_FIT_RADIUS"); + if (status) { psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_DATA_F32, "aperture for psf objects", PSF_FIT_RADIUS); + } } Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtry.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtry.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtry.c (revision 37403) @@ -29,7 +29,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtryFitEXT.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtryFitEXT.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtryFitEXT.c (revision 37403) @@ -29,7 +29,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtryFitPSF.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtryFitPSF.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtryFitPSF.c (revision 37403) @@ -27,7 +27,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtryMakePSF.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtryMakePSF.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtryMakePSF.c (revision 37403) @@ -28,7 +28,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -212,5 +212,5 @@ assert (source->modelEXT); // all unmasked sources should have modelEXT - bool useReff = pmModelUseReff (source->modelEXT->type); + bool useReff = source->modelEXT->class->useReff; psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32, useReff); @@ -218,4 +218,26 @@ e1->data.F32[i] = pol.e1; e2->data.F32[i] = pol.e2; + } + + // weed out extreme e0 outliers here: find the median and exclude points not in the + // range MEDIAN / 5 < e0 < 5 * MEDIAN + { + psStats *e0stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + if (psVectorStats (e0stats, e0, NULL, srcMask, 0xff)) { + float e0med = e0stats->sampleMedian; + + for (int i = 0; i < sources->n; i++) { + // skip any masked sources (failed to fit one of the model steps or get a magnitude) + if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; + + if (e0->data.F32[i] < 0.2*e0med) { + srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_OUTLIER; + } + if (e0->data.F32[i] > 5.0*e0med) { + srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_OUTLIER; + } + } + } + psFree (e0stats); } Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtryMetric.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtryMetric.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtryMetric.c (revision 37403) @@ -28,7 +28,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtryModel.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtryModel.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPSFtryModel.c (revision 37403) @@ -29,7 +29,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPhotObj.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPhotObj.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPhotObj.c (revision 37403) @@ -29,7 +29,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSource.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSource.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSource.c (revision 37403) @@ -33,7 +33,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -146,4 +146,5 @@ source->apMagRaw = NAN; source->apRadius = NAN; + source->apNpixels = 0; source->apFlux = NAN; source->apFluxErr = NAN; @@ -162,4 +163,5 @@ source->sky = NAN; source->skyErr = NAN; + source->extSN = NAN; source->region = psRegionSet(NAN, NAN, NAN, NAN); @@ -174,4 +176,7 @@ source->parent = NULL; source->tmpPtr = NULL; + source->chipNum = -1; + source->chipX = -1000; + source->chipY = -1000; source->imageID = -1; source->nFrames = 0; @@ -236,4 +241,5 @@ source->apMagRaw = in->apMagRaw; source->apRadius = in->apRadius; + source->apNpixels = in->apNpixels; source->apFlux = in->apFlux; source->apFluxErr = in->apFluxErr; @@ -986,5 +992,5 @@ if (!isfinite(oldI0)) return false; - bool useReff = pmModelUseReff (model->type); + bool useReff = model->class->useReff; pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSource.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSource.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSource.h (revision 37403) @@ -97,4 +97,5 @@ float apMagRaw; ///< raw mag in given aperture float apRadius; ///< radius for aperture magnitude + int apNpixels; ///< number of unmasked pixels in aperture float apFlux; ///< apFlux corresponding to psfMag or extMag (depending on type) float apFluxErr; ///< apFluxErr corresponding to psfMag or extMag (depending on type) @@ -113,4 +114,5 @@ float sky; ///< The sky at the center of the object float skyErr; ///< The sky error at the center of the object + float extSN; ///< for externally supplied source the kron signal to noise (used by full force) psRegion region; ///< area on image covered by selected pixels @@ -119,5 +121,5 @@ pmSourceExtendedPars *extpars; ///< extended source parameters pmSourceDiffStats *diffStats; ///< extra parameters for difference detections - pmSourceGalaxyFits *galaxyFits; ///< fits to galaxy models (psphotFullForce only) + psArray *galaxyFits; ///< fits to galaxy models (psphotFullForce only) pmSourceLensing *lensingOBJ; ///< lensing moments parameters (per object) pmSourceLensing *lensingPSF; ///< lensing moments parameters (psf, interpolated) @@ -125,4 +127,7 @@ pmSource *parent; ///< reference to the master source from which this is derived psPtr *tmpPtr; ///< pointer that may be used to store data in a particular module. e.g. psphotKronIterate. + short chipNum; ///< camera dependent of chip suppling pixels for fullforce source + short chipX; ///< chip space X coord of fullforce source + short chipY; ///< chip space Y coord of fullforce source int imageID; psU16 nFrames; Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceContour.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceContour.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceContour.c (revision 37403) @@ -33,7 +33,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceExtendedPars.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceExtendedPars.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceExtendedPars.h (revision 37403) @@ -83,8 +83,16 @@ typedef struct { + int modelType; psVector *Flux; psVector *dFlux; psVector *chisq; - int nPix; + int nPix; + bool reducedTrials; + float fRmajorMin; + float fRmajorMax; + float fRmajorDel; + float fRminorMin; + float fRminorMax; + float fRminorDel; } pmSourceGalaxyFits; Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceFitModel.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceFitModel.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceFitModel.c (revision 37403) @@ -33,7 +33,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -170,5 +170,5 @@ psMinConstraint *constraint = psMinConstraintAlloc(); constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_VECTOR_MASK); - constraint->checkLimits = model->modelLimits; + constraint->checkLimits = model->class->modelLimits; // set parameter mask based on fitting mode @@ -233,6 +233,6 @@ // force the floating parameters to fall within the contraint ranges for (int i = 0; i < params->n; i++) { - model->modelLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL); - model->modelLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL); + model->class->modelLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL); + model->class->modelLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL); } @@ -254,5 +254,5 @@ psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); - fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, model->modelFunc); + fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, model->class->modelFunc); for (int i = 0; i < dparams->n; i++) { if ((constraint->paramMask != NULL) && constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) @@ -306,5 +306,5 @@ altmask->data.PS_TYPE_VECTOR_MASK_DATA[i] = (constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) ? 0 : 1; } - psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, model->modelFunc); + psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, model->class->modelFunc); for (int i = 0; i < dparams->n; i++) { Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceFitPCM.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceFitPCM.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceFitPCM.c (revision 37403) @@ -31,7 +31,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -63,6 +63,6 @@ // force the floating parameters to fall within the contraint ranges for (int i = 0; i < params->n; i++) { - pcm->modelConv->modelLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL); - pcm->modelConv->modelLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL); + pcm->modelConv->class->modelLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL); + pcm->modelConv->class->modelLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL); } @@ -165,5 +165,5 @@ bool pmSourceModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) { - if (!pcm->modelConv->modelGuess(pcm->modelConv, source, maskVal, markVal)) { + if (!pcm->modelConv->class->modelGuess(pcm->modelConv, source, maskVal, markVal)) { return false; } Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceFitSet.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceFitSet.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceFitSet.c (revision 37403) @@ -32,7 +32,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -223,5 +223,5 @@ float *paramOne = params + nParamBase; float *betaOne = betas + nParamBase; - bool status = model->modelLimits (mode, nParamOne, paramOne, betaOne); + bool status = model->class->modelLimits (mode, nParamOne, paramOne, betaOne); return status; } @@ -388,5 +388,5 @@ psVector *derivOne = thisSet->derivSet->data[i]; - chisqOne = model->modelFunc (derivOne, paramOne, x); + chisqOne = model->class->modelFunc (derivOne, paramOne, x); chisqSum += chisqOne; } Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceGroups.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceGroups.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceGroups.c (revision 37403) @@ -18,7 +18,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO.c (revision 37403) @@ -40,7 +40,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -365,5 +365,5 @@ # define PM_SOURCES_WRITE(NAME,TYPE) \ if (!strcmp (exttype, NAME)) { \ - status &= pmSourcesWrite_##TYPE(file->fits, readout, sources, file->header, outhead, dataname, recipe); \ + status = pmSourcesWrite_##TYPE(file->fits, readout, sources, file->header, outhead, dataname, recipe); \ if (xsrcname) { \ status &= pmSourcesWrite_##TYPE##_XSRC(file->fits, readout, sources, file->header, xsrcname, recipe); \ @@ -589,5 +589,5 @@ // these are case-sensitive since the EXTYPE is case-sensitive - status = true; + status = false; PM_SOURCES_WRITE("SMPDATA", SMPDATA); PM_SOURCES_WRITE("PS1_DEV_0", PS1_DEV_0); @@ -601,7 +601,10 @@ PM_SOURCES_WRITE("PS1_SV1", CMF_PS1_SV1); PM_SOURCES_WRITE("PS1_SV2", CMF_PS1_SV2); + PM_SOURCES_WRITE("PS1_SV3", CMF_PS1_SV3); + PM_SOURCES_WRITE("PS1_SV4", CMF_PS1_SV4); PM_SOURCES_WRITE("PS1_DV1", CMF_PS1_DV1); PM_SOURCES_WRITE("PS1_DV2", CMF_PS1_DV2); PM_SOURCES_WRITE("PS1_DV3", CMF_PS1_DV3); + PM_SOURCES_WRITE("PS1_DV4", CMF_PS1_DV4); psFree (outhead); @@ -1121,7 +1124,15 @@ PM_SOURCES_READ_PSF("PS1_SV1", CMF_PS1_SV1); PM_SOURCES_READ_PSF("PS1_SV2", CMF_PS1_SV2); + PM_SOURCES_READ_PSF("PS1_SV3", CMF_PS1_SV3); + PM_SOURCES_READ_PSF("PS1_SV4", CMF_PS1_SV4); PM_SOURCES_READ_PSF("PS1_DV1", CMF_PS1_DV1); PM_SOURCES_READ_PSF("PS1_DV2", CMF_PS1_DV2); PM_SOURCES_READ_PSF("PS1_DV3", CMF_PS1_DV3); + PM_SOURCES_READ_PSF("PS1_DV4", CMF_PS1_DV4); + + if (!sources) { + psError(PS_ERR_IO, false, "reading CMF data from %s with format %s\n", file->filename, exttype); + return false; + } long *sourceIndex = NULL; @@ -1202,5 +1213,11 @@ break; - case PM_FPA_FILE_CFF: + case PM_FPA_FILE_CFF: { + // determine the output table format + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, "PSPHOT"); + if (!status) { + psError(PS_ERR_UNKNOWN, true, "missing recipe PSPHOT in config data"); + return false; + } // read in header, if not yet loaded hdu = pmFPAviewThisHDU (view, file->fpa); @@ -1244,5 +1261,5 @@ } - sources = pmSourcesRead_CFF(file->fits, hdu->header); + sources = pmSourcesRead_CFF(file->fits, hdu->header, recipe); psTrace("psModules.objects", 6, "read CMF table from %s : %s : %s", file->filename, headname, dataname); @@ -1250,4 +1267,5 @@ psFree (dataname); psFree (tableHeader); + } break; @@ -1395,7 +1413,10 @@ PM_SOURCES_READ_XSRC("PS1_SV1", CMF_PS1_SV1); PM_SOURCES_READ_XSRC("PS1_SV2", CMF_PS1_SV2); + PM_SOURCES_READ_XSRC("PS1_SV3", CMF_PS1_SV3); + PM_SOURCES_READ_XSRC("PS1_SV4", CMF_PS1_SV4); PM_SOURCES_READ_XSRC("PS1_DV1", CMF_PS1_DV1); PM_SOURCES_READ_XSRC("PS1_DV2", CMF_PS1_DV2); PM_SOURCES_READ_XSRC("PS1_DV3", CMF_PS1_DV3); + PM_SOURCES_READ_XSRC("PS1_DV4", CMF_PS1_DV4); } psFree(tableHeader); @@ -1435,7 +1456,10 @@ PM_SOURCES_READ_XFIT("PS1_SV1", CMF_PS1_SV1); PM_SOURCES_READ_XFIT("PS1_SV2", CMF_PS1_SV2); + PM_SOURCES_READ_XFIT("PS1_SV3", CMF_PS1_SV3); + PM_SOURCES_READ_XFIT("PS1_SV4", CMF_PS1_SV4); PM_SOURCES_READ_XFIT("PS1_DV1", CMF_PS1_DV1); PM_SOURCES_READ_XFIT("PS1_DV2", CMF_PS1_DV2); PM_SOURCES_READ_XFIT("PS1_DV3", CMF_PS1_DV3); + PM_SOURCES_READ_XFIT("PS1_DV4", CMF_PS1_DV4); } psFree(tableHeader); @@ -1474,7 +1498,10 @@ PM_SOURCES_READ_XRAD("PS1_SV1", CMF_PS1_SV1); PM_SOURCES_READ_XRAD("PS1_SV2", CMF_PS1_SV2); + PM_SOURCES_READ_XRAD("PS1_SV3", CMF_PS1_SV3); + PM_SOURCES_READ_XRAD("PS1_SV4", CMF_PS1_SV4); PM_SOURCES_READ_XRAD("PS1_DV1", CMF_PS1_DV1); PM_SOURCES_READ_XRAD("PS1_DV2", CMF_PS1_DV2); PM_SOURCES_READ_XRAD("PS1_DV3", CMF_PS1_DV3); + PM_SOURCES_READ_XRAD("PS1_DV4", CMF_PS1_DV4); } psFree(tableHeader); @@ -1513,7 +1540,10 @@ PM_SOURCES_READ_XGAL("PS1_SV1", CMF_PS1_SV1); PM_SOURCES_READ_XGAL("PS1_SV2", CMF_PS1_SV2); + PM_SOURCES_READ_XGAL("PS1_SV3", CMF_PS1_SV3); + PM_SOURCES_READ_XGAL("PS1_SV4", CMF_PS1_SV4); PM_SOURCES_READ_XGAL("PS1_DV1", CMF_PS1_DV1); PM_SOURCES_READ_XGAL("PS1_DV2", CMF_PS1_DV2); PM_SOURCES_READ_XGAL("PS1_DV3", CMF_PS1_DV3); + PM_SOURCES_READ_XGAL("PS1_DV4", CMF_PS1_DV4); } psFree(tableHeader); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO.h (revision 37403) @@ -40,7 +40,10 @@ MK_PROTO(CMF_PS1_SV1); MK_PROTO(CMF_PS1_SV2); +MK_PROTO(CMF_PS1_SV3); +MK_PROTO(CMF_PS1_SV4); MK_PROTO(CMF_PS1_DV1); MK_PROTO(CMF_PS1_DV2); MK_PROTO(CMF_PS1_DV3); +MK_PROTO(CMF_PS1_DV4); int pmSourceGetDophotType (pmSource *source); @@ -55,5 +58,5 @@ psArray *pmSourcesReadCMP (char *filename, psMetadata *header); -psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header); +psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header, psMetadata *recipe); bool pmSourcesWrite_CFF (pmReadout *readout, psFits *fits, psArray *sources, psMetadata *header, psMetadata *recipe); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_CFF.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_CFF.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_CFF.c (revision 37403) @@ -37,7 +37,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -53,6 +53,6 @@ #include "pmSourceOutputs.h" -// read in a readout from the fits file -psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header) +// read in sources readout from a cff fits file +psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header, psMetadata *recipe) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -75,12 +75,42 @@ PS_ASSERT_INT_NONNEGATIVE(modelType, NULL); + pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC"); + pmModelType devModelType = pmModelClassGetType("PS_MODEL_DEV"); + + psString modelForce = psMetadataLookupStr(&status, recipe, "EXT_MODEL_TYPE_FORCE"); + psF32 forceDevSersicMin = NAN; + bool forceAll = false; + int modelTypeForce = 0; + if (!strcmp(modelForce, "ALL")) { + forceAll = true; + } else if (!strcmp(modelForce, "PS_MODEL_SERSIC")) { + modelTypeForce = sersicModelType; + forceDevSersicMin = psMetadataLookupF32(&status, recipe, "EXT_MODEL_FORCE_DEV_SERSIC_MIN"); + } else { + modelTypeForce = pmModelClassGetType(modelForce); + PS_ASSERT_INT_NONNEGATIVE(modelTypeForce, NULL); + } + + // skip extended model types for likely stars + // max value of KronMag - psfMag to keep ... + psF32 starCut = psMetadataLookupF32(&status, recipe, "EXT_MODEL_FORCE_MAGDIFF_MAX"); + if (!status) { + starCut = 0; + } + // ... unless SN is less than this value + psF32 SNMinForCut = psMetadataLookupF32(&status, recipe, "EXT_MODEL_FORCE_CUT_SN_MIN"); + if (!status) { + SNMinForCut = 10; + } + // We get the size of the table, and allocate the array of sources first because the table // is large and ephemeral --- when the table gets blown away, whatever is allocated after // the table is read blocks the free. In fact, it's better to read the table row by row. - long numSources = psFitsTableSize(fits); // Number of sources in table - psArray *sources = psArrayAlloc(numSources); // Array of sources, to return - - // convert the table to the pmSource entriesa - for (int i = 0; i < numSources; i++) { + long numRows = psFitsTableSize(fits); // Number of rows in table + psArray *sources = psArrayAllocEmpty(numRows); // Array of sources, to return + + // convert the table to the pmSource entries + pmSource *source = NULL; + for (int i = 0; i < numRows; i++) { psMetadata *row = psFitsReadTableRow(fits, i); // Table row if (!row) { @@ -108,4 +138,5 @@ float kronRadius = psMetadataLookupF32 (&status, row, "KRON_RADIUS"); float petRadius = psMetadataLookupF32 (&status, row, "PETRO_RADIUS"); + float SN = psMetadataLookupF32 (&status, row, "SN"); bool fitGalaxy = psMetadataLookupU8 (&status, row, "FIT_GALAXY"); bool psfStar = psMetadataLookupU8 (&status, row, "PSF_STAR"); @@ -114,7 +145,11 @@ float Rminor = psMetadataLookupF32 (&status, row, "R_MINOR"); float theta = psMetadataLookupF32 (&status, row, "THETA"); + float chisq = psMetadataLookupF32 (&status, row, "CHISQ"); + float nDOF = psMetadataLookupF32 (&status, row, "NDOF"); + float magDiff = psMetadataLookupF32 (&status, row, "MAG_DIFF"); + psS16 modelFlags = psMetadataLookupS32 (&status, row, "MODEL_FLAGS"); int galaxyModelType = psMetadataLookupS32(&status, row, "MODEL_TYPE"); - if (status) { + if (status && galaxyModelType >= 0) { galaxyModelType = pmModelClassGetLocalType(galaxyModelType); } else { @@ -123,79 +158,112 @@ float Sindex = psMetadataLookupF32 (&status, row, "INDEX"); // Should this be PAR_07 not sersic index - pmSource *source = pmSourceAlloc (); - pmModel *model = pmModelAlloc (modelType); - source->modelPSF = model; -// RoughClass wants source type to be unknown -// source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags - source->type = PM_SOURCE_TYPE_UNKNOWN; - - // XXX we can set this in general, but for a specific image, we need to weed out SATSTARS and - // stars that are masked - if (psfStar) { - source->tmpFlags |= PM_SOURCE_TMPF_CANDIDATE_PSFSTAR; - } - - // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. - psF32 *PAR = model->params->data.F32; - psF32 *dPAR = model->dparams->data.F32; - - source->seq = ID; - PAR[PM_PAR_XPOS] = X; - PAR[PM_PAR_YPOS] = Y; - - dPAR[PM_PAR_XPOS] = 0.0; - dPAR[PM_PAR_YPOS] = 0.0; - - PAR[PM_PAR_SKY] = 0.0; - dPAR[PM_PAR_SKY] = 0.0; - - PAR[PM_PAR_I0] = 1.0; - dPAR[PM_PAR_I0] = 0.0; - - source->sky = PAR[PM_PAR_SKY]; - source->skyErr = dPAR[PM_PAR_SKY]; - - source->psfMag = 0.0; - source->psfMagErr = 0.0; - source->apMag = 0.0; - source->apRadius = apRadius; - - // we generate a somewhat fake PSF model here -- - // in most (all?) contexts, we will replace this with a measured psf model - // elsewhere - axes.major = 1.0; - axes.minor = 1.0; - axes.theta = 0.0; - pmPSF_AxesToModel (PAR, axes, modelType); - - // peak->detValue, rawFlux, smoothFlux are all set to the flux argument which is counts per second - source->peak = pmPeakAlloc(X, Y, flux, PM_PEAK_LONE); - source->peak->xf = X; // pmPeakAlloc converts X,Y to int, so reset here - source->peak->yf = Y; // pmPeakAlloc converts X,Y to int, so reset here - source->peak->dx = 0.0; - source->peak->dy = 0.0; - - source->moments = pmMomentsAlloc (); - source->moments->Mx = X; - source->moments->My = Y; - source->moments->Mrf = kronRadius * 0.4; // kronRadius is 2.5 * first radial moment - - // Don't mark the moments as measured because that causes many fields to be left blank. - // The moments code knows not to change the position or the Mrf for external sources - // source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED; - - if (isfinite(petRadius)) { - source->extpars = pmSourceExtendedParsAlloc (); - source->extpars->petrosianRadius = petRadius; - } - + if (!source || ID != source->seq) { + if (source) { + psArrayAdd (sources, 1, source); + psFree(source); + } + source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + // RoughClass wants source type to be unknown + // source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags + source->type = PM_SOURCE_TYPE_UNKNOWN; + + // XXX we can set this in general, but for a specific image, we need to weed out SATSTARS and + // stars that are masked + if (psfStar) { + source->tmpFlags |= PM_SOURCE_TMPF_CANDIDATE_PSFSTAR; + } + + // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + source->seq = ID; + + PAR[PM_PAR_XPOS] = X; + PAR[PM_PAR_YPOS] = Y; + + dPAR[PM_PAR_XPOS] = 0.0; + dPAR[PM_PAR_YPOS] = 0.0; + + PAR[PM_PAR_SKY] = 0.0; + dPAR[PM_PAR_SKY] = 0.0; + + PAR[PM_PAR_I0] = 1.0; + dPAR[PM_PAR_I0] = 0.0; + + source->sky = PAR[PM_PAR_SKY]; + source->skyErr = dPAR[PM_PAR_SKY]; + + source->psfMag = 0.0; + source->psfMagErr = 0.0; + source->apMag = 0.0; + source->apRadius = apRadius; + + // we generate a somewhat fake PSF model here -- + // in most (all?) contexts, we will replace this with a measured psf model + // elsewhere + axes.major = 1.0; + axes.minor = 1.0; + axes.theta = 0.0; + pmPSF_AxesToModel (PAR, axes, model->class->useReff); + + // peak->detValue, rawFlux, smoothFlux are all set to the flux argument which is counts per second + source->peak = pmPeakAlloc(X, Y, flux, PM_PEAK_LONE); + source->peak->xf = X; // pmPeakAlloc converts X,Y to int, so reset here + source->peak->yf = Y; // pmPeakAlloc converts X,Y to int, so reset here + source->peak->dx = 0.0; + source->peak->dy = 0.0; + + source->extSN = SN; + + source->moments = pmMomentsAlloc (); + source->moments->Mx = X; + source->moments->My = Y; + source->moments->Mrf = kronRadius * 0.4; // kronRadius is 2.5 * first radial moment + + // Don't mark the moments as measured because that causes many fields to be left blank. + // The moments code knows not to change the position or the Mrf for external sources + // source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED; + + if (isfinite(petRadius)) { + source->extpars = pmSourceExtendedParsAlloc (); + source->extpars->petrosianRadius = petRadius; + } + + } + bool saveExtModelParams = false; if (fitGalaxy && galaxyModelType >= 0) { - source->modelFits = psArrayAllocEmpty (1); - pmModel *model = pmModelAlloc(galaxyModelType); - psF32 *xPAR = model->params->data.F32; - - xPAR[PM_PAR_SKY] = 0.0; - xPAR[PM_PAR_I0] = 1.0; - xPAR[PM_PAR_XPOS] = X; + // skip likely stars + if (magDiff < starCut || SN < SNMinForCut) { + if (forceAll) { + saveExtModelParams = true; + } else if (galaxyModelType == modelTypeForce) { + // This is the model type that we are looking for + // proceed + saveExtModelParams = true; + } else if (modelTypeForce == sersicModelType && galaxyModelType == devModelType) { + // We're doing sersic models, if sersic index is greater than the recipe's minimum index + // do dev model as well + if (isfinite(forceDevSersicMin) || Sindex >= forceDevSersicMin) { + saveExtModelParams = true; + } + } else { + // not interested in this model + } + } + } + + if (saveExtModelParams) { + if (!source->modelFits) { + source->modelFits = psArrayAllocEmpty (1); + } + pmModel *model = pmModelAlloc(galaxyModelType); + psF32 *xPAR = model->params->data.F32; + + xPAR[PM_PAR_SKY] = 0.0; + xPAR[PM_PAR_I0] = 1.0; + xPAR[PM_PAR_XPOS] = X; xPAR[PM_PAR_YPOS] = Y; @@ -205,29 +273,24 @@ galaxyAxes.theta = theta * PS_RAD_DEG; - pmPSF_AxesToModel (xPAR, galaxyAxes, galaxyModelType); + pmPSF_AxesToModel (xPAR, galaxyAxes, model->class->useReff); if (model->params->n > 7) { xPAR[PM_PAR_7] = 0.5 / Sindex; } + model->chisq = chisq; + model->nDOF = nDOF; + model->flags = modelFlags; + psArrayAdd (source->modelFits, 1, model); -#ifdef notyet - // XXX: set source->modelEXT to this model and flag as extended so that we have a better - // shot at subtracting extended sources? - - // This doesn't work right. Need to do some more work to flesh out the model before it can be - // used in subtaction. Better idea might be to do 2 passes in psphotFullForceReadout to - // First find the best extended models and then do everything else. - source->modelEXT = psMemIncrRefCounter(model); - // is this safe? (no see above) - source->type = PM_SOURCE_TYPE_EXTENDED; - source->mode |= PM_SOURCE_MODE_EXTMODEL; -#endif - psFree (model); } - sources->data[i] = source; psFree(row); + } + if (source) { + // close out last source + psArrayAdd (sources, 1, source); + psFree(source); } @@ -250,12 +313,30 @@ pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC"); + pmModelType devModelType = pmModelClassGetType("PS_MODEL_DEV"); + pmModelType selectedModelType = -1; + bool chooseBest = false; + bool chooseAll = false; psString modelToChoose = psMetadataLookupStr(&mdok, recipe, "EXT_MODEL_TYPE_FOR_CFF"); - pmModelType selectedModelType = -1; + if (mdok && modelToChoose != NULL) { - if (strcmp(modelToChoose, "BEST")) { + if (!strcmp(modelToChoose, "BEST")) { + chooseBest = true; + } else if (!strcmp(modelToChoose, "ALL")) { + chooseAll = true; + } else if (strcmp(modelToChoose, "PS_MODEL_SERSIC")) { + // We have selected a model type other than Sersic. + // Save it's type for use below. Sersic is handled specially selectedModelType = pmModelClassGetType(modelToChoose); } } + + // minimum sersic index to force devModel + psF32 sersicMinDev = psMetadataLookupF32(&mdok, recipe, "EXT_MODEL_FORCE_DEV_SERSIC_MIN"); + if (!mdok) { + sersicMinDev = NAN; + } + + sources = psArraySort (sources, pmSourceSortBySeq); for (int i = 0; i < sources->n; i++) { @@ -263,88 +344,195 @@ pmSource *source = thisSource->parent ? thisSource->parent : thisSource; - psF32 xPos, yPos, flux, rMajor, rMinor, theta; - psS32 modelType = 0; + #define MAX_ROWS_PER_SRC 10 + psF32 xPos[MAX_ROWS_PER_SRC], yPos[MAX_ROWS_PER_SRC], flux[MAX_ROWS_PER_SRC]; + psF32 rMajor[MAX_ROWS_PER_SRC], rMinor[MAX_ROWS_PER_SRC], theta[MAX_ROWS_PER_SRC]; + psF32 chisq[MAX_ROWS_PER_SRC], nDOF[MAX_ROWS_PER_SRC]; + psS32 modelFlags[MAX_ROWS_PER_SRC]; + psS32 modelType[MAX_ROWS_PER_SRC]; + psF32 sersicIndex = NAN; bool fitGalaxy = false; bool psfStar = (source->mode & PM_SOURCE_MODE_PSFSTAR) ? true : false; - psF32 sersicIndex = 0; - // For now only perform galaxy fits on extended objects - if (source->modelEXT == NULL) { - pmModel *model = source->modelPSF; - if (model == NULL) continue; - psF32 *PAR = model->params->data.F32; - if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY]) || - !isfinite(source->psfFlux)) { - continue; - } - - xPos = model->params->data.F32[PM_PAR_XPOS]; - yPos = model->params->data.F32[PM_PAR_YPOS]; - flux = source->psfFlux; - rMajor = 0; - rMinor = 0; - theta = 0; - } else { - // Find the model with the selected type. If selected type is -1 choose the one selected ad - // modelEXT which was the best - int iModel = -1; - if (source->modelEXT) { - pmModelType ext_model_type = selectedModelType != -1 ? selectedModelType : source->modelEXT->type; + int n_rows = 0; + + psF32 kronFlux = source->moments->KronFlux; + psF32 SN = NAN; + psF32 magDiff = NAN; + if (isfinite(kronFlux) && isfinite(source->moments->KronFluxErr) && isfinite(source->psfMag)) { + SN = kronFlux/source->moments->KronFluxErr; + // kronMag - psfMag for use as star/glaxy separator + magDiff = -2.5 * log10(kronFlux) - source->psfMag ; + } + + // start with psf model + pmModel *model = source->modelPSF; + if (model == NULL) continue; + psF32 *PAR = model->params->data.F32; + if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY]) || + !isfinite(source->psfFlux)) { + continue; + } + + // save the PSF model parameters for each object as the first entry + xPos[0] = model->params->data.F32[PM_PAR_XPOS]; + yPos[0] = model->params->data.F32[PM_PAR_YPOS]; + flux[0] = source->psfFlux; + rMajor[0] = 0; + rMinor[0] = 0; + theta[0] = 0; + modelType[0] = -1; + sersicIndex = NAN; + chisq[0] = NAN; + nDOF[0] = NAN; + modelFlags[0] = 0; + n_rows ++; + + if (source->modelFits != NULL) { + // figure out which models to use based on recipe paramters + if (chooseAll) { + // Save parameters for all valid extended models + + // but make sure we aren't going to overflow our arrays + assert (source->modelFits->n < MAX_ROWS_PER_SRC); + for (int j=0; jmodelFits->n; j++) { - pmModel *aModel = source->modelFits->data[j]; - if (aModel->type == ext_model_type) { - iModel = j; - break; + pmModel *model = source->modelFits->data[j]; + psF32 *PAR = model->params->data.F32; + + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SYY]) && isfinite(PAR[PM_PAR_SXY]) && + isfinite(model->mag)) { + + xPos[n_rows] = PAR[PM_PAR_XPOS]; + yPos[n_rows] = PAR[PM_PAR_YPOS]; + + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->class->useReff); + rMajor[n_rows] = axes.major; + rMinor[n_rows] = axes.minor; + theta[n_rows] = axes.theta*PS_DEG_RAD; + flux[n_rows] = pow(10.0, -0.4*model->mag); + modelType[n_rows] = model->type; + chisq[n_rows] = model->chisq; + nDOF[n_rows] = model->nDOF; + modelFlags[n_rows] = model->flags; + fitGalaxy = true; + if (model->params->n == 8) { + // this will save the sersic index in all models but that might + // be useful + sersicIndex = 0.5 / PAR[PM_PAR_7]; + } else { + sersicIndex = NAN; + } + + n_rows++; + } + } + } else { + int jModelSersic = -1; + int jModelDev = -1; + int jModelSelected = -1; + psF32 minChisq = NAN; + for (int j=0; jmodelFits->n; j++) { + pmModel *model = source->modelFits->data[j]; + if (chooseBest) { + // choose the model with lowest chisq + if (isfinite(model->chisq) && (!isfinite(minChisq) || model->chisq < minChisq)) { + jModelSelected = j; + minChisq = model->chisq; + } + } else { + // find the index of models of interest + if (model->type == selectedModelType) { + jModelSelected = j; + } else if (model->type == sersicModelType) { + jModelSersic = j; + } else if (model->type == devModelType) { + jModelDev = j; + } + } + } + if (jModelSelected >= 0 || jModelSersic >= 0) { + // If a specific non-sersic model we take paramers from that one. + // Otherwise we do the sersic model. + pmModel *model = jModelSelected >= 0 ? source->modelFits->data[jModelSelected] : + source->modelFits->data[jModelSersic]; + psF32 *PAR = model->params->data.F32; + + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SYY]) && isfinite(PAR[PM_PAR_SXY]) && + isfinite(model->mag)) { + + xPos[0] = PAR[PM_PAR_XPOS]; + yPos[0] = PAR[PM_PAR_YPOS]; + + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->class->useReff); + rMajor[0] = axes.major; + rMinor[0] = axes.minor; + theta[0] = axes.theta*PS_DEG_RAD; + flux[0] = pow(10.0, -0.4*model->mag); + modelType[0] = model->type; + chisq[0] = model->chisq; + nDOF[0] = model->nDOF; + modelFlags[0] = model->flags; + fitGalaxy = true; + if (model->type == sersicModelType) { + PS_ASSERT_FLOAT_LARGER_THAN(PAR[PM_PAR_7], 0.0, false); + sersicIndex = 0.5 / PAR[PM_PAR_7]; + } + + n_rows = 1; + + // Unless a specific non-sersic model type was selected do dev model for sources with + // sersic index above the recipe limit. + if (jModelSelected == -1 && jModelDev >= 0 && isfinite(sersicMinDev) && + sersicIndex > sersicMinDev) { + + model = source->modelFits->data[jModelDev]; + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SYY]) && isfinite(PAR[PM_PAR_SXY]) && + isfinite(model->mag)) { + + xPos[1] = PAR[PM_PAR_XPOS]; + yPos[1] = PAR[PM_PAR_YPOS]; + + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->class->useReff); + rMajor[1] = axes.major; + rMinor[1] = axes.minor; + theta[1] = axes.theta*PS_DEG_RAD; + flux[1] = pow(10.0, -0.4*model->mag); + modelType[1] = model->type; + sersicIndex = NAN; + chisq[1] = model->chisq; + nDOF[1] = model->nDOF; + modelFlags[1] = model->flags; + n_rows = 2; + } + } } } } - if (iModel == -1) { - // Can this happen? perhaps if model type is qgauss (extended but below s/n for model fits? - // XXX: Should this be an assert? - continue; - } - pmModel *model = source->modelFits->data[iModel]; - psF32 *PAR = model->params->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) { - // XXX: do we need to handle this type - continue; - } else { - if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY]) || - !isfinite(model->mag)) { - // bad model - continue; - } - psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type); - rMajor = axes.major; - rMinor = axes.minor; - theta = axes.theta*PS_DEG_RAD; - flux = pow(10.0, -0.4*model->mag); - fitGalaxy = true; - } - modelType = model->type; - if (modelType == sersicModelType) { - PS_ASSERT_FLOAT_LARGER_THAN(PAR[PM_PAR_7], 0.0, false); - sersicIndex = 0.5 / PAR[PM_PAR_7]; - } - } - psMetadata *row = psMetadataAlloc(); - psMetadataAddU32 (row, PS_LIST_TAIL, "ID", 0, "IPP detection identifier", source->seq); - psMetadataAddF32 (row, PS_LIST_TAIL, "X", 0, "x coordinate", xPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y", 0, "y coordinate", yPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "FLUX", 0, "flux per second", flux/exptime); - psMetadataAddF32 (row, PS_LIST_TAIL, "AP_RADIUS", 0, "aperture radius", source->apRadius); - psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_RADIUS", 0, "Kron radius", source->moments->Mrf * 2.5); - psMetadataAddF32 (row, PS_LIST_TAIL, "PETRO_RADIUS", 0, "Petrosian Radius", source->extpars ? source->extpars->petrosianRadius : NAN); - psMetadataAddBool (row, PS_LIST_TAIL, "FIT_GALAXY", 0, "source has xfit", fitGalaxy); - psMetadataAddBool (row, PS_LIST_TAIL, "PSF_STAR", 0, "source was psf star", psfStar); - psMetadataAddF32 (row, PS_LIST_TAIL, "R_MAJOR", 0, "radius of major axis", rMajor); - psMetadataAddF32 (row, PS_LIST_TAIL, "R_MINOR", 0, "radius of minor axis", rMinor); - psMetadataAddF32 (row, PS_LIST_TAIL, "THETA", 0, "theta", theta); - psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "model type", modelType); - psMetadataAddF32 (row, PS_LIST_TAIL, "INDEX", 0, "sersic index", sersicIndex); - - psArrayAdd(table, 100, row); - psFree(row); + } + + for (int j = 0; j < n_rows; j++) { + psMetadata *row = psMetadataAlloc(); + psMetadataAddU32 (row, PS_LIST_TAIL, "ID", 0, "IPP detection identifier", source->seq); + psMetadataAddF32 (row, PS_LIST_TAIL, "X", 0, "x coordinate", xPos[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y", 0, "y coordinate", yPos[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "FLUX", 0, "flux per second", flux[j]/exptime); + psMetadataAddF32 (row, PS_LIST_TAIL, "SN", 0, "kron flux signal to noise", SN); + psMetadataAddF32 (row, PS_LIST_TAIL, "MAG_DIFF", 0, "psf mag - kron mag", magDiff); + psMetadataAddF32 (row, PS_LIST_TAIL, "AP_RADIUS", 0, "aperture radius", source->apRadius); + psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_RADIUS", 0, "Kron radius", source->moments->Mrf * 2.5); + psMetadataAddF32 (row, PS_LIST_TAIL, "PETRO_RADIUS", 0, "Petrosian Radius", source->extpars ? source->extpars->petrosianRadius : NAN); + psMetadataAddBool (row, PS_LIST_TAIL, "FIT_GALAXY", 0, "source has xfit", fitGalaxy); + psMetadataAddBool (row, PS_LIST_TAIL, "PSF_STAR", 0, "source was psf star", psfStar); + psMetadataAddF32 (row, PS_LIST_TAIL, "R_MAJOR", 0, "radius of major axis", rMajor[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "R_MINOR", 0, "radius of minor axis", rMinor[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "THETA", 0, "theta", theta[j]); + psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "model type", modelType[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "INDEX", 0, "sersic index", sersicIndex); + psMetadataAddF32 (row, PS_LIST_TAIL, "CHISQ", 0, "chisq", chisq[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "NDOF", 0, "n degrees of freedom", nDOF[j]); + psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_FLAGS", 0, "model flags", modelFlags[j]); + + psArrayAdd(table, 100, row); + psFree(row); + } } Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_CMF.c.in (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_CMF.c.in (revision 37403) @@ -37,7 +37,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -142,7 +142,8 @@ @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in reported aperture", source->apMagRaw); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", source->apRadius); @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", source->apFlux); @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", source->apFluxErr); + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "AP_NPIX", PS_DATA_S32, "aperture unmasked pixels", source->apNpixels); @extNsigma); + // PSF shape parameters: @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_CORE", PS_DATA_F32, "k term if defined", outputs.psfCore); + + // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_FWHM_MAJ", PS_DATA_F32, "PSF FWHM (major axis)", outputs.psfMajorFWHM); + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_FWHM_MIN", PS_DATA_F32, "PSF FWHM (minor axis)", outputs.psfMinorFWHM); + + // psf data quality @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", source->pixWeightNotPoor); @@ -197,5 +206,5 @@ } - if (source->lensingOBJ && source->lensingPSF->smear) { + if (source->lensingPSF && source->lensingPSF->smear) { @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X11); @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X12); @@ -205,5 +214,5 @@ } - if (source->lensingOBJ && source->lensingPSF->shear) { + if (source->lensingPSF && source->lensingPSF->shear) { @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X11); @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X12); @@ -212,4 +221,12 @@ @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->e2); } + + // if lensing params exist also include the backmapped chipID and chip coordinates + if (source->lensingPSF && source->lensingPSF->shear) { + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM", PS_DATA_S16, "id of warp input chip", source->chipNum); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X", PS_DATA_S16, "x coord in warp input chip", source->chipX); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_Y", PS_DATA_S16, "y coord in warp input chip", source->chipY); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING3", PS_DATA_S16, "more padding", 0); + } @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); @@ -302,12 +319,16 @@ // define PSF model type - // XXX need to carry the extra model parameters int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); - char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); + // if header does not define the model, default to a gaussian + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSFMODEL"); if (PSF_NAME != NULL) { modelType = pmModelClassGetType (PSF_NAME); } assert (modelType > -1); + + // do we expect to find lensing parameters? + bool haveLensOBJ = psMetadataLookupBool (&status, header, "LENS_OBJ"); + bool haveLensPSF = psMetadataLookupBool (&status, header, "LENS_PSF"); // We get the size of the table, and allocate the array of sources first because the table @@ -344,4 +365,8 @@ @ALL@ axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); @ALL@ axes.theta = axes.theta * PS_RAD_DEG; + + @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) { + @>PS1_V4,>PS1_SV2,>PS1_DV3@ PAR[PM_PAR_7] = psMetadataLookupF32 (&status, row, "PSF_CORE"); + @>PS1_V4,>PS1_SV2,>PS1_DV3@ } @ALL@ PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); @@ -355,4 +380,6 @@ @ALL@ source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); @>PS1_V2,PS1_SV?,>PS1_DV1@ source->apMagRaw = psMetadataLookupF32 (&status, row, "AP_MAG_RAW"); + @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux = psMetadataLookupF32 (&status, row, "AP_FLUX"); + @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr = psMetadataLookupF32 (&status, row, "AP_FLUX_SIG"); // XXX use these to determine PAR[PM_PAR_I0] if they exist? @@ -365,5 +392,5 @@ @ALL@ dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - pmPSF_AxesToModel (PAR, axes, modelType); + pmPSF_AxesToModel (PAR, axes, model->class->useReff); @ALL@ float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); @@ -383,5 +410,6 @@ @ALL@ source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); @ALL@ source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); - @ALL@ source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); + @ALL@ source->apRadius = psMetadataLookupF32 (&status, row, "AP_MAG_RADIUS"); + @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels = psMetadataLookupS32 (&status, row, "AP_NPIX"); // note that some older versions used PSF_PROBABILITY: this was not well defined. @@ -410,4 +438,39 @@ @>PS1_V2,PS1_SV?@ source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S"); @>PS1_V2,PS1_SV?@ source->moments->Myyyy = 0.0; + + // Lensing parameters (on read if PS1_V5+) + if (haveLensOBJ) { + source->lensingOBJ = pmSourceLensingAlloc (); + source->lensingOBJ->smear = pmLensingParsAlloc(); + source->lensingOBJ->shear = pmLensingParsAlloc(); + + @>PS1_V4@ source->lensingOBJ->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SH_OBJ"); + } + + if (haveLensPSF) { + source->lensingPSF = pmSourceLensingAlloc (); + source->lensingPSF->smear = pmLensingParsAlloc(); + source->lensingPSF->shear = pmLensingParsAlloc(); + + @>PS1_V4@ source->lensingPSF->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SM_PSF"); + @>PS1_V4@ source->lensingPSF->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SH_PSF"); + } @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrf = psMetadataLookupF32 (&status, row, "MOMENTS_R1"); @@ -482,6 +545,8 @@ } +#ifdef SORT_EXTENSIONS_BY_FLUX // let's write these out in S/N order sources = psArraySort (sources, pmSourceSortByFlux); +#endif table = psArrayAllocEmpty (sources->n); @@ -563,4 +628,5 @@ // 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 + // XXX NOTE EAM 20140806 : PETRO_MAG_ERR is inverted!! oops!! 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); @@ -821,6 +887,8 @@ psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); +#ifdef SORT_EXTENSIONS_BY_FLUX // let's write these out in S/N order sources = psArraySort (sources, pmSourceSortByFlux); +#endif float magOffset; @@ -847,5 +915,5 @@ } - @>PS1_DV2@ pmChip *chip = readout->parent->parent; + @>PS1_DV2,>PS1_SV3@ pmChip *chip = readout->parent->parent; pmModelStatus badModel = PM_MODEL_STATUS_NONE; @@ -906,12 +974,12 @@ } - @>PS1_DV2@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; - @>PS1_DV2@ float posAngle = 0.0; - @>PS1_DV2@ float pltScale = 0.0; - @>PS1_DV2@ pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); - @>PS1_DV2@ double raPos = ptSky.r*PS_DEG_RAD; - @>PS1_DV2@ double decPos = ptSky.d*PS_DEG_RAD; - @>PS1_DV2@ posAngle *= PS_DEG_RAD; - @>PS1_DV2@ pltScale *= PS_DEG_RAD*3600.0; + @>PS1_DV2,>PS1_SV3@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; + @>PS1_DV2,>PS1_SV3@ float posAngle = 0.0; + @>PS1_DV2,>PS1_SV3@ float pltScale = 0.0; + @>PS1_DV2,>PS1_SV3@ pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + @>PS1_DV2,>PS1_SV3@ double raPos = ptSky.r*PS_DEG_RAD; + @>PS1_DV2,>PS1_SV3@ double decPos = ptSky.d*PS_DEG_RAD; + @>PS1_DV2,>PS1_SV3@ posAngle *= PS_DEG_RAD; + @>PS1_DV2,>PS1_SV3@ pltScale *= PS_DEG_RAD*3600.0; float kronFlux = source->moments ? source->moments->KronFlux : NAN; @@ -927,6 +995,6 @@ psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); - @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "RA_EXT", 0, "EXT model ra coordinate", raPos); - @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "DEC_EXT", 0, "EXT model dec coordinate", decPos); + @>PS1_DV2,>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "RA_EXT", 0, "EXT model ra coordinate", raPos); + @>PS1_DV2,>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "DEC_EXT", 0, "EXT model dec coordinate", decPos); @>PS1_DV2@ float instFlux = isfinite(model->mag) ? pow(10.0, -0.4*model->mag) : NAN; @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_FLUX", 0, "EXT fit instrumental counts", instFlux); @@ -942,5 +1010,5 @@ // EAM : adding for PV2 outputs: - @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", source->modelEXT ? source->modelEXT->flags : 0); + @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", model->flags); @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE", 0, "position angle at source (degrees)", posAngle); @@ -978,5 +1046,5 @@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]); } else { - psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type); + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->class->useReff); psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (major axis), length for trail", axes.major); psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (minor axis), sigma for trail", axes.minor); @@ -1109,4 +1177,5 @@ model->chisq = psMetadataLookupF32(&status, row, "EXT_CHISQ"); model->nDOF = psMetadataLookupF32(&status, row, "EXT_NDOF"); + model->flags = psMetadataLookupS16(&status, row, "EXT_FLAGS"); // EXT_MODEL_TYPE gives the model chosen by psphot as the best. @@ -1127,5 +1196,5 @@ axes.minor = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN"); axes.theta = psMetadataLookupF32(&status, row, "EXT_THETA"); - if (!pmPSF_AxesToModel(PAR, axes, modelType)) { + if (!pmPSF_AxesToModel(PAR, axes, model->class->useReff)) { // Do we need to fail here or can this happen? psError(PS_ERR_UNKNOWN, false, "Failed to convert psf axes to model"); @@ -1225,6 +1294,8 @@ psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES"); +#ifdef SORT_EXTENSIONS_BY_FLUX // let's write these out in S/N order sources = psArraySort (sources, pmSourceSortByFlux); +#endif table = psArrayAllocEmpty (sources->n); @@ -1439,8 +1510,15 @@ psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "galaxy table extension", extname); - psMetadataAddStr (outhead, PS_LIST_TAIL, "HI", PS_META_REPLACE, "does this get through?", "THERE"); - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); + psF32 Q = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_Q"); + psMetadataAddF32 (outhead, PS_LIST_TAIL, "GALAXY_SHAPES_Q", PS_META_REPLACE, "", Q); + + psF32 NSigma = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_NSIGMA"); + psMetadataAddF32 (outhead, PS_LIST_TAIL, "GALAXY_SHAPES_NSIGMA", PS_META_REPLACE, "", NSigma); + + psF32 clampSN = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_CLAMP_SN"); + psMetadataAddF32 (outhead, PS_LIST_TAIL, "GALAXY_SHAPES_CLAMP_SN", PS_META_REPLACE, "", clampSN); + + // They are probably already in this order but ... + sources = psArraySort (sources, pmSourceSortBySeq); psArray *table = psArrayAllocEmpty (sources->n); @@ -1459,28 +1537,48 @@ if (source->galaxyFits == NULL) continue; - pmModel *model = source->modelFits->data[0]; - if (!model) return false; - - // X,Y coordinates are stored with the model parameters - psF32 *PAR = model->params->data.F32; - - psMetadata *row = psMetadataAlloc (); - - // we write out the x,y positions so people can link to the psf either way (position or ID) - psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_FIT", 0, "model x coordinate", PAR[PM_PAR_XPOS]); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_FIT", 0, "model y coordinate", PAR[PM_PAR_YPOS]); - psMetadataAddF32 (row, PS_LIST_TAIL, "NPIX", 0, "number of pixels for fits", source->galaxyFits->nPix); - - psVector *Flux = source->galaxyFits->Flux; - psVector *dFlux = source->galaxyFits->dFlux; - psVector *chisq = source->galaxyFits->chisq; - - psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX", PS_META_REPLACE, "normalization for galaxy flux", Flux); - psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX_ERR", PS_META_REPLACE, "error on normalization", dFlux); - psMetadataAddVector (row, PS_LIST_TAIL, "GAL_CHISQ", PS_META_REPLACE, "galaxy fit chisq", chisq); - - psArrayAdd (table, 100, row); - psFree (row); + for (int iModel = 0; iModel < source->modelFits->n; iModel++) { + pmModel *model = source->modelFits->data[iModel]; + if (!model) continue; + + pmSourceGalaxyFits *galaxyFits = NULL; + for (int iFit = 0; iFit < source->galaxyFits->n; iFit++) { + galaxyFits = source->galaxyFits->data[iFit]; + if (model->type == galaxyFits->modelType) break; + galaxyFits = NULL; + } + + if (!galaxyFits) continue; + + // X,Y coordinates are stored with the model parameters + psF32 *PAR = model->params->data.F32; + + psMetadata *row = psMetadataAlloc (); + + // we write out the x,y positions so people can link to the psf either way (position or ID) + psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); + psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "model type", galaxyFits->modelType); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_FIT", 0, "model x coordinate", PAR[PM_PAR_XPOS]); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_FIT", 0, "model y coordinate", PAR[PM_PAR_YPOS]); + psMetadataAddF32 (row, PS_LIST_TAIL, "NPIX", 0, "number of pixels for fits", galaxyFits->nPix); + // psMetadataAddS32 (row, PS_LIST_TAIL, "REDUCED_TRIALS", 0, "source has reduced number of trials", galaxyFits->reducedTrials); + + psVector *Flux = galaxyFits->Flux; + psVector *dFlux = galaxyFits->dFlux; + psVector *chisq = galaxyFits->chisq; + + psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX", PS_META_REPLACE, "normalization for galaxy flux", Flux); + psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX_ERR", PS_META_REPLACE, "error on normalization", dFlux); + psMetadataAddVector (row, PS_LIST_TAIL, "GAL_CHISQ", PS_META_REPLACE, "galaxy fit chisq", chisq); + + psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MAJOR_MIN", 0, "fractional major axis min", galaxyFits->fRmajorMin); + psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MAJOR_MAX", 0, "fractional major axis max", galaxyFits->fRmajorMax); + psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MAJOR_DEL", 0, "fractional major axis max", galaxyFits->fRmajorDel); + psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MINOR_MIN", 0, "fractional minor axis min", galaxyFits->fRminorMin); + psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MINOR_MAX", 0, "fractional minor axis max", galaxyFits->fRminorMax); + psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MINOR_DEL", 0, "fractional minor axis max", galaxyFits->fRminorDel); + + psArrayAdd (table, 100, row); + psFree (row); + } } @@ -1520,4 +1618,20 @@ return false; } + if (!readout->analysis) { + readout->analysis = psMetadataAlloc(); + } + psF32 Q = psMetadataLookupF32(&status, tableHeader, "GALAXY_SHAPES_Q"); + // XXX: turn this into an assert once we're done + if (status) { + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "GALAXY_SHAPES_Q", PS_META_REPLACE, "", Q); + + psF32 NSigma = psMetadataLookupF32(&status, tableHeader, "GALAXY_SHAPES_NSIGMA"); + psAssert(status, "missing GALAXY_SHAPES_NSIGMA"); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "GALAXY_SHAPES_NSIGMA", PS_META_REPLACE, "", NSigma); + + psF32 clampSN = psMetadataLookupF32(&status, tableHeader, "GALAXY_SHAPES_CLAMP_SN"); + psAssert(status, "missing GALAXY_SHAPES_CLAMP_SN"); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "GALAXY_SHAPES_CLAMP_SN", PS_META_REPLACE, "", clampSN); + } for (long i = 0; i < numSources; i++) { @@ -1541,4 +1655,5 @@ } + int modelType = psMetadataLookupS32 (&status, row, "MODEL_TYPE"); psVector *Flux = psMetadataLookupVector(&status, row, "GAL_FLUX"); psVector *dFlux = psMetadataLookupVector(&status, row, "GAL_FLUX_ERR"); @@ -1546,14 +1661,29 @@ if (Flux && Flux->n > 0) { - psFree(source->galaxyFits); - source->galaxyFits = pmSourceGalaxyFitsAlloc(); - source->galaxyFits->nPix = psMetadataLookupF32(&status, row, "NPIX"); - - psFree(source->galaxyFits->Flux); - source->galaxyFits->Flux = psMemIncrRefCounter(Flux); - psFree(source->galaxyFits->dFlux); - source->galaxyFits->dFlux = psMemIncrRefCounter(dFlux); - psFree(source->galaxyFits->chisq); - source->galaxyFits->chisq = psMemIncrRefCounter(chisq); + if (!source->galaxyFits) { + source->galaxyFits = psArrayAllocEmpty(1); + } + + pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc(); + + psArrayAdd(source->galaxyFits, 1, galaxyFits); + + psFree(galaxyFits); + galaxyFits->modelType = modelType; + galaxyFits->nPix = psMetadataLookupF32(&status, row, "NPIX"); + + galaxyFits->fRmajorMin = psMetadataLookupF32(&status, row, "FR_MAJOR_MIN"); + galaxyFits->fRmajorMax = psMetadataLookupF32(&status, row, "FR_MAJOR_MAX"); + galaxyFits->fRmajorDel = psMetadataLookupF32(&status, row, "FR_MAJOR_DEL"); + galaxyFits->fRminorMin = psMetadataLookupF32(&status, row, "FR_MINOR_MIN"); + galaxyFits->fRminorMax = psMetadataLookupF32(&status, row, "FR_MINOR_MAX"); + galaxyFits->fRminorDel = psMetadataLookupF32(&status, row, "FR_MINOR_DEL"); + + psFree(galaxyFits->Flux); + galaxyFits->Flux = psMemIncrRefCounter(Flux); + psFree(galaxyFits->dFlux); + galaxyFits->dFlux = psMemIncrRefCounter(dFlux); + psFree(galaxyFits->chisq); + galaxyFits->chisq = psMemIncrRefCounter(chisq); } Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_CMP.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_CMP.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_CMP.c (revision 37403) @@ -37,7 +37,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -136,5 +136,5 @@ lsky = (source->sky < 1.0) ? 0.0 : log10(source->sky); - axes = pmPSF_ModelToAxes (PAR, model->type); + axes = pmPSF_ModelToAxes (PAR, model->class->useReff); float psfMagErr = isfinite(source->psfMagErr) ? source->psfMagErr : 999; @@ -293,5 +293,5 @@ goto skip_source; - pmPSF_AxesToModel (PAR, axes, modelType); + pmPSF_AxesToModel (PAR, axes, source->modelPSF->class->useReff); psArrayAdd (sources, 100, source); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_MatchedRefs.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_MatchedRefs.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_MatchedRefs.c (revision 37403) @@ -37,7 +37,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -99,22 +99,11 @@ // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS.SUBSET"); if (rawstars == NULL) continue; // select the raw objects for this readout - psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET"); if (refstars == NULL) continue; psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n); - -# if (0) - // XXX test - FILE *outfile = fopen ("refstars.dat", "w"); - assert (outfile); - for (int nn = 0; nn < refstars->n; nn++) { - pmAstromObj *ref = refstars->data[nn]; - fprintf (outfile, "%lf %lf\n", ref->sky->r*PS_DEG_RAD, ref->sky->d*PS_DEG_RAD); - } - fclose (outfile); -# endif psArray *matches = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH"); @@ -128,16 +117,18 @@ psMetadata *row = psMetadataAlloc (); - psMetadataAdd (row, PS_LIST_TAIL, "RA", PS_DATA_F64, "right ascension (deg, J2000)", PM_DEG_RAD*ref->sky->r); - psMetadataAdd (row, PS_LIST_TAIL, "DEC", PS_DATA_F64, "declination (deg, J2000)", PM_DEG_RAD*ref->sky->d); - psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP", PS_DATA_F32, "x coord on chip", raw->chip->x); - psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP", PS_DATA_F32, "y coord on chip", raw->chip->y); - psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP_FIT",PS_DATA_F32, "x fitted coord on chip", ref->chip->x); - psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP_FIT",PS_DATA_F32, "y fitted coord on chip", ref->chip->y); - psMetadataAdd (row, PS_LIST_TAIL, "X_FPA", PS_DATA_F32, "x coord on focal plane", raw->FP->x); - psMetadataAdd (row, PS_LIST_TAIL, "Y_FPA", PS_DATA_F32, "y coord on focal plane", raw->FP->y); - psMetadataAdd (row, PS_LIST_TAIL, "MAG_INST", PS_DATA_F32, "instrumental magnitude", raw->Mag); - psMetadataAdd (row, PS_LIST_TAIL, "MAG_REF", PS_DATA_F32, "reference star magnitude", ref->Mag); - psMetadataAdd (row, PS_LIST_TAIL, "COLOR_REF",PS_DATA_F32, "reference star color", ref->Color); - psMetadataAdd (row, PS_LIST_TAIL, "CHIP_ID", PS_DATA_STRING, "chip identifier", chipName); + psMetadataAdd (row, PS_LIST_TAIL, "RA_REF", PS_DATA_F64, "right ascension (deg, J2000)", PM_DEG_RAD*ref->sky->r); + psMetadataAdd (row, PS_LIST_TAIL, "DEC_REF", PS_DATA_F64, "declination (deg, J2000)", PM_DEG_RAD*ref->sky->d); + psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP_REF", PS_DATA_F32, "x fitted coord on chip", ref->chip->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP_REF", PS_DATA_F32, "y fitted coord on chip", ref->chip->y); + psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP_RAW", PS_DATA_F32, "x coord on chip", raw->chip->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP_RAW", PS_DATA_F32, "y coord on chip", raw->chip->y); + psMetadataAdd (row, PS_LIST_TAIL, "X_FPA_RAW", PS_DATA_F32, "x coord on focal plane", raw->FP->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_FPA_RAW", PS_DATA_F32, "y coord on focal plane", raw->FP->y); + psMetadataAdd (row, PS_LIST_TAIL, "X_TPA_RAW", PS_DATA_F32, "x coord on focal plane", raw->TP->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_TPA_RAW", PS_DATA_F32, "y coord on focal plane", raw->TP->y); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_INST", PS_DATA_F32, "instrumental magnitude", raw->Mag); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_REF", PS_DATA_F32, "reference star magnitude", ref->Mag); + psMetadataAdd (row, PS_LIST_TAIL, "COLOR_REF", PS_DATA_F32, "reference star color", ref->Color); + psMetadataAdd (row, PS_LIST_TAIL, "CHIP_ID", PS_DATA_STRING, "chip identifier", chipName); // XXX need to add the reference color, but this needs getstar / dvo.photcodes for the reference to be refined. Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_OBJ.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_OBJ.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_OBJ.c (revision 37403) @@ -37,7 +37,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -93,5 +93,5 @@ } - axes = pmPSF_ModelToAxes (PAR, model->type); + axes = pmPSF_ModelToAxes (PAR, model->class->useReff); psLineInit (line); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_PS1_CAL_0.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 37403) @@ -37,7 +37,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -115,5 +115,5 @@ yErr = dPAR[PM_PAR_YPOS]; if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, model->type); + axes = pmPSF_ModelToAxes (PAR, model->class->useReff); } else { axes.major = NAN; @@ -290,5 +290,5 @@ dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - pmPSF_AxesToModel (PAR, axes, modelType); + pmPSF_AxesToModel (PAR, axes, model->class->useReff); float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); @@ -624,5 +624,5 @@ yErr = dPAR[PM_PAR_YPOS]; - axes = pmPSF_ModelToAxes (PAR, model->type); + axes = pmPSF_ModelToAxes (PAR, model->class->useReff); // generate RA,DEC Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_PS1_DEV_0.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 37403) @@ -37,7 +37,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -91,5 +91,5 @@ yErr = dPAR[PM_PAR_YPOS]; - axes = pmPSF_ModelToAxes (PAR, model->type); + axes = pmPSF_ModelToAxes (PAR, model->class->useReff); } else { // XXX: This code seg faults if source->peak is NULL. @@ -216,5 +216,5 @@ source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); - pmPSF_AxesToModel (PAR, axes, modelType); + pmPSF_AxesToModel (PAR, axes, model->class->useReff); float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 37403) @@ -37,7 +37,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -97,5 +97,5 @@ yErr = dPAR[PM_PAR_YPOS]; if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, model->type); + axes = pmPSF_ModelToAxes (PAR, model->class->useReff); } else { axes.major = NAN; @@ -259,5 +259,5 @@ dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - pmPSF_AxesToModel (PAR, axes, modelType); + pmPSF_AxesToModel (PAR, axes, model->class->useReff); float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); @@ -524,5 +524,5 @@ yErr = dPAR[PM_PAR_YPOS]; - axes = pmPSF_ModelToAxes (PAR, model->type); + axes = pmPSF_ModelToAxes (PAR, model->class->useReff); row = psMetadataAlloc (); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_RAW.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_RAW.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_RAW.c (revision 37403) @@ -37,7 +37,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_SMPDATA.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 37403) @@ -37,7 +37,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -93,5 +93,5 @@ lsky = (source->sky < 1.0) ? 0.0 : log10(source->sky); - axes = pmPSF_ModelToAxes (PAR, model->type); + axes = pmPSF_ModelToAxes (PAR, model->class->useReff); } else { @@ -190,5 +190,5 @@ axes.theta = psMetadataLookupF32 (&status, row, "THETA"); - pmPSF_AxesToModel (PAR, axes, modelType); + pmPSF_AxesToModel (PAR, axes, model->class->useReff); source->psfMag = psMetadataLookupF32 (&status, row, "MAG_RAW") - ZERO_POINT; Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_SX.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_SX.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_SX.c (revision 37403) @@ -37,7 +37,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -82,5 +82,5 @@ // pmSourceSextractType (source, &type, &flags); - axes = pmPSF_ModelToAxes (PAR, model->type); + axes = pmPSF_ModelToAxes (PAR, model->class->useReff); psLineInit (line); Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceMatch.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceMatch.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceMatch.c (revision 37403) @@ -18,7 +18,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceMoments.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceMoments.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceMoments.c (revision 37403) @@ -35,7 +35,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceOutputs.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceOutputs.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceOutputs.c (revision 37403) @@ -26,7 +26,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -108,17 +108,29 @@ } if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXY]) && isfinite(PAR[PM_PAR_SYY])) { - axes = pmPSF_ModelToAxes (PAR, model->type); + axes = pmPSF_ModelToAxes (PAR, model->class->useReff); outputs->psfMajor = axes.major; outputs->psfMinor = axes.minor; outputs->psfTheta = axes.theta*PS_DEG_RAD; + + // some models (PS1_V1, QGAUSS) have an extra 'core' parameter + outputs->psfCore = NAN; + if (model->type == pmModelClassGetType ("PS_MODEL_PS1_V1")) { + outputs->psfCore = PAR[PM_PAR_7]; + } + if (model->type == pmModelClassGetType ("PS_MODEL_QGAUSS")) { + outputs->psfCore = PAR[PM_PAR_7]; + } + + outputs->psfMajorFWHM = model->class->modelSetFWHM(model->params, axes.major); + outputs->psfMinorFWHM = model->class->modelSetFWHM(model->params, axes.minor); } else { outputs->psfMajor = NAN; outputs->psfMinor = NAN; outputs->psfTheta = NAN; + outputs->psfCore = NAN; } outputs->chisq = model->chisq; outputs->nDOF = model->nDOF; outputs->nPix = model->nPix; - outputs->apRadius = source->apRadius; } else { bool useMoments = pmSourcePositionUseMoments(source); @@ -138,8 +150,8 @@ outputs->psfMinor = NAN; outputs->psfTheta = NAN; + outputs->psfCore = NAN; outputs->chisq = NAN; outputs->nDOF = 0; outputs->nPix = 0; - outputs->apRadius = NAN; } Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceOutputs.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceOutputs.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceOutputs.h (revision 37403) @@ -26,6 +26,8 @@ float psfMinor; float psfTheta; + float psfCore; + float psfMajorFWHM; + float psfMinorFWHM; float chisq; - float apRadius; int nPix; int nDOF; Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourcePhotometry.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourcePhotometry.c (revision 37403) @@ -33,7 +33,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -330,5 +330,5 @@ // measure fitMag - flux = model->modelFlux (model->params); + flux = model->class->modelFlux (model->params); if (flux > 0) { mag = -2.5*log10(flux); @@ -359,6 +359,5 @@ bool status; - int nPix = 0; - status = pmSourcePhotometryAper(&nPix, &source->apMagRaw, &source->apFlux, &source->apFluxErr, model, image, variance, mask, maskVal); + status = pmSourcePhotometryAper(&source->apNpixels, &source->apMagRaw, &source->apFlux, &source->apFluxErr, model, image, variance, mask, maskVal); if (status) { source->mode |= PM_SOURCE_MODE_AP_MAGS; @@ -490,5 +489,5 @@ // for the full model, add all points - value = fabs(model->modelFunc (NULL, params, coord) - sky); + value = fabs(model->class->modelFunc (NULL, params, coord) - sky); modelSum += value; @@ -893,5 +892,5 @@ // for the full model, add all points - float value = model->modelFunc (NULL, params, coord); + float value = model->class->modelFunc (NULL, params, coord); // fprintf (stderr, "%d, %d : %f, %f : %f - %f : %f\n", Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourcePlotApResid.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourcePlotApResid.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourcePlotApResid.c (revision 37403) @@ -35,7 +35,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourcePlotMoments.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourcePlotMoments.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourcePlotMoments.c (revision 37403) @@ -38,7 +38,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourcePlotPSFModel.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourcePlotPSFModel.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourcePlotPSFModel.c (revision 37403) @@ -39,7 +39,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -147,5 +147,5 @@ // force the axis ratio to be < 20.0 psEllipseAxes axes_mnt = psEllipseMomentsToAxes (moments, 20.0); - psEllipseAxes axes_psf = pmPSF_ModelToAxes (PAR, model->type); + psEllipseAxes axes_psf = pmPSF_ModelToAxes (PAR, model->class->useReff); // moments major axis Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceSky.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceSky.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceSky.c (revision 37403) @@ -34,7 +34,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceUtils.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceUtils.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceUtils.c (revision 37403) @@ -34,7 +34,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -65,5 +65,5 @@ pmModel *model = pmModelAlloc(modelType); - if (!model->modelGuess(model, source, maskVal, markVal)) { + if (!model->class->modelGuess(model, source, maskVal, markVal)) { psFree (model); return NULL; Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceVisual.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceVisual.c (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceVisual.c (revision 37403) @@ -16,7 +16,7 @@ #include "pmMoments.h" #include "pmModelFuncs.h" +#include "pmModelClass.h" #include "pmModel.h" #include "pmModelUtils.h" -#include "pmModelClass.h" #include "pmSourceMasks.h" #include "pmSourceExtendedPars.h" @@ -545,7 +545,31 @@ psFree (model); + bool dumpData = false; + // pause and wait for user input: // continue, save (provide name), ?? - pmVisualAskUser(&plotPSF); +retry: + pmVisualAskUserOrDump(&plotPSF, &dumpData); + if (dumpData) { + char name[128]; + fprintf (stderr, "filename: "); + int status = fscanf (stdin, "%127s", name); + if (status != 1) { + fprintf (stderr, "odd response\n"); + goto retry; + } + + FILE *f = fopen (name, "w"); + if (!f) { + fprintf (stderr, "cannot open %s for output\n", name); + goto retry; + } + for (int i = 0; i < x->n; i++) { + float vModel = pmTrend2DEval (trend, x->data.F32[i], y->data.F32[i]); + fprintf (f, "%f %f %f %f %d\n", x->data.F32[i], y->data.F32[i], param->data.F32[i], vModel, mask->data.PS_TYPE_VECTOR_MASK_DATA[i]); + } + fclose (f); + goto retry; + } return true; Index: branches/eam_branches/ps2-tc3-20130727/psModules/src/psmodules.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psModules/src/psmodules.h (revision 37399) +++ branches/eam_branches/ps2-tc3-20130727/psModules/src/psmodules.h (revision 37403) @@ -95,4 +95,5 @@ #include #include +#include // the following headers are from psModule:imcombine @@ -128,4 +129,5 @@ #include +#include #include #include @@ -148,5 +150,4 @@ #include #include -#include #include #include