叹息线上的密度分布积分

2021-04-11 14:20:48 发布

您现在位置:Python中文网/ 问答频道 /正文

我的问题是这样的:

我知道密度是球体半径的函数。假设密度rho(1000)和半径(1000)已经通过数值计算。我想找到视线上密度的积分,如下图2D所示,尽管这是一个3D问题: enter image description here

这条视线可以从中心移动到边界。我知道我们需要先沿视线插值密度,然后加起来得到视线上的密度积分。但是谁能给我一些快速插值的方法吗?谢谢您。在

1条回答
网友
1楼 ·

我有下面的实现(假设密度分布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)

相关问题