Index: trunk/psphot/src/models/pmModel_STRAIL.c =================================================================== --- trunk/psphot/src/models/pmModel_STRAIL.c (revision 35688) +++ trunk/psphot/src/models/pmModel_STRAIL.c (revision 35769) @@ -563,5 +563,6 @@ // convert to shape terms (SXX,SYY,SXY) // XXX user-defined value for limit? - if (!pmPSF_FitToModel (PAR, 0.1)) { + bool useReff = pmModelUseReff (model->type); + if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); return false; Index: trunk/psphot/src/models/pmModel_TEST1.c =================================================================== --- trunk/psphot/src/models/pmModel_TEST1.c (revision 35688) +++ trunk/psphot/src/models/pmModel_TEST1.c (revision 35769) @@ -264,5 +264,6 @@ // convert to shape terms (SXX,SYY,SXY) // XXX user-defined value for limit? - if (!pmPSF_FitToModel (PAR, 0.1)) { + bool useReff = pmModelUseReff (model->type); + if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); return false; Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 35688) +++ trunk/psphot/src/psphot.h (revision 35769) @@ -67,4 +67,5 @@ bool psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index); +bool psphotModelBackground_Threaded (psThreadJob *job); bool psphotMaskBackground (pmConfig *config, const pmFPAview *view, const char *filerule); @@ -119,4 +120,5 @@ bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, const char *filerule, bool add); bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add); +bool psphotAddOrSubNoise_Threaded (psThreadJob *job); bool psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule); @@ -232,5 +234,5 @@ bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal); bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal); -pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal); +pmModel *psphotFitEXT (pmModel *guessModel, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal); psArray *psphotFitDBL (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal); @@ -243,4 +245,6 @@ bool psphotFitInit (int nThreads); bool psphotFitSummary (void); +bool psphotFitSummaryExtended (void); +bool psphotFitInitExtended (void); bool psphotLoadPSF (pmConfig *config, const pmFPAview *view, const char *filerule); Index: trunk/psphot/src/psphotAddNoise.c =================================================================== --- trunk/psphot/src/psphotAddNoise.c (revision 35688) +++ trunk/psphot/src/psphotAddNoise.c (revision 35769) @@ -1,4 +1,5 @@ # include "psphotInternal.h" +bool psphotAddOrSubNoise_Threaded (psThreadJob *job); bool psphotMaskSource(pmSource *source, bool add, psImageMaskType maskVal); @@ -32,10 +33,10 @@ } -static int Nmasked = 0; - // the return state indicates if any sources were actually replaced bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add) { bool status = false; + + psTimerStart ("psphot.noise"); // find the currently selected readout @@ -55,5 +56,9 @@ if (!sources) return false; - psTimerStart ("psphot.noise"); + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + if (!status) { + nThreads = 0; + } // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -84,4 +89,79 @@ psphotVisualShowImage (readout); + + // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) + int Cx = 1, Cy = 1; + psphotChooseCellSizes (&Cx, &Cy, readout, nThreads); + + psArray *cellGroups = psphotAssignSources (Cx, Cy, sources); + + for (int i = 0; i < cellGroups->n; i++) { + + psArray *cells = cellGroups->data[i]; + + for (int j = 0; j < cells->n; j++) { + + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_ADD_NOISE"); + psArrayAdd(job->args, 1, cells->data[j]); // sources + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, FACTOR, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, SIZE, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, add, PS_TYPE_U8); + +# if (1) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add/sub noise."); + return false; + } +# else + if (!psphotAddOrSubNoise_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add/sub noise."); + return false; + } + psFree(job); +# endif + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false, true)) { + psFree(cellGroups); + psError(PS_ERR_UNKNOWN, false, "Unable to add/sub noise."); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + // we have no returned data from this operation + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } + psFree(job); + } + } + + if (add) { + psLogMsg ("psphot.noise", PS_LOG_WARN, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise")); + } else { + psLogMsg ("psphot.noise", PS_LOG_WARN, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise")); + } + + psFree (cellGroups); + + psphotVisualShowImage (readout); + + return true; +} + +bool psphotAddOrSubNoise_Threaded (psThreadJob *job) { + + psArray *sources = job->args->data[0]; + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[1], PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[2], PS_TYPE_IMAGE_MASK_DATA); + psF32 FACTOR = PS_SCALAR_VALUE(job->args->data[3], F32); + psF32 SIZE = PS_SCALAR_VALUE(job->args->data[4], F32); + psBool add = PS_SCALAR_VALUE(job->args->data[5], U8); // loop over all source @@ -104,13 +184,4 @@ psphotMaskSource (source, add, markVal); } - if (add) { - psLogMsg ("psphot.noise", PS_LOG_WARN, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise")); - } else { - psLogMsg ("psphot.noise", PS_LOG_WARN, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise")); - } - fprintf (stderr, "masked %d objects\n", Nmasked); - - psphotVisualShowImage (readout); - return true; } @@ -142,7 +213,5 @@ } } - Nmasked ++; - - return true; -} - + return true; +} + Index: trunk/psphot/src/psphotBlendFit.c =================================================================== --- trunk/psphot/src/psphotBlendFit.c (revision 35688) +++ trunk/psphot/src/psphotBlendFit.c (revision 35769) @@ -86,4 +86,16 @@ if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) { fitMaxTol = 1.0; + } + + bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance + if (!status) { + // default to the old method (chisqConvergence) + chisqConvergence = true; + } + + int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance + if (!status) { + // default to the old method (chisqConvergence) + gainFactorMode = 0; } @@ -119,4 +131,7 @@ fitOptions->mode = PM_SOURCE_FIT_PSF; fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix + + fitOptions->gainFactorMode = gainFactorMode; + fitOptions->chisqConvergence = chisqConvergence; psphotInitLimitsPSF (recipe, readout); @@ -211,4 +226,5 @@ psLogMsg ("psphot.psphotBlendFit", PS_LOG_WARN, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot.fit.nonlinear"), Nfit, Npsf, Next, Nfail, sources->n - Nfit); + psphotFitSummary (); psphotVisualShowResidualImage (readout, false); @@ -271,32 +287,32 @@ // skip non-astronomical objects (very likely defects) - if (source->mode & PM_SOURCE_MODE_BLEND) continue; - if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; - if (source->type == PM_SOURCE_TYPE_DEFECT) continue; - if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + if (source->mode & PM_SOURCE_MODE_BLEND) goto skip_blend; + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) goto skip_cr; + if (source->type == PM_SOURCE_TYPE_DEFECT) goto skip_defect; + if (source->type == PM_SOURCE_TYPE_SATURATED) goto skip_sat; // skip saturated stars modeled with a radial profile - if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) goto skip_sat; // skip DBL second sources (ie, added by psphotFitBlob - if (source->mode & PM_SOURCE_MODE_PAIR) continue; + if (source->mode & PM_SOURCE_MODE_PAIR) goto skip_blend; // do not include MOMENTS_FAILURES in the fit - if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) goto skip_generic; // limit selection to some SN limit if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) continue; + if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) goto skip_generic; } else { - if (sqrt(source->peak->detValue) < FIT_SN_LIM) continue; + if (sqrt(source->peak->detValue) < FIT_SN_LIM) goto skip_generic; } // exclude sources outside optional analysis region - if (source->peak->xf < AnalysisRegion.x0) continue; - if (source->peak->yf < AnalysisRegion.y0) continue; - if (source->peak->xf > AnalysisRegion.x1) continue; - if (source->peak->yf > AnalysisRegion.y1) continue; + if (source->peak->xf < AnalysisRegion.x0) goto skip_generic; + if (source->peak->yf < AnalysisRegion.y0) goto skip_generic; + if (source->peak->xf > AnalysisRegion.x1) goto skip_generic; + if (source->peak->yf > AnalysisRegion.y1) goto skip_generic; // if model is NULL, we don't have a starting guess - if (source->modelPSF == NULL) continue; + if (source->modelPSF == NULL) goto skip_generic; // skip sources which are insignificant flux? @@ -307,5 +323,5 @@ source->modelPSF->params->data.F32[2], source->modelPSF->params->data.F32[3]); - continue; + goto skip_generic; } @@ -357,4 +373,21 @@ pmSourceCacheModel (source, maskVal); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + continue; + + skip_blend: + psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf); + continue; + skip_cr: + psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf); + continue; + skip_defect: + psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf); + continue; + skip_sat: + psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf); + continue; + skip_generic: + psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf); + continue; } Index: trunk/psphot/src/psphotChoosePSF.c =================================================================== --- trunk/psphot/src/psphotChoosePSF.c (revision 35688) +++ trunk/psphot/src/psphotChoosePSF.c (revision 35769) @@ -158,4 +158,16 @@ } float maxChisqDOF = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MAX_CHISQ"); // Fit tolerance + + bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance + if (!status) { + // default to the old method (chisqConvergence) + chisqConvergence = true; + } + + int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance + if (!status) { + // default to the old method (chisqConvergence) + gainFactorMode = 0; + } // options which modify the behavior of the model fitting @@ -169,6 +181,8 @@ options->fitOptions->mode = PM_SOURCE_FIT_PSF; options->fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix - + options->fitOptions->gainFactorMode = gainFactorMode; + options->fitOptions->chisqConvergence = chisqConvergence; + psArray *stars = psArrayAllocEmpty (sources->n); Index: trunk/psphot/src/psphotExtendedSourceFits.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceFits.c (revision 35688) +++ trunk/psphot/src/psphotExtendedSourceFits.c (revision 35769) @@ -1,3 +1,4 @@ # include "psphotInternal.h" +bool psphotSetWindowTrail (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float newRadius); // for now, let's store the detections on the readout->analysis for each readout @@ -52,4 +53,6 @@ psTimerStart ("psphot.extended"); + psphotFitInitExtended(); + // find the currently selected readout pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest @@ -70,4 +73,5 @@ if (!sources->n) { psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size"); + psphotSersicModelClassCleanup(); return true; } @@ -80,5 +84,5 @@ // do not thread if we are trying to study the fitting process if (psTraceGetLevel ("psphot.psphotFitEXT") >= 6) { - nThreads = 0; + nThreads = 0; } @@ -107,17 +111,93 @@ } + bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance + if (!status) { + // default to the old method (chisqConvergence) + chisqConvergence = true; + } + + int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance + if (!status) { + // default to the old method (chisqConvergence) + gainFactorMode = 0; + } + + // Define source fitting parameters for extended source fits + pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); + fitOptions->mode = PM_SOURCE_FIT_EXT; + fitOptions->saveCovariance = true; // XXX make this a user option? + fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix + fitOptions->nIter = fitIter; + fitOptions->minTol = fitMinTol; + fitOptions->maxTol = fitMaxTol; + + fitOptions->gainFactorMode = gainFactorMode; + fitOptions->chisqConvergence = chisqConvergence; + // maskVal is used to test for rejected pixels, and must include markVal maskVal |= markVal; // select the collection of desired models - psMetadata *models = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS"); + psMetadata *allModels = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS"); if (!status) { psWarning ("extended source model fits requested but model model is missing (EXTENDED_SOURCE_MODELS)\n"); + psFree (fitOptions); + psphotSersicModelClassCleanup(); return true; } + if (allModels->list->n == 0) { + psWarning ("extended source model fits requested but no models are specified\n"); + psFree (fitOptions); + psphotSersicModelClassCleanup(); + return true; + } + + psMetadata *models = NULL; + char *modelSelection = psMetadataLookupStr (&status, recipe, "EXTENDED_SOURCE_MODELS_SELECTION"); + if (!modelSelection || !status) { + models = psMemIncrRefCounter(allModels); + goto got_models; + } + + if (!strcasecmp(modelSelection, "ALL")) { + models = psMemIncrRefCounter(allModels); + goto got_models; + } + + if (!strcasecmp(modelSelection, "NONE")) { + psWarning ("extended source model fits requested but no models are selected (EXTENDED_SOURCE_MODELS_SELECTION = NONE)\n"); + psFree (fitOptions); + psphotSersicModelClassCleanup(); + return true; + } + + psArray *selection = psStringSplitArray (modelSelection, ",", false); + if (selection->n == 0) { + psWarning ("extended source model fits requested but model selection string is empty (EXTENDED_SOURCE_MODELS_SELECTION)\n"); + psFree (fitOptions); + psFree (selection); + psphotSersicModelClassCleanup(); + return true; + } + models = psMetadataAlloc(); + for (int i = 0; i < selection->n; i++) { + psMetadata *model = psMetadataLookupMetadata (&status, allModels, selection->data[i]); + if (!model) { + psWarning ("extended source model selection string includes an invalid name %s\n", (char *) selection->data[i]); + continue; + } + psMetadataAddMetadata (models, PS_LIST_TAIL, selection->data[i], PS_META_REPLACE, "extended model", model); + } + psFree (selection); + if (models->list->n == 0) { - psWarning ("extended source model fits requested but no models are specified\n"); - return true; - } + psWarning ("extended source model fits requested but no valid models in selection string (EXTENDED_SOURCE_MODELS_SELECTION = %s)\n", modelSelection); + psFree (fitOptions); + psFree (models); + psphotSersicModelClassCleanup(); + return true; + } + +got_models: psphotInitRadiusEXT (recipe, readout); @@ -128,34 +208,40 @@ while ((item = psMetadataGetAndIncrement (iter)) != NULL) { - if (item->type != PS_DATA_METADATA) { - // XXX we could cull the bad entries or build a validated model folder - psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name); - } - - psMetadata *model = (psMetadata *) item->data.md; - - // check on the model type - char *modelName = psMetadataLookupStr (&status, model, "MODEL"); - int modelType = pmModelClassGetType (modelName); - if (modelType < 0) { - psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName); - } - psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType); - - // check on the SNLIM, set a float value - char *SNword = psMetadataLookupStr (&status, model, "SNLIM"); - if (!status) { - psAbort("SNLIM not defined for extended source model %s\n", item->name); - } - float SNlim = atof (SNword); - psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim); - - // check on the PSF-Convolution status - char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED"); - if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) { - psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name); - } - bool convolved = !strcasecmp (convolvedWord, "true"); - psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved); + if (item->type != PS_DATA_METADATA) { + // XXX we could cull the bad entries or build a validated model folder + psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name); + } + + psMetadata *model = (psMetadata *) item->data.md; + + // check on the model type + char *modelName = psMetadataLookupStr (&status, model, "MODEL"); + int modelType = pmModelClassGetType (modelName); + if (modelType < 0) { + psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName); + } + psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType); + + // check on the SNLIM, set a float value + char *SNword = psMetadataLookupStr (&status, model, "SNLIM"); + if (!status) { + psAbort("SNLIM not defined for extended source model %s\n", item->name); + } + float SNlim = atof (SNword); + psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim); + + // check on the PSF-Convolution status + char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED"); + if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) { + psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name); + } + bool convolved = !strcasecmp (convolvedWord, "true"); + psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved); + + if (convolved) { + psLogMsg ("psphot", PS_LOG_INFO, "using convolved model class %s (%s)", modelName, item->name); + } else { + psLogMsg ("psphot", PS_LOG_INFO, "using simple model class %s (%s)", modelName, item->name); + } } psFree (iter); @@ -195,11 +281,8 @@ psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL); psArrayAdd(job->args, 1, iter); - psArrayAdd(job->args, 1, AnalysisRegion); // XXX make a pointer + psArrayAdd(job->args, 1, AnalysisRegion); + psArrayAdd(job->args, 1, fitOptions); PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, fitIter, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, fitMinTol, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR(job->args, fitMaxTol, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); @@ -218,4 +301,7 @@ psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); psFree(AnalysisRegion); + psFree (fitOptions); + psFree (models); + psphotSersicModelClassCleanup(); return false; } @@ -224,20 +310,23 @@ psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); psFree(AnalysisRegion); + psFree (fitOptions); + psFree (models); + psphotSersicModelClassCleanup(); return false; } psScalar *scalar = NULL; - scalar = job->args->data[8]; + scalar = job->args->data[9]; Next += scalar->data.S32; - scalar = job->args->data[9]; + scalar = job->args->data[10]; Nconvolve += scalar->data.S32; - scalar = job->args->data[10]; + scalar = job->args->data[11]; NconvolvePass += scalar->data.S32; - scalar = job->args->data[11]; + scalar = job->args->data[12]; Nplain += scalar->data.S32; - scalar = job->args->data[12]; + scalar = job->args->data[13]; NplainPass += scalar->data.S32; - scalar = job->args->data[13]; + scalar = job->args->data[14]; Nfaint += scalar->data.S32; - scalar = job->args->data[14]; + scalar = job->args->data[15]; Nfail += scalar->data.S32; psFree(job->args->data[3]); // iterator allocated above @@ -250,4 +339,7 @@ psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); psFree(AnalysisRegion); + psFree (fitOptions); + psFree (models); + psphotSersicModelClassCleanup(); return false; } @@ -260,17 +352,17 @@ } else { psScalar *scalar = NULL; - scalar = job->args->data[8]; + scalar = job->args->data[9]; Next += scalar->data.S32; - scalar = job->args->data[9]; + scalar = job->args->data[10]; Nconvolve += scalar->data.S32; - scalar = job->args->data[10]; + scalar = job->args->data[11]; NconvolvePass += scalar->data.S32; - scalar = job->args->data[11]; + scalar = job->args->data[12]; Nplain += scalar->data.S32; - scalar = job->args->data[12]; + scalar = job->args->data[13]; NplainPass += scalar->data.S32; - scalar = job->args->data[13]; + scalar = job->args->data[14]; Nfaint += scalar->data.S32; - scalar = job->args->data[14]; + scalar = job->args->data[15]; Nfail += scalar->data.S32; psFree(job->args->data[3]); // metadata iterator allocated above @@ -281,4 +373,6 @@ psFree (cellGroups); psFree(AnalysisRegion); + psFree (fitOptions); + psFree (models); psphotSersicModelClassCleanup(); @@ -290,4 +384,7 @@ psLogMsg ("psphot", PS_LOG_INFO, " %d plain models (%d passed)\n", Nplain, NplainPass); psLogMsg ("psphot", PS_LOG_INFO, " %d too faint to fit, %d failed\n", Nfaint, Nfail); + + psphotFitSummaryExtended(); + return true; } @@ -308,24 +405,16 @@ // arguments: readout, sources, models, region, psfSize, maskVal, markVal - pmReadout *readout = job->args->data[0]; - psArray *sources = job->args->data[1]; - psMetadata *models = job->args->data[2]; + pmReadout *readout = job->args->data[0]; + psArray *sources = job->args->data[1]; + psMetadata *models = job->args->data[2]; psMetadataIterator *iter = job->args->data[3]; - psRegion *region = job->args->data[4]; - int psfSize = PS_SCALAR_VALUE(job->args->data[5],S32); - int fitIter = PS_SCALAR_VALUE(job->args->data[6],S32); - float fitMinTol = PS_SCALAR_VALUE(job->args->data[7],F32); - float fitMaxTol = PS_SCALAR_VALUE(job->args->data[8],F32); - psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[9],PS_TYPE_IMAGE_MASK_DATA); - psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[10],PS_TYPE_IMAGE_MASK_DATA); - - // Define source fitting parameters for extended source fits - pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); - fitOptions->mode = PM_SOURCE_FIT_EXT; - fitOptions->saveCovariance = true; // XXX make this a user option? - fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix - fitOptions->nIter = fitIter; - fitOptions->minTol = fitMinTol; - fitOptions->maxTol = fitMaxTol; + psRegion *region = job->args->data[4]; + pmSourceFitOptions *fitOptions = job->args->data[5]; + + int psfSize = PS_SCALAR_VALUE(job->args->data[6],S32); + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[8],PS_TYPE_IMAGE_MASK_DATA); + + // psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5); // choose the sources of interest @@ -348,4 +437,10 @@ if (source->peak->y > region->y1) continue; + + // XXX for a test, just do the obvious trail + // XXX if (source->peak->xf < 1100) continue; + // XXX if (source->peak->xf > 1400) continue; + // XXX if (source->peak->yf > 245) continue; + // replace object in image if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { @@ -384,6 +479,4 @@ assert (status); - // limit selection to some SN limit - // assert (source->peak); // how can a source not have a peak? // limit selection to some SN limit bool skipSource = false; @@ -424,35 +517,49 @@ } } else { - psFree (source->modelFlux); - source->modelFlux = NULL; - modelFit = psphotFitEXT (readout, source, fitOptions, modelType, maskVal, markVal); - if (!modelFit) { - psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My); - Nfail ++; - continue; - } - psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", - source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter); - Nplain ++; - if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { - NplainPass ++; - source->mode |= PM_SOURCE_MODE_EXTENDED_FIT; - } + bool doneFits = false; + while (!doneFits) { + psFree (source->modelFlux); + source->modelFlux = NULL; + modelFit = psphotFitEXT (modelFit, readout, source, fitOptions, modelType, maskVal, markVal); + if (!modelFit) { + psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My); + Nfail ++; + doneFits = true; + continue; + } + psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter); + Nplain ++; + if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { + NplainPass ++; + source->mode |= PM_SOURCE_MODE_EXTENDED_FIT; + } + doneFits = true; + + if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) { + // calculate the object end points: + float *PAR = modelFit->params->data.F32; + float Xmin = PAR[PM_PAR_XPOS] - 0.5*PAR[PM_PAR_LENGTH]*sin(PAR[PM_PAR_THETA]); + float Ymin = PAR[PM_PAR_YPOS] - 0.5*PAR[PM_PAR_LENGTH]*cos(PAR[PM_PAR_THETA]); + float Xmax = PAR[PM_PAR_XPOS] + 0.5*PAR[PM_PAR_LENGTH]*sin(PAR[PM_PAR_THETA]); + float Ymax = PAR[PM_PAR_YPOS] + 0.5*PAR[PM_PAR_LENGTH]*cos(PAR[PM_PAR_THETA]); + + if (false && (source->peak->xf > 1100) && + (source->peak->xf < 1400) && + (source->peak->yf < 245)) { + fprintf (stderr, "src vs fit : %d %d - %d %d | %f %f - %f %f\n", + source->pixels->col0, source->pixels->row0, + source->pixels->col0 + source->pixels->numCols, source->pixels->row0 + source->pixels->numRows, + Xmin, Ymin, Xmax, Ymax); + } + if (PAR[PM_PAR_LENGTH] > 0.9*fitRadius) { + doneFits = false; + fprintf (stderr, "update window : %f %f : %f -> %f\n", source->peak->xf, source->peak->yf, fitRadius, 2*fitRadius); + psphotSetWindowTrail (&fitRadius, &windowRadius, readout, source, markVal, fitRadius*2.0); + } + } + } } - - // XXX deprecate? - // XXX pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc(); - // XXX extFitPars->Mxx = source->moments->Mxx; - // XXX extFitPars->Mxy = source->moments->Mxy; - // XXX extFitPars->Myy = source->moments->Myy; - // XXX extFitPars->Mrf = source->moments->Mrf; - // XXX extFitPars->Mrh = source->moments->Mrh; - // XXX extFitPars->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN; - - // save kron mag, but assign apMag & psfMag on output (not yet calculated) - // XXX extFitPars->krMag = -2.5*log10(source->moments->KronFlux); - - // psArrayAdd (source->extFitPars, 4, extFitPars); - // psFree (extFitPars); + // XXX really need to do this in a cleaner way: + if (!modelFit) continue; // test for fit quality / result @@ -519,28 +626,54 @@ psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); } - psFree (fitOptions); // change the value of a scalar on the array (wrap this and put it in psArray.h) + scalar = job->args->data[9]; + scalar->data.S32 = Next; + + scalar = job->args->data[10]; + scalar->data.S32 = Nconvolve; + scalar = job->args->data[11]; - scalar->data.S32 = Next; + scalar->data.S32 = NconvolvePass; scalar = job->args->data[12]; - scalar->data.S32 = Nconvolve; + scalar->data.S32 = Nplain; scalar = job->args->data[13]; - scalar->data.S32 = NconvolvePass; + scalar->data.S32 = NplainPass; scalar = job->args->data[14]; - scalar->data.S32 = Nplain; + scalar->data.S32 = Nfaint; scalar = job->args->data[15]; - scalar->data.S32 = NplainPass; - - scalar = job->args->data[16]; - scalar->data.S32 = Nfaint; - - scalar = job->args->data[17]; scalar->data.S32 = Nfail; return true; } + +# define PAD_WINDOW 3.0 +bool psphotSetWindowTrail (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float newRadius) { + + psRegion newRegion; + + psAssert (source, "source not defined??"); + psAssert (source->moments, "moments not defined??"); + + *fitRadius = newRadius; + *windowRadius = *fitRadius + PAD_WINDOW; + + // check to see if new region is completely contained within old region + newRegion = psRegionForSquare (source->peak->xf, source->peak->yf, *windowRadius); + newRegion = psRegionForImage (readout->image, newRegion); + + // redefine the pixels to match + pmSourceRedefinePixelsByRegion (source, readout, newRegion); + + // clear the circular mask + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + + // set the mask to flag the excluded pixels + psImageKeepCircle (source->maskObj, source->peak->xf, source->peak->yf, *fitRadius, "OR", markVal); + + return true; +} Index: trunk/psphot/src/psphotGuessModels.c =================================================================== --- trunk/psphot/src/psphotGuessModels.c (revision 35688) +++ trunk/psphot/src/psphotGuessModels.c (revision 35769) @@ -29,4 +29,6 @@ } +int NpixTotal = 0; + // construct an initial PSF model for each object (new sources only) bool psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index) { @@ -86,4 +88,6 @@ psArray *cellGroups = psphotAssignSources (Cx, Cy, sources); + + NpixTotal = 0; for (int i = 0; i < cellGroups->n; i++) { @@ -103,9 +107,18 @@ PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); +# if (1) if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); return false; } - } +# else + if (!psphotGuessModel_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } + psFree(job); +# endif + } + // wait for the threads to finish and manage results @@ -143,5 +156,5 @@ psFree (cellGroups); - psLogMsg ("psphot.models", PS_LOG_WARN, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.models")); + psLogMsg ("psphot.models", PS_LOG_WARN, "built models for %ld objects: %f sec (%d pixels)\n", sources->n, psTimerMark ("psphot.models"), NpixTotal); return true; } @@ -240,4 +253,5 @@ pmSourceCacheModel (source, maskVal); // ALLOC x14 (!) + NpixTotal += source->pixels->numCols * source->pixels->numRows; } Index: trunk/psphot/src/psphotKronIterate.c =================================================================== --- trunk/psphot/src/psphotKronIterate.c (revision 35688) +++ trunk/psphot/src/psphotKronIterate.c (revision 35769) @@ -166,18 +166,31 @@ } + fprintf (stderr, "--- starting KRON ---\n"); + // We measure the Kron Radius on a smoothed copy of the readout image psImage *smoothedImage = NULL; if (KRON_SMOOTH) { + psTimerStart ("psphot.kron.smooth"); + // Build the smoothed source image // Replace the subtracted sources psphotReplaceAllSourcesReadout(config, view, filerule, index, recipe, false); + psLogMsg ("psphot.kron", PS_LOG_INFO, "replaced %ld raw sources %f sec\n", sources->n, psTimerMark ("psphot.kron.smooth")); + + // smoothedImage = psImageCopy(NULL, readout->image, PS_TYPE_F32); + // psImageSmooth(smoothedImage, KRON_SMOOTH_SIGMA, KRON_SMOOTH_NSIGMA); + // Copy the image and smooth psTimerStart ("psphot.kron.smooth"); - smoothedImage = psImageCopy(NULL, readout->image, PS_TYPE_F32); - psImageSmooth(smoothedImage, KRON_SMOOTH_SIGMA, KRON_SMOOTH_NSIGMA); + bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading in psImageConvolve + smoothedImage = psImageSmoothNoMask_Threaded (NULL, readout->image, KRON_SMOOTH_SIGMA, KRON_SMOOTH_NSIGMA, 0.25); + psImageConvolveSetThreads(oldThreads); psLogMsg ("psphot.kron", PS_LOG_INFO, "smoothed image %f sec\n", psTimerMark ("psphot.kron.smooth")); // remove the sources + psTimerStart ("psphot.kron.smooth"); psphotRemoveAllSourcesReadout( config, view, filerule, index, recipe, false ); + psLogMsg ("psphot.kron", PS_LOG_INFO, "removed %ld raw sources %f sec\n", sources->n, psTimerMark ("psphot.kron.smooth")); + // Now subtract smooth versions of the sources from the smoothed image psTimerStart ("psphot.kron.smooth.sources"); Index: trunk/psphot/src/psphotModelBackground.c =================================================================== --- trunk/psphot/src/psphotModelBackground.c (revision 35688) +++ trunk/psphot/src/psphotModelBackground.c (revision 35769) @@ -28,4 +28,6 @@ } +// track background cell failures for log reporting purposes +static int nFailures = 0; // Generate the background model @@ -33,5 +35,5 @@ // linear interpolation to generate full-scale model // -// NOTE that the 'analysis' metedata pass in here is used to store the binning information. +// NOTE that the 'analysis' metedata passed in here is used to store the binning information. // This may be the analysis for this readout, but it may be the analysis for the pmFPAfile // corresponding to the model. Other information about the background model is saved on the @@ -109,38 +111,33 @@ const psStatsOptions statsOption = statsOptionLocation | statsOptionWidth; - psStats *stats = psStatsAlloc (statsOption); + psStats *statsDefaults = psStatsAlloc (statsOption); // set range for old-version of sky statistic if (statsOptionLocation & PS_STAT_ROBUST_QUARTILE) { - stats->min = 0.25; - stats->max = 0.75; + statsDefaults->min = 0.25; + statsDefaults->max = 0.75; } // set user-option for number of pixels per region - stats->nSubsample = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX"); + statsDefaults->nSubsample = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX"); if (!status) { - stats->nSubsample = 1000; + statsDefaults->nSubsample = 1000; } // optionally set the binsize - stats->binsize = psMetadataLookupF32 (&status, recipe, "SKY_HISTOGRAM_BINS"); + statsDefaults->binsize = psMetadataLookupF32 (&status, recipe, "SKY_HISTOGRAM_BINS"); if (status) { - stats->options |= PS_STAT_USE_BINSIZE; + statsDefaults->options |= PS_STAT_USE_BINSIZE; } // optionally set the sigma clipping - stats->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA"); + statsDefaults->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA"); if (!status) { - if (stats->options & PS_STAT_FITTED_MEAN) { - stats->clipSigma = 1.0; + if (statsDefaults->options & PS_STAT_FITTED_MEAN) { + statsDefaults->clipSigma = 1.0; } else { - stats->clipSigma = 3.0; + statsDefaults->clipSigma = 3.0; } } - - // stats is not initialized by psStats??? use this to save the input options - // XXX re-work this to use the new psStatsInit function - psStats *statsDefaults = psStatsAlloc (statsOption); - *statsDefaults = *stats; // we save the binning structure for use in psphotMagnitudes @@ -154,7 +151,6 @@ // XXXX we can thread this here by running blocks in parallel - - int nFailures = 0; psImageBackgroundInit(); + nFailures = 0; // reset for this pass // we have Nx * Ny model points, but we can use a window which is larger (or smaller) than @@ -163,95 +159,71 @@ // measure clipped median for subimages - 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++) { - - // convert the ruff grid cell to the equivalent fine grid cell - ruffRegion = psRegionSet (ix + 0.5 - 0.5*dXsample, ix + 0.5 + 0.5*dXsample, - iy + 0.5 - 0.5*dYsample, iy + 0.5 + 0.5*dYsample); - fineRegion = psImageBinningSetFineRegion (binning, ruffRegion); - fineRegion = psRegionForImage (image, fineRegion); - - psImage *subset = psImageSubset (image, fineRegion); - if (!subset->numCols || !subset->numRows) { - psFree (subset); - continue; + + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_MODEL_BACKGROUND"); + + psArrayAdd(job->args, 1, image); + psArrayAdd(job->args, 1, mask); + psArrayAdd(job->args, 1, binning); + psArrayAdd(job->args, 1, rng); + psArrayAdd(job->args, 1, statsDefaults); + + psArrayAdd(job->args, 1, modelData); + psArrayAdd(job->args, 1, modelStdevData); + + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + + PS_ARRAY_ADD_SCALAR(job->args, ix, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, iy, PS_TYPE_S32); + + PS_ARRAY_ADD_SCALAR(job->args, dXsample, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, dYsample, PS_TYPE_F32); + + PS_ARRAY_ADD_SCALAR(job->args, statsOptionLocation, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, statsOptionWidth, PS_TYPE_S32); + + PS_ARRAY_ADD_SCALAR(job->args, NAN, PS_TYPE_F32); // this is used as a return value for dQvalue + +# if (1) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return NULL; } - psImage *submask = psImageSubset (mask, fineRegion); - - // reset the default values - *stats = *statsDefaults; - - // Use the selected background statistic for the first pass - // If it fails, fall back on the "ROBUST_MEDIAN" version - // If both fail, set the pixel to NAN and (below) interpolate - // XXX psImageBackground will probably be renamed psImageStats - // XXX don't bother trying if there are no valid pixels... - - psVector *sample = NULL; - - // turn on stats tracing in desired cells - # if (0) - psMetadata *plots = psMetadataLookupPtr (&status, recipe, "DIAGNOSTIC.PLOTS"); - assert (plots); - - int xPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.X"); - assert (status); - int yPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.Y"); - assert (status); - - bool gotX = (xPlot < 0) || (xPlot == ix); - bool gotY = (yPlot < 0) || (yPlot == iy); - - if (gotX && gotY) { - (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 6); - (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 6); - } else { - (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 0); - (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 0); +# else + if (!psphotModelBackground_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Failure to model background."); + return NULL; } - # endif - - if (psImageBackground(stats, &sample, subset, submask, maskVal, rng)) { - if (stats->options & PS_STAT_ROBUST_QUARTILE) { - modelData[iy][ix] = stats->robustMedian; - } else { - modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation); - } - modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth); - // fprintf (stderr, "dQ : %f - %f - %f = %f\n", stats->robustLQ, stats->robustMedian, stats->robustUQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian); - psVectorAppend (dQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian); - - // supply sample to plotting routing - psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample); - } else { - psStatsOptions currentOptions = stats->options; - stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV; - if (!psImageBackground(stats, &sample, subset, submask, maskVal, rng)) { - if ((nFailures < 3) || (nFailures % 100 == 0)) { - psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for " - "(%dx%d, (row0,col0) = (%d,%d)", - subset->numRows, subset->numCols, subset->row0, subset->col0); - } - nFailures ++; - modelData[iy][ix] = modelStdevData[iy][ix] = NAN; - } else { - modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN); - modelStdevData[iy][ix] = psStatsGetValue(stats, PS_STAT_ROBUST_STDEV); - - // supply sample to plotting routing - psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample); - } - // drop errors caused by psImageBackground failures - // XXX we probably should trap and exit on serious failures - psErrorClear(); - stats->options = currentOptions; - } - psFree(sample); - modelData[iy][ix] += SKY_BIAS; - psFree (subset); - psFree (submask); + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = job->args->data[14]; + float dQvalue = scalar->data.F32; + psVectorAppend (dQ, dQvalue); + } + psFree(job); +# endif } + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return NULL; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = job->args->data[14]; + float dQvalue = scalar->data.F32; + psVectorAppend (dQ, dQvalue); + } + psFree(job); } @@ -302,5 +274,4 @@ if (Count == 0) { psError (PSPHOT_ERR_DATA, true, "failed to build background image"); - psFree(stats); psFree(statsDefaults); psFree(binning); @@ -318,4 +289,11 @@ modelData[iy][ix] = Value; modelStdevData[iy][ix] = ValueStdev; + } + } + + // apply artificial offset if desired + for (int iy = 0; iy < model->numRows; iy++) { + for (int ix = 0; ix < model->numCols; ix++) { + modelData[iy][ix] += SKY_BIAS; } } @@ -364,5 +342,4 @@ psFree(dQ); - psFree(stats); psFree(statsDefaults); psFree(binning); @@ -433,2 +410,127 @@ return true; } + +bool psphotModelBackground_Threaded (psThreadJob *job) { + + psImage *image = job->args->data[0]; + psImage *mask = job->args->data[1]; + + psImageBinning *binning = job->args->data[2]; + + psRandom *rng = job->args->data[3]; + psStats *statsDefaults = job->args->data[4]; + + psF32 **modelData = job->args->data[5]; + psF32 **modelStdevData = job->args->data[6]; + + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA); + + int ix = PS_SCALAR_VALUE(job->args->data[8],S32); + int iy = PS_SCALAR_VALUE(job->args->data[9],S32); + + float dXsample = PS_SCALAR_VALUE(job->args->data[10],F32); + float dYsample = PS_SCALAR_VALUE(job->args->data[11],F32); + + psStatsOptions statsOptionLocation = PS_SCALAR_VALUE(job->args->data[12],S32); + psStatsOptions statsOptionWidth = PS_SCALAR_VALUE(job->args->data[13],S32); + + // convert the ruff grid cell to the equivalent fine grid cell + psRegion ruffRegion = psRegionSet (ix + 0.5 - 0.5*dXsample, ix + 0.5 + 0.5*dXsample, iy + 0.5 - 0.5*dYsample, iy + 0.5 + 0.5*dYsample); + psRegion fineRegion = psImageBinningSetFineRegion (binning, ruffRegion); + fineRegion = psRegionForImage (image, fineRegion); + + psImage *subset = psImageSubset (image, fineRegion); + if (!subset->numCols || !subset->numRows) { + psFree (subset); + return false; // XXX do we / should we fail on this? + } + psImage *submask = psImageSubset (mask, fineRegion); + + psStats *stats = psStatsAlloc (PS_STAT_NONE); + *stats = *statsDefaults; + + // Use the selected background statistic for the first pass + // If it fails, fall back on the "ROBUST_MEDIAN" version + // If both fail, set the pixel to NAN and (later) interpolate + + psVector *sample = NULL; + float dQvalue = NAN; + + if (psImageBackground(stats, &sample, subset, submask, maskVal, rng)) { + if (stats->options & PS_STAT_ROBUST_QUARTILE) { + modelData[iy][ix] = stats->robustMedian; + } else { + modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation); + } + modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth); + + // fprintf (stderr, "dQ : %f - %f - %f = %f\n", stats->robustLQ, stats->robustMedian, stats->robustUQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian); + // XXX this operation is not thread safe -- move out somehow... + // psVectorAppend (dQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian); + dQvalue = stats->robustUQ + stats->robustLQ - 2*stats->robustMedian; + // return dQvalue to main thread + + // supply sample to plotting routing + // only allow in a non-threaded context -- can plot more than one cell + // psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample); + } else { + // psStatsOptions currentOptions = stats->options; + stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV; + if (!psImageBackground(stats, &sample, subset, submask, maskVal, rng)) { + if ((nFailures < 3) || (nFailures % 100 == 0)) { + psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for " + "(%dx%d, (row0,col0) = (%d,%d)", + subset->numRows, subset->numCols, subset->row0, subset->col0); + } + nFailures ++; // static + modelData[iy][ix] = modelStdevData[iy][ix] = NAN; + } else { + modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN); + modelStdevData[iy][ix] = psStatsGetValue(stats, PS_STAT_ROBUST_STDEV); + + // supply sample to plotting routing + // only allow in a non-threaded context -- can plot more than one cell + // psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample); + } + // drop errors caused by psImageBackground failures + // NOTE : psStats raises errors when it cannot find a solution; this + // probably should not be errors but exit stats info only, but c'est la code + // XXX we probably should trap and exit on serious failures + psErrorClear(); + // stats->options = currentOptions; // this is not needed (set by init above) + } + psFree (stats); + psFree (sample); + psFree (subset); + psFree (submask); + + // return the dQvalue to the calling thread + psScalar *scalar = job->args->data[14]; + scalar->data.F32 = dQvalue; + + return true; +} + +// code for in-line plotting from above +// turn on stats tracing in desired cells +// XXX this code may need to be re-worked for a threaded context +# if (0) + psMetadata *plots = psMetadataLookupPtr (&status, recipe, "DIAGNOSTIC.PLOTS"); + assert (plots); + + int xPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.X"); + assert (status); + int yPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.Y"); + assert (status); + + bool gotX = (xPlot < 0) || (xPlot == ix); + bool gotY = (yPlot < 0) || (yPlot == iy); + + if (gotX && gotY) { + (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 6); + (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 6); + } else { + (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 0); + (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 0); + } +# endif Index: trunk/psphot/src/psphotModelGroupInit.c =================================================================== --- trunk/psphot/src/psphotModelGroupInit.c (revision 35688) +++ trunk/psphot/src/psphotModelGroupInit.c (revision 35769) @@ -4,10 +4,10 @@ // psModule/src/objects/models -# include "models/pmModel_TEST1.c" -# include "models/pmModel_STRAIL.c" +// # include "models/pmModel_TEST1.c" +// # include "models/pmModel_STRAIL.c" static pmModelClass userModels[] = { - {"PS_MODEL_TEST1", 7, pmModelFunc_TEST1, pmModelFlux_TEST1, pmModelRadius_TEST1, pmModelLimits_TEST1, pmModelGuess_TEST1, pmModelFromPSF_TEST1, pmModelParamsFromPSF_TEST1, pmModelFitStatus_TEST1, NULL}, - {"PS_MODEL_STRAIL", 9, pmModelFunc_STRAIL, pmModelFlux_STRAIL, pmModelRadius_STRAIL, pmModelLimits_STRAIL, pmModelGuess_STRAIL, pmModelFromPSF_STRAIL, pmModelParamsFromPSF_STRAIL, pmModelFitStatus_STRAIL, NULL}, + // {"PS_MODEL_TEST1", 7, pmModelFunc_TEST1, pmModelFlux_TEST1, pmModelRadius_TEST1, pmModelLimits_TEST1, pmModelGuess_TEST1, pmModelFromPSF_TEST1, pmModelParamsFromPSF_TEST1, pmModelFitStatus_TEST1, NULL}, + // {"PS_MODEL_STRAIL", 9, pmModelFunc_STRAIL, pmModelFlux_STRAIL, pmModelRadius_STRAIL, pmModelLimits_STRAIL, pmModelGuess_STRAIL, pmModelFromPSF_STRAIL, pmModelParamsFromPSF_STRAIL, pmModelFitStatus_STRAIL, NULL}, }; Index: trunk/psphot/src/psphotModelTestReadout.c =================================================================== --- trunk/psphot/src/psphotModelTestReadout.c (revision 35688) +++ trunk/psphot/src/psphotModelTestReadout.c (revision 35769) @@ -161,4 +161,9 @@ // get the initial model parameter guess pmModel *model = pmSourceModelGuess (source, modelType, maskVal, markVal); + if (!model) { + fprintf (stderr, "failed to generate model guess\n"); + exit (2); + } + source->modelEXT = model; @@ -186,11 +191,7 @@ fprintf (stderr, "guess: %f @ (%f, %f)\n", params[6]*180/M_PI, params[4], params[5]); } else { - // list model input shape - psEllipseShape shape; - shape.sx = 1.4 / model->params->data.F32[4]; - shape.sy = 1.4 / model->params->data.F32[5]; - shape.sxy = model->params->data.F32[6]; - axes = psEllipseShapeToAxes (shape, 20.0); - fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor); + bool useReff = pmModelUseReff (modelType); + pmModelParamsToAxes (&axes, params[PM_PAR_SXX], params[PM_PAR_SXY], params[PM_PAR_SYY], useReff); + fprintf (stderr, "guess: %f @ (%f, %f) : %f\n", axes.theta*180/M_PI, axes.major, axes.minor, params[PM_PAR_SXY]); } @@ -217,4 +218,12 @@ fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix + bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance + if (!status) chisqConvergence = true; + fitOptions->chisqConvergence = chisqConvergence; + + int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance + if (!status) gainFactorMode = 0; + fitOptions->gainFactorMode = gainFactorMode; + if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { fitOptions->mode = PM_SOURCE_FIT_NO_INDEX; @@ -231,4 +240,18 @@ pmSourceModelGuessPCM (pcm, source, maskVal, markVal); } + + // if we provide a test guess value we want to use that value! + for (int i = 0; i < nParams; i++) { + if (i == PM_PAR_XPOS) continue; + if (i == PM_PAR_YPOS) continue; + + char name[32]; + sprintf (name, "TEST_FIT_PAR%d", i); + float value = psMetadataLookupF32 (&status, recipe, name); + if (status && isfinite (value)) { + params[i] = value; + } + } + pmPCMupdate(pcm, source, fitOptions, model); pmSourceFitPCM (pcm, source, fitOptions, maskVal, markVal, psfSize); @@ -258,4 +281,12 @@ } + if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) { + fprintf (stderr, "result: %f @ (%f, %f)\n", params[6]*180/M_PI, params[4], params[5]); + } else { + bool useReff = pmModelUseReff (modelType); + pmModelParamsToAxes (&axes, params[PM_PAR_SXX], params[PM_PAR_SXY], params[PM_PAR_SYY], useReff); + fprintf (stderr, "result: %f @ (%f, %f) : %f\n", axes.theta*180/M_PI, axes.major, axes.minor, params[PM_PAR_SXY]); + } + // write out psphotSaveImage (NULL, source->pixels, "resid.fits"); Index: trunk/psphot/src/psphotOutput.c =================================================================== --- trunk/psphot/src/psphotOutput.c (revision 35688) +++ trunk/psphot/src/psphotOutput.c (revision 35769) @@ -379,17 +379,15 @@ psF32 *outPar = source->modelEXT->params->data.F32; - psEllipseShape shape; - - shape.sx = outPar[PM_PAR_SXX] / M_SQRT2; - shape.sy = outPar[PM_PAR_SYY] / M_SQRT2; - shape.sxy = outPar[PM_PAR_SXY]; - - psEllipsePol pol = pmPSF_ModelToFit (outPar); + bool useReff = pmModelUseReff (source->modelEXT->type); + + psEllipseAxes axes1; + pmModelParamsToAxes (&axes1, outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY], useReff); + + psEllipsePol pol = pmPSF_ModelToFit (outPar, useReff); inPar[PM_PAR_E0] = pol.e0; inPar[PM_PAR_E1] = pol.e1; inPar[PM_PAR_E2] = pol.e2; - pmPSF_FitToModel (inPar, 0.1); - - psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0); + pmPSF_FitToModel (inPar, 0.1, useReff); + psEllipseAxes axes2 = psEllipsePolToAxes(pol, 0.1); Index: trunk/psphot/src/psphotRadialPlot.c =================================================================== --- trunk/psphot/src/psphotRadialPlot.c (revision 35688) +++ trunk/psphot/src/psphotRadialPlot.c (revision 35769) @@ -15,5 +15,5 @@ if (nCount > 00) { if (*kapa != 0) { - KiiClose (*kapa); + KapaClose (*kapa); *kapa = 0; } Index: trunk/psphot/src/psphotReplaceUnfit.c =================================================================== --- trunk/psphot/src/psphotReplaceUnfit.c (revision 35688) +++ trunk/psphot/src/psphotReplaceUnfit.c (revision 35769) @@ -77,4 +77,6 @@ // maskVal is used to test for rejected pixels, and must include markVal maskVal |= markVal; + + int NpixTotal = 0; for (int i = 0; i < sources->n; i++) { @@ -99,8 +101,9 @@ pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + NpixTotal += source->pixels->numCols * source->pixels->numRows; } psphotVisualShowImage(readout); - psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace")); + psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec (%d pixels)\n", sources->n, psTimerMark ("psphot.replace"), NpixTotal); return true; } Index: trunk/psphot/src/psphotSetThreads.c =================================================================== --- trunk/psphot/src/psphotSetThreads.c (revision 35688) +++ trunk/psphot/src/psphotSetThreads.c (revision 35769) @@ -5,6 +5,18 @@ psThreadTask *task = NULL; + pmPSFThreads (); + + task = psThreadTaskAlloc("PSPHOT_MODEL_BACKGROUND", 15); + task->function = &psphotModelBackground_Threaded; + psThreadTaskAdd(task); + psFree(task); + task = psThreadTaskAlloc("PSPHOT_GUESS_MODEL", 5); task->function = &psphotGuessModel_Threaded; + psThreadTaskAdd(task); + psFree(task); + + task = psThreadTaskAlloc("PSPHOT_ADD_NOISE", 6); + task->function = &psphotAddOrSubNoise_Threaded; psThreadTaskAdd(task); psFree(task); @@ -40,5 +52,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 18); + task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 16); task->function = &psphotExtendedSourceFits_Threaded; psThreadTaskAdd(task); Index: trunk/psphot/src/psphotSourceFits.c =================================================================== --- trunk/psphot/src/psphotSourceFits.c (revision 35688) +++ trunk/psphot/src/psphotSourceFits.c (revision 35769) @@ -3,9 +3,20 @@ // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful +static int NfitBlend = 0; + static int NfitPSF = 0; -static int NfitBlend = 0; +static int NfitIterPSF = 0; + static int NfitDBL = 0; +static int NfitIterDBL = 0; +static int NfitPixDBL = 0; + static int NfitEXT = 0; +static int NfitIterEXT = 0; +static int NfitPixEXT = 0; + static int NfitPCM = 0; +static int NfitIterPCM = 0; +static int NfitPixPCM = 0; bool psphotFitInit (int nThreads) { @@ -17,6 +28,23 @@ bool psphotFitSummary () { - psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n", - NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits")); + psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d psf (%d iter), %5d blend: %5d ext (%d iter, %d pix), %5d dbl (%d iter, %d pix) : %6.2f sec\n", + NfitPSF, NfitIterPSF, NfitBlend, NfitEXT, NfitIterEXT, NfitPixEXT, NfitDBL, NfitIterDBL, NfitPixDBL, psTimerMark ("psphot.fits")); + return true; +} + +bool psphotFitSummaryExtended () { + + psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d ext (%d iter, %d pix), %5d pcm (%d iter, %d pix)\n", + NfitEXT, NfitIterEXT, NfitPixEXT, NfitPCM, NfitIterPCM, NfitPixPCM); + return true; +} + +bool psphotFitInitExtended () { + NfitEXT = 0; + NfitIterEXT = 0; + NfitPixEXT = 0; + NfitPCM = 0; + NfitIterPCM = 0; + NfitPixPCM = 0; return true; } @@ -169,4 +197,5 @@ options.mode = PM_SOURCE_FIT_PSF; pmSourceFitModel (source, PSF, &options, maskVal); + NfitIterPSF += PSF->nIter; if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); @@ -268,5 +297,5 @@ { // XXX need to handle failures better here - EXT = psphotFitEXT (readout, source, fitOptions, modelTypeEXT, maskVal, markVal); + EXT = psphotFitEXT (NULL, readout, source, fitOptions, modelTypeEXT, maskVal, markVal); if (!EXT) goto escape; if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) goto escape; @@ -440,8 +469,11 @@ options.mode = PM_SOURCE_FIT_PSF; pmSourceFitSet (source, modelSet, &options, maskVal); + NfitIterDBL += DBL->nIter; + NfitPixDBL += DBL->nDOF; + return (modelSet); } -pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) { +pmModel *psphotFitEXT (pmModel *guessModel, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) { if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) { @@ -457,5 +489,5 @@ // use the source moments, etc to guess basic model parameters - pmModel *model = pmSourceModelGuess (source, modelType, maskVal, markVal); + pmModel *model = guessModel ? guessModel : pmSourceModelGuess (source, modelType, maskVal, markVal); if (!model) { psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy); @@ -507,4 +539,6 @@ pmSourceFitModel (source, model, &options, maskVal); + NfitIterEXT += model->nIter; + NfitPixEXT += model->nDOF; # if (PS_TRACE_ON) @@ -587,4 +621,6 @@ // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); // NOTE : 1687 allocs in here + NfitIterPCM += pcm->modelConv->nIter; + NfitPixPCM += pcm->modelConv->nDOF; if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); } Index: trunk/psphot/src/psphotVisual.c =================================================================== --- trunk/psphot/src/psphotVisual.c (revision 35688) +++ trunk/psphot/src/psphotVisual.c (revision 35769) @@ -29,7 +29,7 @@ bool psphotVisualClose(void) { - if(kapa1 != -1) KiiClose(kapa1); - if(kapa2 != -1) KiiClose(kapa2); - if(kapa3 != -1) KiiClose(kapa3); + if(kapa1 != -1) KapaClose(kapa1); + if(kapa2 != -1) KapaClose(kapa2); + if(kapa3 != -1) KapaClose(kapa3); return true; } Index: trunk/psphot/test/tap_psphot_galaxies.pro =================================================================== --- trunk/psphot/test/tap_psphot_galaxies.pro (revision 35769) +++ trunk/psphot/test/tap_psphot_galaxies.pro (revision 35769) @@ -0,0 +1,600 @@ +#!/usr/bin/env mana +# -*-sh-*- + +# 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 T + +# basic options for the these images (filter, location, obstype) +$BaseOptions = -type OBJECT +$BaseOptions = $BaseOptions -filter r +$BaseOptions = $BaseOptions -skymags 20.86 +$BaseOptions = $BaseOptions -ra 270.70 +$BaseOptions = $BaseOptions -dec -23.70 +$BaseOptions = $BaseOptions -pa 0.0 +# $BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0 + +# options for the reference image +$RefOptions = $BaseOptions +$RefOptions = $RefOptions -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 + +# options for the simulated images (using the refimage for the stars) +$FakeOptions = $BaseOptions +$FakeOptions = $FakeOptions -exptime 30.0 +$FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS + +# sample alternate options: +# $FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1 +# $FakeOptions = $FakeOptions -Df PSF.ARATIO 1.2 +# $FakeOptions = $FakeOptions -Df PSF.THETA +30.0 + +# create an image with fake sources (these are then inserted into the catdir) +$RefConfig = -camera SIMTEST +$RefConfig = $RefConfig -recipe PPSIM STACKTEST.MAKE + +# basic config for ppSim with randomly distributed stars and gridded galaxies +$FakeConfig = -camera SIMTEST +$FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN +$FakeConfig = $FakeConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref +$FakeConfig = $FakeConfig -Db STARS.FAKE F ; # only use stars from catdir.ref +$FakeConfig = $FakeConfig -Db MATCH.DENSITY F +$FakeConfig = $FakeConfig -Db PSF.CONVOLVE T +$FakeConfig = $FakeConfig -Db GALAXY.FAKE T ; # generate a "realistic" distribution of galaxies +$FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0 +$FakeConfig = $FakeConfig -Db GALAXY.GRID T ; # generate a grid of galaxies (constant mag) +#$FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_GAUSS +#$FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_EXP +#$FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_SERSIC +#$FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_DEV +$FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN 10.0 +$FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX 10.0 +$FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN 0.25 +$FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX 0.25 +$FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0 +$FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180 +$FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.66 +$FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.66 +$FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 120 +$FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 120 + +list fwhm + 1.0 + 1.1 + 1.2 + 1.5 +end + +macro go + mkexp test.exp 1.0 EXP + fitexp test.exp test.exp.fit EXP_CONV + + mkexp test.ser 1.0 SERSIC + fitexp test.ser test.ser.fit SER_CONV + + mkexp test.dev 1.0 DEV + fitexp test.dev test.dev.fit DEV_CONV + + mkexp test.gau 1.0 GAUSS + fitexp test.gau test.gau.fit GAU_CONV + + mkexp test.pg 1.0 PGAUSS + fitexp test.pg test.pg.fit PGA_CONV + + mkexp test.qga 1.0 QGAUSS + fitexp test.qga test.qga.fit QGA_CONV + + mkexp test.p1 1.0 PS1_V1 + fitexp test.p1 test.p1.fit PS1_CONV +end + +macro go.ckgalaxy + foreach type exp ser dev gau pg qga p1 + ckgalaxy test.$type.dat test.$type.fit.cmf + wait $type + end +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 + + echo ppSim $RefOptions $RefConfig refimage + exec ppSim $RefOptions $RefConfig refimage + + 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 -quick-airmass + 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 != 4) + echo "USAGE: mkexp basename (fwhm) (model)" + break + end + + local fwhm basename + $basename = $1 + $fwhm = $2 + + # create the raw image + $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_$3 + echo ppSim $FakeOptions $FakeConfig $basename -seeing $fwhm + exec ppSim $FakeOptions $FakeConfig $basename -seeing $fwhm + 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 + +# create a realistic distribution of fake stars, GAUSS PSF +macro fitexp + if ($0 != 4) + echo "USAGE: fitexp basename outname (fitModel)" + break + end + + local basename fitModel outname + $basename = $1 + $outname = $2 + $fitModel = $3 + + $psphotConfig = + $psphotConfig = $psphotConfig -recipe PSPHOT GALAXY_TEST + $psphotConfig = $psphotConfig -threads 4 + $psphotConfig = $psphotConfig -Db PSPHOT:LMM_FIT_CHISQ_CONVERGENCE F + $psphotConfig = $psphotConfig -Df PSPHOT:EXT_FIT_MIN_TOL 0.1 + $psphotConfig = $psphotConfig -Di PSPHOT:LMM_FIT_GAIN_FACTOR_MODE 2 + $psphotConfig = $psphotConfig -Db PSPHOT:SAVE.RESID T + $psphotConfig = $psphotConfig -D PSPHOT:EXTENDED_SOURCE_MODELS_SELECTION $fitModel + + # ppImage / psphot on the output + echo ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname + exec ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname +end + +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 ckgalaxy + if ($0 != 3) + echo "USAGE: ckgalaxy (dat) (cmf)" + break + end + + data $1 + read Xin_all 1 Yin_all 2 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9 + + subset Xin = Xin_all if (Type == 1) + subset Yin = Yin_all if (Type == 1) + subset Min = Min_all if (Type == 1) + subset Tin = ThetaIn_all if (Type == 1) + + subset RmajIn = RmajIn_all if (Type == 1) + subset RminIn = RminIn_all if (Type == 1) + + data $2 + read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA + set EXT_THETA_ALT = EXT_THETA * (EXT_THETA >= 0.0) + (EXT_THETA + 3.14159265) * (EXT_THETA < 0.0) + set EXT_THETA = EXT_THETA_ALT + + match2d X_EXT Y_EXT Xin Yin 1.0 -index1 index1 -index2 index2 + + reindex Xot_m = X_EXT using index1 + reindex Yot_m = Y_EXT using index1 + + reindex Xin_m = Xin using index2 + reindex Yin_m = Yin using index2 + + set dX = Xin_m - Xot_m + set dY = Yin_m - Yot_m + # lim -n 0 dX dY; clear; box; plot dX dY + + reindex Mot_m = EXT_INST_MAG using index1 + reindex Tot_m = EXT_THETA using index1 + + reindex Min_m = Min using index2 + reindex Tin_m = Tin using index2 + + reindex Rot_m = EXT_WIDTH_MAJ using index1 + reindex rot_m = EXT_WIDTH_MIN using index1 + + reindex Rin_m = RmajIn using index2 + reindex rin_m = RminIn using index2 + + set n = ramp(Xin_m) + + set dM = Min_m - Mot_m + set dT = (Tin_m - Tot_m) * 180 / 3.14159265 + set Tin_deg = Tin_m * 180 / 3.14159265 + + lim -n 0 Tin_deg -5 5; clear; box; plot Tin_deg dT; label -y "delta theta (deg)" + lim -n 1 n dM; clear; box; plot n dM; label -y "delta mag" + lim -n 2 n -5 5; clear; box; plot n dT; label -y "delta theta (deg)" + lim -n 3 -5 5 dM; clear; box; plot dT dM; label -x "delta theta (deg)" -y "delta mag" + + set dR = Rin_m - Rot_m + set dr = rin_m - rot_m + lim -n 4 n dR; clear; box; plot n dR; label -y "delta Rmaj" + lim -n 5 n dr; clear; box; plot n dr; label -y "delta Rmin" +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 + +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 plot.angles + if ($0 != 2) + echo "USAGE: plot.angles (file.dat)" + break + end + + data $1 + read x 1 y 2 type 4 trad 9 + set t = trad * 180 / 3.14159265 + subset xs = x if (type == 1) + subset ys = y if (type == 1) + subset ts = t if (type == 1) + + set cs = dcos(ts) + set sn = dsin(ts) + + delete xp yp + erase red + for i 0 xs[] + point red LINE xs[$i] ys[$i] {10*cs[$i]} {10*sn[$i]} + end + plot xp yp -pt 100 +end + +macro check.flux + if ($0 != 2) + echo "USAGE: plot.angles (file.dat)" + break + end + + data $1 + read f 3 type 4 + subset F = f if (type == 1) + set M = -2.5*log(F) + vstat M +end + +if ($SCRIPT) + echo "no default action defined" + exit 0 +end Index: trunk/psphot/test/tap_psphot_galaxygrid.pro =================================================================== --- trunk/psphot/test/tap_psphot_galaxygrid.pro (revision 35769) +++ trunk/psphot/test/tap_psphot_galaxygrid.pro (revision 35769) @@ -0,0 +1,701 @@ +#!/usr/bin/env mana +# -*-sh-*- + +# 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 T + +# basic options for the these images (filter, location, obstype) +$BaseOptions = -type OBJECT +$BaseOptions = $BaseOptions -filter r +$BaseOptions = $BaseOptions -skymags 20.86 +$BaseOptions = $BaseOptions -ra 270.70 +$BaseOptions = $BaseOptions -dec -23.70 +$BaseOptions = $BaseOptions -pa 0.0 +# $BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0 + +# options for the reference image +$RefOptions = $BaseOptions +$RefOptions = $RefOptions -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 + +# options for the simulated images (using the refimage for the stars) +$FakeOptions = $BaseOptions +$FakeOptions = $FakeOptions -exptime 30.0 +$FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS + +# sample alternate options: +# $FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1 +# $FakeOptions = $FakeOptions -Df PSF.ARATIO 1.2 +# $FakeOptions = $FakeOptions -Df PSF.THETA +30.0 + +# create an image with fake sources (these are then inserted into the catdir) +$RefConfig = -camera SIMTEST +$RefConfig = $RefConfig -recipe PPSIM STACKTEST.MAKE + +# basic config for ppSim with randomly distributed stars and gridded galaxies +$FakeConfig = -camera SIMTEST +$FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN +$FakeConfig = $FakeConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref +$FakeConfig = $FakeConfig -Db STARS.FAKE F ; # only use stars from catdir.ref +$FakeConfig = $FakeConfig -Db MATCH.DENSITY F +$FakeConfig = $FakeConfig -Db PSF.CONVOLVE T +$FakeConfig = $FakeConfig -Db GALAXY.FAKE T ; # generate a "realistic" distribution of galaxies +$FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0 +$FakeConfig = $FakeConfig -Db GALAXY.GRID T ; # generate a grid of galaxies (constant mag) +#$FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_GAUSS +#$FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_EXP +#$FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_SERSIC +#$FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_DEV +$FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN 10.0 +$FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX 10.0 +$FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN 0.25 +$FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX 0.25 +$FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0 +$FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180 +$FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.66 +$FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.66 +$FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 120 +$FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 120 + +list fwhm + 1.0 + 1.1 + 1.2 + 1.5 +end + +macro go.grid.set.devexp + $FakeConfig = -camera SIMTEST + $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN + $FakeConfig = $FakeConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref + $FakeConfig = $FakeConfig -Db STARS.FAKE F ; # only use stars from catdir.ref + $FakeConfig = $FakeConfig -Db MATCH.DENSITY F + $FakeConfig = $FakeConfig -Db PSF.CONVOLVE T + $FakeConfig = $FakeConfig -Db GALAXY.FAKE T ; # generate a "realistic" distribution of galaxies + $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0 + $FakeConfig = $FakeConfig -Db GALAXY.GRID T ; # generate a grid of galaxies (constant mag) + $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0 + $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180 + $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 120 + $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 120 + $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0 + $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0 + + # $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_GAUSS + # $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_EXP + # $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_SERSIC + # $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_DEV + # $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN 10.0 + # $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX 10.0 + # $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN 0.25 + # $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX 0.25 + # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.66 + # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.66 + + $Nseq = 0 + foreach type EXP DEV + foreach Rmajor 3 10 30 + foreach Aratio 0.25 0.5 1.0 + foreach fwhm 0.8 1.0 1.5 + $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor + $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor + $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN $Aratio + $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX $Aratio + + sprint name "sample.%02d" $Nseq + mkexp $name $fwhm $type + fitexp $name $name.fit $type\_CONV + $Nseq ++ + end + end + end + end +end + +macro go.grid.set.sersic + $FakeConfig = -camera SIMTEST + $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN + $FakeConfig = $FakeConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref + $FakeConfig = $FakeConfig -Db STARS.FAKE F ; # only use stars from catdir.ref + $FakeConfig = $FakeConfig -Db MATCH.DENSITY F + $FakeConfig = $FakeConfig -Db PSF.CONVOLVE T + $FakeConfig = $FakeConfig -Db GALAXY.FAKE T ; # generate a "realistic" distribution of galaxies + $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0 + $FakeConfig = $FakeConfig -Db GALAXY.GRID T ; # generate a grid of galaxies (constant mag) + $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0 + $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180 + $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 120 + $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 120 + $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_SERSIC + + # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0 + # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0 + # $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN 10.0 + # $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX 10.0 + # $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN 0.25 + # $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX 0.25 + # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.66 + # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.66 + + $Nseq = 0 + foreach index 1 2 3 4 + foreach Rmajor 3 10 30 + foreach Aratio 0.25 0.5 1.0 + foreach fwhm 0.8 1.0 1.5 + $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor + $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor + $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN $Aratio + $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX $Aratio + $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN $index + $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX $index + + sprint name "sersic.%02d" $Nseq + mkexp $name $fwhm SERSIC + fitexp $name $name.fit SER\_CONV + $Nseq ++ + end + end + end + end +end + +macro go.grid.check.devexp + delete -q Xin_s Yin_s Min_s Tin_s Rin_s rin_s + delete -q Xot_s Yot_s Mot_s Tot_s Rot_s rot_s + + $Nseq = 0 + foreach type EXP DEV + foreach Rmajor 3 10 30 + foreach Aratio 0.25 0.5 1.0 + foreach fwhm 0.8 1.0 1.5 + sprint name "sample.%02d" $Nseq + ckgalaxy.load $name.dat $name.fit.cmf + $Nseq ++ + end + end + end + end +end + +macro go + mkexp test.exp 1.0 EXP + fitexp test.exp test.exp.fit EXP_CONV + + mkexp test.ser 1.0 SERSIC + fitexp test.ser test.ser.fit SER_CONV + + mkexp test.dev 1.0 DEV + fitexp test.dev test.dev.fit DEV_CONV + + mkexp test.gau 1.0 GAUSS + fitexp test.gau test.gau.fit GAU_CONV + + mkexp test.pg 1.0 PGAUSS + fitexp test.pg test.pg.fit PGA_CONV + + mkexp test.qga 1.0 QGAUSS + fitexp test.qga test.qga.fit QGA_CONV + + mkexp test.p1 1.0 PS1_V1 + fitexp test.p1 test.p1.fit PS1_CONV +end + +macro go.ckgalaxy + foreach type exp ser dev gau pg qga p1 + ckgalaxy test.$type.dat test.$type.fit.cmf + wait $type + end +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 + + echo ppSim $RefOptions $RefConfig refimage + exec ppSim $RefOptions $RefConfig refimage + + 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 -quick-airmass + 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 != 4) + echo "USAGE: mkexp basename (fwhm) (model)" + break + end + + local fwhm basename + $basename = $1 + $fwhm = $2 + + # create the raw image + $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_$3 + echo ppSim $FakeOptions $FakeConfig $basename -seeing $fwhm + exec ppSim $FakeOptions $FakeConfig $basename -seeing $fwhm + 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 + +# create a realistic distribution of fake stars, GAUSS PSF +macro fitexp + if ($0 != 4) + echo "USAGE: fitexp basename outname (fitModel)" + break + end + + local basename fitModel outname + $basename = $1 + $outname = $2 + $fitModel = $3 + + $psphotConfig = + $psphotConfig = $psphotConfig -recipe PSPHOT GALAXY_TEST + $psphotConfig = $psphotConfig -threads 4 + $psphotConfig = $psphotConfig -Db PSPHOT:LMM_FIT_CHISQ_CONVERGENCE F + $psphotConfig = $psphotConfig -Df PSPHOT:EXT_FIT_MIN_TOL 0.1 + $psphotConfig = $psphotConfig -Di PSPHOT:LMM_FIT_GAIN_FACTOR_MODE 2 + $psphotConfig = $psphotConfig -Db PSPHOT:SAVE.RESID T + $psphotConfig = $psphotConfig -D PSPHOT:EXTENDED_SOURCE_MODELS_SELECTION $fitModel + + # ppImage / psphot on the output + echo ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname + exec ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname +end + +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 ckgalaxy.load + if ($0 != 3) + echo "USAGE: ckgalaxy (dat) (cmf)" + break + end + + data $1 + read Xin_all 1 Yin_all 2 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9 + + subset Xin = Xin_all if (Type == 1) + subset Yin = Yin_all if (Type == 1) + subset Min = Min_all if (Type == 1) + + subset Tin_rad = ThetaIn_all if (Type == 1) + set Tin = Tin_rad * 180 / 3.14159265 + + subset RmajIn = RmajIn_all if (Type == 1) + subset RminIn = RminIn_all if (Type == 1) + + data $2 + read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA + set EXT_THETA_ALT = EXT_THETA * (EXT_THETA >= 0.0) + (EXT_THETA + 3.14159265) * (EXT_THETA < 0.0) + set EXT_THETA = EXT_THETA_ALT * 180 / 3.14159265 + + match2d X_EXT Y_EXT Xin Yin 1.0 -index1 index1 -index2 index2 + + reindex Xot_m = X_EXT using index1 + reindex Yot_m = Y_EXT using index1 + + reindex Xin_m = Xin using index2 + reindex Yin_m = Yin using index2 + + reindex Mot_m = EXT_INST_MAG using index1 + reindex Tot_m = EXT_THETA using index1 + + reindex Min_m = Min using index2 + reindex Tin_m = Tin using index2 + + reindex Rot_m = EXT_WIDTH_MAJ using index1 + reindex rot_m = EXT_WIDTH_MIN using index1 + + reindex Rin_m = RmajIn using index2 + reindex rin_m = RminIn using index2 + + foreach field X Y M T R r + foreach set in ot + concat $field\$set\_m $field\$set\_s + end + 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 + +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 plot.angles + if ($0 != 2) + echo "USAGE: plot.angles (file.dat)" + break + end + + data $1 + read x 1 y 2 type 4 trad 9 + set t = trad * 180 / 3.14159265 + subset xs = x if (type == 1) + subset ys = y if (type == 1) + subset ts = t if (type == 1) + + set cs = dcos(ts) + set sn = dsin(ts) + + delete xp yp + erase red + for i 0 xs[] + point red LINE xs[$i] ys[$i] {10*cs[$i]} {10*sn[$i]} + end + plot xp yp -pt 100 +end + +macro check.flux + if ($0 != 2) + echo "USAGE: plot.angles (file.dat)" + break + end + + data $1 + read f 3 type 4 + subset F = f if (type == 1) + set M = -2.5*log(F) + vstat M +end + +if ($SCRIPT) + echo "no default action defined" + exit 0 +end