Density profile integral over line of sight

Question

My question is like this:

I know the density as a function of radius for a sphere numerically. Say density rho(1000) and radius(1000) are already calculated numerically. I want to find the integration of the density over a line of sight, as shown below in 2D, although it is a 3D problem:

This line of sight can move from center to the boundary. I know we need to interpolate the density along the line of sight first, then add together to get the integral of density over the line of sight. But anyone can offer me some idea how to do the interpolation fast? Thank you.

Answer 1

I have the implementation below (assume density profile rho = exp(1-log(1+r/rs)/(r/rs))):

The first approach is much faster because it does not need to deal with the singularity from r/np.sqrt(r**2-r_p**2).

import numpy as np
from scipy import integrate as integrate

### From the definition of the LOS integral
def LOS_integration(rs,r_vir,r_p):  #### radius in kpc
    rho = lambda l: np.exp(1 - np.log(1+np.sqrt(l**2 + r_p**2)/rs)/(np.sqrt(l**2 + r_p**2)/rs))
    result = integrate.quad(rho,0,np.sqrt(r_vir**2-r_p**2),epsabs=1.49e-08, epsrel=1.49e-08)
    return result[0]

integration_vec = np.vectorize(LOS_integration)   ### vectorize the function


###  convert LOS integration to radius integration
def LOS_integration1(rs,r_vir,r_p):  #### radius in kpc
    rho = lambda r: np.exp(1 - np.log(1+r/rs)/(r/rs)) * r/np.sqrt(r**2-r_p**2)  
    ### r/np.sqrt(r**2-r_p**2) is the factor convert from LOS integration to radius integration
    result = integrate.quad(rho,r_p,r_vir,epsabs=1.49e-08, epsrel=1.49e-08)
    return result[0]

integration1_vec = np.vectorize(LOS_integration1)