refactor: don't use COSMO from config

docs/tutorials.rst

@@ -9,3 +9,4 @@ The following introductory tutorials will help you get started with ``21cmSense`
  tutorials/getting_started
  tutorials/understanding_21cmsense
  tutorials/cli_tutorial
+ tutorials/reproducing_pober_2015

py21cmsense/__init__.py

@@ -1,5 +1,12 @@
 """A package for calculate sensitivies of 21-cm interferometers."""
-__version__ = "2.0.0.beta"
+from pkg_resources import DistributionNotFound, get_distribution
+
+try:
+ __version__ = get_distribution("21cmSense").version
+except DistributionNotFound:
+ __version__ = "unknown"
+finally:
+ del get_distribution, DistributionNotFound
 
 from . import yaml
 from .antpos import hera

py21cmsense/antpos.py

@@ -19,6 +19,7 @@ def hera(
  separation: tp.Length = 14 * un.m,
  split_core: bool = False,
  outriggers: bool = False,
+ row_separation: tp.Length | None = None,
 ) -> tp.Meters:
  """
  Produce a simple regular hexagonal array.
@@ -45,18 +46,23 @@ def hera(
  """
  sep = separation.to_value("m")
 
+ if row_separation is None:
+ row_sep = sep * np.sqrt(3) / 2
+ else:
+ row_sep = row_separation.to_value("m")
+
  # construct the main hexagon
  positions = []
  for row in range(hex_num - 1, -hex_num + split_core, -1):
  # adding split_core deletes a row if it's true
  for col in range(2 * hex_num - abs(row) - 1):
  x_pos = sep * ((2 - (2 * hex_num - abs(row))) / 2 + col)
- y_pos = row * sep * np.sqrt(3) / 2
+ y_pos = row * row_sep
  positions.append([x_pos, y_pos, 0])
 
  # basis vectors (normalized to sep)
- up_right = sep * np.asarray([0.5, np.sqrt(3) / 2, 0])
- up_left = sep * np.asarray([-0.5, np.sqrt(3) / 2, 0])
+ up_right = sep * np.asarray([0.5, row_sep / sep, 0])
+ up_left = sep * np.asarray([-0.5, row_sep / sep, 0])
 
  # split the core if desired
  if split_core:
@@ -90,7 +96,7 @@ def hera(
  * sep
  * (hex_num - 1)
  )
- y_pos = row * sep * (hex_num - 1) * np.sqrt(3) / 2
+ y_pos = row * (hex_num - 1) * row_sep
  theta = np.arctan2(y_pos, x_pos)
  if np.sqrt(x_pos**2 + y_pos**2) > sep * (hex_num + 1):
  if 0 < theta <= 2 * np.pi / 3 + 0.01:

py21cmsense/config.py

@@ -1,5 +1,2 @@
 """Some global configuration options for 21cmSense."""
-from astropy.cosmology import Planck15 as _Planck15
-
 PROGRESS = True # whether to display progress bars for some calculations.
-COSMO = _Planck15

py21cmsense/conversions.py

@@ -51,7 +51,9 @@ def z2f(z: Union[float, np.array]) -> un.Quantity[un.GHz]:
 
 
 def dL_dth(
- z: Union[float, np.array], cosmo: FLRW = Planck15
+ z: Union[float, np.array],
+ cosmo: FLRW = Planck15,
+ approximate=False,
 ) -> un.Quantity[un.Mpc / un.rad / littleh]:
  """
  Return the factor to convert radians to transverse distance at redshift z.
@@ -70,11 +72,20 @@ def dL_dth(
  -----
  From Furlanetto et al. (2006)
  """
- return cosmo.h * cosmo.comoving_transverse_distance(z) / un.rad / littleh
+ if approximate:
+ return (
+ (1.9 * (1.0 / un.arcmin) * ((1 + z) / 10.0) ** 0.2).to(1 / un.rad)
+ * un.Mpc
+ / littleh
+ )
+ else:
+ return cosmo.h * cosmo.comoving_transverse_distance(z) / un.rad / littleh
 
 
 def dL_df(
- z: Union[float, np.array], cosmo: FLRW = Planck15
+ z: Union[float, np.array],
+ cosmo: FLRW = Planck15,
+ approximate=False,
 ) -> un.Quantity[un.Mpc / un.MHz / littleh]:
  """
  Get the factor to convert bandwidth to line-of-sight distance in Mpc/h.
@@ -84,13 +95,25 @@ def dL_df(
  z : float
  The redshift
  """
- return (cosmo.h * cnst.c * (1 + z) / (z2f(z) * cosmo.H(z) * littleh)).to(
- "Mpc/(MHz*littleh)"
- )
+ if approximate:
+ return (
+ (1.7 / 0.1)
+ * ((1 + z) / 10.0) ** 0.5
+ * (cosmo.Om0 / 0.15) ** -0.5
+ * un.Mpc
+ / littleh
+ / un.MHz
+ )
+ else:
+ return (cosmo.h * cnst.c * (1 + z) / (z2f(z) * cosmo.H(z) * littleh)).to(
+ "Mpc/(MHz*littleh)"
+ )
 
 
 def dk_du(
- z: Union[float, np.array], cosmo: FLRW = Planck15
+ z: Union[float, np.array],
+ cosmo: FLRW = Planck15,
+ approximate=False,
 ) -> un.Quantity[littleh / un.Mpc]:
  """
  Get factor converting bl length in wavelengths to h/Mpc.
@@ -105,11 +128,13 @@ def dk_du(
  Valid for u >> 1
  """
  # from du = 1/dth, which derives from du = d(sin(th)) using the small-angle approx
- return 2 * np.pi / dL_dth(z, cosmo) / un.rad
+ return 2 * np.pi / dL_dth(z, cosmo, approximate=approximate) / un.rad
 
 
 def dk_deta(
- z: Union[float, np.array], cosmo: FLRW = Planck15
+ z: Union[float, np.array],
+ cosmo: FLRW = Planck15,
+ approximate=False,
 ) -> un.Quantity[un.MHz * littleh / un.Mpc]:
  """
  Get gactor converting inverse frequency to inverse distance.
@@ -119,11 +144,13 @@ def dk_deta(
  z: float
  Redshift
  """
- return 2 * np.pi / dL_df(z, cosmo)
+ return 2 * np.pi / dL_df(z, cosmo, approximate=approximate)
 
 
 def X2Y(
- z: Union[float, np.array], cosmo: FLRW = Planck15
+ z: Union[float, np.array],
+ cosmo: FLRW = Planck15,
+ approximate=False,
 ) -> un.Quantity[un.Mpc**3 / littleh**3 / un.steradian / un.GHz]:
  """
  Obtain the conversion factor between observing co-ordinates and cosmological volume.
@@ -139,4 +166,6 @@ def X2Y(
  -------
  astropy.Quantity: the conversion factor. Units are Mpc^3/h^3 / (sr MHz).
  """
- return dL_dth(z, cosmo) ** 2 * dL_df(z, cosmo)
+ return dL_dth(z, cosmo, approximate=approximate) ** 2 * dL_df(
+ z, cosmo, approximate=approximate
+ )

py21cmsense/observation.py

@@ -5,6 +5,7 @@
 import collections
 import numpy as np
 from astropy import units as un
+from astropy.cosmology import LambdaCDM, Planck15
 from astropy.io.misc import yaml
 from attr import validators as vld
 from collections import defaultdict
@@ -79,6 +80,10 @@ class Observation:
  tsky_ref_freq : float or Quantity
  Frequency at which the foreground model is equal to `tsky_amplitude`.
  See `spectral_index`. Default assumed to be in MHz.
+ use_approximate_cosmo : bool
+ Whether to use approximate cosmological conversion factors. Doing so will give
+ the same results as the original 21cmSense code, but non-approximate versions
+ that use astropy are preferred.
  """
 
  observatory: obs.Observatory = attr.ib(validator=vld.instance_of(obs.Observatory))
@@ -124,6 +129,8 @@ class Observation:
  validator=ut.nonnegative,
  )
  tsky_ref_freq: tp.Frequency = attr.ib(default=150 * un.MHz, validator=ut.positive)
+ use_approximate_cosmo: bool = attr.ib(default=False, converter=bool)
+ cosmo: LambdaCDM = attr.ib(default=Planck15)
 
  @classmethod
  def from_yaml(cls, yaml_file):
@@ -283,7 +290,12 @@ def kparallel(self) -> un.Quantity[un.littleh / un.Mpc]:
 
  Order of the values is the same as `fftfreq` (i.e. zero-first)
  """
- return conv.dk_deta(self.redshift, config.COSMO) * self.eta
+ return (
+ conv.dk_deta(
+ self.redshift, self.cosmo, approximate=self.use_approximate_cosmo
+ )
+ * self.eta
+ )
 
  @cached_property
  def total_integration_time(self) -> un.Quantity[un.s]:

py21cmsense/sensitivity.py

@@ -17,6 +17,7 @@
 import numpy as np
 import tqdm
 from astropy import units as un
+from astropy.cosmology import LambdaCDM
 from astropy.cosmology.units import littleh, with_H0
 from astropy.io.misc import yaml
 from attr import validators as vld
@@ -35,16 +36,6 @@
 from . import types as tp
 from .theory import _ALL_THEORY_POWER_SPECTRA, EOS2021, TheoryModel
 
-
-def _kconverter(val, allow_unitless=False):
- if hasattr(val, "unit"):
- return val.to(littleh / un.Mpc, with_H0(config.COSMO.H0))
- if not allow_unitless:
- raise ValueError("no units supplied!")
- # Assume it has the 1/Mpc units (default from 21cmFAST)
- return (val / un.Mpc).to(littleh / un.Mpc, with_H0(config.COSMO.H0))
-
-
 logger = logging.getLogger(__name__)
 
 
@@ -107,6 +98,11 @@ def clone(self, **kwargs):
  """Clone the object with new parameters."""
  return attr.evolve(self, **kwargs)
 
+ @property
+ def cosmo(self) -> LambdaCDM:
+ """The cosmology to use in the sensitivity calculations."""
+ return self.observation.cosmo
+
 
 @attr.s(kw_only=True)
 class PowerSpectrum(Sensitivity):
@@ -138,23 +134,22 @@ class PowerSpectrum(Sensitivity):
  A function that takes a single kperp and an array of kpar, and returns a boolean
  array specifying which of the k's are useable after accounting for systematics.
  that is, it returns False for k's affected by systematics.
+
  """
 
- horizon_buffer: tp.Wavenumber = attr.ib(
- default=0.1 * littleh / un.Mpc,
- validator=(
- tp.vld_unit(littleh / un.Mpc, with_H0(config.COSMO.H0)),
- ut.nonnegative,
- ),
- converter=_kconverter,
- )
+ horizon_buffer: tp.Wavenumber = attr.ib(default=0.1 * littleh / un.Mpc)
  foreground_model: str = attr.ib(
  default="moderate", validator=vld.in_(["moderate", "optimistic"])
  )
  theory_model: TheoryModel = attr.ib()
 
  systematics_mask: Callable | None = attr.ib(None)
 
+ @horizon_buffer.validator
+ def _horizon_buffer_validator(self, att, val):
+ tp.vld_unit(littleh / un.Mpc, with_H0(self.cosmo.H0))(self, att, val)
+ ut.nonnegative(self, att, val)
+
  @classmethod
  def from_yaml(cls, yaml_file) -> Sensitivity:
  """
@@ -202,7 +197,11 @@ def _theory_model_validator(self, att, val):
  def k1d(self) -> tp.Wavenumber:
  """1D array of wavenumbers for which sensitivities will be generated."""
  delta = (
- conv.dk_deta(self.observation.redshift, config.COSMO)
+ conv.dk_deta(
+ self.observation.redshift,
+ self.cosmo,
+ approximate=self.observation.use_approximate_cosmo,
+ )
  / self.observation.bandwidth
  )
  dv = delta.value
@@ -211,7 +210,10 @@ def k1d(self) -> tp.Wavenumber:
  @cached_property
  def X2Y(self) -> un.Quantity[un.Mpc**3 / littleh**3 / un.steradian / un.GHz]:
  """Cosmological scaling factor X^2*Y (eg. Parsons 2012)."""
- return conv.X2Y(self.observation.redshift)
+ return conv.X2Y(
+ self.observation.redshift,
+ approximate=self.observation.use_approximate_cosmo,
+ )
 
  @cached_property
  def uv_coverage(self) -> np.ndarray:
@@ -232,7 +234,8 @@ def uv_coverage(self) -> np.ndarray:
 
  def power_normalisation(self, k: tp.Wavenumber) -> float:
  """Normalisation constant for power spectrum."""
- k = _kconverter(k)
+ assert hasattr(k, "unit")
+ assert k.unit.is_equivalent(littleh / un.Mpc)
 
  return (
  self.X2Y
@@ -254,7 +257,7 @@ def sample_noise(self, k_par: tp.Wavenumber, k_perp: tp.Wavenumber) -> tp.Delta:
  """Sample variance contribution at a particular k mode."""
  k = np.sqrt(k_par**2 + k_perp**2).to_value(
  littleh / un.Mpc if self.theory_model.use_littleh else un.Mpc**-1,
- with_H0(config.COSMO.H0),
+ with_H0(self.cosmo.H0),
  )
  return self.theory_model.delta_squared(self.observation.redshift, k)
 
@@ -282,7 +285,11 @@ def _nsamples_2d(
  continue
 
  umag = np.sqrt(u**2 + v**2)
- k_perp = umag * conv.dk_du(self.observation.redshift, config.COSMO)
+ k_perp = umag * conv.dk_du(
+ self.observation.redshift,
+ self.cosmo,
+ approximate=self.observation.use_approximate_cosmo,
+ )
 
  hor = self.horizon_limit(umag)
 
@@ -437,7 +444,11 @@ def horizon_limit(self, umag: float) -> tp.Wavenumber:
  Horizon limit, in h/Mpc.
  """
  horizon = (
- conv.dk_deta(self.observation.redshift, config.COSMO)
+ conv.dk_deta(
+ self.observation.redshift,
+ self.cosmo,
+ approximate=self.observation.use_approximate_cosmo,
+ )
  * umag
  / self.observation.frequency
  )
@@ -506,7 +517,7 @@ def delta_squared(self) -> tp.Delta:
  """The fiducial 21cm power spectrum evaluated at :attr:`k1d`."""
  k = self.k1d.to_value(
  littleh / un.Mpc if self.theory_model.use_littleh else un.Mpc**-1,
- with_H0(config.COSMO.H0),
+ with_H0(self.cosmo.H0),
  )
  return self.theory_model.delta_squared(self.observation.redshift, k)
 
@@ -600,7 +611,7 @@ def write(
  fl[k] = v
  fl[k.replace("noise", "snr")] = self.delta_squared / v
 
- fl["k"] = self.k1d.to("1/Mpc", with_H0(config.COSMO.H0)).value
+ fl["k"] = self.k1d.to("1/Mpc", with_H0(self.cosmo.H0)).value
  fl["delta_squared"] = self.delta_squared
 
  fl.attrs["total_snr"] = self.calculate_significance()
@@ -622,8 +633,8 @@ def plot_sense_1d(self, sample: bool = True, thermal: bool = True):
  plt.plot(self.k1d, value, label=key)
  plt.xscale("log")
  plt.yscale("log")
- plt.xlabel("k [1/Mpc]")
- plt.ylabel(r"$\Delta^2_N \ [{\rm mK}^2/{\rm Mpc}^3$")
+ plt.xlabel("k [h/Mpc]")
+ plt.ylabel(r"$\Delta^2_N \ [{\rm mK}^2$")
  plt.legend()
  plt.title(f"z={conv.f2z(self.observation.frequency):.2f}")