Refactoring

zodipy/_contour.py

@@ -15,7 +15,7 @@
 
 
 def tabulate_density(
-    grid: npt.NDArray[np.floating] | Sequence[npt.NDArray[np.floating]],
+    grid: npt.NDArray[np.float64] | Sequence[npt.NDArray[np.float64]],
     model: str = "DIRBE",
     earth_pos: u.Quantity[u.AU] = DEFAULT_EARTH_POS,
 ) -> npt.NDArray[np.float64]:

zodipy/blackbody.py

@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import numpy as np
 import numpy.typing as npt
 from astropy import units
@@ -28,27 +30,20 @@ def get_dust_grain_temperature(
     return T_0 * R**-delta
 
 
-def tabulate_center_wavelength_bnu(wavelength: units.Quantity) -> npt.NDArray[np.float64]:
-    """Tabulate blackbody specific intensity for a center wavelength."""
-    return np.asarray(
-        [
-            temperatures.to_value(units.K),
-            blackbody(wavelength).to_value(units.MJy / units.sr),
-        ]
-    )
-
-
-def tabulate_bandpass_integrated_bnu(
-    wavelengths: units.Quantity, normalized_weights: units.Quantity
+def tabulate_blackbody_emission(
+    wavelengths: units.Quantity, weights: units.Quantity | None
 ) -> npt.NDArray[np.float64]:
     """Tabulate bandpass integrated blackbody specific intensity for a range of temperatures."""
-    blackbody_emission = blackbody(wavelengths[:, np.newaxis])
-    integrated_blackbody_emission = integrate.trapezoid(
-        normalized_weights * blackbody_emission.transpose(), wavelengths
-    )
+    if weights is None:
+        tabulated_blackbody_emission = blackbody(wavelengths)
+    else:
+        tabulated_blackbody_emission = integrate.trapezoid(
+            weights * blackbody(wavelengths[:, np.newaxis]).transpose(), wavelengths
+        )
+
     return np.asarray(
         [
             temperatures.to_value(units.K),
-            integrated_blackbody_emission.to_value(units.MJy / units.sr),
+            tabulated_blackbody_emission.to_value(units.MJy / units.sr),
         ]
     )

zodipy/bodies.py

@@ -0,0 +1,53 @@
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+import astropy.coordinates as coords
+import numpy as np
+from astropy import time, units
+
+from zodipy._constants import DISTANCE_FROM_EARTH_TO_SEMB_L2
+
+if TYPE_CHECKING:
+    import numpy.typing as npt
+
+
+def get_sun_earth_moon_barycenter(
+    earthpos: npt.NDArray[np.float64],
+) -> npt.NDArray[np.float64]:
+    """Return a SkyCoord of the heliocentric position of the SEMB-L2 point.
+
+    Note that this is an approximate position, as the SEMB-L2 point is not included in
+    any of the current available ephemerides. We assume that SEMB-L2 is at all times
+    located at a fixed distance from Earth along the vector pointing to Earth from the Sun.
+    """
+    earth_distance = np.linalg.norm(earthpos)
+    SEMB_L2_distance = earth_distance + DISTANCE_FROM_EARTH_TO_SEMB_L2
+    earth_unit_vector = earthpos / earth_distance
+
+    return earth_unit_vector * SEMB_L2_distance
+
+
+def get_earthpos(obs_time: time.Time, ephemeris: str) -> npt.NDArray[np.float64]:
+    """Return the sky coordinates of the Earth in the heliocentric frame."""
+    return (
+        coords.get_body("earth", obs_time, ephemeris=ephemeris)
+        .transform_to(coords.HeliocentricMeanEcliptic)
+        .cartesian.xyz.to_value(units.AU)
+    )
+
+
+def get_obspos(
+    obs: str,
+    obstime: time.Time,
+    earthpos: npt.NDArray[np.float64],
+    ephemeris: str,
+) -> npt.NDArray[np.float64]:
+    """Return the sky coordinates of the observer in the heliocentric frame."""
+    if obs.lower() == "semb-l2":
+        return get_sun_earth_moon_barycenter(earthpos)
+    return (
+        coords.get_body(obs, obstime, ephemeris=ephemeris)
+        .transform_to(coords.HeliocentricMeanEcliptic)
+        .cartesian.xyz.to_value(units.AU)
+    )

zodipy/brightness.py

@@ -30,9 +30,9 @@ def kelsall_brightness_at_step(
     delta: float,
     emissivity: np.float64,
     albedo: np.float64,
-    C1: float,
-    C2: float,
-    C3: float,
+    C1: np.float64,
+    C2: np.float64,
+    C3: np.float64,
     solar_irradiance: np.float64,
 ) -> npt.NDArray[np.float64]:
     """Kelsall uses common line of sight grid from obs to 5.2 AU."""
@@ -78,9 +78,3 @@ def rrm_brightness_at_step(
     blackbody_emission = np.interp(temperature, *bp_interpolation_table)
 
     return blackbody_emission * get_density_function(X_helio) * calibration
-
-
-# EMISSION_MAPPING: dict[type[InterplanetaryDustModel], GetCompEmissionAtStepFn] = {
-#     Kelsall: kelsall,
-#     RRM: rrm,
-# }

zodipy/line_of_sight.py

@@ -11,7 +11,7 @@
 if TYPE_CHECKING:
     import numpy.typing as npt
 
-DIRBE_CUTOFFS: dict[ComponentLabel, tuple[float | np.floating, float | np.floating]] = {
+DIRBE_CUTOFFS: dict[ComponentLabel, tuple[float | np.float64, float | np.float64]] = {
     ComponentLabel.CLOUD: (R_0, R_JUPITER),
     ComponentLabel.BAND1: (R_0, R_JUPITER),
     ComponentLabel.BAND2: (R_0, R_JUPITER),
@@ -20,7 +20,7 @@
     ComponentLabel.FEATURE: (R_EARTH - 0.2, R_EARTH + 0.2),
 }
 
-RRM_CUTOFFS: dict[ComponentLabel, tuple[float | np.floating, float | np.floating]] = {
+RRM_CUTOFFS: dict[ComponentLabel, tuple[float | np.float64, float | np.float64]] = {
     ComponentLabel.FAN: (R_0, RRM[ComponentLabel.FAN].R_outer),  # type: ignore
     ComponentLabel.INNER_NARROW_BAND: (
         RRM[ComponentLabel.INNER_NARROW_BAND].R_inner,  # type: ignore

zodipy/number_density.py

@@ -468,7 +468,6 @@ def construct_density_partials_comps(
     Return a tuple of the density expressions above which has been prepopulated with
     model and configuration parameters, leaving only the `X_helio` argument to be supplied.
     Raises exception for incorrectly defined components or component density functions.
-
     """
     partial_density_funcs: dict[ComponentLabel, ComponentNumberDensityCallable] = {}
     for comp_label, comp in comps.items():

zodipy/scattering.py

@@ -33,7 +33,7 @@ def get_scattering_angle(
 
 
 def get_phase_function(
-    Theta: npt.NDArray[np.float64], C1: float, C2: float, C3: float
+    Theta: npt.NDArray[np.float64], C1: np.float64, C2: np.float64, C3: np.float64
 ) -> npt.NDArray[np.float64]:
     """Return the phase function.
 

zodipy/interpolate.py → zodipy/unpack_model.py

@@ -12,13 +12,16 @@
 
 CompParamDict = dict[ComponentLabel, dict[str, Any]]
 CommonParamDict = dict[str, Any]
+UnpackedModelDicts = tuple[CompParamDict, CommonParamDict]
+T = TypeVar("T", bound=ZodiacalLightModel)
+UnpackModelCallable = Callable[[units.Quantity, Union[units.Quantity, None], T], UnpackedModelDicts]
 
 
-def kelsall_params_to_dicts(
+def unpack_kelsall(
     wavelengths: units.Quantity,
     weights: units.Quantity | None,
     model: Kelsall,
-) -> tuple[CompParamDict, CommonParamDict]:
+) -> UnpackedModelDicts:
     """InterplantaryDustModelToDicts implementation for Kelsall model."""
     model_spectrum = model.spectrum.to(wavelengths.unit, equivalencies=units.spectral())
 
@@ -82,11 +85,11 @@ def kelsall_params_to_dicts(
     return comp_params, common_params
 
 
-def rrm_params_to_dicts(
+def unpack_rrm(
     wavelengths: units.Quantity,
     weights: units.Quantity | None,
     model: RRM,
-) -> tuple[CompParamDict, CommonParamDict]:
+) -> UnpackedModelDicts:
     """InterplantaryDustModelToDicts implementation for Kelsall model."""
     model_spectrum = model.spectrum.to(wavelengths.unit, equivalencies=units.spectral())
 
@@ -121,20 +124,14 @@ def interpolate_spectral_parameter(
     return interpolated_parameter
 
 
-T = TypeVar("T", contravariant=True, bound=ZodiacalLightModel)
-CallableModelToDicts = Callable[
-    [units.Quantity, Union[units.Quantity, None], T], tuple[CompParamDict, CommonParamDict]
-]
-
-
-MODEL_INTERPOLATION_MAPPING: dict[type[ZodiacalLightModel], CallableModelToDicts] = {
-    Kelsall: kelsall_params_to_dicts,
-    RRM: rrm_params_to_dicts,
+model_unpack_mapping: dict[type[ZodiacalLightModel], UnpackModelCallable] = {
+    Kelsall: unpack_kelsall,
+    RRM: unpack_rrm,
 }
 
 
 def get_model_to_dicts_callable(
     model: ZodiacalLightModel,
-) -> CallableModelToDicts:
+) -> UnpackModelCallable:
     """Get the appropriate parameter unpacker for the model."""
-    return MODEL_INTERPOLATION_MAPPING[type(model)]
+    return model_unpack_mapping[type(model)]

zodipy/zodiacal_component.py

@@ -1,6 +1,6 @@
 from __future__ import annotations
 
-from abc import ABC
+import abc
 from dataclasses import dataclass, field
 from enum import Enum
 from typing import TYPE_CHECKING
@@ -12,7 +12,7 @@
 
 
 @dataclass
-class ZodiacalComponent(ABC):
+class ZodiacalComponent(abc.ABC):
     """Base class for storing common model parameters for zodiacal components.
 
     Args:

zodipy/zodiacal_light_model.py

@@ -100,7 +100,7 @@ def brightness_at_step_callable(cls) -> BrightnessAtStepCallable:
 
 
 @dataclass
-class InterplanetaryDustModelRegistry:
+class ModelRegistry:
     """Container for registered models."""
 
     _registry: dict[str, ZodiacalLightModel] = field(init=False, default_factory=dict)
@@ -133,4 +133,4 @@ def get_model(self, name: str) -> ZodiacalLightModel:
         return self._registry[name]
 
 
-model_registry = InterplanetaryDustModelRegistry()
+model_registry = ModelRegistry()