Support passing None for HWP angles to pointing matrix calculation. (#…

…345)

* Support passing None for HWP angles to pointing matrix calculation.

Rather than always passing an array of zeros if there is no HWP,
support passing None value to the underlying C++ code.  Fix an
error in the docstring.  Add a unit test to verify that no HWP
and HWP angles of all zero give matching results.

* Add check for hwp buffer size

src/libtoast/src/toast_tod_pointing.cpp

@@ -92,7 +92,11 @@ void toast::pointing_matrix_healpix(toast::HealpixPixels const & hpix,
         // FIXME:  Switch back to fast version after unit tests improved.
         toast::vatan2(n, by, bx, detang.data());
 
-        if (hwpang != NULL) {
+        if (hwpang == NULL) {
+            for (size_t i = 0; i < n; ++i) {
+                detang[i] *= 2.0;
+            }
+        } else {
             for (size_t i = 0; i < n; ++i) {
                 detang[i] += 2.0 * hwpang[i];
                 detang[i] *= 2.0;

src/toast/_libtoast_tod_pointing.cpp

@@ -9,15 +9,13 @@
 void init_tod_pointing(py::module & m) {
     m.def("pointing_matrix_healpix",
           [](toast::HealpixPixels const & hpix, bool nest, double eps, double cal,
-             std::string const & mode, py::buffer pdata, py::buffer hwpang,
+             std::string const & mode, py::buffer pdata, py::object hwpang,
              py::buffer flags, py::buffer pixels, py::buffer weights) {
               pybuffer_check_1D <double> (pdata);
-              pybuffer_check_1D <double> (hwpang);
               pybuffer_check_1D <uint8_t> (flags);
               pybuffer_check_1D <double> (weights);
               pybuffer_check_1D <int64_t> (pixels);
               py::buffer_info info_pdata = pdata.request();
-              py::buffer_info info_hwpang = hwpang.request();
               py::buffer_info info_flags = flags.request();
               py::buffer_info info_pixels = pixels.request();
               py::buffer_info info_weights = weights.request();
@@ -26,8 +24,7 @@ void init_tod_pointing(py::module & m) {
               if (mode.compare("IQU") == 0) {
                   nw = (size_t)(info_weights.size / 3);
               }
-              if ((info_hwpang.size != n) ||
-                  (info_flags.size != n) ||
+              if ((info_flags.size != n) ||
                   (info_pixels.size != n) || (nw != n)) {
                   auto log = toast::Logger::get();
                   std::ostringstream o;
@@ -36,16 +33,30 @@ void init_tod_pointing(py::module & m) {
                   throw std::runtime_error(o.str().c_str());
               }
               double * rawpdata = reinterpret_cast <double *> (info_pdata.ptr);
-              double * rawhwpang = reinterpret_cast <double *> (info_hwpang.ptr);
               uint8_t * rawflags = reinterpret_cast <uint8_t *> (info_flags.ptr);
               double * rawweights = reinterpret_cast <double *> (info_weights.ptr);
               int64_t * rawpixels = reinterpret_cast <int64_t *> (info_pixels.ptr);
+              double * rawhwpang = NULL;
+              if (!hwpang.is_none()) {
+                  auto hwpbuf = py::cast <py::buffer> (hwpang);
+                  pybuffer_check_1D <double> (hwpbuf);
+                  py::buffer_info info_hwpang = hwpbuf.request();
+                  if (info_hwpang.size != n) {
+                      auto log = toast::Logger::get();
+                      std::ostringstream o;
+                      o << "HWP buffer size is not consistent.";
+                      log.error(o.str().c_str());
+                      throw std::runtime_error(o.str().c_str());
+                  }
+                  rawhwpang = reinterpret_cast <double *> (info_hwpang.ptr);
+              }
               toast::pointing_matrix_healpix(hpix, nest, eps, cal, mode, n,
                                              rawpdata, rawhwpang, rawflags, rawpixels,
                                              rawweights);
               return;
           }, py::arg("hpix"), py::arg("nest"), py::arg("eps"), py::arg("cal"),
-          py::arg("mode"), py::arg("pdata"), py::arg("hwpang"), py::arg("flags"),
+          py::arg("mode"), py::arg("pdata"), py::arg("hwpang").none(true),
+          py::arg("flags"),
           py::arg("pixels"), py::arg(
               "weights"), R"(
         Compute the healpix pixel indices and weights for one detector.

src/toast/tests/ops_pmat.py

@@ -165,6 +165,95 @@ def test_pointing_matrix_healpix(self):
         self.assertFalse(failed)
         return
 
+    def test_pointing_matrix_healpix_hwp(self):
+        nside = 64
+        hpix = HealpixPixels(64)
+        nest = True
+        psivec = np.radians([-180, -135, -90, -45, 0, 45, 90, 135, 180])
+        nsamp = len(psivec)
+        eps = 0.0
+        cal = 1.0
+        mode = "IQU"
+        nnz = 3
+        flags = np.zeros(nsamp, dtype=np.uint8)
+        pix = 49103
+        theta, phi = hp.pix2ang(nside, pix, nest=nest)
+        xaxis, yaxis, zaxis = np.eye(3)
+        thetarot = qa.rotation(yaxis, theta)
+        phirot = qa.rotation(zaxis, phi)
+        pixrot = qa.mult(phirot, thetarot)
+        quats = []
+        for psi in psivec:
+            psirot = qa.rotation(zaxis, psi)
+            quats.append(qa.mult(pixrot, psirot))
+        quats = np.vstack(quats)
+
+        # First with HWP angle == 0.0
+        hwpang = np.zeros(nsamp)
+        pixels_zero = np.zeros(nsamp, dtype=np.int64)
+        weights_zero = np.zeros([nsamp, nnz], dtype=np.float64)
+        pointing_matrix_healpix(
+            hpix,
+            nest,
+            eps,
+            cal,
+            mode,
+            quats.reshape(-1),
+            hwpang,
+            flags,
+            pixels_zero,
+            weights_zero.reshape(-1),
+        )
+
+        # Now passing hwpang == None
+        pixels_none = np.zeros(nsamp, dtype=np.int64)
+        weights_none = np.zeros([nsamp, nnz], dtype=np.float64)
+        pointing_matrix_healpix(
+            hpix,
+            nest,
+            eps,
+            cal,
+            mode,
+            quats.reshape(-1),
+            None,
+            flags,
+            pixels_none,
+            weights_none.reshape(-1),
+        )
+        # print("")
+        # for i in range(nsamp):
+        #     print(
+        #         "HWP zero:  {} {} | {} {} {}".format(
+        #             psivec[i],
+        #             pixels_zero[i],
+        #             weights_zero[i][0],
+        #             weights_zero[i][1],
+        #             weights_zero[i][2],
+        #         )
+        #     )
+        #     print(
+        #         "    none:  {} {} | {} {} {}".format(
+        #             psivec[i],
+        #             pixels_none[i],
+        #             weights_none[i][0],
+        #             weights_none[i][1],
+        #             weights_none[i][2],
+        #         )
+        #     )
+        failed = False
+        if not np.all(np.equal(pixels_zero, pixels_none)):
+            print("HWP pixels do not agree {} != {}".format(pixels_zero, pixels_none))
+            failed = True
+
+        if not np.allclose(weights_zero, weights_none):
+            print(
+                "HWP weights do not agree {} != {}".format(weights_zero, weights_none)
+            )
+            failed = True
+
+        self.assertFalse(failed)
+        return
+
     def test_hpix_simple(self):
         rank = 0
         if self.comm is not None:

src/toast/todmap/pointing.py

@@ -33,7 +33,7 @@ class OpPointingHpix(Operator):
     the total response is:
 
     .. math::
-        d = cal \\left[\\frac{(1+eps)}{2} I + \\frac{(1-eps)}{2} \\left[Q \\cos{4(a+w)} + U \\sin{4(a+w)}\\right]\\right]
+        d = cal \\left[\\frac{(1+eps)}{2} I + \\frac{(1-eps)}{2} \\left[Q \\cos{2a+4w} + U \\sin{2a+4w}\\right]\\right]
 
     Args:
         pixels (str): write pixels to the cache with name <pixels>_<detector>.
@@ -164,9 +164,7 @@ def exec(self, data):
             try:
                 hwpang = tod.local_hwp_angle()
             except:
-                pass
-            if hwpang is None:
-                hwpang = np.zeros(nsamp, dtype=np.float64)
+                hwpang = None
 
             # read the common flags and apply bitmask
 
@@ -227,7 +225,9 @@ def exec(self, data):
                         # Use cached version
                         detp = pdata[bslice, :]
 
-                    hslice = hwpang[bslice].reshape(-1)
+                    hslice = None
+                    if hwpang is not None:
+                        hslice = hwpang[bslice].reshape(-1)
                     fslice = common[bslice].reshape(-1)
 
                     pointing_matrix_healpix(