diff --git a/zodipy/brightness.py b/zodipy/brightness.py
index d35ca3f..b824d18 100644
--- a/zodipy/brightness.py
+++ b/zodipy/brightness.py
@@ -22,8 +22,6 @@ def kelsall_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],
@@ -39,9 +37,7 @@ def kelsall_brightness_at_step(
 ) -> npt.NDArray[np.float64]:
     """Kelsall uses common line of sight grid from obs to 5.2 AU."""
     # 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)
+    R_los = 0.5 * (stop - start) * r + 0.5 * (stop + start)
 
     X_los = R_los * u_los
     X_helio = X_los + X_obs
@@ -64,8 +60,6 @@ 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],
@@ -76,9 +70,7 @@ def rrm_brightness_at_step(
 ) -> npt.NDArray[np.float64]:
     """RRM is implented with component specific line-of-sight grids."""
     # 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)
+    R_los = 0.5 * (stop - start) * r + 0.5 * (stop + start)
 
     X_los = R_los * u_los
     X_helio = X_los + X_obs
diff --git a/zodipy/line_of_sight.py b/zodipy/line_of_sight.py
index 78c8a5b..4114407 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_laguerre(
+def integrate_leggauss(
     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-Laguerre quadrature."""
-    return np.squeeze(sum(np.exp(x) * fn(x) * w for x, w in zip(points, weights)))
+    return np.squeeze(sum(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 4c70158..f05b5f8 100644
--- a/zodipy/zodipy.py
+++ b/zodipy/zodipy.py
@@ -16,7 +16,7 @@
 from zodipy.bodies import get_earthpos_xyz, get_obspos_xyz
 from zodipy.line_of_sight import (
     get_line_of_sight_range,
-    integrate_gauss_laguerre,
+    integrate_leggauss,
 )
 from zodipy.model_registry import model_registry
 from zodipy.number_density import populate_number_density_with_model
@@ -99,7 +99,7 @@ def __init__(
         self._comp_parameters, self._common_parameters = brightness_callable_dicts
 
         self._ephemeris = ephemeris
-        self._gauss_x, self.gauss_w = np.polynomial.laguerre.laggauss(gauss_quad_degree)
+        self._leggauss_points_and_weights = np.polynomial.legendre.leggauss(gauss_quad_degree)
 
     def evaluate(
         self,
@@ -210,8 +210,6 @@ 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],
                         )
@@ -220,17 +218,15 @@ def evaluate(
 
                     proc_chunks = [
                         pool.apply_async(
-                            integrate_gauss_laguerre,
-                            args=(comp_integrand, self._gauss_x, self.gauss_w),
+                            integrate_leggauss,
+                            args=(comp_integrand, *self._leggauss_points_and_weights),
                         )
                         for comp_integrand in comp_integrands
                     ]
                     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]
-                    )
+                    emission[idx] *= 0.5 * (stop[comp_label] - start[comp_label])
 
         else:
             for idx, comp_label in enumerate(self._interplanetary_dust_model.comps.keys()):
@@ -239,19 +235,15 @@ def evaluate(
                     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],
                 )
-                emission[idx] = integrate_gauss_laguerre(
-                    comp_integrand, self._gauss_x, self.gauss_w
+                emission[idx] = integrate_leggauss(
+                    comp_integrand, *self._leggauss_points_and_weights
                 )
-
                 # Correct for change of integral limits
-                emission[idx] *= (stop[comp_label] - start[comp_label]) / (
-                    self._gauss_x[-1] - self._gauss_x[0]
-                )
+                emission[idx] *= 0.5 * (stop[comp_label] - start[comp_label])
+
         if contains_duplicates:
             emission = emission[:, inverse]