diff --git a/src/libtoast/src/toast_tod_pointing.cpp b/src/libtoast/src/toast_tod_pointing.cpp index 0d7501382..89de35255 100644 --- a/src/libtoast/src/toast_tod_pointing.cpp +++ b/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; diff --git a/src/toast/_libtoast_tod_pointing.cpp b/src/toast/_libtoast_tod_pointing.cpp index 4906be7ca..65267dc45 100644 --- a/src/toast/_libtoast_tod_pointing.cpp +++ b/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 (pdata); - pybuffer_check_1D (hwpang); pybuffer_check_1D (flags); pybuffer_check_1D (weights); pybuffer_check_1D (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 (info_pdata.ptr); - double * rawhwpang = reinterpret_cast (info_hwpang.ptr); uint8_t * rawflags = reinterpret_cast (info_flags.ptr); double * rawweights = reinterpret_cast (info_weights.ptr); int64_t * rawpixels = reinterpret_cast (info_pixels.ptr); + double * rawhwpang = NULL; + if (!hwpang.is_none()) { + auto hwpbuf = py::cast (hwpang); + pybuffer_check_1D (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 (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. diff --git a/src/toast/tests/ops_pmat.py b/src/toast/tests/ops_pmat.py index 1672ed832..00b80bb90 100644 --- a/src/toast/tests/ops_pmat.py +++ b/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: diff --git a/src/toast/todmap/pointing.py b/src/toast/todmap/pointing.py index 08647ed7f..c9c3d0f25 100644 --- a/src/toast/todmap/pointing.py +++ b/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 _. @@ -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(