test: reproduce Pober+2014 results

.github/workflows/run-docs-code.yaml

@@ -14,7 +14,7 @@ jobs:
  strategy:
  fail-fast: false
  matrix:
- demo: ["getting_started", "understanding_21cmsense"]
+ demo: ["getting_started", "understanding_21cmsense", "reproducing_pober_2014"]
  steps:
  - uses: actions/checkout@master
  with:

py21cmsense/beam.py

@@ -94,6 +94,11 @@ class GaussianBeam(PrimaryBeam):
  validator=(tp.vld_physical_type("length"), ut.positive)
  )
 
+ @property
+ def wavelength(self) -> un.Quantity[un.m]:
+ """The wavelength of the observation."""
+ return (cnst.c / self.frequency).to("m")
+
  @property
  def dish_size_in_lambda(self) -> float:
  """The dish size in units of wavelengths."""

py21cmsense/observation.py

@@ -13,6 +13,7 @@
 from os import path
 
 from . import _utils as ut
+from . import config
 from . import conversions as conv
 from . import observatory as obs
 from . import types as tp
@@ -282,7 +283,7 @@ 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) * self.eta
+ return conv.dk_deta(self.redshift, config.COSMO) * self.eta
 
  @cached_property
  def total_integration_time(self) -> un.Quantity[un.s]:

py21cmsense/sensitivity.py

@@ -201,7 +201,10 @@ def _theory_model_validator(self, att, val):
  @cached_property
  def k1d(self) -> tp.Wavenumber:
  """1D array of wavenumbers for which sensitivities will be generated."""
- delta = conv.dk_deta(self.observation.redshift) / self.observation.bandwidth
+ delta = (
+ conv.dk_deta(self.observation.redshift, config.COSMO)
+ / self.observation.bandwidth
+ )
  dv = delta.value
  return np.arange(dv, dv * self.observation.n_channels, dv) * delta.unit
 
@@ -279,7 +282,7 @@ def _nsamples_2d(
  continue
 
  umag = np.sqrt(u**2 + v**2)
- k_perp = umag * conv.dk_du(self.observation.redshift)
+ k_perp = umag * conv.dk_du(self.observation.redshift, config.COSMO)
 
  hor = self.horizon_limit(umag)
 
@@ -405,9 +408,7 @@ def calculate_sensitivity_2d_grid(
  # Get the kperp bin it's in.
  kperp_indx = np.where(k_perp >= kperp_edges)[0][-1]
 
- k = np.sqrt(self.observation.kparallel**2 + k_perp**2)
-
- kpar_indx = np.digitize(k, kpar_edges) - 1
+ kpar_indx = np.digitize(self.observation.kparallel, kpar_edges) - 1
  good_ks = kpar_indx >= 0
  good_ks &= kpar_indx < len(kpar_edges) - 1
 
@@ -436,7 +437,9 @@ def horizon_limit(self, umag: float) -> tp.Wavenumber:
  Horizon limit, in h/Mpc.
  """
  horizon = (
- conv.dk_deta(self.observation.redshift) * umag / self.observation.frequency
+ conv.dk_deta(self.observation.redshift, config.COSMO)
+ * umag
+ / self.observation.frequency
  )
  # calculate horizon limit for baseline of length umag
  if self.foreground_model in ["moderate", "pessimistic"]: