added TOD documentation

docs/source/usage.rst

@@ -9,4 +9,5 @@ Usage
    usage/sites.rst
    usage/plans.rst
    usage/simulations.rst
+   usage/tods.rst
    usage/mapping.rst

docs/source/usage/mapping.rst

@@ -1,36 +1,8 @@
-#######
+=======
 Mapping
-#######
+=======
 
-We can map a ``TOD`` (or several of them) with a ``Mapper``. The simplest possible mapper just bins the data:
+.. toctree::
+   :maxdepth: 2
 
-.. code-block:: python
-
-    from maria.map.mappers import BinMapper
-
-    mapper = BinMapper(center=(150, 10),
-                       frame="ra_dec",
-                       width=1e0,
-                       height=1e0,
-                       resolution=5e-3,
-                       tod_preprocessing={
-                            "window": {"tukey": {"alpha": 0.1}},
-                            "remove_modes": {"n": 1},
-                            "filter": {"f_lower": 0.01},
-                            "despline": {"knot_spacing": 5},
-                                        },
-                        map_postprocessing={
-                            "gaussian_filter": {"sigma": 1},
-                            "median_filter": {"size": 1},
-                                        },
-                       units="K_RJ",
-                       tods=[tod],
-                       )
-
-    output_map = mapper.run()
-
-where we define the preprocessing to be done on the ``TOD``. We can see the output with
-
-.. code-block:: python
-
-    output_map.plot()
+   mapping/index.rst

docs/source/usage/mapping/index.rst

@@ -0,0 +1,36 @@
+#######
+Mapping
+#######
+
+We can map a ``TOD`` (or several of them) with a ``Mapper``. The simplest possible mapper just bins the data:
+
+.. code-block:: python
+
+    from maria.map.mappers import BinMapper
+
+    mapper = BinMapper(center=(150, 10),
+                       frame="ra_dec",
+                       width=1e0,
+                       height=1e0,
+                       resolution=5e-3,
+                       tod_preprocessing={
+                            "window": {"tukey": {"alpha": 0.1}},
+                            "remove_modes": {"n": 1},
+                            "filter": {"f_lower": 0.01},
+                            "despline": {"knot_spacing": 5},
+                                        },
+                        map_postprocessing={
+                            "gaussian_filter": {"sigma": 1},
+                            "median_filter": {"size": 1},
+                                        },
+                       units="K_RJ",
+                       tods=[tod],
+                       )
+
+    output_map = mapper.run()
+
+where we define the preprocessing to be done on the ``TOD``. We can see the output with
+
+.. code-block:: python
+
+    output_map.plot()

docs/source/usage/simulations.rst

@@ -6,7 +6,6 @@ Simulations
    :maxdepth: 2
 
    simulations/index.rst
-   simulations/tods.rst
    simulations/atmosphere.rst
    simulations/cmb.rst
    simulations/maps.rst

docs/source/usage/tods.rst

@@ -0,0 +1,8 @@
+====
+TODs
+====
+
+.. toctree::
+   :maxdepth: 2
+
+   tods/index.rst

docs/source/usage/tods/index.rst

@@ -0,0 +1,47 @@
+####
+TODs
+####
+
+The result of a simulation is a ``TOD``, which encapsulates the generated time-ordered data (and all metadata, like the pointing coordinates).
+The overall signal in the
+
+
+.. code-block:: python
+
+    tod.plot()
+
+It has a few useful features. We can see what our data looks like with
+
+.. code-block:: python
+
+    tod.plot()
+
+==========
+Components
+==========
+
+The total signal in each detector can be accessed as
+
+.. code-block:: python
+
+    tod.signal # returns an array
+
+which is the sum of all of the simulated fields (e.g. noise, atmosphere, CMB) separately, contributing to the incident power. For convenience, we can also access the individual fields as
+
+.. code-block:: python
+
+    tod.get_field("atmosphere") # returns an array
+
+We can see all the available fields with ``tod.fields``.
+
+
+=====
+Units
+=====
+
+TODs are by default in units of picowatts, but we can convert to any unit that is a combination of an SI prefix and a base unit (one of `K_RJ`, `K_CMB`, or `W`).
+
+.. code-block:: python
+
+    tod_in_rj_units = tod.to(units="mK_RJ")
+    tod_in_cmb_units = tod.to(units="uK_CMB")

maria/base/__init__.py

@@ -145,9 +145,12 @@ def run(self):
 
         for k, data in self.data.items():
             # scaling floats doesn't make them more accurate, unless they're huge or tiny
-            offset = data.mean(axis=-1)[..., None]
+            offset = data.mean(axis=-1)
             # scale = data.std(axis=-1)[..., None]
-            tod_data[k] = {"data": (data - offset).astype(self.dtype), "offset": offset}
+            tod_data[k] = {
+                "data": (data - offset[..., None]).astype(self.dtype),
+                "offset": offset,
+            }
 
         tod = TOD(
             data=tod_data,

maria/coords/coordinates.py

@@ -217,6 +217,12 @@ def compute_transforms(self):
             setattr(self, frames[frame]["phi"], frame_phi)
             setattr(self, frames[frame]["theta"], frame_theta)
 
+    @property
+    def timestep(self):
+        if len(self.time):
+            return np.mean(np.gradient(self.time))
+        return None
+
     def downsample(self, timestep: float = None, factor: int = None):
         if timestep is None and factor is None:
             raise ValueError("You must supply either 'timestep' or 'factor'.")

maria/functions.py

@@ -1,7 +1,7 @@
 import numpy as np
 import scipy as sp
 
-from maria.constants import c, h, k_B
+from maria.units.constants import c, h, k_B
 
 
 def rayleigh_jeans_spectrum(nu: float, T: float):

maria/instrument/band/__init__.py

@@ -1,19 +1,25 @@
 import glob
 import os
+import re
 from collections.abc import Mapping
 from typing import Sequence, Union
 
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
 
-from ...constants import T_CMB, c, k_B
 from ...functions import planck_spectrum
 from ...io import flatten_config, read_yaml
 from ...spectrum import AtmosphericSpectrum
+from ...units import BASE_TOD_UNITS, prefixes
+from ...units.constants import T_CMB, c, k_B
 
 here, this_filename = os.path.split(__file__)
 
+prefix_phrase = "|".join(prefixes.index)
+units_phrase = "|".join(BASE_TOD_UNITS)
+units_pattern = re.compile(rf"(?P<prefix>({prefix_phrase}))?(?P<base>.+)")
+
 FIELD_FORMATS = pd.read_csv(f"{here}/format.csv", index_col=0)
 
 BAND_CONFIGS = {}
@@ -66,18 +72,26 @@ def parse_bands(bands):
     return band_list
 
 
-def parse_calibration_signature(s):
-    valid_units = ["pW", "K_RJ", "K_CMB"]
-    for sep in ["/", "->"]:
+def parse_tod_calibration_signature(s):
+    res = {}
+    for sep in ["->"]:
         if s.count(sep) == 1:
             if sep is not None:
                 items = [u.strip() for u in s.split(sep)]
                 if len(items) == 2:
-                    for u in items:
-                        if u not in valid_units:
-                            raise ValueError(f"Invalid units '{u}'.")
-                    return items
-
+                    for io, u in zip(["in", "out"], items):
+                        match = units_pattern.search(u)
+                        if match is None:
+                            raise ValueError(
+                                f"Invalid units '{u}'. Valid units are a combination of an SI prefix "
+                                f"(one of {prefixes.index}) and a base units (one of {BASE_TOD_UNITS})."
+                            )
+                        prefix = match["prefix"]
+                        res[f"{io}_factor"] = (
+                            prefixes.loc[prefix].factor if prefix else 1e0
+                        )
+                        res[f"{io}_units"] = match["base"]
+        return res
     raise ValueError("Calibration must have signature 'units1 -> units2'.")
 
 
@@ -189,10 +203,10 @@ def set_sensitivity(
         )
 
         if kind.lower() == "rayleigh-jeans":
-            self.NEP = self.dP_dTRJ * value / transmission
+            self.NEP = 1e12 * self.dP_dTRJ * value / transmission
 
         elif kind.lower() == "cmb":
-            self.NEP = self.dP_dTCMB * value / transmission
+            self.NEP = 1e12 * self.dP_dTCMB * value / transmission
 
         self._sensitivity = value
         self._sensitivity_kind = kind
@@ -263,7 +277,7 @@ def passband(self, nu):
     @property
     def dP_dTRJ(self) -> float:
         """
-        In picowatts per kelvin Rayleigh-Jeans, assuming perfect transmission.
+        In watts per kelvin Rayleigh-Jeans, assuming perfect transmission.
         """
 
         nu = np.linspace(self.nu_min, self.nu_max, 256)
@@ -272,15 +286,15 @@ def dP_dTRJ(self) -> float:
         # dP_dTRJ = np.trapezoid(dI_dTRJ * self.passband(nu), 1e9 * nu)
         dP_dTRJ = k_B * np.trapezoid(self.passband(nu), 1e9 * nu)
 
-        return 1e12 * self.efficiency * dP_dTRJ
+        return self.efficiency * dP_dTRJ
 
     @property
     def dP_dTCMB(self) -> float:
         """
-        In picowatts per kelvin CMB, assuming perfect transmission.
+        In watts per kelvin CMB, assuming perfect transmission.
         """
 
-        eps = 1e-2
+        eps = 1e-3
         delta_T = np.array([-eps / 2, eps / 2])
 
         nu = np.linspace(self.nu_min, self.nu_max, 256)
@@ -291,38 +305,37 @@ def dP_dTCMB(self) -> float:
         )
 
         return (
-            1e12
-            * self.efficiency
+            self.efficiency
             * k_B
             * np.diff(np.trapezoid(TRJ * self.passband(nu), 1e9 * nu))[0]
             / eps
         )
 
     def cal(self, signature: str) -> float:
         """
-        We compute this as
-
-
+        Remember that:
         d(out units) / d(in units) = (d(out units) / d(pW)) / (d(in units) / d(pW))
         """
 
-        in_units, out_units = parse_calibration_signature(signature)
+        res = parse_tod_calibration_signature(signature)
 
-        if in_units == "K_RJ":
-            d_in_d_pW = 1 / self.dP_dTRJ
-        elif in_units == "K_CMB":
-            d_in_d_pW = 1 / self.dP_dTCMB
+        if res["in_units"] == "K_RJ":
+            d_in_d_W = 1 / self.dP_dTRJ
+        elif res["in_units"] == "K_CMB":
+            d_in_d_W = 1 / self.dP_dTCMB
         else:
-            d_in_d_pW = 1
+            d_in_d_W = 1
 
-        if out_units == "K_RJ":
-            d_out_d_pW = 1 / self.dP_dTRJ
-        elif out_units == "K_CMB":
-            d_out_d_pW = 1 / self.dP_dTCMB
+        if res["out_units"] == "K_RJ":
+            d_out_d_W = 1 / self.dP_dTRJ
+        elif res["out_units"] == "K_CMB":
+            d_out_d_W = 1 / self.dP_dTCMB
         else:
-            d_out_d_pW = 1
+            d_out_d_W = 1
+
+        overall_factor = res["in_factor"] / res["out_factor"]
 
-        return d_out_d_pW / d_in_d_pW
+        return overall_factor * d_out_d_W / d_in_d_W
 
     @property
     def wavelength(self):

maria/sim/atmosphere.py

@@ -6,7 +6,7 @@
 from tqdm import tqdm
 
 from ..atmosphere.turbulent_layer import TurbulentLayer
-from ..constants import k_B
+from ..units.constants import k_B
 
 here, this_filename = os.path.split(__file__)
 

maria/sim/cmb.py

@@ -4,8 +4,8 @@
 import scipy as sp
 from tqdm import tqdm
 
-from ..constants import T_CMB
 from ..functions import planck_spectrum
+from ..units.constants import T_CMB
 
 
 class CMBMixin:

maria/sim/map.py

@@ -6,8 +6,8 @@
 import scipy as sp
 from tqdm import tqdm
 
-from ..constants import k_B
 from ..instrument import beam
+from ..units.constants import k_B
 
 here, this_filename = os.path.split(__file__)
 

maria/sim/simulation.py

@@ -111,10 +111,9 @@ def __init__(
             )
 
             # give it the simulation, so that it knows about pointing, site, etc.
+            # kind of cursed
             self.atmosphere.initialize(self)
 
-            # self._initialize_2d_atmosphere(**atmosphere_kwargs)
-
             duration = ttime.monotonic() - ref_time
             logger.info(f"Initialized atmosphere in {int(1e3 * duration)} ms.")
             ref_time = ttime.monotonic()