I am trying to compute the surface density profile, given the spherical density profile in 3D for different parameters in order to interpolate and have it as a function of them to use later on for some gravitational lensing modelling. However, looping through all parameters and evaluating the integral over the 2D Radius array R is extremely time consuming (many hours for complicated profiles).
The code that follows describe the procedure I am following in order to do that. Does anyone have a better suggestion to improve the performance of the following code in 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])
quad
and how doesdensity_profile
look like? Also is this with Python 2.7 or Python 3.x? \$\endgroup\$