FAQ¶
How do I use my own theoretical model?¶
The theoretical model is required to compute sample (or cosmic) variance.
By default, 21cmSense
uses a state-of-the-art model run with 21cmFAST
from https://arxiv.org/abs/2110.13919. This model provides power spectra at multiple
redshifts and wavenumbers, and should be sufficient for most sensitivity calculations.
If you would like to provide your own theoretical model, you should provide your own
implementation of the py21cmsense.theory.TheoryModel
class. This class
has two required definitions: a class variable use_littleh
, which determines the
units of k
that are passed to it, and the method delta_squared
. A simple
power-law theoretical model may be implemented like this:
from py21cmsense.theory import TheoryModel
from astropy import units as un
class PowerLaw(TheoryModel):
use_littleh = False
def delta_squared(self, z: float, k: np.ndarray) -> un.Quantity[un.mK**2]:
return ((1 + z)**2 * (k / 100.0)**-2) << un.mK**2
Note that the delta_squared
method must exactly match this signature, as it is
called from within the py21cmsense.sensitivity.PowerSpectrum
class.
The default py21cmsense.theory.EOS2021
model loads a 2D array of values
in redshift and k and interpolates them. You could easily set up something similar to
this.
How do I use my custom theory function from the CLI?¶
If you have defined a custom theory model as in the previous question, you can use it
from the command line by adding a couple of lines to your sensitivity YAML
configuration. Imagine you wrote a package called my_theory
, within which there
was a module called module
containing the PowerLaw
class defined above (in
practice, whatever module PowerLaw
resides in just has to be importable on the
PYTHONPATH
of your environment). You would then add the following lines to your
sensitivity YAML configuration:
plugins:
my_theory.module
theory_model:
PowerLaw
How do I change cosmology?¶
You can change the cosmology in the Obervation
class:
>>> from py21cmsense import Observation
>>> from astropy.cosmology import WMAP9
>>> obs = Observation(cosmo=WMAP9, ...)
In practice, the cosmology makes little difference in the calculcations.
How do I define a new PrimaryBeam
?¶
First, you should consider whether you really need one. 21cmSense
only takes into
account the frequency, beam width and UV resolution – not the precise shape.
The beam width specifies the beam-crossing time, which determines how long one can
integrate coherently. The UV resolution is assumed to be a single number, by default
based on the inverse of the beam width. Unless these specific quantities are significantly
different for your desired beam, there is no benefit to being more precise than this.
If you do want to specify your own beam, you must sub-class from PrimaryBeam
, you
must define the same methods as defined by the following definition of the
GaussianBeam
:
@attr.s(frozen=True)
class GaussianBeam(PrimaryBeam):
"""
A simple Gaussian Primary beam.
Parameters
----------
frequency
The fiducial frequency at which the beam operates, assumed to be in MHz
unless otherwise defined.
dish_size
The size of the (assumed circular) dish, assumed to be in meters unless
otherwise defined. This generates the beam size.
"""
dish_size: tp.Length = attr.ib(
validator=(tp.vld_physical_type("length"), ut.positive)
)
@property
def dish_size_in_lambda(self) -> float:
"""The dish size in units of wavelengths."""
return (self.dish_size / (cnst.c / self.frequency)).to("").value
@property
def uv_resolution(self) -> un.Quantity[1 / un.radian]:
"""The appropriate resolution of a UV cell given the beam size."""
return self.dish_size_in_lambda
@property
def area(self) -> un.Quantity[un.steradian]:
"""The integral of the beam over angle, in sr."""
return 1.13 * self.fwhm**2
@property
def width(self) -> un.Quantity[un.radian]:
"""The width of the beam (i.e. sigma), in radians."""
return un.rad * 0.45 / self.dish_size_in_lambda
@property
def fwhm(self) -> un.Quantity[un.radians]:
"""The full-width half maximum of the beam."""
return 2.35 * self.width
@property
def sq_area(self) -> un.Quantity[un.steradian]:
"""The integral of the squared beam, in sr."""
return self.area / 2
@property
def first_null(self) -> un.Quantity[un.radians]:
"""The angle of the first null of the beam."""
return un.rad * 1.22 / self.dish_size_in_lambda