GalSim-developers/GalSim — pull request #433 by rmjarvis (2013-06-11 23:36:18+00:00)
Title: #426 Improve behavior of Sersic class

There are a number of improvements to the Sersic class on this branch triggered by David noticing some ripples in images that were higher amplitude than they should have been. But then I pulled in some other improvements from other branches in order to get them onto master sooner, so they would be used by the various testing efforts that are ongoing.

- The proximate cause of the feature David was seeing was the lookup table for the Sersic Hankel transform being a linear interpolation, so discontinuous second derivative. Switching to a spline helped a lot.
- It improved even more when I switched to interpolating F(logk) * k^2, rather than just F(logk). This was inspired by David pointing out that the asymptotic expansion of F is a/k^2 + b/k^3 +...
- I added code to estimate a,b as it is building up the FT values for the lookup table. It now uses this approximation for the high-k limit (`k > _maxk`) once the approximation is found to be accurate at the `kvalue_accuracy` level. This speeds up the code a lot, since it can stop building the lookup table much earlier than it used to.
- From #325, I took the additional photon shooting GSParams and the use of these in the Interpolants, OneDimensionalDeviate, etc. (This also accidentally pulled in some versions of Barney's `test_sersic_highn` scripts. You can ignore them for this PR I suppose.)
- This also includes a switch to using GSParamsPtr, since the old way the GSParams were used in LRUCache actually had a bug. It turned up during testing here, and I think it's the same bug that Barney reported in #422.
- From #343, I took the improvement of adding explicit splits at the zeros of j0(kr) for the integrals. This helps improve their accuracy and also speed up the convergence.
- Since I was doing all this work on the Sersic class, I also refactored the back end (the SersicInfo class) to simplify things a lot. Mostly I switched it to use r0 as its basic unit rather than re. I also broke up parts of what was the constructor into more discrete units, which meant that its use in SersicImpl was able to be simplified considerably.
- I added `table_spacing` as a new GSParams item. It represents the spacing to use in lookup tables relative to the normal calculated value (by whatever calculation is done in each case). So it has a default value of 1. If it is <1, the lookup table will be more accurate at the expense of high time to calculate the values.
- I added `stepk_minimum_hlr` as another GSParams item. For the stepk calculation, we typically start with a calculation based on `alias_threshold`. But then, each class imposed some minimum value, e.g. 4*sigma for Gaussian. These seemed rather arbitrary to me, so I regularized them with a multiple of the half-light-radius for all profiles that have this feature. The default is 5, which roughly comports with the typical values that the various classes had been using.
- I added unit tests that `Sersic(n=0.5)` and `Sersic(n=1)` agree with `Gaussian` and `Exponential` respectively both in xValues and kValues. (They didn't used to!)
- I changed the `test_rescale` unit test from using `wmult` to using `alias_threshold` and `maxk_threshold`. `wmult` really only affects the effective maxk, not stepk. So it is not generally sufficient to improve the flux values. (i.e. The accuracy starts to level out above a certain value of `wmult`, since the relevant parameter that needs to drop becomes `alias_threshold`.)