Implemented the factory method pattern to models

zodipy/_coordinates.py

@@ -28,14 +28,14 @@ def get_target_coordinates(
             Either a list of epochs in JD or MJD format or a dictionary
             defining a range of times and dates; the range dictionary has to
             be of the form {``'start'``:'YYYY-MM-DD [HH:MM:SS]',
-            ``'stop'``:'YYYY-MM-DD [HH:MM:SS]', ``'step'``:'n[y|d|m|s]'}.
-            Epoch timescales depend on the type of query performed: UTC for
-            ephemerides queries, TDB for element queries, CT for vector queries.
+            ``'stop'``:'YYYY-MM-DD [HH:MM:SS]', ``'step'``:'n[y|d|h|m|s]'}.
             If no epochs are provided, the current time is used.
 
     Returns
     -------
-        Heliocentric cartesian coordinates of the target.
+    coordinates
+        Heliocentric cartesian coordinates of the target. The shape is 
+        (`n_observations`, 3)
     """
 
     if target.lower() in TARGET_ALIASES:

zodipy/_functions.py

@@ -1,12 +1,11 @@
-from typing import Union
-
 import astropy.constants as const
 import numpy as np
 
+DELTA = 0.4668626
+T_0 = 286
+
 
-def blackbody_emission(
-    T: Union[float, np.ndarray], freq: float
-) -> Union[float, np.ndarray]:    
+def blackbody_emission(T: np.ndarray, freq: float) -> np.ndarray:    
     """Returns the blackbody emission.
     
     Assumes the frequency to be in units of GHz.
@@ -31,8 +30,8 @@ def blackbody_emission(
 
 
 def interplanetary_temperature(
-    R: Union[float, np.ndarray], T_0: float = 286, delta: float = 0.4668626
-) -> Union[float, np.ndarray]:
+    R: np.ndarray, T_0: float = T_0, delta: float = DELTA
+) -> np.ndarray:
     """Returns the interplanetary temperature as a function of radius.
     
     Parameters

zodipy/_integration.py

@@ -1,11 +1,11 @@
-from dataclasses import dataclass
 from typing import Callable
 
 import numpy as np
-import scipy
+
+ZERO = np.finfo(float).eps
+RADIAL_CUTOFF = 6
 
 
-@dataclass
 class IntegrationConfig:
     """Config object for the LOS/shell integration in the simulation.
     
@@ -20,42 +20,33 @@ class IntegrationConfig:
         Function which performs the integration.
     """
 
-    R_max : float
-    n : int
-    integrator : Callable[[float, float, int], np.ndarray] = np.trapz
-
-    @property
-    def R(self) -> np.ndarray:
-        """Linearly spaced grid of distances to shells around the observer."""
-
-        ZERO = np.finfo(float).eps
-        return np.expand_dims(np.linspace(ZERO, self.R_max, self.n), axis=1)
-
-    @property
-    def dR(self) -> np.ndarray:
-        """Distance between grid points in R"""
-
-        return np.diff(self.R)
-
+    def __init__(
+        self, 
+        R_max: float, 
+        n:int, 
+        integrator: Callable[[float, float, int], np.ndarray] = np.trapz
+    ):
+        self.R = np.expand_dims(np.linspace(ZERO, R_max, n), axis=1)   
+        self.dR = np.diff(self.R)
+        self.integrator = integrator
 
-_CUT_OFF = 6
 
 DEFAULT = {
-    'cloud': IntegrationConfig(R_max=_CUT_OFF, n=250),
-    'band1': IntegrationConfig(R_max=_CUT_OFF, n=50),
-    'band2': IntegrationConfig(R_max=_CUT_OFF, n=50),
-    'band3': IntegrationConfig(R_max=_CUT_OFF, n=50),
-    'ring': IntegrationConfig(R_max=2.25, n=50),
-    'feature': IntegrationConfig(R_max=1, n=50),
+    'cloud': IntegrationConfig(RADIAL_CUTOFF, 250),
+    'band1': IntegrationConfig(RADIAL_CUTOFF, 50),
+    'band2': IntegrationConfig(RADIAL_CUTOFF, 50),
+    'band3': IntegrationConfig(RADIAL_CUTOFF, 50),
+    'ring': IntegrationConfig(2.25, 50),
+    'feature': IntegrationConfig(1, 50),
 }
 
 HIGH = {
-    'cloud': IntegrationConfig(R_max=_CUT_OFF, n=500),
-    'band1': IntegrationConfig(R_max=_CUT_OFF, n=500),
-    'band2': IntegrationConfig(R_max=_CUT_OFF, n=500),
-    'band3': IntegrationConfig(R_max=_CUT_OFF, n=500),
-    'ring': IntegrationConfig(R_max=2.25, n=200),
-    'feature': IntegrationConfig(R_max=1, n=200),
+    'cloud': IntegrationConfig(RADIAL_CUTOFF, 500),
+    'band1': IntegrationConfig(RADIAL_CUTOFF, 500),
+    'band2': IntegrationConfig(RADIAL_CUTOFF, 500),
+    'band3': IntegrationConfig(RADIAL_CUTOFF, 500),
+    'ring': IntegrationConfig(2.25, 200),
+    'feature': IntegrationConfig(1, 200),
 }
 
 

zodipy/_model.py

@@ -1,56 +1,57 @@
+from dataclasses import dataclass
 from typing import Iterable, Dict
 
 from zodipy._emissivity import Emissivity
-from zodipy import components
+from zodipy.components import BaseComponent, Cloud, Band, Ring, Feature
 
 
+class ModelFactory:
+    """Class that is responsible for registring and book-keeping models."""
+
+    def __init__(self):
+        self._models = {}
+
+    def register_model(self, name, components, parameters, emissivities):
+        model = init_model(components, parameters, emissivities)
+        self._models[name] = model
+
+    def get_model(self, name):
+        model = self._models.get(name)
+        if model is None:
+            raise ValueError(
+                f'model {name} is not registered. Available models are'
+                f'{list(self._models.keys())}'
+            )
+
+        return model
+
+
+@dataclass
 class InterplanetaryDustModel:
-    """Class that represents a model of the IPD.
-    
-    Attributes
-    ----------
-    components : dict
-        Dictionary containing initialized `zodipy.components.BaseComponent`
-        objects.
-    emissivities : `zodipy._emissivity.Emissivity`
-        Emissivity object.
-    """
-
-    def __init__(
-        self,
-        components: Iterable[str], 
-        parameters: Dict[str, Dict[str, float]], 
-        emissivities: Emissivity
-    ) -> None: 
-        """Initilizes a Model object.
-        
-        Parameters
-        ----------
-        components
-            Iterable containing component labels as strings.
-        parameters
-            Dictionary containing component parameters.
-        emissivities
-            Emissivity object.
-        """
-
-        self.components = self._init_components(components, parameters)
-        self.emissivities = emissivities
-
-    def _init_components(self, comp_labels: Iterable, parameters: dict) -> dict:
-        """Initialize component dictionary."""
-
-        comps = {}
-        for label in comp_labels:
-            if label.startswith('cloud'):
-                comp_type = components.Cloud
-            elif label.startswith('band'):
-                comp_type = components.Band
-            elif label.startswith('ring'):
-                comp_type = components.Ring
-            elif label.startswith('feature'):
-                comp_type = components.Feature
-
-            comps[label] = comp_type(**parameters[label])
-
-        return comps
+    """Data class representing an Interplanetary dust model."""
+
+    components: Dict[str, BaseComponent]
+    emissivities: Emissivity
+
+
+def init_model(
+    components: Iterable[str], 
+    parameters: Dict[str, Dict[str, float]], 
+    emissivities: Emissivity
+) -> InterplanetaryDustModel:
+
+    initialized_components = {}
+    for comp in components:
+        if comp.startswith('cloud'):
+            comp_type = Cloud
+        elif comp.startswith('band'):
+            comp_type = Band
+        elif comp.startswith('ring'):
+            comp_type = Ring
+        elif comp.startswith('feature'):
+            comp_type = Feature
+        initialized_components[comp] = comp_type(**parameters[comp])
+    return InterplanetaryDustModel(
+        components=initialized_components,
+        emissivities=emissivities
+    )

zodipy/components.py

@@ -6,7 +6,7 @@
 
 import numpy as np
 
-import zodipy._functions as F
+from zodipy._functions import interplanetary_temperature, blackbody_emission
 
 
 @dataclass
@@ -172,10 +172,10 @@ def get_emission(
 
         coords, R_helio = self.get_coordinates(X_observer, X_earth, X_unit, R)
         density = self.get_density(*coords)
-        temperature = F.interplanetary_temperature(R_helio)
-        blackbody_emission = F.blackbody_emission(temperature, freq)
+        temperature = interplanetary_temperature(R_helio)
+        emission = blackbody_emission(temperature, freq)
 
-        return blackbody_emission * density
+        return emission * density
 
 
 @dataclass

zodipy/models.py

@@ -1,29 +1,25 @@
+from zodipy._model import ModelFactory
 from zodipy import emissivities
 from zodipy import parameters
-from zodipy._model import InterplanetaryDustModel
 
 
-PLANCK_2013 = InterplanetaryDustModel(
+models = ModelFactory()
+
+models.register_model(
+    name='planck 2013',
     components=('cloud', 'band1', 'band2', 'band3', 'ring', 'feature'),
     parameters=parameters.K98,
     emissivities=emissivities.PLANCK_2013
 )
-
-PLANCK_2015 = InterplanetaryDustModel(
-    components=('cloud', 'band1', 'band2', 'band3'),
+models.register_model(
+    name='planck 2015',
+    components=('cloud', 'band1', 'band2', 'band3',),
     parameters=parameters.K98,
     emissivities=emissivities.PLANCK_2015
 )
-
-PLANCK_2018 = InterplanetaryDustModel(
-    components=('cloud', 'band1', 'band2', 'band3'),
+models.register_model(
+    name='planck 2018',
+    components=('cloud', 'band1', 'band2', 'band3',),
     parameters=parameters.K98,
     emissivities=emissivities.PLANCK_2018
 )
-
-
-MODELS = {
-    'planck 2013' : PLANCK_2013,
-    'planck 2015' : PLANCK_2015,
-    'planck 2018' : PLANCK_2018
-}

zodipy/simulation.py

@@ -1,6 +1,5 @@
 from abc import ABC, abstractmethod
 from dataclasses import dataclass
-from typing import Iterable
 import warnings
 
 import healpy as hp
@@ -32,8 +31,8 @@ class SimulationStrategy(ABC):
 
     model: InterplanetaryDustModel
     integration_config: IntegrationConfig
-    observer_locations: Iterable
-    earth_locations: Iterable
+    observer_locations: np.ndarray
+    earth_locations: np.ndarray
     hit_maps: np.ndarray
 
 
@@ -106,8 +105,8 @@ def simulate(self, nside: int, freq: float) -> np.ndarray:
 class TimeOrderedStrategy(SimulationStrategy):
     """Simulation strategy for time-ordered emission.
     
-    This returns the pixel weighted average of multiple simultaneous 
-    observations.
+    This strategy simulates the sky at multiple different times and returns
+    the pixel weighted average of all observations.
     """
 
     def simulate(self, nside: int, freq: float) -> np.ndarray:

zodipy/zodi.py

@@ -5,7 +5,7 @@
 
 from zodipy._coordinates import get_target_coordinates, change_coordinate_system
 from zodipy._integration import INTEGRATION_CONFIGS
-from zodipy.models import MODELS
+from zodipy.models import models
 from zodipy.simulation import InstantaneousStrategy, TimeOrderedStrategy
 
 
@@ -31,8 +31,6 @@ def __init__(
             defining a range of times and dates; the range dictionary has to
             be of the form {``'start'``:'YYYY-MM-DD [HH:MM:SS]',
             ``'stop'``:'YYYY-MM-DD [HH:MM:SS]', ``'step'``:'n[y|d|h|m|s]'}.
-            Epoch timescales depend on the type of query performed: UTC for
-            ephemerides queries, TDB for element queries, CT for vector queries.
             If no epochs are provided, the current time is used.
         hit_maps
             The number of times each pixel is observed for a given observation.
@@ -46,17 +44,11 @@ def __init__(
             'high'. Defaults to 'default'.
         """
 
-        try:
-            model = MODELS[model.lower()]
-        except KeyError:
-            raise KeyError(
-                f"Model {model!r} not found. Available models are: "
-                f"{list(MODELS.keys())}"
-        )
-        try:
-            integration_config = INTEGRATION_CONFIGS[integration_config.lower()]
-        except KeyError:
-            raise KeyError(
+        model = models.get_model(model)
+
+        integration_config = INTEGRATION_CONFIGS.get(integration_config.lower())
+        if integration_config is None:
+            raise ValueError(
                 f"Config {integration_config!r} not found. Available configs "
                 f"are: {list(INTEGRATION_CONFIGS.keys())}"
         )