diff --git a/drift/core/beamtransfer.py b/drift/core/beamtransfer.py
index 2fac5a75..0d744fee 100644
--- a/drift/core/beamtransfer.py
+++ b/drift/core/beamtransfer.py
@@ -538,7 +538,7 @@ def _generate_mfiles(self, regen=False):
         # Calculate the number of baselines to deal with at any one time. Aim
         # to have a maximum of "mem_chunk" GB in memory at any one time
         fbsize = self.telescope.num_pol_sky * nl * 2 * nm * 16.0
-        nodemem = self.mem_chunk * 2 ** 30.0
+        nodemem = self.mem_chunk * 2**30.0
 
         num_fb_per_node = int(nodemem / fbsize)
         num_fb_per_chunk = num_fb_per_node * mpiutil.size
diff --git a/drift/core/manager.py b/drift/core/manager.py
index dc0a7649..55791c05 100644
--- a/drift/core/manager.py
+++ b/drift/core/manager.py
@@ -14,6 +14,7 @@
     focalplane,
     restrictedcylinder,
     exotic_cylinder,
+    external_beam,
 )
 from drift.core import beamtransfer
 
@@ -35,6 +36,7 @@
     "RestrictedExtra": restrictedcylinder.RestrictedExtra,
     "GradientCylinder": exotic_cylinder.GradientCylinder,
     "PertCylinder": exotic_cylinder.CylinderPerturbed,
+    "PolarisedCylinderExternalBeam": external_beam.PolarisedCylinderTelescopeExternalBeam,
 }
 
 
diff --git a/drift/core/psestimation.py b/drift/core/psestimation.py
index 86815607..52c2aad9 100644
--- a/drift/core/psestimation.py
+++ b/drift/core/psestimation.py
@@ -47,7 +47,7 @@ def bandfunc_2d_cart(kpar_s, kpar_e, kperp_s, kperp_e):
     def band(k, mu):
 
         kpar = k * mu
-        kperp = k * (1.0 - mu ** 2) ** 0.5
+        kperp = k * (1.0 - mu**2) ** 0.5
 
         parb = (kpar >= kpar_s) * (kpar <= kpar_e)
         perpb = (kperp >= kperp_s) * (kperp < kperp_e)
@@ -327,7 +327,7 @@ def genbands(self):
                 zip(self.kpar_start, self.kpar_end, self.kperp_start, self.kperp_end)
             )
 
-            self.k_center = (self.kpar_center ** 2 + self.kperp_center ** 2) ** 0.5
+            self.k_center = (self.kpar_center**2 + self.kperp_center**2) ** 0.5
 
             # Make a list of functions of the band window functions
             self.band_func = [bandfunc_2d_cart(*bound) for bound in bounds]
@@ -649,7 +649,7 @@ def q_estimator(self, mi, vec1, vec2=None, noise=False):
                     lyvec.conj()
                     * np.dot(self.clarray[bi][li].astype(np.complex128), lxvec),
                     axis=0,
-                ).astype(
+                ).real.astype(
                     np.float64
                 )  # TT only.
 
@@ -817,11 +817,11 @@ def _work_fisher_bias_m(self, mi):
 
         for ia in range(self.nbands):
             c_a = self.getproj(mi, ia)
-            fisher[ia, ia] = np.sum(c_a * c_a.T * ci ** 2)
+            fisher[ia, ia] = np.sum(c_a * c_a.T * ci**2)
 
             for ib in range(ia):
                 c_b = self.getproj(mi, ib)
-                fisher[ia, ib] = np.sum(c_a * c_b.T * ci ** 2)
+                fisher[ia, ib] = np.sum(c_a * c_b.T * ci**2)
                 fisher[ib, ia] = np.conj(fisher[ia, ib])
 
         self.delproj(mi)
diff --git a/drift/core/psmc.py b/drift/core/psmc.py
index 59e1d1e6..03aebcd7 100644
--- a/drift/core/psmc.py
+++ b/drift/core/psmc.py
@@ -224,7 +224,7 @@ def sim_skyvec(trans, n):
 
     gaussvars = (
         np.random.standard_normal(matshape) + 1.0j * np.random.standard_normal(matshape)
-    ) / 2.0 ** 0.5
+    ) / 2.0**0.5
 
     for i in range(lside):
         gaussvars[i] = np.dot(trans[i], gaussvars[i])
diff --git a/drift/core/telescope.py b/drift/core/telescope.py
index a134f448..ee372a7c 100644
--- a/drift/core/telescope.py
+++ b/drift/core/telescope.py
@@ -118,7 +118,7 @@ def max_lm(baselines, wavelengths, uwidth, vwidth=0.0):
     vmax = (np.abs(baselines[:, 1]) + vwidth) / wavelengths
 
     mmax = np.ceil(2 * np.pi * umax).astype(np.int64)
-    lmax = np.ceil((mmax ** 2 + (2 * np.pi * vmax) ** 2) ** 0.5).astype(np.int64)
+    lmax = np.ceil((mmax**2 + (2 * np.pi * vmax) ** 2) ** 0.5).astype(np.int64)
 
     return lmax, mmax
 
@@ -546,7 +546,7 @@ def _unique_baselines(self):
         bl2 = np.around(bl1[..., 0] + 1.0j * bl1[..., 1], self._bl_tol)
 
         # Construct array of baseline lengths
-        blen = np.sum(bl1 ** 2, axis=-1) ** 0.5
+        blen = np.sum(bl1**2, axis=-1) ** 0.5
 
         # Create mask of included baselines
         mask = np.logical_and(blen >= self.minlength, blen <= self.maxlength)
@@ -791,7 +791,7 @@ def transfer_matrices(self, bl_indices, f_indices, global_lmax=True):
         tshape = bl_indices.shape + (self.num_pol_sky, lside + 1, 2 * lside + 1)
         logger.info(
             "Size: %i elements. Memory %f GB."
-            % (np.prod(tshape), 2 * np.prod(tshape) * 8.0 / 2 ** 30)
+            % (np.prod(tshape), 2 * np.prod(tshape) * 8.0 / 2**30)
         )
         tarray = np.zeros(tshape, dtype=np.complex128)
 
diff --git a/drift/pipeline/timestream.py b/drift/pipeline/timestream.py
index f757748e..68b2f784 100644
--- a/drift/pipeline/timestream.py
+++ b/drift/pipeline/timestream.py
@@ -631,7 +631,7 @@ def _q_estimate(mi):
     fisher, bias = ps.fisher_bias()
 
     # Subtract bias and reshape into new array
-    qtotal = (qtotal - bias).reshape(nstream ** 2, ps.nbands).T
+    qtotal = (qtotal - bias).reshape(nstream**2, ps.nbands).T
 
     powerspectrum = np.dot(la.inv(fisher), qtotal)
     powerspectrum = powerspectrum.T.reshape(nstream, nstream, ps.nbands)
diff --git a/drift/telescope/cylbeam.py b/drift/telescope/cylbeam.py
index 2b55c0b9..9614a012 100644
--- a/drift/telescope/cylbeam.py
+++ b/drift/telescope/cylbeam.py
@@ -108,7 +108,7 @@ def fraunhofer_cylinder(antenna_func, width, res=1.0):
 
     ua = -1.0 * np.linspace(-1.0, 1.0, num, endpoint=False)[::-1]
 
-    ax = antenna_func(2 * ua / (1 + ua ** 2))
+    ax = antenna_func(2 * ua / (1 + ua**2))
 
     axe = np.zeros(res * num)
 
diff --git a/drift/telescope/exotic_cylinder.py b/drift/telescope/exotic_cylinder.py
index 0225685a..e3960195 100644
--- a/drift/telescope/exotic_cylinder.py
+++ b/drift/telescope/exotic_cylinder.py
@@ -51,7 +51,7 @@ def feed_positions_cylinder(self, cylinder_index):
         i = np.arange(nf)
 
         pos[:, 0] = cylinder_index * self.cylinder_spacing
-        pos[:, 1] = a * i + 0.5 * b * i ** 2
+        pos[:, 1] = a * i + 0.5 * b * i**2
 
         return pos
 
diff --git a/drift/telescope/external_beam.py b/drift/telescope/external_beam.py
new file mode 100644
index 00000000..a8d412d2
--- /dev/null
+++ b/drift/telescope/external_beam.py
@@ -0,0 +1,399 @@
+"""Telescope classes and beam transfer routines related to external beam models."""
+
+import abc
+import logging
+
+import numpy as np
+import healpy
+from scipy.interpolate import RectBivariateSpline
+
+from caput import cache
+from caput import config
+
+from cora.util import coord, hputil, units
+
+from drift.core import telescope
+from drift.telescope import cylinder, cylbeam
+
+from draco.core.containers import ContainerBase, GridBeam, HEALPixBeam
+
+
+# Create logger object
+logger = logging.getLogger(__name__)
+
+
+class PolarisedTelescopeExternalBeam(telescope.PolarisedTelescope):
+    """Base class for polarised telescope with external beam model.
+
+    Beam model is read in from a file.
+
+    Attributes
+    ----------
+    primary_beam_filename : str
+        Path to the file containing the primary beam. Can either be a Healpix beam or a
+        GridBeam.
+    freq_interp_beam : bool, optional
+        Interpolate between neighbouring frequencies if we don't have a beam for every
+        frequency channel. Default: False.
+    force_real_beam : bool, optional
+        Ensure the output beam is real, regardless of what the datatype of the beam file
+        is. This can help save memory if the saved beam is complex but you know the
+        imaginary part is zero. Default: False.
+    """
+
+    primary_beam_filename = config.Property(proptype=str)
+    freq_interp_beam = config.Property(proptype=bool, default=False)
+    force_real_beam = config.Property(proptype=bool, default=False)
+
+    def _finalise_config(self):
+        """Get the beam file object."""
+        (
+            self._primary_beam,
+            self._is_grid_beam,
+            self._beam_freq,
+            self._beam_nside,
+            self._beam_pol_map,
+            self._output_dtype,
+        ) = self._load_external_beam(self.primary_beam_filename)
+
+        self._x_grid = None
+        self._y_grid = None
+        self._x_tel = None
+        self._pix_mask = None
+
+        self._pvec_x = None
+        self._pvec_y = None
+
+    def _load_external_beam(self, filename):
+        """Load beam from file, and return container and metadata."""
+        logger.debug("Reading beam model from {}...".format(filename))
+        beam = ContainerBase.from_file(
+            filename, mode="r", distributed=False, ondisk=True
+        )
+
+        is_grid_beam = isinstance(beam, GridBeam)
+        is_healpix_beam = isinstance(beam, HEALPixBeam)
+
+        # cache axes
+        beam_freq = beam.freq[:]
+        beam_nside = None if is_grid_beam else beam.nside
+
+        # TODO must use bytestring here because conversion doesn't work with ondisk=True
+        if is_grid_beam:
+            beam_pol_map = {
+                "X": list(beam.pol[:]).index(b"XX"),
+                "Y": list(beam.pol[:]).index(b"YY"),
+            }
+        else:
+            beam_pol_map = {
+                "X": list(beam.pol[:]).index(b"X"),
+                "Y": list(beam.pol[:]).index(b"Y"),
+            }
+
+        if len(beam.input) > 1:
+            raise ValueError("Per-feed beam model not supported for now.")
+
+        # If a HEALPixBeam, must check types of theta and phi fields
+        if is_healpix_beam:
+            hpb_types = [v[0].type for v in beam.beam.dtype.fields.values()]
+
+            complex_beam = np.all(
+                [np.issubclass_(hpbt, np.complexfloating) for hpbt in hpb_types]
+            )
+        else:
+            complex_beam = np.issubclass_(beam.beam.dtype.type, np.complexfloating)
+
+        output_dtype = (
+            np.complex128 if complex_beam and not self.force_real_beam else np.float64
+        )
+
+        return beam, is_grid_beam, beam_freq, beam_nside, beam_pol_map, output_dtype
+
+    def beam(self, feed, freq_id):
+        """Compute the beam pattern.
+
+        Parameters
+        ----------
+        feed : int
+            Feed index.
+        freq_id : int
+            Frequency ID.
+
+        Returns
+        -------
+        beam : np.ndarray[pixel, pol]
+            Return the vector beam response at each point in the Healpix grid. This
+            array is of type `np.float64` if the input beam pattern is real, or if
+            `force_real_beam` is set; otherwise it is of type `np.complex128`.
+        """
+        if self._is_grid_beam:
+            # Either we haven't set up interpolation coords yet, or the nside has
+            # changed
+            if (self._beam_nside is None) or (self._beam_nside != self._nside):
+                self._beam_nside = self._nside
+                (
+                    self._x_grid,
+                    self._y_grid,
+                    self._x_tel,
+                    self._pix_mask,
+                ) = self._setup_gridbeam_interpolation(self._primary_beam)
+                self._setup_polpattern()
+
+        map_out = self._evaluate_external_beam(
+            feed,
+            freq_id,
+            self._primary_beam,
+            self._beam_pol_map,
+            self._beam_freq,
+            self._is_grid_beam,
+            self._beam_nside,
+            self._output_dtype,
+            self._x_grid,
+            self._y_grid,
+            self._x_tel,
+            self._pix_mask,
+        )
+
+        return map_out
+
+    def _evaluate_external_beam(
+        self,
+        feed,
+        freq_id,
+        primary_beam,
+        beam_pol_map,
+        beam_freq,
+        is_grid_beam,
+        beam_nside,
+        output_dtype,
+        x_grid=None,
+        y_grid=None,
+        x_tel=None,
+        pix_mask=None,
+    ):
+        tel_freq = self.frequencies
+        nside = self._nside
+        npix = healpy.nside2npix(nside)
+
+        # Get beam model polarization index corresponding to pol of requested feed
+        if self.polarisation[feed] == "X":
+            pol_ind = beam_pol_map["X"]
+        elif self.polarisation[feed] == "Y":
+            pol_ind = beam_pol_map["Y"]
+        else:
+            raise ValueError(
+                f"Unexpected polarisation ({self.polarisation[feed]} for feed {feed}!)"
+            )
+
+        # Find nearest frequency
+        freq_sel = _nearest_freq(
+            tel_freq, beam_freq, freq_id, single=(not self.freq_interp_beam)
+        )
+        # Raise an error if we can't find any suitable frequency
+        if len(freq_sel) == 0:
+            raise ValueError(f"No beam model spans frequency {tel_freq[freq_id]}.")
+
+        if is_grid_beam:
+            # Interpolate gridbeam onto Healpix
+            beam_map = self._interpolate_gridbeam(
+                freq_sel,
+                pol_ind,
+                primary_beam,
+                beam_pol_map,
+                x_grid,
+                y_grid,
+                x_tel,
+                pix_mask,
+            )
+
+        else:  # Healpix input beam - just need to change to the required resolution
+            beam_map = primary_beam.beam[freq_sel, pol_ind, 0, :]
+
+            # Check resolution and resample to a better resolution if needed
+            if nside != beam_nside:
+                if nside > beam_nside:
+                    logger.warning(
+                        f"Requested nside={nside} higher than that of "
+                        f"beam {beam_nside}"
+                    )
+
+                logger.debug(
+                    "Resampling external beam from nside {:d} to {:d}".format(
+                        beam_nside,
+                        nside,
+                    )
+                )
+                beam_map_new = np.zeros((len(freq_sel), npix), dtype=beam_map.dtype)
+                beam_map_new["Et"] = healpy.ud_grade(beam_map["Et"], nside)
+                beam_map_new["Ep"] = healpy.ud_grade(beam_map["Ep"], nside)
+                beam_map = beam_map_new
+
+        map_out = np.empty((npix, 2), dtype=output_dtype)
+
+        # Pull out the real part of the beam if we are forcing a conversion. This should
+        # do nothing if the array is already real
+        def _conv_real(x):
+            if self.force_real_beam:
+                x = x.real
+            return x
+
+        if len(freq_sel) == 1:
+            # Exact match
+            map_out[:, 0] = _conv_real(beam_map["Et"][0])
+            map_out[:, 1] = _conv_real(beam_map["Ep"][0])
+        else:
+            # Interpolate between pair of frequencies
+            freq_high = beam_freq[freq_sel[1]]
+            freq_low = beam_freq[freq_sel[0]]
+            freq_int = tel_freq[freq_id]
+
+            alpha = (freq_high - freq_int) / (freq_high - freq_low)
+            beta = (freq_int - freq_low) / (freq_high - freq_low)
+
+            map_out[:, 0] = _conv_real(
+                beam_map["Et"][0] * alpha + beam_map["Et"][1] * beta
+            )
+            map_out[:, 1] = _conv_real(
+                beam_map["Ep"][0] * alpha + beam_map["Ep"][1] * beta
+            )
+
+        return map_out
+
+    def _setup_gridbeam_interpolation(self, primary_beam):
+        # Grid beam coordinates
+        x_grid = primary_beam.phi[:]
+        y_grid = primary_beam.theta[:]
+
+        # Celestial coordinates
+        angpos = hputil.ang_positions(self._nside)
+        x_cel = coord.sph_to_cart(angpos).T
+
+        # Rotate to telescope coords:
+        # first align y with N, then polar axis with NCP
+        x_tel = cylbeam.rotate_ypr(
+            (1.5 * np.pi, np.radians(90.0 - self.latitude), 0), *x_cel
+        )
+
+        # Mask any pixels outside grid
+        x_t, y_t, z_t = x_tel
+        pix_mask = (
+            (z_t > 0)
+            & (np.abs(x_t) < np.abs(x_grid.max()))
+            & (np.abs(y_t) < np.abs(y_grid.max()))
+        )
+
+        return x_grid, y_grid, x_tel, pix_mask
+
+    def _setup_polpattern(self):
+        # Pre-compute polarisation pattern.
+        # Taken from driftscan
+        zenith = np.array([np.pi / 2.0 - np.radians(self.latitude), 0.0])
+        that, phat = coord.thetaphi_plane_cart(zenith)
+        xhat, yhat, zhat = cylbeam.rotate_ypr(
+            [0.0, 0.0, 0.0], phat, -that, coord.sph_to_cart(zenith)
+        )
+
+        angpos = hputil.ang_positions(self._nside)
+        self._pvec_x = cylbeam.polpattern(angpos, xhat)
+        self._pvec_y = cylbeam.polpattern(angpos, yhat)
+
+    def _interpolate_gridbeam(
+        self, f_sel, p_ind, primary_beam, beam_pol_map, x_grid, y_grid, x_tel, pix_mask
+    ):
+        x, y = x_grid, y_grid
+        x_t, y_t, z_t = x_tel
+        mask = pix_mask
+
+        # Interpolation routine requires increasing axes
+        reverse_x = (np.diff(x) < 0).any()
+        if reverse_x:
+            x = x[::-1]
+
+        npix = healpy.nside2npix(self._nside)
+        beam_out = np.zeros(
+            (len(f_sel), npix), dtype=HEALPixBeam._dataset_spec["beam"]["dtype"]
+        )
+        for i, fi in enumerate(f_sel):
+            # For now we just use the magnitude. Assumes input is power beam
+            beam = primary_beam.beam[fi, p_ind, 0]
+            if reverse_x:
+                beam = beam[:, ::-1]
+            beam_spline = RectBivariateSpline(y, x, np.sqrt(np.abs(beam)))
+
+            # Beam amplitude
+            amp = np.zeros(npix, dtype=beam.real.dtype)
+            amp[mask] = beam_spline(y_t[mask], x_t[mask], grid=False)
+
+            # Polarisation projection
+            pvec = self._pvec_x if beam_pol_map["X"] == p_ind else self._pvec_y
+
+            beam_out[i]["Et"] = amp * pvec[:, 0]
+            beam_out[i]["Ep"] = amp * pvec[:, 1]
+
+        return beam_out
+
+
+def _nearest_freq(tel_freq, map_freq, freq_id, single=False):
+    """Find nearest neighbor frequencies. Assumes map frequencies
+    are uniformly spaced.
+
+    Parameters
+    ----------
+    tel_freq : float
+        frequencies from telescope object.
+    map_freq : float
+        frequencies from beam map file.
+    freq_id : int
+        frequency selection.
+    single : bool
+        Only return the single nearest neighbour.
+
+    Returns
+    -------
+    freq_ind : list of neighboring map frequencies matched to tel_freq.
+
+    """
+
+    diff_freq = abs(map_freq - tel_freq[freq_id])
+    if single:
+        return np.array([np.argmin(diff_freq)])
+
+    map_freq_width = abs(map_freq[1] - map_freq[0])
+    match_mask = diff_freq < map_freq_width
+
+    freq_ind = np.nonzero(match_mask)[0]
+
+    return freq_ind
+
+
+class PolarisedCylinderTelescopeExternalBeam(
+    cylinder.CylinderTelescope, PolarisedTelescopeExternalBeam
+):
+    """Class for polarised cylinder telescope with external beam model.
+
+    Repeats some code from SimplePolarisedTelescope.
+    """
+
+    @property
+    def polarisation(self):
+        """
+        Polarisation map.
+
+        Returns
+        -------
+        pol : np.ndarray
+            One-dimensional array with the polarization for each feed ('X' or 'Y').
+        """
+        return np.asarray(
+            ["X" if feed % 2 == 0 else "Y" for feed in self.beamclass], dtype=np.str
+        )
+
+    @property
+    def beamclass(self):
+        """Simple beam mode of dual polarisation feeds."""
+        nsfeed = self._single_feedpositions.shape[0]
+        return np.concatenate((np.zeros(nsfeed), np.ones(nsfeed))).astype(np.int64)
+
+    @property
+    def feedpositions(self):
+        return np.concatenate((self._single_feedpositions, self._single_feedpositions))
diff --git a/drift/telescope/focalplane.py b/drift/telescope/focalplane.py
index 1912a485..134d487c 100644
--- a/drift/telescope/focalplane.py
+++ b/drift/telescope/focalplane.py
@@ -38,7 +38,7 @@ def beam_circular(angpos, zenith, uv_diameter):
 def gaussian_beam(angpos, pointing, fwhm):
 
     sigma = np.radians(fwhm) / (8.0 * np.log(2.0)) ** 0.5
-    x2 = (1.0 - coord.sph_dot(angpos, pointing) ** 2) / (4 * sigma ** 2)
+    x2 = (1.0 - coord.sph_dot(angpos, pointing) ** 2) / (4 * sigma**2)
 
     return np.exp(-x2)
 
@@ -94,16 +94,13 @@ def beam_square(self, feed, freq):
 
         pointing = self.beam_pointings[feed]
         bdist = self._angpos - pointing[np.newaxis, :]
-        bdist = (
-            np.abs(
-                np.where(
-                    (bdist[:, 1] < np.pi)[:, np.newaxis],
-                    bdist,
-                    bdist - np.array([0, 2 * np.pi])[np.newaxis, :],
-                )
+        bdist = np.abs(
+            np.where(
+                (bdist[:, 1] < np.pi)[:, np.newaxis],
+                bdist,
+                bdist - np.array([0, 2 * np.pi])[np.newaxis, :],
             )
-            / np.radians(self.beam_size)
-        )
+        ) / np.radians(self.beam_size)
         # bdist = np.abs(np.where((bdist[:, 1] < np.pi)[:, np.newaxis], bdist, bdist - np.array([0, 2*np.pi])[np.newaxis, :])) / np.radians(self.beam_size)
         beam = np.logical_and(bdist[:, 0] < 0.5, bdist[:, 1] < 0.5).astype(np.float64)
 
diff --git a/drift/telescope/gmrt.py b/drift/telescope/gmrt.py
index d2640583..c1a59b2e 100644
--- a/drift/telescope/gmrt.py
+++ b/drift/telescope/gmrt.py
@@ -113,7 +113,7 @@ def beam(self, feed, freq):
                 [np.pi / 2.0 - np.radians(self.pointing), self.zenith[1]]
             )
 
-            x2 = (1.0 - coord.sph_dot(self._angpos, pointing) ** 2) / (4 * sigma ** 2)
+            x2 = (1.0 - coord.sph_dot(self._angpos, pointing) ** 2) / (4 * sigma**2)
             self._bc_map = np.exp(-x2)
 
             self._bc_freq = freq
diff --git a/drift/telescope/restrictedcylinder.py b/drift/telescope/restrictedcylinder.py
index 0c6132b4..76cb7c22 100644
--- a/drift/telescope/restrictedcylinder.py
+++ b/drift/telescope/restrictedcylinder.py
@@ -8,7 +8,7 @@
 def gaussian_fwhm(x, fwhm):
 
     sigma = fwhm / (8.0 * np.log(2.0)) ** 0.5
-    x2 = x ** 2 / (2 * sigma ** 2)
+    x2 = x**2 / (2 * sigma**2)
 
     return np.exp(-x2)
 
diff --git a/drift/util/plotutil.py b/drift/util/plotutil.py
index f18314be..7f7a048a 100644
--- a/drift/util/plotutil.py
+++ b/drift/util/plotutil.py
@@ -27,7 +27,7 @@ def regrid_polar(polar_img, r_bins, theta_bins, res=1024):
     rpar = ra[:, np.newaxis]
     rperp = ra[np.newaxis, :]
 
-    r = (rpar ** 2 + rperp ** 2) ** 0.5
+    r = (rpar**2 + rperp**2) ** 0.5
     th = np.arccos(rpar / r)
     th[0, 0] = 0.0