implement Sersic / double-Sersic models

[pmelchior]

because: why not

[pmelchior]

@AstroJacobLi is interested in a stable parametric radial profile.

I suggest to implement the Spergel Bessel function class (described here) as a numerically more stable alternative to the Sersic. However, unlike the epsilon, phi (strength and orientation) parameterization for elliptical profiles in that paper, I recommend that we use the epsilon1, epsilon2 (2-component shear vector) parameterization (see e.g. equations B.6 and following here)

[pmelchior]

A good starting place is the PointSource morphology implementation here

[AstroJacobLi]

Hi all, I'm trying to implement the Spergel profile in scarlet. I implement the Spergel profile directly in get_model: https://github.com/AstroJacobLi/scarlet/blob/a168cb79292f9e529e8f9d1f6c457f6315d51042/scarlet/profile.py#L184. I use np.meshgrid to evaluate the profile because it’s more intuitive to shear the profile. However, it seems autograd cannot work on np.meshgrid. I got an error from autograd: VJP of meshgrid wrt argnums (0, 1) not defined.

Also really have a hard time making the autograd recognize the parameters of Spergel profile...

[AstroJacobLi]

@pmelchior I noticed that you added normalization to the Spergel morphology (https://github.com/pmelchior/scarlet/blob/53cf193654c5753f4de8da74ee385409a62a9d9e/scarlet/morphology.py#L279). This actually makes it hard to calculate the intrinsic luminosity via $$\Sigma{\nu}(r)=\frac{c{\nu}^{2} L{0}}{2 \pi r{0}^{2}} f{\nu}\left(\frac{c{\nu} r}{r_{0}}\right)$$, since the normalization wipes out the amplitude information. It would be useful to add an attribute such as self.amplitude before normalizing the profile? Or, just don't normalize it?

[pmelchior]

I'm not sure I agree. With normalization $L_0=1$, right?

The normalization is important because it clarifies whether SED or morphology carry luminosity information. With this normalization, it's the former. Without it, it's degenerate and will make the optimizer unhappy.

[AstroJacobLi]

@pmelchior I reply here such that you can see the latex equations. The Spergel profile in surface brightness is defined as $$\Sigma\nu(r) = \frac{c\nu^2 L_0}{2\pi r0^2} f\nu\left(\frac{c_\nu r}{r0}\right)$$ (notice that there's a 2pi that Spergel omitted in eqn. (3) of his paper). In your implementation https://github.com/pmelchior/scarlet/blob/53cf193654c5753f4de8da74ee385409a62a9d9e/scarlet/morphology.py#L443 and before any normalization, you effectively implement $$\Sigma\nu(r) = f\nu\left(\frac{c\nu r}{r0}\right)$$ Therefore, this implicitly set $\frac{c\nu^2 L_0}{2\pi r_0^2} = 1$, thus we can express the total luminosity $L_0 = \frac{2\pi r0^2}{c\nu^2}$. This $L_0$ is different from morpho.sum() because it is a result by integrating the profile to infinity.

Let's look at an example. If we don't normalize, then a Spergel object will be like SED = [1.20914725, 3.00661596, 4.7323393 , 6.28362644] and morph.sum() = 102.39, and L_0 = \frac{2\pi r_0^2}{c_\nu^2} = 157.97. Therefore, the total luminosity in each band will be L_0 * SED = [191.00899103, 474.95512287, 747.56763971, 992.62446796].

If we normalize it, then SED = [123.80529874, 307.84917798, 484.54700727, 643.38420999] and morph.sum() = 1.0 and still L_0 = \frac{2\pi r_0^2}{c_\nu^2} = 157.97. Now we don't know how to calculate the total luminosity. However, if we record the M = morph.sum()=102.39 before normalization, we can still work out the total luminosity to be L_0 * SED / M. So I suggest a "silent" attribute to SpergelMorphology.

[pmelchior]

I see, it's total luminosity you're after. The problem with the implicit normalization is that we can't set the initial SED. If I know that the peak pixel of the model has intensity 1, then I can read off the initial SED from the observed peak pixel. Similarly, if the model morphology sums to one, I can set the SED with point-source photometry. I can't set these if I don't know the normalization of the source. However, we do know the normalization of this morphology (and likely of any morphology based on a radial profile). So, the best way is to use that in the source's __init__ routine. I'll look into it.