Merge pull request #248 from HERA-Team/mutual_coupling_fix

Mutual coupling fix

.flake8

@@ -1,20 +1,28 @@
 [flake8]
 ignore =
-    E203  # No space before colon
-    W503
-    C408  # Ignore using dict() function.
-    D107  # Missing docstring in __init__ (we do it in the class)
-    D401  # "First line should be in imperative mood" -- this doesn't work for properties, see https://github.com/PyCQA/pydocstyle/issues/301
-    A003  # allow method names to be the same as python builtins
-    RST210  # inline strong start-string without end-string. This is OK in the case of **kwargs in parameters.
-    G004  # Logging statement uses f-string.
+    # No space before colon
+    E203,
+    W503,
+    # Ignore using dict() function.
+    C408,
+    # Missing docstring in __init__ (we do it in the class)
+    D107,
+    # "First line should be in imperative mood" -- this doesn't work for properties, see https://github.com/PyCQA/pydocstyle/issues/301
+    D401,
+    # allow method names to be the same as python builtins
+    A003,
+    # inline strong start-string without end-string. This is OK in the case of **kwargs in parameters.
+    RST210,
+    # Logging statement uses f-string.
+    G004,
 max-line-length = 88
 # Should be 18.
 max-complexity = 35
 exclude =
     development/*
 per-file-ignores =
-    hera_sim/tests/*:D,RST,T001,T201   # print statements allowed in tests
+    # print statements allowed in tests
+    hera_sim/tests/*:D,RST,T001,T201
     docs/conf.py:D,A
     */__init__.py:F401
     scripts/*:T001,T201

.pre-commit-config.yaml

@@ -2,7 +2,7 @@ exclude: 'config.yaml|config_examples/.*.yaml|hera_sim/config/H1C.yaml|hera_sim/
 
 repos:
 - repo: https://github.com/pre-commit/pre-commit-hooks
-  rev: v4.3.0
+  rev: v4.4.0
   hooks:
   - id: trailing-whitespace
   - id: check-added-large-files
@@ -18,7 +18,7 @@ repos:
   - id: mixed-line-ending
     args: ['--fix=no']
 -   repo: https://github.com/PyCQA/flake8
-    rev: 5.0.4  # pick a git hash / tag to point to
+    rev: 6.0.0  # pick a git hash / tag to point to
     hooks:
     -   id: flake8
         additional_dependencies:
@@ -33,26 +33,26 @@ repos:
           - flake8-comprehensions
           - flake8-print
 -   repo: https://github.com/psf/black
-    rev: 22.10.0
+    rev: 23.1.0
     hooks:
     - id: black
 -   repo: https://github.com/pre-commit/pygrep-hooks
-    rev: v1.9.0
+    rev: v1.10.0
     hooks:
       - id: rst-backticks
 
 - repo: https://github.com/PyCQA/isort
-  rev: 5.10.1
+  rev: 5.12.0
   hooks:
   - id: isort
 
 - repo: https://github.com/asottile/pyupgrade
-  rev: v3.1.0
+  rev: v3.3.1
   hooks:
   - id: pyupgrade
     args: [--py38-plus]
 
 - repo: https://github.com/asottile/setup-cfg-fmt
-  rev: v2.1.0
+  rev: v2.2.0
   hooks:
   - id: setup-cfg-fmt

CHANGELOG.rst

@@ -5,6 +5,14 @@ Changelog
 dev-version
 ===========
 
+v3.1.1 [2023.02.23]
+===================
+
+Changed
+-------
+- Coupling matrix calculation in :class:`~.sigchain.MutualCoupling` has been updated
+  to correctly calculate the coupling coefficients from the provided E-field beam.
+
 v3.1.0 [2023.01.17]
 ===================
 

hera_sim/__yaml_constructors.py

@@ -41,7 +41,10 @@ def astropy_unit_constructor(loader, node):
     if value is None:
         return None
     elif units is None:
-        warnings.warn("You have not specified the units for this item. Returning None.")
+        warnings.warn(
+            "You have not specified the units for this item. Returning None.",
+            stacklevel=2,
+        )
         return None
     else:
         try:

hera_sim/beams.py

@@ -288,7 +288,6 @@ def __init__(
         ref_freq=1e8,
         polarized=False,
     ):
-
         self.ref_freq = ref_freq
         self.spectral_index = spectral_index
         self.polarized = polarized
@@ -486,7 +485,6 @@ def __init__(
         perturb_zeropoint=None,
         **kwargs
     ):
-
         # Initialize base class
         super().__init__(beam_coeffs=beam_coeffs, **kwargs)
 

hera_sim/cli_utils.py

@@ -137,7 +137,8 @@ def write_calfits(
         if sim_x_orientation is None:
             warnings.warn(
                 "x_orientation not specified in simulation object."
-                "Assuming that the x-direction points north."
+                "Assuming that the x-direction points north.",
+                stacklevel=1,
             )
         else:
             x_orientation = sim_x_orientation

hera_sim/defaults.py

@@ -295,7 +295,7 @@ def _check_config(self):
                 "Please check your configuration, as only the last "
                 "value specified for each parameter will be used."
             )
-            warnings.warn(warning)
+            warnings.warn(warning, stacklevel=1)
 
     def _handler(self, func, *args, **kwargs):
         """Decorator for applying new function parameter defaults."""

hera_sim/foregrounds.py

@@ -209,7 +209,6 @@ def __init__(
         spectral_index_std=0.5,
         reference_freq=0.15,
     ):
-
         super().__init__(
             nsrcs=nsrcs,
             Smin=Smin,

hera_sim/interpolators.py

@@ -181,7 +181,8 @@ def _interpolator(self):
             warnings.warn(
                 "The provided LSTs do not sufficiently cover [0, 2*pi). "
                 "The interpolated sky temperature may have unexpected behavior "
-                "near 0 and 2*pi."
+                "near 0 and 2*pi.",
+                stacklevel=1,
             )
 
         lsts = np.concatenate(

hera_sim/io.py

@@ -69,7 +69,8 @@ def empty_uvdata(
             "deprecated and will be removed in the future. Please "
             "update your code to use the Nfreqs, Ntimes, and "
             "array_layout parameters instead.",
-            DeprecationWarning,
+            category=DeprecationWarning,
+            stacklevel=2,
         )
 
     # for backwards compatability

hera_sim/noise.py

@@ -243,5 +243,9 @@ def white_noise(*args, **kwargs):
 
     Deprecated. Use ``utils.gen_white_noise`` instead.
     """
-    warnings.warn("white_noise is being deprecated. Use utils.gen_white_noise instead.")
+    warnings.warn(
+        "white_noise is being deprecated. Use utils.gen_white_noise instead.",
+        category=DeprecationWarning,
+        stacklevel=2,
+    )
     return utils.gen_white_noise(*args, **kwargs)

hera_sim/rfi.py

@@ -53,7 +53,6 @@ def __init__(
         std: float = 10.0,
         timescale: float = 100.0,
     ):
-
         self.f0 = f0
         self.duty_cycle = duty_cycle
         self.strength = strength
@@ -165,7 +164,10 @@ def __call__(self, lsts, freqs, **kwargs):
         rfi = np.zeros((lsts.size, freqs.size), dtype=complex)
 
         if stations is None:
-            warnings.warn("You did not specify any stations to simulate.")
+            warnings.warn(
+                "You did not specify any stations to simulate.",
+                stacklevel=2,
+            )
             return rfi
         elif isinstance(stations, (str, Path)):
             # assume that it's a path to a npy file
@@ -411,7 +413,8 @@ def __call__(self, lsts, freqs, **kwargs):
             warnings.warn(
                 "The DTV band does not overlap with any of the passed "
                 "frequencies. Please ensure that you are passing the "
-                "correct set of parameters."
+                "correct set of parameters.",
+                stacklevel=2,
             )
 
         # define an iterator, just to keep things neat

hera_sim/sigchain.py

@@ -5,9 +5,11 @@
 """
 from __future__ import annotations
 
+import astropy_healpix as aph
+import copy
 import numpy as np
 import warnings
-from astropy import constants
+from astropy import constants, units
 from pathlib import Path
 from pyuvdata import UVBeam
 from pyuvsim import AnalyticBeam
@@ -231,7 +233,8 @@ def _type_check(arr):
                     if arr.shape[0] == Nfreqs:
                         warnings.warn(
                             "The input array had lengths Nfreqs "
-                            "and is being reshaped as (Ntimes,1)."
+                            "and is being reshaped as (Ntimes,1).",
+                            stacklevel=1,
                         )
                 elif arr.ndim > 1:
                     assert arr.shape[1] in (1, Nfreqs), (
@@ -672,17 +675,18 @@ class MutualCoupling(Crosstalk):
     .. math::
 
         {\bf X}_{jk} \equiv \frac{i\eta_0}{4\lambda} \frac{\Gamma_k}{R_k}
-        \frac{e^{i2\pi\nu\tau_{jk}}}{b_{jk}} {\bf J}_j (\hat{\bf b}_{jk})
-        {\bf J}_k(\hat{\bf b}_{kj})^\dagger,
+        \frac{e^{-i2\pi\nu\tau_{jk}}}{b_{jk}} {\bf J}_j (\hat{\bf b}_{jk})
+        {\bf J}_k(\hat{\bf b}_{kj})^\dagger h_0^2,
 
     where :math:`\Gamma` is the reflection coefficient, :math:`R` is the real
     part of the impedance, :math:`\eta_0` is the impedance of free space,
     :math:`\lambda` is the wavelength of the radiation, :math:`\nu` is the
     frequency of the radiation, :math:`\tau=b/c` is the delay of the baseline,
     :math:`b` is the baseline length, :math:`\hat{\bf b}_{ij}` is a unit
-    vector pointing from antenna :math:`i` to antenna :math:`j`, and
-    :math:`{\bf J}` is the Jones matrix describing the fully-polarized
-    response of the antenna to incident radiation.
+    vector pointing from antenna :math:`i` to antenna :math:`j`, :math:`{\bf J}`
+    is the Jones matrix describing the antenna's peak-normalized far-field
+    radiation pattern, and :math:`h_0` is the amplitude of the antenna's
+    effective height.
 
     The boldfaced variables without any overhead decorations indicate 2x2
     matrices:
@@ -693,10 +697,25 @@ class MutualCoupling(Crosstalk):
             V_{XX} & V_{XY} \\ V_{YX} & V_{YY}
         \end{pmatrix},
         \quad
-        {\bf J} = \begin{pmatrix}
+        {\bf J} = \frac{1}{h_0} \begin{pmatrix}
             h_{X\theta} & h_{X\phi} \\ h_{Y\theta} & h_{Y\phi}
         \end{pmatrix}
 
+    The effective height can be rewritten as
+
+    .. math::
+
+        h_0^2 = \frac{4\lambda^2 R}{\eta_0 \Omega_p}
+
+    where :math:`\Omega_p` is the beam area (i.e. integral of the peak-normalized
+    power beam). Substituting this in to the previous expression for the coupling
+    coefficient and taking antennas to be identical gives
+
+    .. math::
+
+        {\bf X}_{jk} = \frac{i\Gamma}{\Omega_p} \frac{e^{-i2\pi\nu\tau_{jk}}}
+        {b_{jk}/\lambda} {\bf J}(\hat{\bf b}_{jk}) {\bf J}(\hat{\bf b}_{kj})^\dagger.
+
     In order to efficiently simulate the mutual coupling, the antenna and
     polarization axes of the visibilities and coupling matrix are combined
     into a single "antenna-polarization" axis, and the problem is recast as a
@@ -715,10 +734,9 @@ class MutualCoupling(Crosstalk):
         Should be either a :class:`np.ndarray` or an interpolation object that
         gives the reflection coefficient as a function of frequency (in GHz).
         Not required if providing a pre-calculated coupling matrix.
-    resistance
-        The real part of the impedance to use for calculating the coupling
-        matrix. Has the same requirements as the ``reflection`` parameter,
-        and the resistance should be measured in Ohms.
+    omega_p
+        The integral of the peak-normalized power beam as a function of frequency
+        (in GHz). Not required if providing a pre-calculated coupling matrix.
     ant_1_array
         Array of integers specifying the number of the first antenna in each
         visibility. Required for calculating the coupling matrix and the
@@ -769,7 +787,7 @@ def __init__(
         self,
         uvbeam: UVBeam | str | Path | None,
         reflection: np.ndarray | Callable | None = None,
-        resistance: np.ndarray | Callable | None = None,
+        omega_p: np.ndarray | Callable | None = None,
         ant_1_array: np.ndarray | None = None,
         ant_2_array: np.ndarray | None = None,
         pol_array: np.ndarray | None = None,
@@ -783,7 +801,7 @@ def __init__(
         super().__init__(
             uvbeam=uvbeam,
             reflection=reflection,
-            resistance=resistance,
+            omega_p=omega_p,
             ant_1_array=ant_1_array,
             ant_2_array=ant_2_array,
             pol_array=pol_array,
@@ -831,7 +849,7 @@ def __call__(
         (
             uvbeam,
             reflection,
-            resistance,
+            omega_p,
             ant_1_array,
             ant_2_array,
             pol_array,
@@ -891,7 +909,7 @@ def __call__(
                 array_layout=array_layout,
                 uvbeam=uvbeam,
                 reflection=reflection,
-                resistance=resistance,
+                omega_p=omega_p,
                 pixel_interp=pixel_interp,
                 freq_interp=freq_interp,
                 **beam_kwargs,
@@ -927,7 +945,7 @@ def build_coupling_matrix(
         array_layout: dict,
         uvbeam: UVBeam | str,
         reflection: np.ndarray | Callable | None = None,
-        resistance: np.ndarray | Callable | None = None,
+        omega_p: np.ndarray | Callable | None = None,
         pixel_interp: str | None = "az_za_simple",
         freq_interp: str | None = "cubic",
         **beam_kwargs,
@@ -962,11 +980,10 @@ def build_coupling_matrix(
             The reflection coefficient to use for calculating the coupling matrix.
             Should be either a :class:`np.ndarray` or an interpolation object that
             gives the reflection coefficient as a function of frequency (in GHz).
-            Not required if providing a pre-calculated coupling matrix.
-        resistance
-            The real part of the impedance to use for calculating the coupling
-            matrix. Has the same requirements as the `reflection` parameter,
-            and the resistance should be measured in Ohms.
+        omega_p
+            The integral of the peak-normalized power beam as a function of frequency
+            (in GHz). If this is not provided, then it will be calculated from the
+            provided beam model.
         pixel_interp
             The name of the spatial interpolation method used for the beam. Not
             required if using an analytic beam or if providing a pre-computed
@@ -999,17 +1016,56 @@ def build_coupling_matrix(
         if reflection.size != freqs.size:
             raise ValueError("Reflection coefficients have the wrong shape.")
 
-        if resistance is None:
-            resistance = np.ones_like(freqs)
-        elif callable(resistance):
-            resistance = resistance(freqs)
-        if resistance.size != freqs.size:
-            raise ValueError("Resistances have the wrong shape.")
+        if omega_p is None:
+            warnings.warn(
+                "Calculating the power beam integral; this may take a while.",
+                stacklevel=1,
+            )
+            # Since AnalyticBeam doesn't have a method for calculating the
+            # beam integral, we need to do it manually.
+            if isinstance(uvbeam, AnalyticBeam):
+                power_beam = copy.deepcopy(uvbeam)
+                power_beam.efield_to_power()
+                nside = 128
+                npix = aph.nside_to_npix(nside)
+                pix_area = aph.nside_to_pixel_area(nside).to("sr").value
+                pix_inds = np.arange(npix)
+                lon, lat = aph.healpix_to_lonlat(pix_inds, nside)
+                above_horizon = lat.value > 0
+                beam_vals = power_beam.interp(
+                    az_array=lon.to("rad").value,
+                    za_array=np.pi / 2 - lat.to("rad").value,
+                    freq_array=freqs * units.GHz.to("Hz"),
+                )[0][
+                    0, 0, 0
+                ]  # Just take the XX polarization
+                beam_vals[:, ~above_horizon] = 0  # Apply horizon cut
+                omega_p = beam_vals.sum(axis=1).real * pix_area
+            else:
+                power_beam = uvbeam.copy()
+                power_beam.efield_to_power()
+                if power_beam.interpolation_function is None:
+                    power_beam.interpolation_function = pixel_interp
+                if power_beam.freq_interp_kind is None:
+                    power_beam.freq_interp_kind = freq_interp
+                power_beam = power_beam.interp(
+                    freq_array=freqs * units.GHz.to("Hz"), new_object=True
+                )  # Interpolate to the desired frequencies
+                power_beam.to_healpix()
+                power_beam.peak_normalize()
+                omega_p = power_beam.get_beam_area(pol="xx").real
+            del power_beam
+        elif callable(omega_p):
+            omega_p = omega_p(freqs)
+        if omega_p.size != freqs.size:
+            raise ValueError("Beam integral has the wrong shape.")
 
         # Check the beam is OK and make it smaller if it's too big.
         uvbeam = MutualCoupling._handle_beam(uvbeam, **beam_kwargs)
         MutualCoupling._check_beam_is_ok(uvbeam)
         if isinstance(uvbeam, UVBeam):
+            uvbeam = uvbeam.copy()
+            uvbeam.peak_normalize()
             if uvbeam.Naxes2 > 5:
                 # We only need two points on either side of the horizon.
                 za_array = uvbeam.axis2_array
@@ -1039,10 +1095,10 @@ def build_coupling_matrix(
             uvbeam.interp(
                 az_array=unique_angles,
                 za_array=np.ones_like(unique_angles) * np.pi / 2,
-                freq_array=freqs * 1e9,  # Since we default to GHz...
+                freq_array=freqs * units.GHz.to("Hz"),
             )[0]
             .squeeze()
-            .transpose(3, 2, 0, 1)
+            .transpose(3, 2, 1, 0)
         )
         jones_matrices = {
             angle: jones_matrices[i] for i, angle in enumerate(unique_angles)
@@ -1083,14 +1139,10 @@ def build_coupling_matrix(
                 coupling_matrix[0, :, 1::2, 1::2][:, j, i] = coupling[:, 1, 1]
 
         # Now let's tack on the prefactor
-        eta0 = np.sqrt(constants.mu0.value / constants.eps0.value)
-        coupling_matrix *= (
-            1j
-            * eta0
-            * reflection
-            * freqs
-            / (4 * constants.c.to("m/ns").value * resistance)
-        ).reshape(1, -1, 1, 1)
+        wavelengths = constants.c.si.value / (freqs * units.GHz.to("Hz"))
+        coupling_matrix *= (1j * reflection * wavelengths / omega_p).reshape(
+            1, -1, 1, 1
+        )
         return coupling_matrix
 
     @staticmethod

hera_sim/simulate.py

@@ -289,7 +289,8 @@ def add(
         if seed is None and add_vis:
             warnings.warn(
                 "You have not specified how to seed the random state. "
-                "This effect might not be exactly recoverable."
+                "This effect might not be exactly recoverable.",
+                stacklevel=2,
             )
 
         # Simulate the effect by iterating over baselines and polarizations.
@@ -1007,7 +1008,8 @@ def _iteratively_apply(
             warnings.warn(
                 "You have chosen to neither add nor return the effect "
                 "you are trying to simulate, so nothing will be "
-                f"computed. This warning was raised for the model: {model}"
+                f"computed. This warning was raised for the model: {model}",
+                stacklevel=2,
             )
             return
 
@@ -1030,7 +1032,8 @@ def _iteratively_apply(
                 "You are attempting to compute a component but have "
                 "not specified an ``is_multiplicative`` attribute for "
                 "the component. The component will be added under "
-                "the assumption that it is *not* multiplicative."
+                "the assumption that it is *not* multiplicative.",
+                stacklevel=2,
             )
             is_multiplicative = False
 
@@ -1487,14 +1490,14 @@ def _sanity_check(self, model):
         if is_multiplicative and not has_data:
             warnings.warn(
                 "You are trying to compute a multiplicative "
-                "effect, but no visibilities have been "
-                "simulated yet."
+                "effect, but no visibilities have been simulated yet.",
+                stacklevel=1,
             )
         elif not is_multiplicative and contains_multiplicative_effect:
             warnings.warn(
                 "You are adding visibilities to a data array "
-                "*after* multiplicative effects have been "
-                "introduced."
+                "*after* multiplicative effects have been introduced.",
+                stacklevel=1,
             )
 
     def _update_history(self, model, **kwargs):

hera_sim/tests/conftest.py

@@ -6,7 +6,6 @@
 
 @pytest.fixture(autouse=True, scope="session")
 def setup_and_teardown_package():
-
     # Try to download the latest IERS table. If the download succeeds, run a
     # computation that requires the values, so they are cached for all future
     # tests. If it fails, turn off auto downloading for the tests and turn it

hera_sim/tests/test_beams.py

@@ -129,7 +129,7 @@ def run_sim(
         uvdata=uvdata,
         beams=beams,
         sky_model=SkyModel(
-            freq_array=freqs,
+            freq_array=freqs * units.Hz,
             ra=Longitude(ra_dec[:, 0] * units.rad),
             dec=Latitude(ra_dec[:, 1] * units.rad),
             spectral_type="full",
@@ -214,7 +214,6 @@ def get_perturbed_beams(
         return beams
 
     def test_rotations(self, antennas, sources):
-
         # Rotate the beam from 0 to 180 degrees, and check that autocorrelation
         # of antenna 0 has approximately the same value when pixel beams are
         # used, and when pixel beams not used (direct beam calculation).

hera_sim/tests/test_compare_pyuvsim.py

@@ -3,6 +3,7 @@
 
 import copy
 import numpy as np
+from astropy import units
 from astropy.coordinates import Latitude, Longitude
 from astropy.time import Time
 from astropy.units import Quantity
@@ -89,7 +90,7 @@ def get_sky_model(uvdata, nsource):
         dec=Latitude(ra_dec[:, 1], "rad"),
         spectral_type="spectral_index",
         spectral_index=sources[:, 3],
-        stokes=stokes,
+        stokes=stokes * units.Jy,
         reference_frequency=Quantity(reference_frequency, "Hz"),
     )
 

hera_sim/tests/test_io.py

@@ -64,7 +64,7 @@ def test_conj_convention(time_params, freq_params, array_layout):
         **time_params,
         **freq_params,
     )
-    for (ai, aj) in uvd.get_antpairs():
+    for ai, aj in uvd.get_antpairs():
         assert ai <= aj
 
 

hera_sim/tests/test_sigchain.py

@@ -406,9 +406,7 @@ def test_mutual_coupling(use_numba):
         vis_amps[(aj, ai)] = vis_amps[(ai, aj)]
 
     # Set the coupling parameters so that it's simple to check the amplitudes.
-    refl_amp = 1e-2
-    eta0 = np.sqrt(constants.mu0 / constants.eps0).to("ohm").value
-    resistance = eta0 * freqs / (4 * constants.c.si.value)
+    refl_amp = 1
 
     # Actually simulate the coupling.
     mutual_coupling = sigchain.MutualCoupling(
@@ -418,7 +416,7 @@ def test_mutual_coupling(use_numba):
         pol_array=uvdata.polarization_array,
         array_layout=dict(zip(*uvdata.get_ENU_antpos()[::-1])),
         reflection=np.ones(uvdata.Nfreqs) * refl_amp,
-        resistance=resistance,
+        omega_p=constants.c.si.value / uvdata.freq_array[0],
     )
     uvdata.data_array += mutual_coupling(
         freqs=freqs / 1e9,
@@ -542,33 +540,35 @@ def uvbeam(tmp_path):
 
 
 def test_mutual_coupling_with_uvbeam(uvbeam, sample_uvdata, sample_coupling):
+    np.random.seed(0)
     data = np.random.normal(size=sample_uvdata.data_array.shape) + 0j
     sample_uvdata.data_array[...] = data[...]
     sample_uvdata.data_array += sample_coupling(
         freqs=sample_uvdata.freq_array.squeeze() / 1e9,
         visibilities=sample_uvdata.data_array,
+        reflection=np.ones(sample_uvdata.Nfreqs) * 10,
         uvbeam=uvbeam,
     )
     assert not np.any(np.isclose(data, sample_uvdata.data_array))
 
 
-@pytest.mark.parametrize("resistance", [None, "callable"])
+@pytest.mark.parametrize("omega_p", [None, "callable"])
 @pytest.mark.parametrize("reflection", [None, "callable"])
 def test_mutual_coupling_input_types(
-    resistance, reflection, sample_uvdata, sample_coupling
+    omega_p, reflection, sample_uvdata, sample_coupling
 ):
     def func(freqs):
         return np.ones_like(freqs)
 
-    resistance = func if resistance == "callable" else None
+    omega_p = func if omega_p == "callable" else None
     reflection = func if reflection == "callable" else None
     data = np.random.normal(size=sample_uvdata.data_array.shape) + 0j
     sample_uvdata.data_array += data
     sample_uvdata.data_array += sample_coupling(
         freqs=sample_uvdata.freq_array.squeeze() / 1e9,
         visibilities=sample_uvdata.data_array,
-        resistance=resistance,
         reflection=reflection,
+        omega_p=omega_p,
     )
     assert not np.any(np.isclose(data, sample_uvdata.data_array))
 
@@ -584,11 +584,11 @@ def test_mutual_coupling_bad_ants(sample_uvdata, sample_coupling):
         )
 
 
-@pytest.mark.parametrize("isbad", ["reflection", "resistance"])
+@pytest.mark.parametrize("isbad", ["reflection", "omega_p"])
 def test_mutual_coupling_bad_feed_params(isbad, sample_uvdata, sample_coupling):
-    kwargs = {"reflection": None, "resistance": None}
+    kwargs = {"omega_p": None, "reflection": None}
     kwargs[isbad] = np.ones(sample_uvdata.Nfreqs + 1)
-    with pytest.raises(ValueError, match="have the wrong shape"):
+    with pytest.raises(ValueError, match="the wrong shape"):
         _ = sample_coupling(
             freqs=sample_uvdata.freq_array.squeeze() / 1e9,
             visibilities=sample_uvdata.data_array,

hera_sim/tests/test_vis.py

@@ -248,10 +248,11 @@ def make_point_sky(uvdata, ra: np.ndarray, dec: np.ndarray, align=True):
                 np.zeros((len(freqs), len(ra))),
                 np.zeros((len(freqs), len(ra))),
             ]
-        ),
+        )
+        * units.Jy,
         name=np.array(["derp"] * len(ra)),
         spectral_type="full",
-        freq_array=freqs,
+        freq_array=freqs * units.Hz,
     )
 
 
@@ -322,7 +323,7 @@ def create_uniform_sky(freq, nbase=4, scale=1) -> SkyModel:
         * units.Jy
         / units.sr,
         spectral_type="full",
-        freq_array=freq,
+        freq_array=freq * units.Hz,
         name=np.array([str(i) for i in range(npix)]),
     )
 
@@ -516,7 +517,6 @@ def test_comparison(simulator, uvdata2, sky_model, beam_model):
 @pytest.mark.parametrize("order", ["time", "baseline", "ant1", "ant2"])
 @pytest.mark.parametrize("conj", ["ant1<ant2", "ant2<ant1"])
 def test_ordering(uvdata_linear, simulator, order, conj):
-
     uvdata_linear.reorder_blts(order=order, conj_convention=conj)
 
     sky_model = make_point_sky(

hera_sim/utils.py

@@ -471,8 +471,8 @@ def Jy2T(freqs, omega_p):
     Deprecated in v1.0.0. Will be removed in v1.1.0
     """
     warnings.warn(
-        "The function Jy2T has been renamed 'jansky_to_kelvin'. It will be removed in "
-        "v1.1."
+        "This function has been deprecated. Please use `jansky_to_kelvin` instead.",
+        stacklevel=1,
     )
     return jansky_to_kelvin(freqs, omega_p)
 

hera_sim/visibilities/healvis_wrapper.py

@@ -46,11 +46,10 @@ def __init__(self, fov=180, nprocesses=1, sky_ref_chan=0):
             raise ImportError("to use the healvis wrapper, you must install healvis!")
 
         warnings.warn(
-            (
-                "The healvis package is deprecated. Please use pyuvsim instead. "
-                "The healvis wrapper will be removed from hera_sim in version 4",
-            ),
+            "The healvis package is deprecated. Please use pyuvsim instead. "
+            "The healvis wrapper will be removed from hera_sim in version 4",
             category=DeprecationWarning,
+            stacklevel=2,
         )
 
         self.fov = fov
@@ -72,7 +71,8 @@ def validate(self, model_data: ModelData):
             warnings.warn(
                 "Using pyuvsim.AnalyticBeam for healvis is not really supported. "
                 "model_data.beams is being automatically modified to be a single "
-                "healvis.AnalyticBeam of the same type."
+                "healvis.AnalyticBeam of the same type.",
+                stacklevel=1,
             )
             old_args = model_data.beams[0].__dict__
 

hera_sim/visibilities/pyuvsim_wrapper.py

@@ -34,7 +34,8 @@ def simulate(self, data_model: ModelData):
             data_model.sky_model.name = np.array(data_model.sky_model.name)
 
         warnings.warn(
-            "UVSim requires time-ordered data. Ensuring that order in UVData..."
+            "UVSim requires time-ordered data. Ensuring that order in UVData...",
+            stacklevel=1,
         )
         data_model.uvdata.reorder_blts("time")
 
@@ -44,7 +45,8 @@ def simulate(self, data_model: ModelData):
         # at least check whether reordering is necessary once uvdata has that ability.
         if np.any(data_model.uvdata.polarization_array != np.array([-5, -6, -7, -8])):
             warnings.warn(
-                "In UVSim, polarization array must be in AIPS order. Reordering..."
+                "In UVSim, polarization array must be in AIPS order. Reordering...",
+                stacklevel=1,
             )
             data_model.uvdata.reorder_pols("AIPS")
 

hera_sim/visibilities/simulators.py

@@ -74,7 +74,6 @@ def __init__(
         beams: BeamListType | None = None,
         normalize_beams: bool = False,
     ):
-
         self.uvdata = self._process_uvdata(uvdata)
 
         # NOT Nants because we only want ants with data
@@ -92,7 +91,6 @@ def __init__(
         self._validate()
 
     def _process_uvdata(self, uvdata: UVData | str | Path):
-
         if isinstance(uvdata, (str, Path)):
             out = UVData()
             out.read(str(uvdata))