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)
print 'Half radius: ', f_int(integral[-1]/2.)


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(&lt; 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.)