Index: /branches/eam_branches/ipp-20120405/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphot.h (revision 33945) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphot.h (revision 33946) @@ -99,5 +99,5 @@ bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final); -bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final); +bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode); bool psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize); Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotDefineFiles.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotDefineFiles.c (revision 33945) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotDefineFiles.c (revision 33946) @@ -85,4 +85,26 @@ output->save = true; } +// // optionally save the smoothed variance model (small FITS image) +// if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKMDL")) { +// int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN"); +// int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); +// pmFPAfile *output = pmFPAfileDefineFromFile (config, input, DX, DY, "PSPHOT.VARMDL"); +// if (!output) { +// psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.VARMDL"); +// return false; +// } +// output->save = true; +// } +// // optionally save the smoothed variance model's standard deviation (small FITS image) +// if (psMetadataLookupBool(NULL, recipe, "SAVE.VARMDL.STDEV")) { +// int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN"); +// int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); +// pmFPAfile *output = pmFPAfileDefineFromFile (config, input, DX, DY, "PSPHOT.VARMDL.STDEV"); +// if (!output) { +// psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.VARMDL.STDEV"); +// return false; +// } +// output->save = true; +// } // optionally save the PSF Model if (psMetadataLookupBool(NULL, recipe, "SAVE.PSF")) { Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotEfficiency.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotEfficiency.c (revision 33945) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotEfficiency.c (revision 33946) @@ -420,5 +420,5 @@ // psphotFitSourcesLinearReadout subtracts the model fits - if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) { + if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST)) { psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources."); psFree(fakeSources); Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotFitSourcesLinear.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotFitSourcesLinear.c (revision 33945) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotFitSourcesLinear.c (revision 33946) @@ -12,4 +12,8 @@ static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal); +bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources); +pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe); +bool psphotFreeModelVariance (pmReadout *readout, psArray *sources); + // for now, let's store the detections on the readout->analysis for each readout bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final) @@ -24,4 +28,14 @@ assert (recipe); + pmSourceFitVarMode fitVarMode = psphotGetFitVarMode (recipe); + if (!fitVarMode) { + psError (PSPHOT_ERR_CONFIG, true, "failed to get LINEAR_FIT_VARIANCE_MODE"); + return false; + } + // MODEL_VAR requires 2 passes -- in the first, we get the rough fluxes; in the second, we + // use the flux to define the model variance before fitting the objects. Other modes only + // do a single pass. + pmSourceFitVarMode fitVarModePass1 = (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) ? PM_SOURCE_PHOTFIT_CONST : fitVarMode; + int num = psphotFileruleCount(config, filerule); @@ -50,8 +64,34 @@ psAssert (psf, "missing psf?"); - if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) { + if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1)) { psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i); return false; } + + // the MODEL_VAR weighting scheme requires knowledge of the model fluxes to generate the variance + // after we have determined the initial set of fits, then we can generate the variance image and + // re-run the fit against that variance. + if (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) { + // generate the model variance image & source pointers + if (!psphotGenerateModelVariance (config, view, file, i, recipe, readout, sources)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i); + return false; + } + + // replace all sources (use TMPF_SUBTRACTED as test) + psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe, false); + + // rerun fit with correct fitVarMode + if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarMode)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i); + return false; + } + + // free the model variance image & source pointers + if (!psphotFreeModelVariance (readout, sources)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i); + return false; + } + } psphotVisualShowResidualImage (readout, (num > 0)); @@ -62,5 +102,34 @@ } -bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final) { +// look up the fit variance mode from the recipe; older recipes do not have the value +// 'LINEAR_FIT_VARIANCE_MODE'; in those cases, look for 'CONSTANT_PHOTOMETRIC_WEIGHTS' as a boolean and +// set the value to either CONST or IMAGE_VAR +pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe) { + + bool status = false; + + char *fitVarModeString = psMetadataLookupStr(&status, recipe, "LINEAR_FIT_VARIANCE_MODE"); + if (!status) { + bool CONSTANT_PHOTOMETRIC_WEIGHTS = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); + if (!status) { + psAbort("You must provide a value for LINEAR_FIT_VARIANCE_MODE or CONSTANT_PHOTOMETRIC_WEIGHTS"); + } + pmSourceFitVarMode fitVarMode = CONSTANT_PHOTOMETRIC_WEIGHTS ? PM_SOURCE_PHOTFIT_CONST : PM_SOURCE_PHOTFIT_IMAGE_VAR; + return fitVarMode; + } + if (!strcasecmp(fitVarModeString, "CONSTANT") || !strcasecmp(fitVarModeString, "CONST")) { + return PM_SOURCE_PHOTFIT_CONST; + } + if (!strcasecmp(fitVarModeString, "IMAGE") || !strcasecmp(fitVarModeString, "IMAGE_VAR")) { + return PM_SOURCE_PHOTFIT_IMAGE_VAR; + } + if (!strcasecmp(fitVarModeString, "MODEL") || !strcasecmp(fitVarModeString, "MODEL_VAR")) { + return PM_SOURCE_PHOTFIT_MODEL_VAR; + } + psError (PSPHOT_ERR_CONFIG, false, "Invalid value for LINEAR_FIT_VARIANCE_MODE (%s)", fitVarModeString); + return PM_SOURCE_PHOTFIT_NONE; +} + +bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode) { bool status; @@ -99,9 +168,4 @@ if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); - bool CONSTANT_PHOTOMETRIC_WEIGHTS = - psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); - if (!status) { - psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS"); - } int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER"); if (!status) { @@ -231,4 +295,12 @@ psSparseBorder *border = psSparseBorderAlloc (sparse, nBorder); + // if fitVarMode is MODEL_VAR, then we need to generate the model image variance + // XXX we have two possibilities here: + + // 1) do 2 passes, where in the first case we use the CONST weighting, and in the second + // use the fitted model values to define the model + + // 2) do a single pass, and use the model guess to define the model variance (but do I trust the Model Guess?) + // fill out the sparse matrix elements and border elements (B) // SRCi is the current source of interest @@ -238,10 +310,10 @@ // diagonal elements of the sparse matrix (auto-cross-product) - f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); + f = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal); psSparseMatrixElement (sparse, i, i, f); - // the formal error depends on the weighting scheme - if (CONSTANT_PHOTOMETRIC_WEIGHTS) { - float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal); + // if we have used CONSTANT errors, then we need to calculate the value of the parameter error + if (fitVarMode != PM_SOURCE_PHOTFIT_IMAGE_VAR) { + float var = pmSourceModelDotModel (SRCi, SRCi, PM_SOURCE_PHOTFIT_IMAGE_VAR, covarFactor, maskVal); errors->data.F32[i] = 1.0 / sqrt(var); } else { @@ -251,5 +323,5 @@ // find the image x model value - f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); + f = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal); psSparseVectorElement (sparse, i, f); @@ -257,11 +329,11 @@ switch (SKY_FIT_ORDER) { case 1: - f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); + f = pmSourceModelWeight (SRCi, 1, fitVarMode, covarFactor, maskVal); psSparseBorderElementB (border, i, 1, f); - f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); + f = pmSourceModelWeight (SRCi, 2, fitVarMode, covarFactor, maskVal); psSparseBorderElementB (border, i, 2, f); case 0: - f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); + f = pmSourceModelWeight (SRCi, 0, fitVarMode, covarFactor, maskVal); psSparseBorderElementB (border, i, 0, f); break; @@ -283,5 +355,5 @@ // got an overlap; calculate cross-product and add to output array - f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); + f = pmSourceModelDotModel (SRCi, SRCj, fitVarMode, covarFactor, maskVal); psSparseMatrixElement (sparse, j, i, f); } @@ -321,5 +393,5 @@ // set the sky, sky_x, sky_y components of border matrix - SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal); + SetBorderMatrixElements (border, readout, fitSources, (fitVarMode == PM_SOURCE_PHOTFIT_CONST), SKY_FIT_ORDER, markVal); psSparseConstraint constraint; @@ -479,2 +551,112 @@ return true; } + +bool psphotModelBackgroundReadout(psImage *model, // Model image + psImage *modelStdev, // Model stdev image + psMetadata *analysis, // Analysis metadata for outputs + pmReadout *readout, // Readout for which to generate a background model + psImageBinning *binning, // Binning parameters + const pmConfig *config,// Configuration + bool useVarianceImage + ); + +bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) { + + bool status = false; + psRegion fullRegion = psRegionSet (0, 0, 0, 0); + + // bit-masks to test for good/bad pixels + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + assert (maskVal); + + // bit-mask to mark pixels not used in analysis + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); + assert (markVal); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // create a model variance image + psImage *modelVar = psImageCopy (NULL, readout->variance, PS_TYPE_F32); + + // find the binning information + psImageBinning *backBinning = psphotBackgroundBinning (modelVar, config); + assert (backBinning); + + psImage *varModel = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model + psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model + + if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate background model"); + psFree (varModel); + psFree (varModelStdev); + return false; + } + + // linear interpolation to full-scale + if (!psImageUnbin (modelVar, varModel, backBinning)) { + psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning"); + psFree (varModel); + psFree (varModelStdev); + return false; + } + + psFree (varModel); + psFree (varModelStdev); + + // XXX for a test: + psphotSaveImage (NULL, modelVar, "model.bck.fits"); + + // insert all of the source models + for (int i = 0; i < sources->n; i++) { + + // source of interest + pmSource *source = sources->data[i]; + + // skip sources which were not fitted already + if (!(source->mode & PM_SOURCE_MODE_LINEAR_FIT)) continue; + + // pixel region appropriate for the source + psRegion region = psRegionForImage (source->pixels, fullRegion); + + // define the source->modelVar pixels (view on modelVar image) + psAssert (!source->modelVar, "programming error : modelVar should be NULL here"); + psAssert (source->modelFlux, "programming error : modelFlux should not be NULL here"); + psAssert (source->modelFlux->data.F32, "programming error : modelFlux should not be NULL here"); + source->modelVar = psImageSubset(modelVar, region); + + // add the source model to the model variance image + pmSourceAdd (source, PM_MODEL_OP_MODELVAR, maskVal); + } + + // XXX for a test: + psphotSaveImage (NULL, modelVar, "model.var.fits"); + psphotSaveImage (NULL, readout->variance, "image.var.fits"); + + // we save the model variance for future reference + psMetadataAddImage(readout->analysis, PS_LIST_TAIL, "PSPHOT.MODEL.VAR", PS_META_REPLACE, "model variance", modelVar); + psFree (modelVar); + + return true; +} + +bool psphotFreeModelVariance (pmReadout *readout, psArray *sources) { + + bool status = false; + + // find the binning information + psImage *modelVar = psMetadataLookupPtr(&status, readout->analysis, "PSPHOT.MODEL.VAR"); + assert (modelVar); + + psMetadataRemoveKey (readout->analysis, "PSPHOT.MODEL.VAR"); + + // clear modelVar pointers for all of the source models + for (int i = 0; i < sources->n; i++) { + + // source of interest + pmSource *source = sources->data[i]; + psFree (source->modelVar); + } + + return true; +} Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotKronIterate.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotKronIterate.c (revision 33945) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotKronIterate.c (revision 33946) @@ -83,4 +83,9 @@ if (!status) { MIN_KRON_RADIUS = 0.25*RADIUS; + } + + int KRON_ITERATIONS = psMetadataLookupS32 (&status, recipe, "KRON_ITERATIONS"); + if (!status) { + KRON_ITERATIONS = 2; } @@ -140,7 +145,8 @@ PS_ARRAY_ADD_SCALAR(job->args, RADIUS, PS_TYPE_F32); PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, KRON_ITERATIONS, PS_TYPE_S32); // set this to 0 to run without threading -# if (1) +# if (0) if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); @@ -150,5 +156,4 @@ if (!psphotKronIterate_Threaded(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - // psFree(AnalysisRegion); return false; } @@ -188,7 +193,8 @@ float RADIUS = PS_SCALAR_VALUE(job->args->data[5],F32); float MIN_KRON_RADIUS = PS_SCALAR_VALUE(job->args->data[6],F32); + int KRON_ITERATIONS = PS_SCALAR_VALUE(job->args->data[7],S32); // XXX TEST : set iteration to 1 - for (int j = 0; j < 1; j++) { + for (int j = 0; j < KRON_ITERATIONS; j++) { for (int i = 0; i < sources->n; i++) { @@ -217,8 +223,18 @@ // XXX float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), maxWindow); + // Sextractor apparently takes something like the skyRadius window, measures the moments inside that window, + // then uses an elliptical outline based on the 2nd moments, but 6x the size + // XXX TEST : use a window based on the radial profile numbers: max is skyRadius, min is RADIUS // if we lack the skyRadius (eg MATCHED sources), go to the default value - float maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS; - float windowRadius = PS_MAX(RADIUS, maxWindow); + float maxWindow, windowRadius; + + if (j == 0) { + maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS; + windowRadius = PS_MAX(RADIUS, maxWindow); + } else { + maxWindow = isfinite(source->moments->Mrf) ? 6.0*source->moments->Mrf : RADIUS; + windowRadius = PS_MAX(RADIUS, maxWindow); + } #ifdef notdef Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotModelBackground.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotModelBackground.c (revision 33945) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotModelBackground.c (revision 33946) @@ -37,10 +37,11 @@ // corresponding to the model. Other information about the background model is saved on the // readout->analysis -static bool psphotModelBackgroundReadout(psImage *model, // Model image - psImage *modelStdev, // Model stdev image - psMetadata *analysis, // Analysis metadata for outputs - pmReadout *readout, // Readout for which to generate a background model - psImageBinning *binning, // Binning parameters - const pmConfig *config // Configuration +bool psphotModelBackgroundReadout(psImage *model, // Model image + psImage *modelStdev, // Model stdev image + psMetadata *analysis, // Analysis metadata for outputs + pmReadout *readout, // Readout for which to generate a background model + psImageBinning *binning, // Binning parameters + const pmConfig *config,// Configuration + bool useVarianceImage ) { @@ -49,5 +50,6 @@ bool status = true; - psImage *image = readout->image, *mask = readout->mask; // Image and mask for readout + psImage *image = useVarianceImage ? readout->variance : readout->image; + psImage *mask = readout->mask; // Image and mask for readout // select the appropriate recipe information @@ -143,5 +145,7 @@ // we save the binning structure for use in psphotMagnitudes - psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning); + if (!useVarianceImage) { + psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning); + } psVector *dQ = psVectorAllocEmpty (100, PS_TYPE_F32); @@ -159,6 +163,6 @@ // measure clipped median for subimages - psRegion ruffRegion = {0,0,0,0}; - psRegion fineRegion = {0,0,0,0}; + psRegion ruffRegion = psRegionSet (0,0,0,0); + psRegion fineRegion = psRegionSet (0,0,0,0); for (int iy = 0; iy < model->numRows; iy++) { for (int ix = 0; ix < model->numCols; ix++) { @@ -330,25 +334,30 @@ psVectorStats (statsDQ, dQ, NULL, NULL, 0); - psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value); - psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev); - psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_DQ", PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian); - psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian); - - // measure image and background stats and save for later output - psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN | - PS_STAT_SAMPLE_STDEV | - PS_STAT_MIN | - PS_STAT_MAX); - psImageStats (statsBck, model, NULL, 0); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN", PS_META_REPLACE, "sky model mean", statsBck->sampleMean); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_SIG", PS_META_REPLACE, "sky model stdev", statsBck->sampleStdev); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DEV", PS_META_REPLACE, "sky stdev", ValueStdev); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DQ", PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MAX", PS_META_REPLACE, "sky model maximum value", statsBck->max); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MIN", PS_META_REPLACE, "sky model minimum value", statsBck->min); - psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NX", PS_META_REPLACE, "sky model size (x)", model->numCols); - psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NY", PS_META_REPLACE, "sky model size (y)", model->numRows); - psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f", - statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev); + if (!useVarianceImage) { + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_DQ", PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian); + psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian); + + // measure image and background stats and save for later output + psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN | + PS_STAT_SAMPLE_STDEV | + PS_STAT_MIN | + PS_STAT_MAX); + psImageStats (statsBck, model, NULL, 0); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN", PS_META_REPLACE, "sky model mean", statsBck->sampleMean); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_SIG", PS_META_REPLACE, "sky model stdev", statsBck->sampleStdev); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DEV", PS_META_REPLACE, "sky stdev", ValueStdev); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DQ", PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MAX", PS_META_REPLACE, "sky model maximum value", statsBck->max); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MIN", PS_META_REPLACE, "sky model minimum value", statsBck->min); + psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NX", PS_META_REPLACE, "sky model size (x)", model->numCols); + psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NY", PS_META_REPLACE, "sky model size (y)", model->numRows); + psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f", + statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev); + psFree(statsBck); + } else { + psLogMsg ("psphot", PS_LOG_INFO, "variance data : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian); + } psFree(statsDQ); @@ -356,5 +365,4 @@ psFree(stats); - psFree(statsBck); psFree(statsDefaults); psFree(binning); @@ -375,5 +383,5 @@ psImage *modelStdev = psImageAlloc(binning->nXruff, binning->nYruff, PS_TYPE_F32); // Standard deviation - if (!psphotModelBackgroundReadout(model, modelStdev, readout->analysis, readout, binning, config)) { + if (!psphotModelBackgroundReadout(model, modelStdev, readout->analysis, readout, binning, config, false)) { psFree(model); psFree(modelStdev); @@ -400,5 +408,5 @@ pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning, index); - if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config)) { + if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config, false)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate background model"); return false; Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotSetThreads.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotSetThreads.c (revision 33945) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotSetThreads.c (revision 33946) @@ -30,5 +30,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 7); + task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 8); task->function = &psphotKronIterate_Threaded; psThreadTaskAdd(task); Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotTest.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotTest.c (revision 33945) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotTest.c (revision 33946) @@ -103,5 +103,5 @@ # if (0) - psRegion region = {0,0,0,0}; // a region representing the entire array +psRegion region = psRegionSet (0,0,0,0); // a region representing the entire array psphotTestArguments (&argc, argv); Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotTestSourceOutput.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotTestSourceOutput.c (revision 33945) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotTestSourceOutput.c (revision 33946) @@ -136,5 +136,5 @@ // generate working image for this source - psRegion region = {ix - dx, ix + dx, iy - dy, iy + dy}; + psRegion region = psRegionSet(ix - dx, ix + dx, iy - dy, iy + dy); psImage *vM = psImageSubset (imMo, region); Index: /branches/eam_branches/ipp-20120405/psphot/test/tap_psphot_varmodel.pro =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/test/tap_psphot_varmodel.pro (revision 33946) +++ /branches/eam_branches/ipp-20120405/psphot/test/tap_psphot_varmodel.pro (revision 33946) @@ -0,0 +1,642 @@ +#!/usr/bin/env mana +# -*-sh-*- + +# $KAPA = kapa -noX + +# PSF.CONVOLVE : if true, we insert delta functions (and optionally +# galaxies) and smooth the image with the psf model +# (uses a GAUSS regardless of the model). Note that +# PSF.CONVOLVE = T is faster than F, but (a) only +# allows Gauss models and (b) only yields quantized +# locations + +# config for ppImage to generate chip, mask, weight +$ppImageConfig = -recipe PPIMAGE PPIMAGE_N +$ppImageConfig = $ppImageConfig -Db BACKGROUND T +$ppImageConfig = $ppImageConfig -Db CHIP.FITS T +$ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T +$ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T +$ppImageConfig = $ppImageConfig -Db BASE.FITS F +$ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T +$ppImageConfig = $ppImageConfig -Db PHOTOM F + +# basic options for the these images (filter, location, obstype) +$BaseOptions = -type OBJECT -filter r -skymags 20.86 -ra 270.70 -dec -23.70 -pa 0.0 +$BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0 + +# create an image with fake sources and insert the resulting cmf file into a dvodb +$RefConfig = -camera SIMTEST +$RefConfig = $RefConfig -recipe PPSIM STACKTEST.MAKE +$RefConfig = $RefConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref +$RefConfig = $RefConfig -Db PSF.CONVOLVE F + +# options for the reference image +$RefOptions = $BaseOptions -exptime 100.0 +$RefOptions = $RefOptions -seeing 1.0 +$RefOptions = $RefOptions -D PSF.MODEL PS_MODEL_GAUSS +$RefOptions = $RefOptions -Df STARS.DENSITY 10.0 +$RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5 + +# basic config for ppSim with randomly distributed stars and NO galaxies +$RealConfig = -camera SIMTEST +$RealConfig = $RealConfig -recipe PPSIM STACKTEST.RUN +$RealConfig = $RealConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref +$RealConfig = $RealConfig -Db STARS.FAKE F +$RealConfig = $RealConfig -Db STARS.REAL T +$RealConfig = $RealConfig -Db MATCH.DENSITY F +$RealConfig = $RealConfig -Db PSF.CONVOLVE F +$RealConfig = $RealConfig -Df STARS.DENSITY 10.0 +$RealConfig = $RealConfig -Df STARS.SIGMA.LIM 2.5 +$RealConfig = $RealConfig -Db GALAXY.FAKE F +$RealConfig = $RealConfig -Db GALAXY.GRID F + +# options for the repeated images +$RealOptions = $BaseOptions -exptime 30.0 + +$ExtraOptions = -D PSF.MODEL PS_MODEL_GAUSS + +# sample alternate options: +# $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1 +# $ppSimOptions = $FakeOptions -Df PSF.ARATIO 1.2 +# $ppSimOptions = $FakeOptions -Df PSF.THETA +30.0 +# $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS + +list fwhm + 1.0 + 1.1 + 1.2 + 1.5 +end + +macro go + + mkdir test + + $ExtraOptions = -D PSF.MODEL PS_MODEL_GAUSS + mkexp test/image.00 1.0 + $ExtraOptions = -D PSF.MODEL PS_MODEL_PGAUSS + mkexp test/image.01 1.0 + $ExtraOptions = -D PSF.MODEL PS_MODEL_PS1_V1 + mkexp test/image.02 1.0 +end + +# create a reference database of fake stars to be used by ppSim below +macro mkref + exec rm -rf catdir.ref + exec rm -f refimage.fits + + exec time ppSim $RefConfig $RefOptions refimage -nx 3000 -ny 3000 + + file synth.photcodes found + if (not($found)) + echo "making photcodes file" + mkphotcodes synth.photcodes + end + + exec addstar -D CAMERA simtest -D CATDIR catdir.ref -accept-astrom -photcode SYNTH.r -D PHOTCODE_FILE synth.photcodes refimage.cmf + exec relphot -averages -D CATDIR catdir.ref -update -region 260 280 -33 -13 +end + +# create a realistic distribution of fake stars, GAUSS PSF +macro mkexp + if ($0 != 3) + echo "USAGE: mkexp basename fwhm" + break + end + + local fwhm basename + $basename = $1 + $fwhm = $2 + + # create the raw image + echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename + exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename + exec /bin/mv -f $basename.cmf $basename.in.cmf + + # create the chip output + echo ppImage $ppImageConfig -file $basename.fits $basename + exec ppImage $ppImageConfig -file $basename.fits $basename +end + +macro runphot + if ($0 != 4) + echo "USAGE: runphot basename outname options" + break + end + + local basename + $basename = $1 + $outname = $2 + $options = $3 + + # create the chip output + echo psphot -threads 4 -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $outname $options + exec psphot -threads 4 -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $outname $options +end + +# compare two cmf files with extname Chip.psf +# things to compare: +# * completeness (which sources in (1) are not detected in (2) +# * positions (X_PSF, Y_PSF) +# * instrumental psf mags +# * position errors (no input errors; use a model?) +# * measured FWHM? +# * kron mags (fluxes) +# * etc, etc +macro ckchip + if ($0 != 5) + echo "USAGE: ckchip (raw) (out) (output) (zpt_off)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 Chip.psf out + + # images generated with convolution will not have the right output positions + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + set M_raw = PSF_INST_MAG_raw + $4 + set K_out = -2.5*log(KRON_FLUX_out) + match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2 + + local i NX NY nx ny N + + device -n compare + resize 1000 1000 + + # plot trends as a function of mag + $NX = 2 + $NY = 5 + $nx = 0 + $ny = 0 + $N = 0 + clear -s + for i 0 $pairs:n + section a$nx\$ny {$nx/$NX} {$ny/$NY} {1/$NX} {1/$NY} + show.pair $i + $ny ++ + if ($ny == $NY) + $ny = 0 + $nx ++ + end + if ($nx == $NX) + png -name $3.$N.png + clear -s + $nx = 0 + $ny = 0 + $N ++ + end + end + + # plot (input - output) vs mag +end + +macro stchip + if ($0 != 5) + echo "USAGE: stchip (raw) (out) (output) (zpt_off)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 Chip.psf out + + # images generated with convolution will not have the right output positions + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + set M_raw = PSF_INST_MAG_raw + $4 + set K_out = -2.5*log(KRON_FLUX_out) + match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2 + + local i + + for i 0 $spairs:n + stats.pair $i $3 + end +end + +# compare chip to warp +macro ckwarp + if ($0 != 5) + echo "USAGE: ckwarp (raw) (out) (output) (zpt_off)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 SkyChip.psf out + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + set M_raw = PSF_INST_MAG_raw + $4 + set K_out = -2.5*log(KRON_FLUX_out) + match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.0 -index1 index1 -index2 index2 + + local i nx ny NX NY N + + device -n compare + resize 1000 1000 + + # plot trends as a function of mag + $NX = 2 + $NY = 5 + $nx = 0 + $ny = 0 + $N = 0 + clear -s + + for i 0 $pairs:n + section a$nx\$ny {$nx/$NX} {$ny/$NY} {1/$NX} {1/$NY} + show.pair $i + $ny ++ + if ($ny == $NY) + $ny = 0 + $nx ++ + end + if ($nx == $NX) + png -name $3.$N.png + clear -s + $nx = 0 + $ny = 0 + $N ++ + end + end +end + +macro stwarp + if ($0 != 5) + echo "USAGE: stwarp (raw) (out) (output) (zpt_off)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 SkyChip.psf out + + # images generated with convolution will not have the right output positions + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + set M_raw = PSF_INST_MAG_raw + $4 + set K_out = -2.5*log(KRON_FLUX_out) + match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2 + + local i + + for i 0 $spairs:n + stats.pair $i $3 + end +end + +macro stats.pair + if ($0 != 3) + echo "USAGE: stats.pair (N) (output)" + break + end + + list word -split $spairs:$1 + if ($word:n != 8) + echo "invalid pair $1" + break + end + + $Nr = $word:3 + + reindex v1 = $word:0 using index1 + reindex v2 = $word:1 using index2 + reindex rv = $word:2 using index$Nr + + set delta = v1 - v2 + subset d1 = delta if ($word:4 < rv) && (rv < $word:5) && (abs(delta) < $word:7) + subset d2 = delta if ($word:5 < rv) && (rv < $word:6) && (abs(delta) < $word:7) + + vstats -q d1 -sigma-clip 3.0 + $M1 = $MEAN + $S1 = $SIGMA + vstats -q d2 -sigma-clip 3.0 + $M2 = $MEAN + $S2 = $SIGMA + + output $2 + fprintf "%-18s %7.4f %7.4f %7.4f %7.4f" $word:0 $M1 $S1 $M2 $S2 + output stdout +end + +macro show.pair + if ($0 != 2) + echo "USAGE: show.pair (N)" + break + end + + list word -split $pairs:$1 + if ($word:n != 7) + echo "invalid pair $1" + break + end + + $Nr = $word:3 + + reindex v1 = $word:0 using index1 + if ("$word:6" == "V") + reindex v2 = $word:1 using index2 + end + if ("$word:6" == "S") + set v2 = $word:1 + zero(index1) + end + reindex rv = $word:2 using index$Nr + + set delta = v1 - v2 + if (("$word:4" == "def") || ("$word:5" == "def")) + lim rv delta; box; plot rv delta + else + lim rv $word:4 $word:5; box; plot rv delta + end + label -y '$word:0' -x '$word:2' +end + +# This list is used to compare a pair of vectors (sans error) or a +# vector and an expected (constant) value. The last field defines a +# vector or constant for the comparison. It is assumed that the +# vector sets have been loaded and matched with match2d to generate +# index vectors 'index1' and index2'. The macro 'show.pair' generates +# a plot of the range vector vs (v1 - v2). The indices for v1, v2 are +# index1 and 2 respectively. The index for the range vector is defined +# by the integer following that vector. the y-limits of the plot are +# given by the last two numbers +list pairs + X_PSF_out X_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + Y_PSF_out Y_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + X_PSF_out X_PSF_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + Y_PSF_out Y_PSF_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + X_PSF_SIG_out X_PSF_SIG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + Y_PSF_SIG_out Y_PSF_SIG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + #PSF_INST_MAG_out PSF_INST_MAG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + PSF_INST_MAG_out M_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + PSF_INST_MAG_SIG_out PSF_INST_MAG_SIG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + #PSF_INST_FLUX_out PSF_INST_FLUX_raw PSF_INST_MAG_raw 2 def def V + #PSF_INST_FLUX_SIG_out PSF_INST_FLUX_SIG_raw PSF_INST_MAG_raw 2 def def V + AP_MAG_out M_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + AP_MAG_RAW_out M_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + AP_MAG_RADIUS_out 0.0 PSF_INST_MAG_raw 2 -0.01 20.1 S + SKY_out 0.0 PSF_INST_MAG_raw 2 def def S + SKY_SIGMA_out 0.0 PSF_INST_MAG_raw 2 def def S + PSF_CHISQ_out 1.0 PSF_INST_MAG_raw 2 def def S + CR_NSIGMA_out 0.0 PSF_INST_MAG_raw 2 def def S + EXT_NSIGMA_out 0.0 PSF_INST_MAG_raw 2 -5.01 5.01 S + PSF_MAJOR_out 0.0 PSF_INST_MAG_raw 2 -0.01 5.01 S + PSF_MINOR_out 0.0 PSF_INST_MAG_raw 2 -0.01 5.01 S + PSF_THETA_out 0.0 PSF_INST_MAG_raw 2 -1.61 1.61 S + PSF_QF_out 0.0 PSF_INST_MAG_raw 2 -0.10 1.10 S + PSF_QF_PERFECT_out 0.0 PSF_INST_MAG_raw 2 -0.10 1.10 S + PSF_NDOF_out 0.0 PSF_INST_MAG_raw 2 def def S + PSF_NPIX_out 0.0 PSF_INST_MAG_raw 2 def def S + MOMENTS_XX_out 0.0 PSF_INST_MAG_raw 2 -0.01 3.01 S + MOMENTS_XY_out 0.0 PSF_INST_MAG_raw 2 -3.01 3.01 S + MOMENTS_YY_out 0.0 PSF_INST_MAG_raw 2 -0.01 3.01 S + MOMENTS_M3C_out 0.0 PSF_INST_MAG_raw 2 -3.01 3.01 S + MOMENTS_M3S_out 0.0 PSF_INST_MAG_raw 2 -3.01 3.01 S + MOMENTS_M4C_out 0.0 PSF_INST_MAG_raw 2 -2.01 2.01 S + MOMENTS_M4S_out 0.0 PSF_INST_MAG_raw 2 -2.01 2.01 S + MOMENTS_R1_out 0.0 PSF_INST_MAG_raw 2 -5.01 5.01 S + MOMENTS_RH_out 0.0 PSF_INST_MAG_raw 2 -5.01 5.01 S + K_out M_raw PSF_INST_MAG_raw 2 def def V + KRON_FLUX_ERR_out 0.0 PSF_INST_MAG_raw 2 def def S + KRON_FLUX_INNER_out 0.0 PSF_INST_MAG_raw 2 def def S + KRON_FLUX_OUTER_out 0.0 PSF_INST_MAG_raw 2 def def S +# CAL_PSF_MAG CAL_PSF_MAG none Mraw 2 -1.01 1.01 V +# CAL_PSF_MAG_SIG CAL_PSF_MAG_SIG none Mraw 2 -1.01 1.01 V +# RA_PSF RA_PSF none Mraw 2 -1.01 1.01 V +# DEC_PSF DEC_PSF none Mraw 2 -1.01 1.01 V +# PEAK_FLUX_AS_MAG PEAK_FLUX_AS_MAG none Mraw 2 -1.01 1.01 V +# FLAGS FLAGS 0.0 Mraw 2 -1.01 1.01 S +# FLAGS2 FLAGS2 0.0 Mraw 2 -1.01 1.01 S +end + +macro load.cmf + if ($0 != 4) + echo "load.cmf (filename) (ext) (label)" + break + end + + data $1 + + # create the list of fields to load + delete -q myFields + for i 0 $fields:n + if ($?myFields) + $myFields = $myFields $fields:$i + else + $myFields = $fields:$i + end + end + + read -fits $2 $myFields + + # rename the loaded vectors appending the supplied lable + for i 0 $fields:n + set $fields:$i\_$3 = $fields:$i + delete $fields:$i + end +end + +# this list defines the fields to be loaded from file +list fields + X_PSF + Y_PSF + X_PSF_SIG + Y_PSF_SIG + PSF_INST_MAG + PSF_INST_MAG_SIG + PSF_INST_FLUX + PSF_INST_FLUX_SIG + AP_MAG + AP_MAG_RAW + AP_MAG_RADIUS + SKY + SKY_SIGMA + PSF_CHISQ + CR_NSIGMA + EXT_NSIGMA + PSF_MAJOR + PSF_MINOR + PSF_THETA + PSF_QF + PSF_QF_PERFECT + PSF_NDOF + PSF_NPIX + MOMENTS_XX + MOMENTS_XY + MOMENTS_YY + MOMENTS_M3C + MOMENTS_M3S + MOMENTS_M4C + MOMENTS_M4S + MOMENTS_R1 + MOMENTS_RH + KRON_FLUX + KRON_FLUX_ERR + KRON_FLUX_INNER + KRON_FLUX_OUTER +# CAL_PSF_MAG +# CAL_PSF_MAG_SIG +# RA_PSF +# DEC_PSF +# PEAK_FLUX_AS_MAG +# FLAGS +# FLAGS2 +end + +# use these cmf entries to measure average stats of the given pairs +list spairs + X_PSF_out X_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 + Y_PSF_out Y_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 + PSF_INST_MAG_out M_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 + AP_MAG_out M_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 +end + +# # XXX this is a hack: the cmf file created by ppSim is not +# # compatible with the image for warp (because one is a chip-mosaic +# # and the other is not) +# exec rm -f fix.hdr +# output fix.hdr +# echo "PSMOSAIC= 'CHIP ' / Mosaicked level" +# output stdout +# fix the header to be compatible with the chip file (to avoid running psphot) +# exec fits_insert $basename.cmf fix.hdr + +macro completeness + if ($0 != 4) + echo "USAGE: completeness (raw) (out) (output)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 Chip.psf out + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + match2d X_PSF_raw Y_PSF_raw X_PSF_out Y_PSF_out 1.5 -index1 index1 -index2 index2 -closest + + histogram PSF_INST_MAG_raw nMag -16.0 -3.0 0.25 -range dMag + set fMag = zero(dMag) + for i 0 {dMag[]-1} + set inrange = (PSF_INST_MAG_raw > dMag[$i]) && (PSF_INST_MAG_raw <= dMag[$i+1]) + subset all = index1 if (inrange) + subset got = index1 if (inrange) && (index1 >= 0) + if (all[] == 0) + fMag[$i] = 0 + else + fMag[$i] = got[] / all[] + end + end + + device -n complete + resize 1000 600 + + clear -s + + section default 0 0 1 1 + lim dMag fMag; clear; box -ypad 5 +ypad 5 -ticks 1110; plot -x 1 dMag fMag + label -x mag_inst -y det_frac + + set found = (index1 >= 0) + plot PSF_INST_MAG_raw found -c red + + set ldmag = log(PSF_INST_MAG_SIG_raw) + section overlay 0 0 1 1; lim dMag -5 1.2; box -ypad 5 +ypad 5 -ticks 1011 -labels 1001; plot PSF_INST_MAG_raw ldmag + label +y log(S/N) + + png -name $3 +end + +macro show.dpair + if ($0 != 2) + echo "USAGE: show.dpair (N)" + break + end + + list word -split $dpairs:$1 + if ($word:n != 7) + echo "invalid dpair $1" + break + end + + $Nr = $word:4 + + reindex v1 = $word:0 using index1 + reindex dv = $word:1 using index1 + reindex v2 = $word:2 using index2 + reindex rv = $word:3 using index$Nr + + set delta = (v1 - v2) / dv + if (("$word:5" == "def") || ("$word:6" == "def")) + lim rv delta; box; plot rv delta + else + lim rv $word:5 $word:6; box; plot rv delta + end + label -y '$word:0' -x '$word:3' +end + +# this list is used to compare a pair of vectors with an error it is +# assumed that the vector sets have been loaded and matched with +# match2d to generate index vectors 'index1' and index2'. the macro +# show.dpair generates a plot of the range vector vs (v1 - v2) / dv1. +# The indices for v1, dv1, and v2 are index1,1, and 2 respectively. The +# index for the range vector is defined by the integer following that +# vector. The y-limits of the plot are given by the last two numbers +# (use 'def') for the full default range of the delta vector +list dpairs + # v1 dv v2 range + X_PSF X_PSF_SIG Xraw Mraw 2 -10.0 10.0 + Y_PSF Y_PSF_SIG Yraw Mraw 2 -10.0 10.0 + PSF_INST_MAG PSF_INST_MAG_SIG Mraw Mraw 2 -10.0 10.0 + PSF_INST_FLUX PSF_INST_FLUX_SIG Fraw Mraw 2 -10.0 10.0 + AP_MAG PSF_INST_MAG_SIG Mraw Mraw 2 -10.0 10.0 + AP_MAG_RAW PSF_INST_MAG_SIG Mraw Mraw 2 -10.0 10.0 +end + +# if we run this test as a stand-alone program somewhere, we may need to create a local copy of the photcode file: +macro mkphotcodes + if ($0 != 2) + echo "USAGE: mkphotcodes (filename)" + break + end + + exec /bin/rm -f $1 + output $1 + echo "# airmass color astrometry mag photom astrom mask photom mask" + echo "# code name type zero slope offset c1 c2 slope zero equiv sys scale scale sys poor bad poor bad" + echo " 1 g_SYNTH sec 0.000 0.000 0.000 1 3 0.0000 0 21 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 2 r_SYNTH sec 0.000 0.000 0.000 2 3 0.0000 0 22 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3 i_SYNTH sec 0.000 0.000 0.000 2 3 0.0000 0 23 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 4 z_SYNTH sec 0.000 0.000 0.000 3 4 0.0000 0 24 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 5 y_SYNTH sec 0.000 0.000 0.000 4 5 0.0000 0 25 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3001 SYNTH.g ref 0.000 0.000 0.000 - - 0.0000 0 1 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3002 SYNTH.r ref 0.000 0.000 0.000 - - 0.0000 0 2 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3003 SYNTH.i ref 0.000 0.000 0.000 - - 0.0000 0 3 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3004 SYNTH.z ref 0.000 0.000 0.000 - - 0.0000 0 4 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3005 SYNTH.y ref 0.000 0.000 0.000 - - 0.0000 0 5 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + output stdout +end + +macro ckradialflux + data $1 + read -fits Chip.xrad X_APER Y_APER PSF_FWHM APER_FLUX APER_FLUX_ERR APER_FLUX_STDEV APER_FILL + read -fits Chip.psf X_PSF Y_PSF PSF_INST_MAG + set mag = PSF_INST_MAG + # XXX include nradii in header + for i 0 11 + set mA$i = -2.5*log(APER_FLUX:$i) + set dM$i = mA$i - PSF_INST_MAG + vstat dM$i + $DM$i = $MEDIAN + end + for i 1 11 + $j = $i - 1 + echo {$DM$i - $DM$j} + end + + lim mag -2 2; clear; box + for i 0 11 + plot mag dM$i -c black + end +end + +if ($SCRIPT) + fulltest 4 + exit 0 +end