视线范围内的密度分布积分
Density profile integral over line of sight
我的问题是这样的:
我知道密度是球体半径的函数。说密度 rho(1000) 和 radius(1000) 已经用数字计算出来了。我想找到一条视线上的密度积分,如下图 2D 所示,尽管它是 3D 问题:
这条视线可以从中心移动到边界。我知道我们需要先沿视线对密度进行插值,然后相加以获得视线上的密度积分。但是任何人都可以给我一些想法如何快速进行插值?谢谢。
我有以下实现(假设密度分布 rho = exp(1-log(1+r/rs)/(r/rs))
):
第一种方法要快得多,因为它不需要处理来自 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)
我的问题是这样的:
我知道密度是球体半径的函数。说密度 rho(1000) 和 radius(1000) 已经用数字计算出来了。我想找到一条视线上的密度积分,如下图 2D 所示,尽管它是 3D 问题:
这条视线可以从中心移动到边界。我知道我们需要先沿视线对密度进行插值,然后相加以获得视线上的密度积分。但是任何人都可以给我一些想法如何快速进行插值?谢谢。
我有以下实现(假设密度分布 rho = exp(1-log(1+r/rs)/(r/rs))
):
第一种方法要快得多,因为它不需要处理来自 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)