# Retrieving input effective radius from Sérsic profile

I'm playing around with Sérsic profiles, and would like to retrieve the effective (half-light) radius that I put in by integration of the profile. I'm not managing to, so perhaps someone can help me out.

I create the profile using astropy:

from astropy.modeling.functional_models import Sersic1D
from scipy.integrate import simps
Sersic = Sersic1D(amplitude=1,r_eff=1,n=1)


I can then make a light profile:

x = np.arange(1000)
I = Sersic(x)


Then we can perform an integral over this I(r):

integral = [simps(I[0:i],x[0:i]) for i in np.arange(len(I)-1)+1]


From which one would finally find the radius at which the enclosed flux is half that of the total:

from scipy.interpolate import interp1d

f_int = interp1d(integral,np.arange(len(I)-1)+1)


However, this does not give me the desired result -- the input radius of 1. Why not? Is this a conceptual error or a numerical one - or both?

The conceptual problem is that the half-light radius is the radius of a circle enclosing half the total light, assuming that the object is circular with a radial intensity profile equal to the Sersic function.

So you need to integrate over circular rings, each of which has an area of $2 \pi r \, dr$ and an intensity (or surface brightness) of $S(r)$: $$L(< R) = \int_{0}^{R} 2 \pi r \, S(r) \, dr$$ where $S(r)$ is the Sersic function evaluated at radius $r$.

Instead, you're just doing the integral of a line of fixed width: $$\int_{0}^{R} S(r) \, dr$$

In Python terms, the (pre-calculated) function you integrate should be

I = 2 * np.pi * x * Sersic(x)


not

I = Sersic(x)


(Alex is also correct that you have almost no resolution near the value of r_eff that you are trying to find.)

There are a number of conceptual and numeric problems with your code. They contribute, but I haven't found the conceptual problem yet. I will edit this answer if / when I find it.

To start, you set the half-light radius for an exponential profile at 1, create a profile out to 1000, and with a resolution of 1. You need more resolution at the centre, and the last 990 or so bins (probably more than that) are useless. I changed your x to

x = np.arange(0,10,0.01) # tried with 100, too -- no change


The half-light radius doesn't change if you go from an upper end of 1000 to 100 or even to 10. It does change quite dramatically, though, if you change the resolution from 1 to 0.1 to 0.01. If we increase it to 0.001, it starts to slowly stabilize.

This brings us to the next problem: These are the results I get when I increase the resolution:

• 1 : 1.518
• 0.1 : 5.139
• 0.01: 42.30
• 0.001:414.0

You can see the problem: the interpolation counts bins and assumes a width of 1.

Next, I did a small test to see what I'd get:

import numpy as np
cs = np.cumsum(I) # Array of length I, with cs[n] = sum(I[:n])
total = cs[-1] # Total flux count
halflight = np.where(cs > (total / 2.)) # First entry where I > total / 2


Result is 412, consistent with yours. So, your integration and interpolation works. Right now, I think the problem might be in the definition of the half-light radius (missing pre-factor?). You should be able to find that in the astropy documentation.

Quick check: Changing r_eff to 10, and integrating out to 100 gives a factor 10 larger half-light radius.

Edit before submitting: it looks like you're missing a factor of $(\pi / 2)^2$.