Telescope class using externally-defined beam model

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

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,
 }
 
 

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)

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])

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)
 

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)

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)
 

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