diff --git a/docs/analytic_beams.rst b/docs/analytic_beams.rst new file mode 100644 index 000000000..9d25a3950 --- /dev/null +++ b/docs/analytic_beams.rst @@ -0,0 +1,77 @@ +Analytic Beams +============== + +pyuvdata defines several analytic primary beams for radio telescopes. While these +are not realistic models for true antennas (like those represented in +:class:`pyuvdata.UVBeam`), they can be useful in simulation because they are +lightweight and fast to evaluate (as opposed to having to interpolate). + +The analytic beams defined in pyuvdata are based on an abstract base class, +:class:`pyuvdata.analytic_beam.AnalyticBeam`, which ensures a standard interface +and can be used to define other analytic beams in a consistent way. To evaluate +analytic beams in particular directions at particular frequencies, use the +:meth:`pyuvdata.analytic_beam.AnalyticBeam.efield_eval` or +:meth:`pyuvdata.analytic_beam.AnalyticBeam.power_eval` methods as appropriate. + +The ``AnalyticBeam`` base class also provides a yaml constructor that can enable +analytic beams to be instantiated directly from yaml files (see +:ref:`yaml_constructors`, similar constructors are also available for UVBeam +objects) and a plugin infrastructure that can automatically include any imported +subclass even if they are defined in other packages. This infrastructure, along +with the :class:`pyuvdata.BeamInterface` class, can simplify the setup for +simulations. + +.. autoclass:: pyuvdata.analytic_beam.AnalyticBeam + :members: + + +.. autoclass:: pyuvdata.AiryBeam + :members: + +.. autoclass:: pyuvdata.GaussianBeam + :members: + +.. autoclass:: pyuvdata.ShortDipoleBeam + :members: + +.. autoclass:: pyuvdata.UniformBeam + :members: + + +.. _yaml_constructors: + +yaml constructors +----------------- + +Analytic beams can be instantiated directly from yaml files using the +``!AnalyticBeam`` tag. The ``class`` parameter must be specified and it can be +set to one of the pyuvdata provided analytic beams or to any AnalyticBeam +subclass. If the subclass is not defined in pyuvdata, either the subclass must +already be imported or it must be specified with the dot-pathed modules included +(i.e. ``my_module.MyAnalyticBeam``). Some analytic beams have required parameters +(e.g. ``diameter`` for AiryBeams), which must also be provided, see the object +definitions for details. + +Some examples are provided below, note that the node key can be anything, it +does not need to be ``beam``: + +A 16 meter diameter Airy beam:: + + beam: !AnalyticBeam + class: AiryBeam + diameter: 16 + +A classical short dipole beam (the dot-pathed module notation is not required +for pyvudata beams but is shown here as an example):: + + beam: !AnalyticBeam + class: pyuvdata.ShortDipoleBeam + +A gaussian shaped beam with an E-Field beam sigma of 0.26 radians that has +width that scales as a power law with frequency:: + + beam: !AnalyticBeam + class: GaussianBeam + reference_frequency: 120000000. + spectral_index: -1.5 + sigma: 0.26 diff --git a/docs/beam_interface.rst b/docs/beam_interface.rst new file mode 100644 index 000000000..1f3af3b82 --- /dev/null +++ b/docs/beam_interface.rst @@ -0,0 +1,13 @@ +Beam Interface +============== + +The BeamInterface object is designed to provide a unified interface for UVBeam +and AnalyticBeam objects to compute beam response values. It can be constructed +with either a UVBeam or AnalyticBeam and it provides the +:meth:`pyuvdata.BeamInterface.compute_response` method, which can be used to +either evaluate the response of an AnalyticBeam or interpolate a UVBeam to get +the beam response. This interface provides a simplified way for simulators to +call either UVBeams or analytic beams. + +.. autoclass:: pyuvdata.BeamInterface + :members: diff --git a/docs/cst_settings_yaml.rst b/docs/cst_settings_yaml.rst index 44811e79b..0894363d5 100644 --- a/docs/cst_settings_yaml.rst +++ b/docs/cst_settings_yaml.rst @@ -7,8 +7,7 @@ via keywords when the files are read in, but it is better for the metadata to be specified once and carried with the data files. To that end, we developed a yaml settings file specification to carry all the metadata. This format is very human readable and writeable and we encourage using such a file as the best way to -ensure the metadata is preserved. Note that reading a yaml settings file into -UVBeam requires that pyyaml is installed. +ensure the metadata is preserved. Required Fields *************** diff --git a/docs/make_index.py b/docs/make_index.py index d089b7219..0ca44337d 100644 --- a/docs/make_index.py +++ b/docs/make_index.py @@ -49,6 +49,8 @@ def write_index_rst(readme_file=None, write_file=None): " uvbeam\n" " uvflag\n" " telescope\n" + " analytic_beams\n" + " beam_interface\n" " fast_uvh5_meta\n" " fast_calh5_meta\n" " utility_functions\n" diff --git a/docs/make_uvbeam.py b/docs/make_uvbeam.py index 673a44f55..752e605db 100644 --- a/docs/make_uvbeam.py +++ b/docs/make_uvbeam.py @@ -21,7 +21,9 @@ def write_uvbeam_rst(write_file=None): "transforming the data (interpolating/regridding, selecting, converting\n" "types) and can be interacted with directly.\n\n" "Note that there are some tricks that can help with reading in CST beam\n" - "simulation files in `CST Settings Files`_.\n\n" + "simulation files in `CST Settings Files`_.\n" + "UVBeam also provides a yaml constructor that can enable beams to be\n" + "instantiated directly from yaml files (see: `UVbeam yaml constructor`_)\n\n" "Attributes\n----------\n" "The attributes on UVBeam hold all of the metadata and data required to\n" "describe primary beam models. Under the hood, the attributes are implemented\n" @@ -67,6 +69,33 @@ def write_uvbeam_rst(write_file=None): out += "Methods\n-------\n.. autoclass:: pyuvdata.UVBeam\n :members:\n\n" + out += """ +UVBeam yaml constructor +----------------------- + +UVbeams can be instantiated directly from yaml files using the +``!UVBeam`` tag. The ``filename`` parameter must be specified and and +any other parameter that can be passed to the :meth:`pyuvdata.UVBeam.read` +method can also be specified. + +Some examples are provided below, note that the node key can be anything, +it does not need to be ``beam``: + +A simple UVBeam specification:: + + beam: !UVBeam + filename: hera.beamfits + +An UVbeam specification with some with some keywords to pass to ``UVBeam.read``:: + + beam: !UVBeam + filename: mwa_full_EE_test.h5 + pixels_per_deg: 1 + freq_range: [100.e+6, 200.e+6] + +\n\n +""" + with open("cst_settings_yaml.rst", encoding="utf-8") as cst_settings_file: cst_setting_text = cst_settings_file.read() diff --git a/docs/uvcal_tutorial.rst b/docs/uvcal_tutorial.rst index 83b61ae50..96676d9d6 100644 --- a/docs/uvcal_tutorial.rst +++ b/docs/uvcal_tutorial.rst @@ -324,6 +324,7 @@ that should be noted: on the object itself). * Regardless of the value of ``undo``, the convention that is inferred for the calibration solutions is determined as follows: + * If neither ``uvc_pol_convention`` nor ``uvcal.pol_convention`` are specified, a a warning is raised (since the resulting calibrated data is not well-determined), and it is *assumed* that the solutions have the same convention as the ``UVData`` @@ -332,7 +333,9 @@ that should be noted: corrections are applied, and the result is ambiguous. * If both ``uvc_pol_convention`` and ``uvcal.pol_convention`` are specified and are different, an error is raised. + * When **calibrating** in :func:`pyuvdata.utils.uvcalibrate` (i.e. ``undo=False``): + * If ``uvdata.pol_convention`` is specified, an error is raised, because you are trying to calibrate already-calibrated data. * The convention applied to the resulting ``UVData`` object is inferred in the @@ -341,7 +344,9 @@ that should be noted: or ``uvcal.pol_convention``, see above), (iii) if still unspecified, no convention will be used and a warning will be raised. This was always the behaviour in earlier versions of ``pyuvdata`` (pre-v3). + * When **un-calibrating** with :func:`pyuvdata.utils.uvcalibrate` (i.e. ``undo=True``): + * If both ``uvd_pol_convention`` and ``uvdata.pol_convention`` are defined and are different, an error is raised. * If neither are set, a warning is raised, since the resulting un-calibrated values @@ -349,6 +354,7 @@ that should be noted: convention could have been used to calibrate originally than is being used to de-calibrate). However, calibration will continue, assuming that the ``UVData`` object has the same convention as the ``UVCal`` object used to de-calibrate. + * It is not supported to have ``pol_convention`` set on ``UVCal``, but *not* ``gain_scale``. A ``pol_convention`` only makes sense in the context of having a scale for the gains. diff --git a/src/pyuvdata/analytic_beam.py b/src/pyuvdata/analytic_beam.py index 5c8fcaed8..c218708a8 100644 --- a/src/pyuvdata/analytic_beam.py +++ b/src/pyuvdata/analytic_beam.py @@ -27,7 +27,7 @@ @dataclass(kw_only=True) class AnalyticBeam(ABC): """ - Define an analytic beam base class. + Analytic beam base class. Attributes ---------- @@ -37,7 +37,8 @@ class AnalyticBeam(ABC): x_orientation : str Physical orientation of the feed for the x feed. Not meaningful for unpolarized analytic beams, but clarifies the orientation of the dipole - for the ShortDipoleBeam and matches with the meaning on UVBeam objects. + for polarized beams like the ShortDipoleBeam and matches with the meaning + on UVBeam objects. Parameters ---------- @@ -47,7 +48,8 @@ class AnalyticBeam(ABC): x_orientation : str Physical orientation of the feed for the x feed. Not meaningful for unpolarized analytic beams, but clarifies the orientation of the dipole - for the ShortDipoleBeam and matches with the meaning on UVBeam objects. + for for polarized beams like the ShortDipoleBeam and matches with the + meaning on UVBeam objects include_cross_pols : bool Option to include the cross polarized beams (e.g. xy and yx or en and ne) for the power beam. @@ -529,6 +531,114 @@ def _analytic_beam_representer(dumper, beam): yaml.add_multi_representer(AnalyticBeam, _analytic_beam_representer, Dumper=yaml.Dumper) +@dataclass(kw_only=True) +class AiryBeam(AnalyticBeam): + """ + A zenith pointed Airy beam. + + Airy beams are the diffraction pattern of a circular aperture, so represent + an idealized dish. Requires a dish diameter in meters and is inherently + chromatic and unpolarized. + + The unpolarized nature leads to some results that may be surprising to radio + astronomers: if two feeds are specified they will have identical responses + and the cross power beam between the two feeds will be identical to the + spower beam for a single feed. + + Attributes + ---------- + diameter : float + Dish diameter in meters. + feed_array : np.ndarray of str + Feeds to define this beam for, e.g. x & y or n & e (for "north" and "east") + or r & l. + x_orientation : str + Physical orientation of the feed for the x feed. Not meaningful for + AiryBeams, which are unpolarized. + + Parameters + ---------- + diameter : float + Dish diameter in meters. + feed_array : np.ndarray of str + Feeds to define this beam for, e.g. n & e or x & y or r & l. + x_orientation : str + Physical orientation of the feed for the x feed. Not meaningful for + AiryBeams, which are unpolarized. + include_cross_pols : bool + Option to include the cross polarized beams (e.g. xy and yx or en and ne) for + the power beam. + + """ + + diameter: float + feed_array: npt.NDArray[str] | None = field(default=None, repr=False, compare=False) + x_orientation: Literal["east", "north"] = field( + default="east", repr=False, compare=False + ) + + include_cross_pols: InitVar[bool] = True + + basis_vector_type = "az_za" + + def _efield_eval( + self, + *, + az_array: npt.NDArray[float], + za_array: npt.NDArray[float], + freq_array: npt.NDArray[float], + ) -> npt.NDArray[float]: + """Evaluate the efield at the given coordinates.""" + data_array = self._get_empty_data_array(az_array.size, freq_array.size) + + za_grid, f_grid = np.meshgrid(za_array, freq_array) + kvals = (2.0 * np.pi) * f_grid / speed_of_light.to("m/s").value + xvals = (self.diameter / 2.0) * np.sin(za_grid) * kvals + values = np.zeros_like(xvals) + nz = xvals != 0.0 + ze = xvals == 0.0 + values[nz] = 2.0 * j1(xvals[nz]) / xvals[nz] + values[ze] = 1.0 + + # This is different than what is in pyuvsim because it means that we have the + # same response to any polarized source in both feeds. We think this is correct + # for an azimuthally symmetric analytic beam but it is a change. + # On pyuvsim we essentially have each basis vector go into only one feed. + for fn in np.arange(self.Nfeeds): + data_array[0, fn, :, :] = values / np.sqrt(2.0) + data_array[1, fn, :, :] = values / np.sqrt(2.0) + + return data_array + + def _power_eval( + self, + *, + az_array: npt.NDArray[float], + za_array: npt.NDArray[float], + freq_array: npt.NDArray[float], + ) -> npt.NDArray[float]: + """Evaluate the power at the given coordinates.""" + data_array = self._get_empty_data_array( + az_array.size, freq_array.size, beam_type="power" + ) + + za_grid, f_grid = np.meshgrid(za_array, freq_array) + kvals = (2.0 * np.pi) * f_grid / speed_of_light.to("m/s").value + xvals = (self.diameter / 2.0) * np.sin(za_grid) * kvals + values = np.zeros_like(xvals) + nz = xvals != 0.0 + ze = xvals == 0.0 + values[nz] = (2.0 * j1(xvals[nz]) / xvals[nz]) ** 2 + values[ze] = 1.0 + + for fn in np.arange(self.Npols): + # For power beams the first axis is shallow because we don't have to worry + # about polarization. + data_array[0, fn, :, :] = values + + return data_array + + def diameter_to_sigma(diameter: float, freq_array: npt.NDArray[float]) -> float: """ Find the sigma that gives a beam width similar to an Airy disk. @@ -563,7 +673,24 @@ def diameter_to_sigma(diameter: float, freq_array: npt.NDArray[float]) -> float: @dataclass(kw_only=True) class GaussianBeam(AnalyticBeam): """ - Define a Gaussian beam, optionally with frequency dependent size. + A circular, zenith pointed Gaussian beam. + + Requires either a dish diameter in meters or a standard deviation sigma in + radians. Gaussian beams specified by a diameter will have their width + matched to an Airy beam at each simulated frequency, so are inherently + chromatic. For Gaussian beams specified with sigma, the sigma_type defines + whether the width specified by sigma specifies the width of the E-Field beam + (default) or power beam in zenith angle. If only sigma is specified, the + beam is achromatic, optionally both the spectral_index and reference_frequency + parameters can be set to generate a chromatic beam with standard deviation + defined by a power law: + + stddev(f) = sigma * (f/ref_freq)**(spectral_index) + + The unpolarized nature leads to some results that may be + surprising to radio astronomers: if two feeds are specified they will have + identical responses and the cross power beam between the two feeds will be + identical to the power beam for a single feed. Attributes ---------- @@ -589,8 +716,7 @@ class GaussianBeam(AnalyticBeam): or r & l. x_orientation : str Physical orientation of the feed for the x feed. Not meaningful for - unpolarized analytic beams, but clarifies the orientation of the dipole - for the ShortDipoleBeam and matches with the meaning on UVBeam objects. + GaussianBeams, which are unpolarized. Parameters ---------- @@ -615,8 +741,7 @@ class GaussianBeam(AnalyticBeam): Feeds to define this beam for, e.g. n & e or x & y or r & l. x_orientation : str Physical orientation of the feed for the x feed. Not meaningful for - unpolarized analytic beams, but clarifies the orientation of the dipole - for the ShortDipoleBeam and matches with the meaning on UVBeam objects. + GaussianBeams, which are unpolarized. include_cross_pols : bool Option to include the cross polarized beams (e.g. xy and yx or en and ne) for the power beam. @@ -743,111 +868,17 @@ def _power_eval( return data_array -@dataclass(kw_only=True) -class AiryBeam(AnalyticBeam): - """ - Define an Airy beam. - - Attributes - ---------- - diameter : float - Dish diameter in meters. - feed_array : np.ndarray of str - Feeds to define this beam for, e.g. x & y or n & e (for "north" and "east") - or r & l. - x_orientation : str - Physical orientation of the feed for the x feed. Not meaningful for - unpolarized analytic beams, but clarifies the orientation of the dipole - for the ShortDipoleBeam and matches with the meaning on UVBeam objects. - - Parameters - ---------- - diameter : float - Dish diameter in meters. - feed_array : np.ndarray of str - Feeds to define this beam for, e.g. n & e or x & y or r & l. - x_orientation : str - Physical orientation of the feed for the x feed. Not meaningful for - unpolarized analytic beams, but clarifies the orientation of the dipole - for the ShortDipoleBeam and matches with the meaning on UVBeam objects. - include_cross_pols : bool - Option to include the cross polarized beams (e.g. xy and yx or en and ne) for - the power beam. - - """ - - diameter: float - feed_array: npt.NDArray[str] | None = field(default=None, repr=False, compare=False) - x_orientation: Literal["east", "north"] = field( - default="east", repr=False, compare=False - ) - - include_cross_pols: InitVar[bool] = True - - basis_vector_type = "az_za" - - def _efield_eval( - self, - *, - az_array: npt.NDArray[float], - za_array: npt.NDArray[float], - freq_array: npt.NDArray[float], - ) -> npt.NDArray[float]: - """Evaluate the efield at the given coordinates.""" - data_array = self._get_empty_data_array(az_array.size, freq_array.size) - - za_grid, f_grid = np.meshgrid(za_array, freq_array) - kvals = (2.0 * np.pi) * f_grid / speed_of_light.to("m/s").value - xvals = (self.diameter / 2.0) * np.sin(za_grid) * kvals - values = np.zeros_like(xvals) - nz = xvals != 0.0 - ze = xvals == 0.0 - values[nz] = 2.0 * j1(xvals[nz]) / xvals[nz] - values[ze] = 1.0 - - # This is different than what is in pyuvsim because it means that we have the - # same response to any polarized source in both feeds. We think this is correct - # for an azimuthally symmetric analytic beam but it is a change. - # On pyuvsim we essentially have each basis vector go into only one feed. - for fn in np.arange(self.Nfeeds): - data_array[0, fn, :, :] = values / np.sqrt(2.0) - data_array[1, fn, :, :] = values / np.sqrt(2.0) - - return data_array - - def _power_eval( - self, - *, - az_array: npt.NDArray[float], - za_array: npt.NDArray[float], - freq_array: npt.NDArray[float], - ) -> npt.NDArray[float]: - """Evaluate the power at the given coordinates.""" - data_array = self._get_empty_data_array( - az_array.size, freq_array.size, beam_type="power" - ) - - za_grid, f_grid = np.meshgrid(za_array, freq_array) - kvals = (2.0 * np.pi) * f_grid / speed_of_light.to("m/s").value - xvals = (self.diameter / 2.0) * np.sin(za_grid) * kvals - values = np.zeros_like(xvals) - nz = xvals != 0.0 - ze = xvals == 0.0 - values[nz] = (2.0 * j1(xvals[nz]) / xvals[nz]) ** 2 - values[ze] = 1.0 - - for fn in np.arange(self.Npols): - # For power beams the first axis is shallow because we don't have to worry - # about polarization. - data_array[0, fn, :, :] = values - - return data_array - - @dataclass(kw_only=True) class ShortDipoleBeam(AnalyticBeam): """ - Define an analytic short (Hertzian) dipole beam for two crossed dipoles. + A zenith pointed analytic short dipole beam with two crossed feeds. + + A classical short (Hertzian) dipole beam with two crossed feeds aligned east + and north. Short dipole beams are intrinsically polarized but achromatic. + Does not require any parameters, but the orientation of the dipole labelled + as "x" can be specified to align "north" or "east" via the x_orientation + parameter (matching the parameter of the same name on UVBeam and UVData + objects). Attributes ---------- @@ -855,9 +886,9 @@ class ShortDipoleBeam(AnalyticBeam): Feeds to define this beam for, e.g. x & y or n & e (for "north" and "east") or r & l. x_orientation : str - Physical orientation of the feed for the x feed. Not meaningful for - unpolarized analytic beams, but clarifies the orientation of the dipole - for the ShortDipoleBeam and matches with the meaning on UVBeam objects. + The orientation of the dipole labeled 'x'. The default ("east") means + that the x dipole is aligned east-west and that the y dipole is aligned + north-south. Parameters ---------- @@ -932,7 +963,16 @@ def _power_eval( @dataclass(kw_only=True) class UniformBeam(AnalyticBeam): """ - Define a uniform beam. + A uniform beam. + + Uniform beams have identical responses in all directions, so are quite + unphysical but can be useful for testing other aspects of simulators. They + are unpolarized and achromatic and do not take any parameters. + + The unpolarized nature leads to some results that may be surprising to radio + astronomers: if two feeds are specified they will have identical responses + and the cross power beam between the two feeds will be identical to the + power beam for a single feed. Attributes ---------- @@ -941,17 +981,15 @@ class UniformBeam(AnalyticBeam): or r & l. x_orientation : str Physical orientation of the feed for the x feed. Not meaningful for - unpolarized analytic beams, but clarifies the orientation of the dipole - for the ShortDipoleBeam and matches with the meaning on UVBeam objects. + UniformBeams, which are unpolarized. Parameters ---------- feed_array : np.ndarray of str Feeds to define this beam for, e.g. n & e or x & y or r & l. x_orientation : str - The orientation of the dipole labeled 'x'. The default ("east") means - that the x dipole is aligned east-west and that the y dipole is aligned - north-south. + Physical orientation of the feed for the x feed. Not meaningful for + UniformBeams, which are unpolarized. include_cross_pols : bool Option to include the cross polarized beams (e.g. xy and yx or en and ne) for the power beam.