Projecting a 3D radial density profile into 2D

Question

In Python, I am using eq. 1.78 in Binney and Tremaine (2008) to project a 3D density radial distribution into 2D -- see also eq. 18 below:

https://www.mpe.mpg.de/~saglia/journals_pdf/denicola2020.pdf

I don't know how to do the integration without fudging the integral limits. I think I need to change variables but I can't see how.

For a 3D density distribution density_stars as a function of 3D radial coordinate r:

#project density profile from rho (3D) to sigma (2D)
R = np.logspace(-3,np.log10(np.max(r)))
sigma = R*0.
rhomod = splrep(r,density_stars)

for i in range(R.size):
    rr = np.logspace(np.log10(R[i]+0.001),np.log10(R[-1]),1000)
    model = 2. * splev(rr,rhomod) * rr / (rr**2. - R[i]**2.)**0.5
    model2 = splrep(rr,model)
    sigma[i] = splint(R[i],R[-1],model2)

But this isn't stable because if I change the 0.001 fudge factor then I get different answers. I think I need to change variables but I can't see how. Thanks so much for your help!