Add support for specifying an obstime and obspos per coordinate.

tests/test_evaluate.py

@@ -109,6 +109,48 @@ def test_output_shape() -> None:
     assert test_model.evaluate(skycoord, return_comps=False).shape == (6,)
 
 
+def test_input_shape() -> None:
+    """Test that shape of the input sky_coord is correct and obspos and obstime are correct."""
+    obstimes = time.Time([TEST_TIME.mjd + i for i in range(4)], format="mjd")
+    obsposes = (
+        coords.get_body("earth", obstimes)
+        .transform_to(coords.HeliocentricMeanEcliptic)
+        .cartesian.xyz.to(units.AU)
+    )
+    test_model.evaluate(SkyCoord([20, 30], [30, 40], unit=units.deg, obstime=TEST_TIME))
+    test_model.evaluate(
+        SkyCoord([20, 30, 40, 50], [30, 40, 30, 20], unit=units.deg, obstime=obstimes)
+    )
+    test_model.evaluate(
+        SkyCoord([20, 30, 40, 50], [30, 40, 30, 20], unit=units.deg, obstime=obstimes),
+        obspos=obsposes,
+    )
+
+    # obstime > sky_coord
+    with pytest.raises(ValueError):
+        test_model.evaluate(SkyCoord([20, 30], [30, 40], unit=units.deg, obstime=obstimes))
+
+    # obspos > sky_coord
+    with pytest.raises(ValueError):
+        test_model.evaluate(SkyCoord(20, 30, unit=units.deg, obstime=TEST_TIME), obspos=obsposes)
+
+    # obspos > obstime
+    with pytest.raises(ValueError):
+        test_model.evaluate(SkyCoord(20, 30, unit=units.deg, obstime=TEST_TIME), obspos=obsposes)
+
+    # obstime > obspos
+    with pytest.raises(ValueError):
+        test_model.evaluate(
+            SkyCoord([20, 30, 40, 50], [30, 40, 30, 20], unit=units.deg, obstime=obstimes),
+            obspos=[[1, -0.4, 0.2]] * units.AU,
+        )
+    with pytest.raises(ValueError):
+        test_model.evaluate(
+            SkyCoord([20, 30, 40, 50], [30, 40, 30, 20], unit=units.deg, obstime=obstimes),
+            obspos=obsposes[:2],
+        )
+
+
 def test_return_comps() -> None:
     """Test that the return_comps function works for valid user input."""
     emission_comps = test_model.evaluate(TEST_SKY_COORD, return_comps=True)

zodipy/__init__.py

@@ -1,6 +1,6 @@
+from zodipy.model import Model
 from zodipy.model_registry import model_registry
 from zodipy.number_density import grid_number_density
-from zodipy.zodipy import Model
 
 __all__ = (
     "Model",

zodipy/bodies.py

@@ -5,6 +5,7 @@
 import astropy.coordinates as coords
 import numpy as np
 from astropy import time, units
+from scipy import interpolate
 
 if TYPE_CHECKING:
     import numpy.typing as npt
@@ -29,11 +30,32 @@ def get_sun_earth_moon_barycenter(
 
 def get_earthpos_xyz(obstime: time.Time, ephemeris: str) -> npt.NDArray[np.float64]:
     """Return the sky coordinates of the Earth in the heliocentric frame."""
-    return (
-        coords.get_body("earth", obstime, ephemeris=ephemeris)
+    if obstime.size == 1:
+        return (
+            coords.get_body("earth", obstime, ephemeris=ephemeris)
+            .transform_to(coords.HeliocentricMeanEcliptic)
+            .cartesian.xyz.to_value(units.AU)
+        )
+    return get_interpolated_body_xyz("earth", obstime, ephemeris)
+
+
+def get_interpolated_body_xyz(
+    body: str,
+    obstimes: time.Time,
+    ephemeris: str,
+) -> npt.NDArray[np.float64]:
+    """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")
+
+    bodypos = (
+        coords.get_body(body, times, ephemeris=ephemeris)
         .transform_to(coords.HeliocentricMeanEcliptic)
         .cartesian.xyz.to_value(units.AU)
     )
+    f = interpolate.interp1d(times.mjd, bodypos, axis=-1)
+    return f(obstimes.mjd)
 
 
 def get_obspos_xyz(
@@ -46,12 +68,17 @@ def get_obspos_xyz(
     if isinstance(obspos, str):
         if obspos.lower() == "semb-l2":
             return get_sun_earth_moon_barycenter(earthpos)
+        if obspos.lower() == "earth":
+            return earthpos
+
         try:
-            return (
-                coords.get_body(obspos, obstime, ephemeris=ephemeris)
-                .transform_to(coords.HeliocentricMeanEcliptic)
-                .cartesian.xyz.to_value(units.AU)
-            )
+            if obstime.size == 1:
+                return (
+                    coords.get_body(obspos, obstime, ephemeris=ephemeris)
+                    .transform_to(coords.HeliocentricMeanEcliptic)
+                    .cartesian.xyz.to_value(units.AU)
+                )
+            return get_interpolated_body_xyz(obspos, obstime, ephemeris)
         except KeyError as error:
             valid_obs = [*coords.solar_system_ephemeris.bodies, "semb-l2"]
             msg = f"Invalid observer string: '{obspos}'. Valid observers are: {valid_obs}"

zodipy/brightness.py

@@ -36,7 +36,8 @@ 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 [0, inf] to the true ecliptic positions
+    # Convert the quadrature range from [-1, 1] to the true ecliptic positions
+    # and back again at the end
     R_los = 0.5 * (stop - start) * r + 0.5 * (stop + start)
 
     X_los = R_los * u_los
@@ -46,14 +47,13 @@ def kelsall_brightness_at_step(
     temperature = get_dust_grain_temperature(R_helio, T_0, delta)
     blackbody_emission = np.interp(temperature, *bp_interpolation_table)
     emission = (1 - albedo) * (emissivity * blackbody_emission)
-
     if albedo != 0:
         solar_flux = solar_irradiance / R_helio**2
         scattering_angle = get_scattering_angle(R_los, R_helio, X_los, X_helio)
         phase_function = get_phase_function(scattering_angle, C1, C2, C3)
         emission += albedo * solar_flux * phase_function
 
-    return emission * get_density_function(X_helio)
+    return emission * get_density_function(X_helio) * 0.5 * (stop - start)
 
 
 def rrm_brightness_at_step(
@@ -69,7 +69,8 @@ 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 [0, inf] to the true ecliptic positions
+    # Convert the quadrature range from [-1, 1] to the true ecliptic positions
+    # and back again at the end
     R_los = 0.5 * (stop - start) * r + 0.5 * (stop + start)
 
     X_los = R_los * u_los
@@ -79,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
+    return blackbody_emission * get_density_function(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_EARTH - 0.2, R_EARTH + 0.2),
-    ComponentLabel.FEATURE: (R_EARTH - 0.2, R_EARTH + 0.2),
+    ComponentLabel.RING: (R_0, R_EARTH + 0.2),
+    ComponentLabel.FEATURE: (R_0, R_EARTH + 0.2),
 }
 
 RRM_CUTOFFS: dict[ComponentLabel, tuple[float | np.float64, float | np.float64]] = {
@@ -67,10 +67,12 @@ 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.flatten()
+    x, y, z = obs_pos
     r_obs = np.sqrt(x**2 + y**2 + z**2)
-    if r_obs > cutoff:
-        return np.asarray([np.finfo(float).eps])
+    if (r_obs > cutoff).any():
+        if obs_pos.ndim == 1:
+            return np.array([np.finfo(float).eps])
+        return np.full(obs_pos.shape[-1], np.finfo(float).eps)
 
     u_x, u_y, u_z = unit_vectors
     lon = np.arctan2(u_y, u_x)