Fix bug with parallel evaluations. Move building of partial functions…

… to the model initialization. Refactor variable and function names.

README.md

@@ -33,23 +33,20 @@ from astropy.time import Time
 import zodipy
 
 # Initialize a zodiacal light model at a wavelength/frequency or over a bandpass
-model = zodipy.Model(25*u.micron, name="DIRBE")
+model = zodipy.Model(25*u.micron)
 
 # Use Astropy's `SkyCoord` to specify coordinate
 lon = [10, 10.1, 10.2] * u.deg
 lat = [90, 89, 88] * u.deg
-skycoord = SkyCoord(
-    lon,
-    lat,
-    obstime=Time("2022-01-01 12:00:00"),
-    frame="galactic",
-)
+obstimes = Time(["2022-01-01 12:00:00", "2022-01-01 12:01:00", "2022-01-01 12:02:00"])
+
+skycoord = SkyCoord(lon, lat, obstime=obstimes, frame="galactic")
 
 # Compute the zodiacal light as seen from Earth
 emission = model.evaluate(skycoord, obspos="earth")
 
 print(emission)
-#> [27.52410841 27.66581351 27.81270207] MJy / sr
+#> [27.52410841 27.66572294 27.81251906] MJy / sr
 ```
 
 ## Related scientific papers

docs/index.md

@@ -26,23 +26,20 @@ from astropy.time import Time
 import zodipy
 
 # Initialize a zodiacal light model at a wavelength/frequency or over a bandpass
-model = zodipy.Model(25*u.micron, name="DIRBE")
+model = zodipy.Model(25*u.micron)
 
 # Use Astropy's `SkyCoord` to specify coordinate
 lon = [10, 10.1, 10.2] * u.deg
 lat = [90, 89, 88] * u.deg
-skycoord = SkyCoord(
-    lon,
-    lat,
-    obstime=Time("2022-01-01 12:00:00"),
-    frame="galactic",
-)
+obstimes = Time(["2022-01-01 12:00:00", "2022-01-01 12:01:00", "2022-01-01 12:02:00"])
+
+skycoord = SkyCoord(lon, lat, obstime=obstimes, frame="galactic")
 
 # Compute the zodiacal light as seen from Earth
 emission = model.evaluate(skycoord, obspos="earth")
 
 print(emission)
-#> [27.52410841 27.66581351 27.81270207] MJy / sr
+#> [27.52410841 27.66572294 27.81251906] MJy / sr
 ```
 
 For more information on using ZodiPy, see [the usage section](usage.md).

tests/test_evaluate.py

@@ -95,7 +95,7 @@ def test_evaluate_invalid_obspos() -> None:
 
 def test_output_shape() -> None:
     """Test that the return_comps function works for valid user input."""
-    n_comps = test_model._interplanetary_dust_model.ncomps
+    n_comps = test_model._ipd_model.ncomps
     assert test_model.evaluate(TEST_SKY_COORD, return_comps=True).shape == (n_comps, 1)
     assert test_model.evaluate(TEST_SKY_COORD, return_comps=False).shape == (1,)
 
@@ -168,7 +168,7 @@ def test_return_comps() -> None:
     assert_array_equal(emission_comps.sum(axis=0), emission)
 
 
-def test_multiprocessing_nproc() -> None:
+def test_multiprocessing_nproc_inst() -> None:
     """Test that the multiprocessing works with n_proc > 1."""
     model = Model(x=20 * units.micron)
 
@@ -198,3 +198,27 @@ def test_multiprocessing_nproc() -> None:
     emission = model.evaluate(skycoord, nprocesses=1)
 
     assert_array_equal(emission_multi, emission)
+
+
+def test_multiprocessing_nproc_time() -> None:
+    """Test that the multiprocessing works with n_proc > 1."""
+    model = Model(x=20 * units.micron)
+
+    lon = np.random.randint(low=0, high=360, size=10000)
+    lat = np.random.randint(low=-90, high=90, size=10000)
+    obstime = np.linspace(0, 300, 10000) + TEST_TIME.mjd
+    skycoord = SkyCoord(
+        lon,
+        lat,
+        unit=units.deg,
+        obstime=time.Time(obstime, format="mjd"),
+    )
+    emission_multi = model.evaluate(skycoord, nprocesses=4)
+    emission = model.evaluate(skycoord, nprocesses=1)
+    assert_array_equal(emission_multi, emission)
+
+    # model = Model(x=75 * units.micron, name="rrm-experimental")
+    # emission_multi = model.evaluate(skycoord, nprocesses=4)
+    # emission = model.evaluate(skycoord, nprocesses=1)
+
+    # assert_array_equal(emission_multi, emission)

zodipy/blackbody.py

@@ -9,14 +9,14 @@
 MIN_TEMP = 40 * units.K
 MAX_TEMP = 550 * units.K
 N_TEMPS = 100
-temperatures = np.linspace(MIN_TEMP, MAX_TEMP, N_TEMPS)
-blackbody = BlackBody(temperatures)
+TEMPERATURES = np.linspace(MIN_TEMP, MAX_TEMP, N_TEMPS)
+blackbody = BlackBody(TEMPERATURES)
 
 
 def get_dust_grain_temperature(
     R: npt.NDArray[np.float64], T_0: float, delta: float
 ) -> npt.NDArray[np.float64]:
-    """Return the dust grain temperature given a radial distance from the Sun.
+    """Return the dust grain temperature given at a radial distance from the Sun.
 
     Args:
         R: Radial distance from the sun in ecliptic heliocentric coordinates [AU / 1AU].
@@ -43,7 +43,7 @@ def tabulate_blackbody_emission(
 
     return np.asarray(
         [
-            temperatures.to_value(units.K),
+            TEMPERATURES.to_value(units.K),
             tabulated_blackbody_emission.to_value(units.MJy / units.sr),
         ]
     )

zodipy/bodies.py

@@ -47,7 +47,7 @@ def get_interpolated_body_xyz(
     """Return interpolated Earth positions in the heliocentric frame."""
     dt = (1 * units.hour).to_value(units.day)
     t0, t1 = obstimes[0].mjd, obstimes[-1].mjd
-    times = time.Time(np.arange(t0, max(t0 + 365, t1), dt), format="mjd")
+    times = time.Time(np.arange(t0, max(t0 + 366, t1) + dt, dt), format="mjd")
 
     bodypos = (
         coords.get_body(body, times, ephemeris=ephemeris)

zodipy/brightness.py

@@ -9,7 +9,7 @@
 from zodipy.scattering import get_phase_function, get_scattering_angle
 
 if TYPE_CHECKING:
-    from zodipy.number_density import ComponentNumberDensityCallable
+    from zodipy.number_density import NumberDensityFunc
 
 """
 Function that return the zodiacal emission at a step along all lines of sight given
@@ -25,7 +25,7 @@ def kelsall_brightness_at_step(
     X_obs: npt.NDArray[np.float64],
     u_los: npt.NDArray[np.float64],
     bp_interpolation_table: npt.NDArray[np.float64],
-    get_density_function: ComponentNumberDensityCallable,
+    number_density_func: NumberDensityFunc,
     T_0: float,
     delta: float,
     emissivity: np.float64,
@@ -53,7 +53,7 @@ def kelsall_brightness_at_step(
         phase_function = get_phase_function(scattering_angle, C1, C2, C3)
         emission += albedo * solar_flux * phase_function
 
-    return emission * get_density_function(X_helio) * 0.5 * (stop - start)
+    return emission * number_density_func(X_helio) * 0.5 * (stop - start)
 
 
 def rrm_brightness_at_step(
@@ -63,7 +63,7 @@ def rrm_brightness_at_step(
     X_obs: npt.NDArray[np.float64],
     u_los: npt.NDArray[np.float64],
     bp_interpolation_table: npt.NDArray[np.float64],
-    get_density_function: ComponentNumberDensityCallable,
+    number_density_func: NumberDensityFunc,
     T_0: float,
     delta: float,
     calibration: np.float64,
@@ -80,4 +80,4 @@ def rrm_brightness_at_step(
     temperature = get_dust_grain_temperature(R_helio, T_0, delta)
     blackbody_emission = np.interp(temperature, *bp_interpolation_table)
 
-    return blackbody_emission * get_density_function(X_helio) * calibration * 0.5 * (stop - start)
+    return blackbody_emission * number_density_func(X_helio) * calibration * 0.5 * (stop - start)

zodipy/line_of_sight.py

@@ -21,8 +21,8 @@
     ComponentLabel.BAND1: (R_0, R_JUPITER),
     ComponentLabel.BAND2: (R_0, R_JUPITER),
     ComponentLabel.BAND3: (R_0, R_JUPITER),
-    ComponentLabel.RING: (R_0, R_EARTH + 0.2),
-    ComponentLabel.FEATURE: (R_0, R_EARTH + 0.2),
+    ComponentLabel.RING: (R_EARTH - 0.2, R_EARTH + 0.2),
+    ComponentLabel.FEATURE: (R_EARTH - 0.3, R_EARTH + 0.3),
 }
 
 RRM_CUTOFFS: dict[ComponentLabel, tuple[float | np.float64, float | np.float64]] = {
@@ -53,12 +53,12 @@
 
 
 def integrate_leggauss(
-    fn: Callable[[float], npt.NDArray[np.float64]],
+    func: 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(fn(x) * w for x, w in zip(points, weights)))
+    return np.squeeze(sum(func(x) * w for x, w in zip(points, weights)))
 
 
 def get_sphere_intersection(
@@ -67,8 +67,8 @@ def get_sphere_intersection(
     cutoff: float | np.float64,
 ) -> npt.NDArray[np.float64]:
     """Returns the distance from the observer to a heliocentric sphere with radius `cutoff`."""
-    x, y, z = obs_pos
-    r_obs = np.sqrt(x**2 + y**2 + z**2)
+    x_0, y_0, z_0 = obs_pos
+    r_obs = np.sqrt(x_0**2 + y_0**2 + z_0**2)
     if (r_obs > cutoff).any():
         if obs_pos.ndim == 1:
             return np.array([np.finfo(float).eps])
@@ -79,15 +79,15 @@ def get_sphere_intersection(
     lat = np.arcsin(u_z)
 
     cos_lat = np.cos(lat)
-    b = 2 * (x * cos_lat * np.cos(lon) + y * cos_lat * np.sin(lon))
+    b = 2 * (x_0 * cos_lat * np.cos(lon) + y_0 * cos_lat * np.sin(lon))
     c = r_obs**2 - cutoff**2
 
     q = -0.5 * b * (1 + np.sqrt(b**2 - 4 * c) / np.abs(b))
 
     return np.maximum(q, c / q)
 
 
-def get_line_of_sight_range(
+def get_line_of_sight_range_dicts(
     components: Iterable[ComponentLabel],
     unit_vectors: npt.NDArray[np.float64],
     obs_pos: npt.NDArray[np.float64],