From 596c1df243e34f3e931c5ccb2944c74179973eaf Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Fri, 13 May 2022 15:28:58 +0200
Subject: [PATCH] add failing test case

---
 test/test_mueller_convolver.py | 53 ++++++++++++++++++++++++++++++++++
 1 file changed, 53 insertions(+)

diff --git a/test/test_mueller_convolver.py b/test/test_mueller_convolver.py
index 5d85b1b0..8ac8e2fe 100644
--- a/test/test_mueller_convolver.py
+++ b/test/test_mueller_convolver.py
@@ -1,6 +1,7 @@
 import ducc0
 import numpy as np
 import pytest
+import healpy
 from litebird_sim import MuellerConvolver
 
 
@@ -62,3 +63,55 @@ def test_trivial_mueller_matrix(fct, lkmax, ncomp):
     )
     ref_sig = ref_conv.interpol(ptg)[0] * fct
     np.testing.assert_allclose(sig, ref_sig, atol=1e-3)
+
+
+# This test case fails for reasons I don't yet understand
+# We are using a polarized beam, which is observing a uniform, completely
+# unpolarized sky through a rotating polarizer.
+# I would expect a signal modulated with exp(i*2*alpha), but the result is
+# actually constant. Any insights are welcome!
+
+@pmp("lmax", [100])
+def test_polarized(lmax):
+    rng = np.random.default_rng(41)
+
+    ncomp = 3
+    kmax = 2
+    nptg = 100
+    epsilon = 1e-4
+    ofactor = 1.5
+    nthreads = 0
+
+    # completely dark sky
+    slm = random_alm(lmax, lmax, ncomp, rng) * 0
+    # add uniform unpolarized emission
+    slm[0,0] = 1
+    # generate a Gaussian beam using healpy
+    blm = healpy.blm_gauss(1.*np.pi/180., lmax=lmax, pol=True)
+
+    ptg = np.empty((nptg, 3))
+    ptg[:, 0] = 0.5 * np.pi
+    ptg[:, 1] = 0.
+    ptg[:, 2] = 0.
+    alpha = rng.random(nptg) * 2 * np.pi
+
+    # Linear polarizer (see last page of https://www.brown.edu/research/labs/mittleman/sites/brown.edu.research.labs.mittleman/files/uploads/lecture17_0.pdf)
+    mueller = np.zeros((4,4))
+    mueller[:2,:2] = 1
+
+    fullconv = MuellerConvolver(
+        lmax,
+        kmax,
+        slm,
+        blm,
+        mueller,
+        single_precision=False,
+        epsilon=epsilon,
+        ofactor=ofactor,
+        nthreads=nthreads,
+    )
+    sig = fullconv.signal(ptg, alpha)
+
+    # I'm testing for near-constness for now, to detect that I'm not getting the
+    # result I expect. This has to improve once we have found the bug.
+    np.testing.assert_array_less(1e-4, np.max(sig)-np.min(sig))