表面密度剖面图

2021-09-27 05:00:45 发布

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

我试图计算表面密度分布,给定不同参数的三维球形密度分布,以便插值,并将其作为它们的函数,稍后用于一些引力透镜建模。然而,循环遍历所有参数并计算二维半径阵列R上的积分非常耗时(对于复杂的剖面,需要花费许多小时)。下面的代码描述了我为实现这一点而遵循的过程。有人有更好的建议来改进python中以下代码的性能吗?非常感谢你!你知道吗

def density_profile(r, par_a, par_b, par_c):
    ...
    return ...

def surface_denstiy_profile(R, density_profile, *args):

    surface_density_array = []

    for R_i in R:

        surface_density_array.append( quad(lambda r : 2.0 * density_profile(r , *args, ) * r / np.sqrt(r**2. - R_i**2.), R_i, np.inf, epsabs=1.49e-03, epsrel=1.49e-03)[0] ) )

    return surface_density_array


surface_density = np.zeros((len(par_a_array), len(par_b_array), len(par_b_array)), dtype = object)

for par_a in range(0, len(par_a_array)) :

    for par_b in range(0, len(par_b_array)):

        for par_c in range(0, len(par_c_array)):

            surface_density[par_a, par_b, par_c] = surface_denstiy_profile(np.logspace(-2., 2., 100), density_profile, par_a_array[par_a], par_b_array[par_b], par_c_array[par_c])