From d14ba5cd3ddb33b059ce15361c832f34430ab603 Mon Sep 17 00:00:00 2001 From: Metin San Date: Tue, 25 Jun 2024 11:47:14 +0200 Subject: [PATCH] Chane from gauss legendre to gauss laguerre quadrature to get better integration of near observer dust where its brightest --- zodipy/brightness.py | 20 +++++++++--- zodipy/line_of_sight.py | 6 ++-- zodipy/zodipy.py | 68 ++++++++++++++++++++++------------------- 3 files changed, 55 insertions(+), 39 deletions(-) diff --git a/zodipy/brightness.py b/zodipy/brightness.py index 0a2e593..d35ca3f 100644 --- a/zodipy/brightness.py +++ b/zodipy/brightness.py @@ -20,8 +20,10 @@ def kelsall_brightness_at_step( r: npt.NDArray[np.float64], - start: np.float64, + start: npt.NDArray[np.float64], stop: npt.NDArray[np.float64], + quad_root_0: float, + quad_root_n: float, X_obs: npt.NDArray[np.float64], u_los: npt.NDArray[np.float64], bp_interpolation_table: npt.NDArray[np.float64], @@ -36,8 +38,11 @@ def kelsall_brightness_at_step( solar_irradiance: np.float64, ) -> npt.NDArray[np.float64]: """Kelsall uses common line of sight grid from obs to 5.2 AU.""" - # Convert the quadrature range from [-1, 1] to the true ecliptic positions - R_los = ((stop - start) / 2) * r + (stop + start) / 2 + # Convert the quadrature range from [0, inf] to the true ecliptic positions + a, b = start, stop + i, j = quad_root_0, quad_root_n + R_los = ((b - a) / (j - i)) * r + (j * a - i * b) / (j - i) + X_los = R_los * u_los X_helio = X_los + X_obs R_helio = np.sqrt(X_helio[0] ** 2 + X_helio[1] ** 2 + X_helio[2] ** 2) @@ -59,6 +64,8 @@ def rrm_brightness_at_step( r: npt.NDArray[np.float64], start: npt.NDArray[np.float64], stop: npt.NDArray[np.float64], + quad_root_0: float, + quad_root_n: float, X_obs: npt.NDArray[np.float64], u_los: npt.NDArray[np.float64], bp_interpolation_table: npt.NDArray[np.float64], @@ -68,8 +75,11 @@ def rrm_brightness_at_step( calibration: np.float64, ) -> npt.NDArray[np.float64]: """RRM is implented with component specific line-of-sight grids.""" - # Convert the quadrature range from [-1, 1] to the true ecliptic positions - R_los = ((stop - start) / 2) * r + (stop + start) / 2 + # Convert the quadrature range from [0, inf] to the true ecliptic positions + a, b = start, stop + i, j = quad_root_0, quad_root_n + R_los = ((b - a) / (j - i)) * r + (j * a - i * b) / (j - i) + X_los = R_los * u_los X_helio = X_los + X_obs R_helio = np.sqrt(X_helio[0] ** 2 + X_helio[1] ** 2 + X_helio[2] ** 2) diff --git a/zodipy/line_of_sight.py b/zodipy/line_of_sight.py index c2d0f76..78c8a5b 100644 --- a/zodipy/line_of_sight.py +++ b/zodipy/line_of_sight.py @@ -52,13 +52,13 @@ COMPONENT_CUTOFFS = {**DIRBE_CUTOFFS, **RRM_CUTOFFS} -def integrate_gauss_legendre( +def integrate_gauss_laguerre( fn: Callable[[float], npt.NDArray[np.float64]], points: npt.NDArray[np.float64], weights: npt.NDArray[np.float64], ) -> npt.NDArray[np.float64]: - """Integrate a function using Gauss-Legendre quadrature.""" - return np.squeeze(sum(fn(x) * w for x, w in zip(points, weights))) + """Integrate a function using Gauss-Laguerre quadrature.""" + return np.squeeze(sum(np.exp(x) * fn(x) * w for x, w in zip(points, weights))) def get_sphere_intersection( diff --git a/zodipy/zodipy.py b/zodipy/zodipy.py index 8fe5407..4c70158 100644 --- a/zodipy/zodipy.py +++ b/zodipy/zodipy.py @@ -4,7 +4,7 @@ import multiprocessing import multiprocessing.pool import platform -from typing import cast +import typing import numpy as np import numpy.typing as npt @@ -14,7 +14,10 @@ from zodipy.blackbody import tabulate_blackbody_emission from zodipy.bodies import get_earthpos_xyz, get_obspos_xyz -from zodipy.line_of_sight import get_line_of_sight_range, integrate_gauss_legendre +from zodipy.line_of_sight import ( + get_line_of_sight_range, + integrate_gauss_laguerre, +) from zodipy.model_registry import model_registry from zodipy.number_density import populate_number_density_with_model from zodipy.unpack_model import get_model_to_dicts_callable @@ -46,7 +49,7 @@ def __init__( (Jy/sr). name: Interplanetary dust model to use. For a list of available models, see https://cosmoglobe.github.io/zodipy/introduction/. Defaults to 'dirbe'. - gauss_quad_degree: Order of the Gaussian-Legendre quadrature used to evaluate the + gauss_quad_degree: Order of the Gaussian-laguerre quadrature used to evaluate the line-of-sight integral in the simulations. Default is 50 points. extrapolate: If `True` all spectral quantities in the selected model are extrapolated to the requested frequencies or wavelengths. If `False`, an exception is raised on @@ -96,7 +99,7 @@ def __init__( self._comp_parameters, self._common_parameters = brightness_callable_dicts self._ephemeris = ephemeris - self._gauss_points_and_weights = np.polynomial.legendre.leggauss(gauss_quad_degree) + self._gauss_x, self.gauss_w = np.polynomial.laguerre.laggauss(gauss_quad_degree) def evaluate( self, @@ -135,7 +138,7 @@ def evaluate( """ try: - obstime = cast(time.Time, skycoord.obstime) + obstime = typing.cast(time.Time, skycoord.obstime) except AttributeError as error: msg = "The input coordinates must be an astropy SkyCoord object." raise TypeError(msg) from error @@ -143,10 +146,9 @@ def evaluate( msg = "The `obstime` attribute of the `SkyCoord` object must be set." raise ValueError(msg) - ncoords_in = skycoord.size if contains_duplicates: _, index, inverse = np.unique( - cast( + typing.cast( list[npt.NDArray[np.float64]], [skycoord.spherical.lon.value, skycoord.spherical.lat.value], ), @@ -154,18 +156,20 @@ def evaluate( return_inverse=True, axis=1, ) - skycoord = cast(coords.SkyCoord, skycoord[index]) # filter out identical coordinates + skycoord = typing.cast( + coords.SkyCoord, skycoord[index] + ) # filter out identical coordinates ncoords = skycoord.size earth_xyz = get_earthpos_xyz(obstime, self._ephemeris) obs_xyz = get_obspos_xyz(obstime, obspos, earth_xyz, self._ephemeris) skycoord = skycoord.transform_to(coords.BarycentricMeanEcliptic) if skycoord.isscalar: - skycoord_xyz = cast( + skycoord_xyz = typing.cast( npt.NDArray[np.float64], skycoord.cartesian.xyz.value[:, np.newaxis] ) else: - skycoord_xyz = cast(npt.NDArray[np.float64], skycoord.cartesian.xyz.value) + skycoord_xyz = typing.cast(npt.NDArray[np.float64], skycoord.cartesian.xyz.value) start, stop = get_line_of_sight_range( components=self._interplanetary_dust_model.comps.keys(), @@ -189,10 +193,10 @@ def evaluate( **self._common_parameters, ) - distribute_to_cores = ncoords > nprocesses and nprocesses > 1 - if distribute_to_cores: - unit_vector_chunks = np.array_split(skycoord_xyz, nprocesses, axis=-1) - integrated_comp_emission = np.zeros((self._interplanetary_dust_model.ncomps, ncoords)) + emission = np.zeros((self._interplanetary_dust_model.ncomps, ncoords)) + dist_to_cores = ncoords > nprocesses and nprocesses > 1 + if dist_to_cores: + skycoord_xyz_splits = np.array_split(skycoord_xyz, nprocesses, axis=-1) with multiprocessing.get_context(_PLATFORM_METHOD).Pool(nprocesses) as pool: for idx, comp_label in enumerate(self._interplanetary_dust_model.comps.keys()): stop_chunks = np.array_split(stop[comp_label], nprocesses, axis=-1) @@ -206,48 +210,50 @@ def evaluate( u_los=np.expand_dims(vec, axis=-1), start=np.expand_dims(start, axis=-1), stop=np.expand_dims(stop, axis=-1), + quad_root_0=self._gauss_x[0], + quad_root_n=self._gauss_x[-1], get_density_function=density_callables[comp_label], **self._comp_parameters[comp_label], ) - for vec, start, stop in zip(unit_vector_chunks, start_chunks, stop_chunks) + for vec, start, stop in zip(skycoord_xyz_splits, start_chunks, stop_chunks) ] proc_chunks = [ pool.apply_async( - integrate_gauss_legendre, - args=(comp_integrand, *self._gauss_points_and_weights), + integrate_gauss_laguerre, + args=(comp_integrand, self._gauss_x, self.gauss_w), ) for comp_integrand in comp_integrands ] - integrated_comp_emission[idx] += ( - np.concatenate([result.get() for result in proc_chunks]) - * 0.5 - * (stop[comp_label] - start[comp_label]) + emission[idx] = np.concatenate([result.get() for result in proc_chunks]) + + # Correct for change of integral limits + emission[idx] *= (stop[comp_label] - start[comp_label]) / ( + self._gauss_x[-1] - self._gauss_x[0] ) else: - integrated_comp_emission = np.zeros((self._interplanetary_dust_model.ncomps, ncoords)) for idx, comp_label in enumerate(self._interplanetary_dust_model.comps.keys()): comp_integrand = functools.partial( common_integrand, u_los=np.expand_dims(skycoord_xyz, axis=-1), start=np.expand_dims(start[comp_label], axis=-1), stop=np.expand_dims(stop[comp_label], axis=-1), + quad_root_0=self._gauss_x[0], + quad_root_n=self._gauss_x[-1], get_density_function=density_callables[comp_label], **self._comp_parameters[comp_label], ) - - integrated_comp_emission[idx] = ( - integrate_gauss_legendre(comp_integrand, *self._gauss_points_and_weights) - * 0.5 - * (stop[comp_label] - start[comp_label]) + emission[idx] = integrate_gauss_laguerre( + comp_integrand, self._gauss_x, self.gauss_w ) + # Correct for change of integral limits + emission[idx] *= (stop[comp_label] - start[comp_label]) / ( + self._gauss_x[-1] - self._gauss_x[0] + ) if contains_duplicates: - emission = np.zeros((self._interplanetary_dust_model.ncomps, ncoords_in)) - emission = integrated_comp_emission[:, inverse] - else: - emission = integrated_comp_emission + emission = emission[:, inverse] emission <<= units.MJy / units.sr return emission if return_comps else emission.sum(axis=0)