Index: /trunk/psphot/src/psphotBlendFit.c =================================================================== --- /trunk/psphot/src/psphotBlendFit.c (revision 36085) +++ /trunk/psphot/src/psphotBlendFit.c (revision 36086) @@ -273,6 +273,6 @@ pmSource *source = sources->data[i]; +# if (0) int TEST_ON = false; -# if (0) # define TEST_X 653 # define TEST_Y 466 @@ -342,7 +342,9 @@ if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { if (psphotFitBlob (readout, source, newSources, psf, fitOptions, maskVal, markVal)) { +# if (PS_TRACE_ON) if (TEST_ON) { psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); } +# endif psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf); Next ++; @@ -351,7 +353,9 @@ } else { if (psphotFitBlend (readout, source, psf, fitOptions, maskVal, markVal)) { +# if (PS_TRACE_ON) if (TEST_ON) { psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); } +# endif source->type = PM_SOURCE_TYPE_STAR; psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf); Index: /trunk/psphot/src/psphotDeblendSatstars.c =================================================================== --- /trunk/psphot/src/psphotDeblendSatstars.c (revision 36085) +++ /trunk/psphot/src/psphotDeblendSatstars.c (revision 36086) @@ -128,4 +128,6 @@ // int display = psphotKapaChannel (1); // psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 0); + + int Nsatstar = 0; // examine sources in decreasing SN order @@ -179,4 +181,6 @@ source->mode2 |= PM_SOURCE_MODE2_SATSTAR_PROFILE; // and we have in fact subtracted the profile + Nsatstar ++; + // XXX visualize, model, and subtract // if (!psphotVisualRadialProfileSatstar (source, maskVal)) { @@ -192,5 +196,5 @@ psFree (index); - psLogMsg ("psphot", PS_LOG_INFO, "deblend satstar: %f sec\n", psTimerMark ("psphot.deblend.sat")); + psLogMsg ("psphot", PS_LOG_INFO, "deblend %d satstars: %f sec\n", Nsatstar, psTimerMark ("psphot.deblend.sat")); return true; } Index: /trunk/psphot/src/psphotExtendedSourceFits.c =================================================================== --- /trunk/psphot/src/psphotExtendedSourceFits.c (revision 36085) +++ /trunk/psphot/src/psphotExtendedSourceFits.c (revision 36086) @@ -65,4 +65,6 @@ psphotVisualShowImage(readout); + // psphotSaveImage (NULL, readout->image, "test.01.fits"); + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); psAssert (detections, "missing detections?"); @@ -123,7 +125,11 @@ } + // perform full extended source non-linear fits? + bool isInteractive = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_INTERACTIVE"); + if (!status) isInteractive = false; + // Define source fitting parameters for extended source fits pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); - fitOptions->mode = PM_SOURCE_FIT_EXT; + fitOptions->mode = PM_SOURCE_FIT_EXT_AND_SKY; fitOptions->saveCovariance = true; // XXX make this a user option? fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix @@ -134,4 +140,5 @@ fitOptions->gainFactorMode = gainFactorMode; fitOptions->chisqConvergence = chisqConvergence; + fitOptions->isInteractive = isInteractive; // maskVal is used to test for rejected pixels, and must include markVal @@ -296,41 +303,42 @@ PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail -// set this to 0 to run without threading -# if (1) - if (!psThreadJobAddPending(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree(AnalysisRegion); - psFree (fitOptions); - psFree (models); - psphotSersicModelClassCleanup(); - return false; - } -# else - if (!psphotExtendedSourceFits_Threaded(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree(AnalysisRegion); - psFree (fitOptions); - psFree (models); - psphotSersicModelClassCleanup(); - return false; +// XXX TEST + if (!isInteractive) { + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + psFree (fitOptions); + psFree (models); + psphotSersicModelClassCleanup(); + return false; + } + } else { + // run without threading + if (!psphotExtendedSourceFits_Threaded(job)) { + 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[9]; + Next += scalar->data.S32; + scalar = job->args->data[10]; + Nconvolve += scalar->data.S32; + scalar = job->args->data[11]; + NconvolvePass += scalar->data.S32; + scalar = job->args->data[12]; + Nplain += scalar->data.S32; + scalar = job->args->data[13]; + NplainPass += scalar->data.S32; + scalar = job->args->data[14]; + Nfaint += scalar->data.S32; + scalar = job->args->data[15]; + Nfail += scalar->data.S32; + psFree(job->args->data[3]); // iterator allocated above + psFree(job); } - psScalar *scalar = NULL; - scalar = job->args->data[9]; - Next += scalar->data.S32; - scalar = job->args->data[10]; - Nconvolve += scalar->data.S32; - scalar = job->args->data[11]; - NconvolvePass += scalar->data.S32; - scalar = job->args->data[12]; - Nplain += scalar->data.S32; - scalar = job->args->data[13]; - NplainPass += scalar->data.S32; - scalar = job->args->data[14]; - Nfaint += scalar->data.S32; - scalar = job->args->data[15]; - Nfail += scalar->data.S32; - psFree(job->args->data[3]); // iterator allocated above - psFree(job); -# endif } @@ -378,4 +386,6 @@ psphotSersicModelClassCleanup(); + // psphotSaveImage (NULL, readout->image, "test.02.fits"); + psphotVisualShowResidualImage (readout, false); @@ -521,4 +531,8 @@ psFree (source->modelFlux); source->modelFlux = NULL; + // if ((fabs(source->peak->xf - 2572) < 20) && (fabs(source->peak->yf - 5874) < 20)) { + // fprintf (stderr, "test object\n"); + // } + modelFit = psphotFitEXT (modelFit, readout, source, fitOptions, modelType, maskVal, markVal); if (!modelFit) { @@ -544,7 +558,5 @@ 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)) { + if ((fabs(source->peak->xf - 2572) < 20) && (fabs(source->peak->yf - 5874) < 20)) { fprintf (stderr, "src vs fit : %d %d - %d %d | %f %f - %f %f\n", source->pixels->col0, source->pixels->row0, @@ -579,5 +591,5 @@ if (model->flags & (PM_MODEL_STATUS_BADARGS)) continue; - if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue; + // if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue; if (model->flags & (PM_MODEL_STATUS_OFFIMAGE)) continue; @@ -598,4 +610,11 @@ pmSourceCacheModel (source, maskVal); + pmModel *model = source->modelFits->data[0]; + int flags = 0xffffffff; + if (model) { + flags = model->flags; + } + + fprintf (stderr, "failed to fit extended source model to object %d @ %f, %f (%x)\n", source->id, source->moments->Mx, source->moments->My, flags); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); @@ -615,6 +634,8 @@ // cache the model flux if (source->modelEXT->isPCM) { + // fprintf (stderr, "subtract PCM extended source model for object %d @ %f, %f\n", source->id, source->moments->Mx, source->moments->My); pmPCMCacheModel (source, maskVal, psfSize); } else { + // fprintf (stderr, "subtract non-PCM extended source model for object %d @ %f, %f\n", source->id, source->moments->Mx, source->moments->My); pmSourceCacheModel (source, maskVal); } Index: /trunk/psphot/src/psphotRadiusChecks.c =================================================================== --- /trunk/psphot/src/psphotRadiusChecks.c (revision 36085) +++ /trunk/psphot/src/psphotRadiusChecks.c (revision 36086) @@ -171,6 +171,9 @@ # undef MIN_WINDOW +// XXX EAM : 20130724 : for a test, double the window size parameters +// # define MIN_WINDOW 10.0 +// # define SCALE1 7.0 # define MIN_WINDOW 5.0 -# define SCALE1 5.0 +# define SCALE1 3.0 # define PAD_WINDOW 3.0 Index: /trunk/psphot/src/psphotSourceFits.c =================================================================== --- /trunk/psphot/src/psphotSourceFits.c (revision 36085) +++ /trunk/psphot/src/psphotSourceFits.c (revision 36086) @@ -559,4 +559,10 @@ # define TIMING 0 +# define EXTRA_VERBOSE 0 + +bool psphotSersicModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize); +bool psphotFitSersicShapeAndIndex (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize); +bool psphotFitSersicShapeAndIndexGrid (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize); +bool psphotFitSersicShapeAndIndexGridAuto (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize); pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { @@ -580,5 +586,16 @@ float t1, t2, t4, t5; + t1 = t2 = t4 = t5 = 0.0; if (TIMING) { psTimerStart ("psphotFitPCM"); } + + if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { + options.mode = PM_SOURCE_FIT_NO_INDEX; // XXX note that there may be a conflict with psphotExtendedSourceFits.c:133 + options.mode = PM_SOURCE_FIT_EXT_AND_SKY; + // if we are ever (in a given psphot implementation) going to fit a parameter, we must set the options here to include + // that parameter (otherwise pmPCMupdate will fail to allocate the dmodelFlux image) + // thus, if the sersic analysis below uses an index fit, need to use this EXT_AND_SKY mode for init + } else { + options.mode = PM_SOURCE_FIT_EXT_AND_SKY; + } pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize); @@ -594,7 +611,18 @@ if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { // use the source moments, etc to guess basic model parameters - if (!psphotSersicModelClassGuessPCM (pcm, source)) { + if (!psphotSersicModelGuessPCM (pcm, source, maskVal, psfSize)) { psFree (pcm); model->flags |= PM_MODEL_STATUS_BADARGS; + return model; + } + if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); } + + // psphotFitSersicShapeAndIndex (pcm, readout, source, fitOptions, maskVal, markVal, psfSize); + options.mode = PM_SOURCE_FIT_NO_INDEX; + if (!psphotFitSersicShapeAndIndexGridAuto (pcm, readout, source, &options, maskVal, markVal, psfSize)) { + psFree (pcm); + model->flags |= PM_MODEL_STATUS_BADARGS; + psError(PS_ERR_UNKNOWN, true, "Failed to find a index & shape"); + psErrorClear (); // clear the polynomial error return model; } @@ -608,13 +636,4 @@ } - if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); } - - if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { - options.mode = PM_SOURCE_FIT_NO_INDEX; - } else { - options.mode = PM_SOURCE_FIT_EXT; - } - // update the pcm elements if we have changed the circumstance (options.mode or source->pixels) - pmPCMupdate(pcm, source, &options, model); if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); } @@ -629,6 +648,10 @@ fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d of %5d, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, nPixBig, t1, t2, t4, t5); } - - // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); + if (EXTRA_VERBOSE && !TIMING) { + int nPixBig = source->pixels->numCols * source->pixels->numRows; + float *PAR = model->params->data.F32; + fprintf (stderr, "source %d : %f - %f %f - %f %f %f - %f | nIter: %2d, radius: %6.1f, npix: %5d of %5d, chisq %f\n", source->id, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1], model->nIter, model->fitRadius, model->nPix, nPixBig, model->chisqNorm); + } + psFree (pcm); @@ -748,2 +771,472 @@ return true; } + +// float indexGuessInv[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083}; +float indexGuessInv[] = {0.5, 0.4, 0.3, 0.25, 0.20, 0.15, 0.125}; +# define N_INDEX_GUESS_INV 7 + +// float reffGuess[] = {3.0, 10.0, 20.0, 30.0, 40.0}; +float reffGuess[] = {0.5, 0.75, 1.0, 1.4, 2.0}; +# define N_REFF_GUESS 5 + +// A sersic model is very sensitive to the index. attempt to find the index first by grid search in just the index +// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations +bool psphotSersicModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize) { + + // we get a reasonable guess from: + // * Reff = Kron R1 + // * Rmajor / Rminor & Theta from moments + // * Io from total Kron flux + + // the guesses are used to fill in PAR: + psF32 *PAR = pcm->modelConv->params->data.F32; + + // convert the moments to Major,Minor,Theta + psEllipseMoments moments; + + if (!isfinite(source->moments->Mrf)) return false; + if (!isfinite(source->moments->Mxx)) return false; + if (!isfinite(source->moments->Mxy)) return false; + if (!isfinite(source->moments->Myy)) return false; + + moments.x2 = source->moments->Mxx; + moments.y2 = source->moments->Myy; + moments.xy = source->moments->Mxy; + + // limit axis ratio < 20.0 + psEllipseAxes momentAxes = psEllipseMomentsToAxes (moments, 20.0); + + if (0) { + psEllipseAxes guessAxes; + guessAxes.major = source->moments->Mrf; + guessAxes.minor = (momentAxes.minor / momentAxes.minor) * guessAxes.major; + guessAxes.theta = momentAxes.theta; + + if (!isfinite(guessAxes.major)) return false; + if (!isfinite(guessAxes.minor)) return false; + if (!isfinite(guessAxes.theta)) return false; + + // convert the major,minor,theta to shape parameters for an Reff-like model + pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true); + } + + // set the model position + if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) { + return false; + } + + // sky is zero (no longer fitted, but not yet deprecated) + PAR[PM_PAR_SKY] = 0.0; + + // for the index loop, use Io = 1.0, use fitted values to determine Io + PAR[PM_PAR_I0] = 1.0; + + float xMin = NAN; + float iMin = NAN; + float sMin = NAN; + float rMin = NAN; + + // loop over index and Reff, keeping the ARatio and Theta constant? + // loop over index guesses and find the best fit + for (int j = 0; j < N_REFF_GUESS; j++) { + + psEllipseAxes guessAxes; + guessAxes.major = reffGuess[j] * source->moments->Mrf; + guessAxes.minor = guessAxes.major * (momentAxes.minor / momentAxes.major); + guessAxes.theta = momentAxes.theta; + + if (!isfinite(guessAxes.major)) return false; + if (!isfinite(guessAxes.minor)) return false; + if (!isfinite(guessAxes.theta)) return false; + + // convert the major,minor,theta to shape parameters for an Reff-like model + pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true); + + for (int i = 0; i < N_INDEX_GUESS_INV; i++) { + PAR[PM_PAR_7] = indexGuessInv[i]; + + // generated the modelFlux + pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize); + + float YY = 0.0; + float YM = 0.0; + float MM = 0.0; + bool usePoisson = false; + + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + // skip masked points + if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) { + continue; + } + // skip zero-variance points + if (source->variance->data.F32[iy][ix] == 0) { + continue; + } + // skip nan value points + if (!isfinite(source->pixels->data.F32[iy][ix])) { + continue; + } + + float fy = source->pixels->data.F32[iy][ix]; + float fm = source->modelFlux->data.F32[iy][ix]; + float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0; + + YY += PS_SQR(fy) * wt; + YM += fm * fy * wt; + MM += PS_SQR(fm) * wt; + } + } + + float Io = YM / MM; + float Chisq = YY - 2 * Io * YM + Io * Io * MM; + if (isnan(xMin) || (Chisq < xMin)) { + xMin = Chisq; + iMin = Io; + sMin = indexGuessInv[i]; + rMin = reffGuess[j]; + } + // fprintf (stderr, "%d | %f %f %f %f | %f %f %f %f", i, indexGuessInv[i], reffGuess[j], Io, Chisq, sMin, rMin, iMin, xMin); + // fprintf (stderr, "\n"); + } + } + + { + psEllipseAxes guessAxes; + guessAxes.major = rMin * source->moments->Mrf; + guessAxes.minor = guessAxes.major * (momentAxes.minor / momentAxes.major); + guessAxes.theta = momentAxes.theta; + + if (!isfinite(guessAxes.major)) return false; + if (!isfinite(guessAxes.minor)) return false; + if (!isfinite(guessAxes.theta)) return false; + + // convert the major,minor,theta to shape parameters for an Reff-like model + pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true); + } + + PAR[PM_PAR_I0] = iMin; + PAR[PM_PAR_7] = sMin; + + return true; +} + +// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX +bool psphotFitSersicShapeAndIndex (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + + pmModel *model = pcm->modelConv; + + assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); + + pmSourceFitOptions options = *fitOptions; + + for (int i = 0; i < 3; i++) { + // fit EXT (not PSF) model (set/unset the pixel mask) + options.mode = PM_SOURCE_FIT_SHAPE; + options.nIter = 2; + + // update the pcm elements if we have changed the circumstance (here, options.mode) + pmPCMupdate(pcm, source, &options, model); + + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + if (EXTRA_VERBOSE) { + float *PAR = model->params->data.F32; + fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]); + } + + // fit EXT (not PSF) model (set/unset the pixel mask) + options.mode = PM_SOURCE_FIT_INDEX; + // options.mode = PM_SOURCE_FIT_EXT_AND_SKY; + options.nIter = 30; + + // update the pcm elements if we have changed the circumstance (here, options.mode) + pmPCMupdate(pcm, source, &options, model); + + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + if (EXTRA_VERBOSE) { + float *PAR = model->params->data.F32; + fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]); + } + } + + // update the pcm elements if we have changed the circumstance (here, options.mode) + pmPCMupdate(pcm, source, fitOptions, model); + + return true; +} + +// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX +bool psphotFitSersicShapeAndIndexGridAuto (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + + pmModel *model = pcm->modelConv; + + assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); + + pmSourceFitOptions options = *fitOptions; + + psF32 *PAR = pcm->modelConv->params->data.F32; + + options.mode = PM_SOURCE_FIT_SHAPE; + options.nIter = 7; + + // update the pcm elements if we have changed the circumstance (here, options.mode) + pmPCMupdate(pcm, source, &options, model); + + // we have been provided a guess at the index (P[7]) from the list of indexGuessInv + + // find the matching indexGuessInv + int nStart = -1; + for (int i = 0; i < N_INDEX_GUESS_INV; i++) { + if (fabs(PAR[PM_PAR_7] - indexGuessInv[i]) < 0.01) { + nStart = i; + break; + } + } + if (nStart == -1) { + fprintf (stderr, "WARNING: could not find start guess %f\n", PAR[PM_PAR_7]); + return false; + } + + psVector *chi2 = psVectorAllocEmpty (16, PS_TYPE_F32); + psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32); + + PAR[PM_PAR_7] = indexGuessInv[nStart]; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]); + psVectorAppend (Sidx, 100*PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + + PAR[PM_PAR_7] = (nStart < N_INDEX_GUESS_INV - 1) ? 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]) : 0.1; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]); + psVectorAppend (Sidx, 100*PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + + PAR[PM_PAR_7] = (nStart > 0) ? 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]) : 0.55; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]); + psVectorAppend (Sidx, 100*PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + + if (chi2->data.F32[1] < chi2->data.F32[2]) { + if (nStart == N_INDEX_GUESS_INV - 1) { + PAR[PM_PAR_7] = 0.11; + } else { + PAR[PM_PAR_7] = indexGuessInv[nStart + 1]; + } + } else { + if (nStart == 0) { + PAR[PM_PAR_7] = 0.52; + } else { + PAR[PM_PAR_7] = indexGuessInv[nStart - 1]; + } + } + + psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2); + if (!psVectorFitPolynomial1D (poly, NULL, 0, chi2, NULL, Sidx)) { + psError(PS_ERR_UNKNOWN, true, "Failed to find a good chisq parabola"); + psFree (chi2); + psFree (Sidx); + psFree (poly); + return false; + } + + // where is the minimum of this polynomial fit? + float Smin = -0.5 * poly->coeff[1] / poly->coeff[2] / 100.0; + + // constrain Smin to be in a valid range (1.0 - 0.1, corresponding to 0.5 (Gauss) to 5.0 (slightly peakier than Dev) + Smin = PS_MAX(PS_MIN(Smin, 1.0), 0.1); + PAR[PM_PAR_7] = Smin; + + // return to the original fitting mode (fitOptions) + pmPCMupdate(pcm, source, fitOptions, model); + + psFree (chi2); + psFree (Sidx); + psFree (poly); + + return true; +} + + +// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX +bool psphotFitSersicShapeAndIndexGridAutoScaled (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + + pmModel *model = pcm->modelConv; + + assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); + + pmSourceFitOptions options = *fitOptions; + + psF32 *PAR = pcm->modelConv->params->data.F32; + + options.mode = PM_SOURCE_FIT_SHAPE; + options.nIter = 5; + + // update the pcm elements if we have changed the circumstance (here, options.mode) + pmPCMupdate(pcm, source, &options, model); + + float parStart[8]; + for (int i = 0; i < 8; i++) parStart[i] = PAR[i]; + + // we start with a guess at the index (P[7]) + + // get chisq for P[7], P[7]*1.1, P[7]*1.25 (or *0.75 depending on the result of 1.1) + + psVector *chi2 = psVectorAllocEmpty (16, PS_TYPE_F32); + psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32); + + PAR[PM_PAR_7] = parStart[7]; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]); + psVectorAppend (Sidx, 100*PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + + float fI = 1.1; + PAR[PM_PAR_7] = parStart[7]*fI; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]); + psVectorAppend (Sidx, 100*PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + + if (chi2->data.F32[1] < chi2->data.F32[0]) { + fI = 1.3; + } else { + fI = 1.0 / 1.3; + } + + PAR[PM_PAR_7] = parStart[7]*fI; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]); + psVectorAppend (Sidx, 100*PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + + // can we fit the 3 pts with a parabola? + int nTry = 0; + psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2); + while (!psVectorFitPolynomial1D (poly, NULL, 0, chi2, NULL, Sidx)) { + psErrorClear (); // clear the polynomial error + if (nTry > 4) { + psError(PS_ERR_UNKNOWN, true, "Failed to find a good chisq parabola"); + psFree (chi2); + psFree (Sidx); + psFree (poly); + return false; + } + fI = (fI < 1.0) ? fI / 1.3 : fI * 1.3; + PAR[PM_PAR_7] = parStart[7]*fI; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]); + psVectorAppend (Sidx, 100*PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + nTry ++; + } + + // where is the minimum of this polynomial fit? + float Smin = -0.5 * poly->coeff[1] / poly->coeff[2] / 100.0; + + // constrain Smin to be in a valid range (1.0 - 0.1, corresponding to 0.5 (Gauss) to 5.0 (slightly peakier than Dev) + Smin = PS_MAX(PS_MIN(Smin, 1.0), 0.1); + PAR[PM_PAR_7] = Smin; + + // pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + // if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]); + + //// for (int i = 0; i < 8; i++) PAR[i] = parStart[i]; + //// + //// for (float fI = 0.0; fI < 0.15; fI += 0.01) { + //// PAR[PM_PAR_7] = parStart[7] - fI; + //// + //// // fit EXT (not PSF) model (set/unset the pixel mask) + //// + //// pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + //// if (TIMING) { + //// float *PAR = model->params->data.F32; + //// fprintf (stderr, "%d %f : %f - %f %f %f - %f\n", model->nIter, model->chisqNorm, PAR[7], PAR[4], PAR[5], PAR[6], PAR[1]); + //// } + //// } + + // return to the original fitting mode (fitOptions) + pmPCMupdate(pcm, source, fitOptions, model); + + psFree (chi2); + psFree (Sidx); + psFree (poly); + + return true; +} + + +// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX +bool psphotFitSersicShapeAndIndexGrid (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + + pmModel *model = pcm->modelConv; + + assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); + + pmSourceFitOptions options = *fitOptions; + + psF32 *PAR = pcm->modelConv->params->data.F32; + + options.mode = PM_SOURCE_FIT_SHAPE; + options.nIter = 10; + + // update the pcm elements if we have changed the circumstance (here, options.mode) + pmPCMupdate(pcm, source, &options, model); + + psVector *chi2 = psVectorAllocEmpty (16, PS_TYPE_F32); + psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32); + + float par7[] = {0.100, 0.125, 0.150, 0.175, 0.200, 0.225, 0.250}; + for (int i = 0; i < 7; i++) { + PAR[PM_PAR_7] = par7[i]; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]); + psVectorAppend (Sidx, PAR[PM_PAR_7]); + psVectorAppend (chi2, model->chisqNorm); + } + + psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2); + if (!psVectorFitPolynomial1D (poly, NULL, 0, chi2, NULL, Sidx)) { + psError(PS_ERR_UNKNOWN, true, "Failed to find a good chisq parabola"); + psFree (chi2); + psFree (Sidx); + psFree (poly); + return false; + } + + // where is the minimum of this polynomial fit? + fprintf (stderr, "fit1d: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); + float Smin = -0.5 * poly->coeff[1] / poly->coeff[2]; + + // constrain Smin to be in a valid range (1.0 - 0.1, corresponding to 0.5 (Gauss) to 5.0 (slightly peakier than Dev) + Smin = PS_MAX(PS_MIN(Smin, 1.0), 0.1); + PAR[PM_PAR_7] = Smin; + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]); + + //// for (int i = 0; i < 8; i++) PAR[i] = parStart[i]; + //// + //// for (float fI = 0.0; fI < 0.15; fI += 0.01) { + //// PAR[PM_PAR_7] = parStart[7] - fI; + //// + //// // fit EXT (not PSF) model (set/unset the pixel mask) + //// + //// pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + //// if (TIMING) { + //// float *PAR = model->params->data.F32; + //// fprintf (stderr, "%d %f : %f - %f %f %f - %f\n", model->nIter, model->chisqNorm, PAR[7], PAR[4], PAR[5], PAR[6], PAR[1]); + //// } + //// } + + // return to the original fitting mode (fitOptions) + pmPCMupdate(pcm, source, fitOptions, model); + + psFree (chi2); + psFree (Sidx); + psFree (poly); + + return true; +} + + Index: /trunk/psphot/test/tap_psphot_galaxygrid.pro =================================================================== --- /trunk/psphot/test/tap_psphot_galaxygrid.pro (revision 36085) +++ /trunk/psphot/test/tap_psphot_galaxygrid.pro (revision 36086) @@ -28,4 +28,5 @@ $RefOptions = $RefOptions -Df STARS.DENSITY 10.0 $RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5 +$RefOptions = $RefOptions -nx 3000 -ny 3000 # options for the simulated images (using the refimage for the stars) @@ -33,4 +34,5 @@ $FakeOptions = $FakeOptions -exptime 30.0 $FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS +$FakeOptions = $FakeOptions -nx 3000 -ny 3000 # sample alternate options: @@ -65,6 +67,6 @@ $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 +$FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300 +$FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300 list fwhm @@ -75,5 +77,6 @@ end -macro go.grid.set.devexp +# generate fake images and run psphot on them +macro normtest.mkexp.devexp $FakeConfig = -camera SIMTEST $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN @@ -84,11 +87,13 @@ $FakeConfig = $FakeConfig -Db GALAXY.FAKE T ; # generate a "realistic" distribution of galaxies $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0 + $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5 $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 -Di GALAXY.GRID.DX 300 + $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300 $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0 $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0 + $BaseConfig = $FakeConfig # $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_GAUSS @@ -103,9 +108,12 @@ # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.66 + mkdir normtest + $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 + foreach fwhm 1.0 + $FakeConfig = $BaseConfig $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor @@ -113,7 +121,6 @@ $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX $Aratio - sprint name "sample.%02d" $Nseq + sprint name "normtest/test.%02d" $Nseq mkexp $name $fwhm $type - fitexp $name $name.fit $type\_CONV $Nseq ++ end @@ -123,5 +130,13 @@ end -macro go.grid.set.sersic +# generate fake images and run psphot on them +macro grid.mkexp.devexp + if ($0 != 2) + echo "USAGE: grid.mkexp.devexp (dir)" + break + end + + mkdir $1 + $FakeConfig = -camera SIMTEST $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN @@ -132,19 +147,155 @@ $FakeConfig = $FakeConfig -Db GALAXY.FAKE T ; # generate a "realistic" distribution of galaxies $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0 + $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5 $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 -Di GALAXY.GRID.DX 300 + $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300 + $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0 + $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0 + $BaseConfig = $FakeConfig + + $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 = $BaseConfig + $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 "$1/sample.%02d" $Nseq + mkexp $name $fwhm $type + $Nseq ++ + end + end + end + end +end + +# generate fake images and run psphot on them +macro grid.fitexp.devexp + if ($0 != 3) + echo "USAGE: grid.fitexp.devexp (srcdir) (outdir)" + break + end + + #$Nseq = 0 + #foreach type EXP DEV + + mkdir $2 + + $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 srcname "$1/sample.%02d" $Nseq + sprint outname "$2/sample.%02d.fit" $Nseq + fitexp $srcname $outname $type\_CONV + $Nseq ++ + end + end + end + end +end + +macro grid.load.devexp + if ($0 != 3) + echo "USAGE: grid.load.devexp (srcdir) (fitdir)" + break + end + + 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 + delete -q min_S Min_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 + cmf.load.concat $1/$name.dat $2/$name.fit.cmf $type + $Nseq ++ + end + end + end + end +end + +macro grid.plots.devexp + if ($0 != 2) + echo "USAGE: grid.plot.devexp (version)" + break + end + # things to examine: theta, Rmajor, AR + + # go.grid.check.devexp + + # check on theta + set dT = Tot_s - Tin_s + set ARin = rin_s / Rin_s + # lim ARin dT; clear; box; plot ARin dT + + subset dTx = dT if (ARin < 0.95) + histogram dTx NdT -8 8 0.1 -range dTi + lim -n 0$1 dTi NdT; clear; box; plot -x 1 dTi NdT + label -x "angle offset (degrees, only non-circular)" -y "number count" + resize 700 320 + + # check on Rmajor + set dR = Rin_s - Rot_s + # lim Rin_s dR; clear; box; plot Rin_s dR + # lim Rot_s dR; clear; box; plot Rot_s dR + # lim Rin_s dR; clear; box; plot Rin_s dR + create n 0 dR[] + set dRf = dR / Rin_s + lim -n 1$1 n -0.5 0.5; clear; box; plot n dRf + label -x sequence -y "1 - R_out| / R_in|" + resize 700 320 + + # check on A.Ratio + set ARot = rot_s / Rot_s + # lim ARin ARot; clear; box; plot ARin ARot + set fAR = ARot / ARin + lim -n 2$1 n 0.5 1.5; clear; box; plot n fAR + label -x sequence -y "AR_out| / AR_in|" + resize 700 320 + + # check on magnitude + set dM = Mot_s - Min_s + lim -n 3$1 n -0.5 0.5; clear; box; plot n dM + label -x sequence -y "M_out| - M_in|" + resize 700 320 +end + +macro grid.mkexp.sersic + if ($0 != 2) + echo "USAGE: grid.mkexp.devexp (dir)" + break + end + + mkdir $1 + + $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 -Df GALAXY.GRID.MAG 14.5 + $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 300 + $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300 $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 + $BaseConfig = $FakeConfig $Nseq = 0 @@ -153,4 +304,5 @@ foreach Aratio 0.25 0.5 1.0 foreach fwhm 0.8 1.0 1.5 + $FakeConfig = $BaseConfig $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor @@ -160,7 +312,7 @@ $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX $index - sprint name "sersic.%02d" $Nseq + sprint name "$1/sersic.%02d" $Nseq mkexp $name $fwhm SERSIC - fitexp $name $name.fit SER\_CONV + echo "$Nseq : $index $Rmajor $Aratio $fwhm" $Nseq ++ end @@ -170,15 +322,19 @@ 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 +macro grid.fitexp.sersic.devexp + if ($0 != 3) + echo "USAGE: grid.fitexp.sersic.devexp (srcdir) (outdir)" + break + end + + mkdir $2 $Nseq = 0 - foreach type EXP DEV + 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 - sprint name "sample.%02d" $Nseq - ckgalaxy.load $name.dat $name.fit.cmf + sprint name "sersic.%02d" $Nseq + fitexp $1/$name $2/$name.fit EXP_CONV,DEV_CONV $Nseq ++ end @@ -188,32 +344,105 @@ 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 +macro grid.fitexp.sersic + if ($0 != 3) + echo "USAGE: grid.fitexp.sersic (srcdir) (outdir)" + break + end + + mkdir $2 + + $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 + sprint name "sersic.%02d" $Nseq + fitexp $1/$name $2/$name.fit SER\_CONV + $Nseq ++ + end + end + end + end +end + +macro grid.load.sersic + if ($0 != 3) + echo "USAGE: grid.load.devexp (srcdir) (fitdir)" + break + end + + delete -q Xin_s Yin_s Min_s Tin_s Rin_s rin_s MTin_s Iin_s + delete -q Xot_s Yot_s Mot_s Tot_s Rot_s rot_s MTot_s Iot_s + + $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 + sprint name "sersic.%02d" $Nseq + cmf.load.concat $1/$name.dat $2/$name.fit.cmf SERSIC + $Nseq ++ + end + end + end + end +end + +# I want to make plots of Iin_s vs Mkron, Mxx, and similar things +# this means I need to be able to join Chip.xfit things against Chip.psf +# and to join Chip.xfit(DEV) to Chip.xfit(EXP) + +# I think I need a generic 'JOIN' function + +macro grid.plots.sersic + if ($0 != 2) + echo "USAGE: grid.plot.devexp (version)" + break + end + # things to examine: theta, Rmajor, AR + + # go.grid.check.devexp + + # check on theta + set dT = Tot_s - Tin_s + set ARin = rin_s / Rin_s + # lim ARin dT; clear; box; plot ARin dT + + subset dTx = dT if (ARin < 0.95) + histogram dTx NdT -8 8 0.1 -range dTi + lim -n 0$1 dTi NdT; clear; box; plot -x 1 dTi NdT + label -x "angle offset (degrees, only non-circular)" -y "number count" + resize 700 320 + + # check on Rmajor + set dR = Rin_s - Rot_s + # lim Rin_s dR; clear; box; plot Rin_s dR + # lim Rot_s dR; clear; box; plot Rot_s dR + # lim Rin_s dR; clear; box; plot Rin_s dR + create n 0 dR[] + set dRf = dR / Rin_s + lim -n 1$1 n -0.5 0.5; clear; box; plot n dRf + label -x sequence -y "1 - R_out| / R_in|" + resize 700 320 + + # check on A.Ratio + set ARot = rot_s / Rot_s + # lim ARin ARot; clear; box; plot ARin ARot + set fAR = ARot / ARin + lim -n 2$1 n 0.5 1.5; clear; box; plot n fAR + label -x sequence -y "AR_out| / AR_in|" + resize 700 320 + + # check on magnitude + set dM = Mot_s - Min_s + lim -n 3$1 n -0.5 0.5; clear; box; plot n dM + label -x sequence -y "M_out| - M_in|" + resize 700 320 + + # check on index + set dI = Iot_s - Iin_s + lim -n 4$1 n -0.5 0.5; clear; box; plot n dI + label -x sequence -y "I_out| - I_in|" + resize 700 320 end @@ -280,68 +509,25 @@ # ppImage / psphot on the output + break -auto off 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)" + exec echo ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname >& $outname.log + exec ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname >>& $outname.log + break -auto on +end + +macro cmf.load.concat + if ($0 != 4) + echo "USAGE: cmf.load.concat (dat) (cmf) (type)" 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 + read Xin_all 1 Yin_all 2 Fin_all 3 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9 IndexIn_all 10 subset Xin = Xin_all if (Type == 1) subset Yin = Yin_all if (Type == 1) subset Min = Min_all if (Type == 1) + subset Fin = Fin_all if (Type == 1) + set min = -2.5*log(Fin) subset Tin_rad = ThetaIn_all if (Type == 1) @@ -351,6 +537,30 @@ subset RminIn = RminIn_all if (Type == 1) + subset IndexIn = IndexIn_all if (Type == 1) + + $TYPE_S = 83 + $TYPE_D = 68 + $TYPE_E = 69 + if ("$3" == "SERSIC") + $InType = $TYPE_S + end + if ("$3" == "DEV") + $InType = $TYPE_D + end + if ("$3" == "EXP") + $InType = $TYPE_E + end + data $2 - read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA + + break -auto off + read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE EXT_PAR_07 + $reread = not($STATUS) + break -auto on + if ($reread) + read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE + set EXT_PAR_07 = (MODEL_TYPE:9 == 68)*4 + (MODEL_TYPE:9 == 69) + end + 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 @@ -364,4 +574,7 @@ reindex Yin_m = Yin using index2 + set MTin_m = $InType + zero(Xin_m) + reindex MTot_m = MODEL_TYPE:9 using index1 + reindex Mot_m = EXT_INST_MAG using index1 reindex Tot_m = EXT_THETA using index1 @@ -376,261 +589,183 @@ reindex rin_m = RminIn using index2 - foreach field X Y M T R r + if ("$3" == "SERSIC") + reindex Iot_m = EXT_PAR_07 using index1 + reindex Iin_m = IndexIn using index2 + $fields = X Y M T R r MT I + else + $fields = X Y M T R r MT + end + + foreach field $fields foreach set in ot concat $field\$set\_m $field\$set\_s end end -end - -macro stats.pair + + concat min min_S + concat Min Min_S +end + +macro grid.load.sersic.test 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 - + echo "USAGE: grid.load.devexp.test (srcdir) (fitdir)" + break + end + + $fields_bt = X Y M T R r I + $fields_ot = Pmag Kmag Amag IDx MTot + + foreach field $fields_bt + foreach set in ot + delete -q $field\$set\_dev_s + delete -q $field\$set\_exp_s + end + end + foreach field $fields_ot + delete -q $field\_dev_s + delete -q $field\_exp_s + end + + $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 + sprint name "sersic.%02d" $Nseq + cmf.load.sersic.test $1/$name.dat $2/$name.fit.cmf + $Nseq ++ + end + end + end + end +end + +# I have run DEV and EXP against input models of type SERSIC +macro cmf.load.sersic.test + if ($0 != 3) + echo "USAGE: cmf.load.sersic.test (dat) (cmf)" + break + end + + # input parameters 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 + read Xin_all 1 Yin_all 2 Fin_all 3 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9 IndexIn_all 10 + + # galaxies only + subset Xin = Xin_all if (Type == 1) + subset Yin = Yin_all if (Type == 1) + subset Min = Min_all if (Type == 1) + subset Fin = Fin_all if (Type == 1) + set min = -2.5*log(Fin) + + 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) + + subset IndexIn = IndexIn_all if (Type == 1) + + $TYPE_S = 83 + $TYPE_D = 68 + $TYPE_E = 69 + + data $2 + + # load measured values from xfit + read -fits Chip.xfit IPP_IDET X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE PSF_INST_MAG AP_MAG KRON_MAG + set EXT_PAR_07 = (MODEL_TYPE:9 == $TYPE_D)*4 + (MODEL_TYPE:9 == $TYPE_E) + + 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 + set IPP_IDET_EXT = IPP_IDET + + 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 + + reindex Pmag_m = PSF_INST_MAG using index1 + reindex Kmag_m = KRON_MAG using index1 + reindex Amag_m = AP_MAG using index1 + + reindex IDx_m = IPP_IDET_EXT using index1 + + reindex MTot_m = MODEL_TYPE:9 using index1 + + reindex Iot_m = EXT_PAR_07 using index1 + reindex Iin_m = IndexIn using index2 + + # load moments and other kron values from Chip.psf + read -fits Chip.psf IPP_IDET X_PSF Y_PSF MOMENTS_XX MOMENTS_XY MOMENTS_YY KRON_FLUX_INNER MOMENTS_R1 MOMENTS_RH + set IPP_IDET_PSF = IPP_IDET + + join -outer IPP_IDET_PSF IDx_m + reindex Xp = X_PSF using index2 + reindex Yp = Y_PSF using index2 + + reindex Mxx_m = MOMENTS_XX using index2 + reindex Mxy_m = MOMENTS_XY using index2 + reindex Myy_m = MOMENTS_YY using index2 + reindex Mr1_m = MOMENTS_R1 using index2 + reindex Mrh_m = MOMENTS_RH using index2 + reindex Kfi_m = KRON_FLUX_INNER using index2 + set Kmi_m = -2.5*log(Kfi_m) + + $fields_bt = X Y M T R r I + $fields_ot = Pmag Kmag Amag IDx MTot Mxx Mxy Myy Mr1 Mrh Kmi + + foreach field $fields_ot + subset $field\_exp_m = $field\_m where (MTot_m == $TYPE_E) + subset $field\_dev_m = $field\_m where (MTot_m == $TYPE_D) + end + foreach field $fields_bt + foreach set in ot + subset $field\$set\_exp_m = $field\$set\_m where (MTot_m == $TYPE_E) + subset $field\$set\_dev_m = $field\$set\_m where (MTot_m == $TYPE_D) + end + end + + join IDx_exp_m IDx_dev_m + foreach field $fields_ot + reindex $field\_exp_mr = $field\_exp_m using index1 + reindex $field\_dev_mr = $field\_dev_m using index2 + end + foreach field $fields_bt + foreach set in ot + reindex $field\$set\_exp_mr = $field\$set\_exp_m using index1 + reindex $field\$set\_dev_mr = $field\$set\_dev_m using index2 + end + end + + # concat + foreach field $fields_bt + foreach set in ot + concat $field\$set\_dev_mr $field\$set\_dev_s + concat $field\$set\_exp_mr $field\$set\_exp_s + end + end + foreach field $fields_ot + concat $field\_dev_mr $field\_dev_s + concat $field\_exp_mr $field\_exp_s + end + + concat min min_S + concat Min Min_S end @@ -659,39 +794,301 @@ end -macro plot.angles +macro sersic.integral if ($0 != 2) - echo "USAGE: plot.angles (file.dat)" + echo "sersic.integral (index)" + break + end + + $index = $1 + create r 0.0 200.0 0.001 + set q = r^(1/$index) + set f = exp(-q) + set rf = f*r + integrate r rf r[0] r[-1] +end + +macro sersic.integral.rmax + if ($0 != 4) + echo "sersic.integral.rmax (index) (kappa) (rmax)" + echo "kappa is a guess for kappa" + break + end + + local index kappa + + $index = $1 + $kappa = $2 + + create r 0.0 $3 0.003 + set q = r^(1/$index) + set f = exp(-$kappa*q) + set rf = f*r + integrate r rf r[0] r[-1] +end + +# S(r) = exp(-kappa*(r/Reff)^(1/index)) +# integrate S(r) r dr [ignores 2pi and change-of-variable factors (Rmaj*Rmin)] +macro sersic.integral.reff.rmax + if ($0 != 5) + echo "sersic.integral.rmax (index) (kappa) (reff) (rmax)" + echo "kappa is a guess for kappa" + break + end + + local index kappa Reff + + $index = $1 + $kappa = $2 + $Reff = $3 + + create r 0.0 $4 0.01 + set q = (r/$Reff)^(1/$index) + set f = exp(-$kappa*q) + set rf = f*r + integrate r rf r[0] r[-1] +end + +# integrate to r = reff (rho = 1), applying kappa +macro sersic.integral.reff + if ($0 != 3) + echo "sersic.integral.reff (index) (kappa)" + break + end + + local index kappa + + $index = $1 + $kappa = $2 + create r 0.0 100.0 0.001 + set q = r^(1/$index) + set f = exp(-$kappa*q) + set rf = f*r + integrate r rf r[0] r[-1] + $S0 = $sum + integrate r rf r[0] 1.0 + $S1 = $sum + echo $S0 $S1 {$S1 / $S0} +end + +macro sersic.integral.find.reff + if ($0 != 3) + echo "sersic.integral.find.reff (index) (kappa)" + echo "kappa is a guess for kappa" + break + end + + local kappa index Rmax i + $index = $1 + $kappa = $1 + + $Rmin = 0; $Vmin = 0 + $Rmax = 4000 + sersic.integral.rmax $index $kappa $Rmax + $Vmax = $sum + + $Rtry = 0.5*($Rmin + $Rmax) + + while (abs($Rmax - $Rmin) > 0.001) + sersic.integral.rmax $index $kappa $Rtry + $Vtry = $sum + + $Rold = $Rtry + if ($Vtry > 0.5*$Vmax) + $Rtry = 0.5*($Rmin + $Rtry) + $Rmax = $Rold + else + $Rtry = 0.5*($Rmax + $Rtry) + $Rmin = $Rold + end + # fprintf "%5.2f %5.2f %5.2f | %5.2f | %6.3f %6.3f | %6.3f" $Rmin $Rold $Rmax $Rtry {$Vtry / $Vmax} $Vmax $Vtry + end + # echo $Rtry + # echo {$kappa * $Rtry ^ (1.0/$index)} + $KappaReal = $kappa * $Rtry ^ (1.0/$index) +end + +macro sersic.integral.index + + local i kappa + + delete -q sersic_sum sersic_norm + # vlist idx 0.5 0.75 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 + # create idx 0.5 5.1 0.1 + vlist idx 1 2 4 + for i 0 idx[] + $kappa = -0.275552 + 1.972625*idx[$i] + 0.003487*idx[$i]*idx[$i] + sersic.integral.rmax idx[$i] $kappa 300 + concat $sum sersic_sum + + $bn = 1.9992*idx[$i] - 0.3271; + $Io = exp($bn); + + # the integral of a Sersic (supposedly) has an analytical form as follows: + $logGamma = lgamma(2.0*idx[$i]); + $bnFactor = $bn^(2.0*idx[$i]); + $norm = idx[$i] * $Io * exp($logGamma) / $bnFactor; + concat $norm sersic_norm + + echo idx[$i] $kappa $sum $logGamma $bn $bnFactor $norm + end + set lsersic_sum = ln(sersic_sum) + break + + $order = 2 + delete -q sersic_fit fit_idx + + # fit to ranges: + delete lsersic_fit fit_idx + vlist bound 0.0 1.0 2.0 3.0 4.0 5.5 + for i 0 {bound[] - 1} + subset idxs = idx if (idx >= bound[$i]) && (idx < bound[$i+1]) + subset sums = lsersic_sum if (idx >= bound[$i]) && (idx < bound[$i+1]) + fit idxs sums $order + applyfit idxs fits + concat fits lsersic_fit + concat idxs fit_idx + end + set sersic_fit = exp(lsersic_fit) + + lim -n 1 idx sersic_sum; clear; box; plot idx sersic_sum; plot -c blue -pt 7 idx sersic_norm; plot -x 0 -c red fit_idx sersic_fit + lim -n 2 idx lsersic_sum; clear; box; plot idx lsersic_sum; plot -x 0 -c red fit_idx lsersic_fit +end + +macro sersic.kappa.index + + local i + + delete -q sersic_kappa + # vlist idx 0.5 0.75 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 + create idx 0.5 5.1 0.1 + for i 0 idx[] + sersic.integral.find.reff idx[$i] 10 + concat $KappaReal sersic_kappa + end + + # $order = 2 + # delete -q sersic_kfit kfit_idx + # + # # fit to ranges: + # vlist bound 0.0 1.0 2.0 3.0 4.0 5.5 + # for i 0 {bound[] - 1} + # subset idxs = idx if (idx >= bound[$i]) && (idx < bound[$i+1]) + # subset sums = lsersic_sum if (idx >= bound[$i]) && (idx < bound[$i+1]) + # fit idxs sums $order + # applyfit idxs fits + # concat fits lsersic_fit + # concat idxs fit_idx + # end + # set sersic_fit = exp(lsersic_fit) + + lim -n 1 idx sersic_kappa; clear; box; plot idx sersic_kappa; # plot -c blue -pt 7 idx sersic_norm; plot -x 0 -c red fit_idx sersic_fit + # lim -n 2 idx lsersic_sum; clear; box; plot idx lsersic_sum; plot -x 0 -c red fit_idx lsersic_fit + + # kappa(index) : + # y = -0.275552 x^0 1.972625 x^1 0.003487 x^2 +end + +macro sersic.central.pixel + if ($0 != 3) + echo "USAGE: sersic.central.pixel (index) (Reff)" + break + end + + local index Reff + + $index = $1 + $Reff = $2 + + $kappa = -0.275552 + 1.972625*$index + 0.003487 * $index^2 + # echo "kappa: $kappa" + + sersic.norm $index + sersic.integral.reff.rmax $index $kappa $Reff {10*$Reff}; set sumFull = $sum; echo $sum {$Reff^2*$myNorm} {$sum / ($Reff^2*$myNorm)} + # sersic.integral.reff.rmax $index $kappa $Reff $Reff ; # echo $sum {$sum / $sumFull} + # sersic.integral.reff.rmax $index $kappa $Reff 0.1 ; # echo $sum {$sum / $sumFull} + # sersic.integral.reff.rmax $index $kappa $Reff 0.2 ; # echo $sum {$sum / $sumFull} + # sersic.integral.reff.rmax $index $kappa $Reff 0.5 ; # echo $sum {$sum / $sumFull} + # sersic.integral.reff.rmax $index $kappa $Reff 1.0 ; # echo $sum {$sum / $sumFull} + # sersic.integral.reff.rmax $index $kappa $Reff 2.0 ; # echo $sum {$sum / $sumFull} + sersic.integral.reff.rmax $index $kappa $Reff 0.564 ; echo $index $Reff {$sum / $sumFull} + +end + +macro sersic.norm + if ($0 != 2) + echo "USAGE: sersic.norm (index)" + break + end + + local index + $index = $1 + + if (($index >= 0.0) && ($index < 1.0)) + $norm = 0.201545 - 0.950965 * $index - 0.315248 * $index^2 + echo $norm {exp($norm)} + $myNorm = exp($norm) + return + end + + if (($index >= 1.0) && ($index < 2.0)) + $norm = 0.402084 - 1.357775 * $index - 0.105102 * $index^2 + echo $norm {exp($norm)} + $myNorm = exp($norm) + return + end + + if (($index >= 2.0) && ($index < 3.0)) + $norm = 0.619093 - 1.591674 * $index - 0.041576 * $index^2 + echo $norm {exp($norm)} + $myNorm = exp($norm) + return + end + + if (($index >= 3.0) && ($index < 4.0)) + $norm = 0.770263 - 1.696421 * $index - 0.023363 * $index^2 + echo $norm {exp($norm)} + $myNorm = exp($norm) + return + end + + if (($index >= 4.0) && ($index < 5.5)) + $norm = 0.885891 - 1.755684 * $index - 0.015753 * $index^2 + echo $norm {exp($norm)} + $myNorm = exp($norm) + return + end +end + +macro load.model.im + if ($0 != 2) + echo "USAGE: load.model.im (N)" + break + end + + rd obj$1 obj.$1.fits + rd cnv$1 cnv.$1.fits + rd var$1 var.$1.fits + rd msk$1 msk.$1.fits + for i 1 7 + rd dpar$i.$1 dpar.$i.$1.fits + end + set dC$1 = (msk$1 == 0) * (obj$1 - cnv$1)^2 / var$1 + tv dC$1 -0.01 3.0 +end + +macro load.normdata + if ($0 != 2) + echo "USAGE: load.normdata (file)" 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 + read f 3 t 4 m 5 + set M = -2.5*log(f) + subset Mg = M if (t == 1) + subset mg = m if (t == 1) + set dm = mg - Mg + set n = ramp(dm) + lim n dm; clear; box; plot n dm end