From 3a03c76e99db465e760f6cce42cbfb4562e4a3a6 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Fri, 26 Nov 2021 09:46:42 +0100
Subject: [PATCH 1/7] temporary, non-working state

---
 beamconv/instrument.py      | 58 +++++++++++++++++++++++++++++++++++--
 tests/test_scan_strategy.py | 25 +++++++++++-----
 2 files changed, 74 insertions(+), 9 deletions(-)

diff --git a/beamconv/instrument.py b/beamconv/instrument.py
index 8a53468..e1f7840 100644
--- a/beamconv/instrument.py
+++ b/beamconv/instrument.py
@@ -9,6 +9,7 @@
 import numpy as np
 import qpoint as qp
 import healpy as hp
+import ducc0.totalconvolve
 
 from . import tools
 from .detector import Beam
@@ -3213,7 +3214,11 @@ def scan(self, beam, add_to_tod=False, interp=False,
             # Convert qpoint output to input healpy (in-place),
             # needed for interpolation later.
             tools.radec2colatlong(ra, dec)
-
+            import matplotlib.pyplot as plt
+            plt.plot(ra)
+            plt.show()
+            plt.plot(dec)
+            plt.show()
         else:
             # In no interpolation is required, we can go straight
             # from quaternion to pix and pa.
@@ -3347,7 +3352,44 @@ def scan(self, beam, add_to_tod=False, interp=False,
             del tod_tmp
                                     
             # Normal unpolarized part.
-            self._scan_modulate_pa(tod, pix, pa,
+            if spinmaps['s0a0']['interpolators'] is not None:
+                print("Interpolator found!")
+                print("nspinmaps: ",  len(spinmaps['s0a0']['maps']))
+                print(spinmaps['s0a0']['s_vals'])
+                pa[()]=0.
+                if interp:
+                    ptg = np.zeros((pix[0].shape[0], 3))
+                    ptg[:, 0] = pix[1]
+                    ptg[:, 1] = pix[0]
+                else:
+                    nside = hp.npix2nside(len(spinmaps['s0a0']['maps'][0]))
+                    print("nside=",nside)
+                    ptg = np.zeros((len(pix), 3))
+                    tptg = hp.pix2ang(nside,pix,nest=False)
+                    ptg[:, 0] = tptg[0]
+                    ptg[:, 1] = tptg[1]
+                ptg[:, 2] = np.pi
+                print(pa)
+                print(np.min(pa), np.max(pa))
+                tmp = spinmaps['s0a0']['interpolators'].interpol(ptg)
+                tod += 2*tmp[0]
+                xtmp = tod.copy()*0
+                self._scan_modulate_pa(xtmp, pix, pa,
+                                   spinmaps['s0a0']['maps'],
+                                   spinmaps['s0a0']['s_vals'],
+                                   reality=True, interp=interp)
+                import matplotlib.pyplot as plt
+  #              xtmp -= np.mean(xtmp)
+  #              tmp[0] -= np.mean(tmp[0])
+  #              tmp[0] *= np.max(xtmp)/np.max(tmp[0])
+                plt.plot(xtmp)
+                plt.plot(2*tmp[0])
+  #              plt.plot(ptg[:,0])
+  #              plt.plot(ptg[:,1])
+                plt.show()
+                print("blaaah",xtmp-2*tmp[0])
+            else:
+                self._scan_modulate_pa(tod, pix, pa,
                                    spinmaps['s0a0']['maps'],
                                    spinmaps['s0a0']['s_vals'],
                                    reality=True, interp=interp)
@@ -3652,6 +3694,9 @@ def _init_spinmaps(alm, blm, max_spin, nside,
         almB = alm[2]
 
 
+        xspin = 0 if symmetric else max_spin
+        nblm = ((xspin+1)*(xspin+2))//2 + (xspin+1)*(lmax-xspin)
+
         if hwp_mueller is not None:
 
             hwp_spin = tools.mueller2spin(hwp_mueller)
@@ -3663,12 +3708,21 @@ def _init_spinmaps(alm, blm, max_spin, nside,
             spinmap_dict['s0a0']['maps'] = ScanStrategy._spinmaps_real(
                 alm[0], blm_s0a0, spin_values_unpol, nside)
 
+
             if input_v:
                 blm_s0a0_v = ScanStrategy.blmxhwp(blm, hwp_spin, 's0a0_v',
                                                     beam_v=beam_v)
                 spinmap_dict['s0a0']['maps'] += ScanStrategy._spinmaps_real(
                     alm[3], blm_s0a0_v, spin_values_unpol, nside)
+                t_slm = np.array([alm[0],alm[3]])
+                t_blm = np.array([blm_s0a0[:nblm],blm_s0a0_v[:nblm]])
                 del blm_s0a0_v
+            else:
+                t_slm = np.array([alm[0]])
+                t_blm = np.array([blm_s0a0[:nblm]])
+
+            spinmap_dict['s0a0']['interpolators'] = ducc0.totalconvolve.Interpolator(
+                t_slm.copy(), t_blm.copy(), False, lmax, xspin, 1e-6, 2., 1)
 
             spinmap_dict['s0a0']['s_vals'] = spin_values_unpol
             del blm_s0a0
diff --git a/tests/test_scan_strategy.py b/tests/test_scan_strategy.py
index 271edb1..ec805b4 100644
--- a/tests/test_scan_strategy.py
+++ b/tests/test_scan_strategy.py
@@ -25,9 +25,12 @@ def rand_alm(lmax):
             alm += 1j * np.random.randn(hp.Alm.getsize(lmax))
             # Make m=0 modes real.
             alm[:lmax+1] = np.real(alm[:lmax+1])
+#            alm[:] = 0
+            alm[0]= 0
             return alm
 
         cls.alm = tuple([rand_alm(lmax) for i in range(3)])
+        cls.alm = (cls.alm[0], cls.alm[1]*0, cls.alm[2]*0)
         cls.lmax = lmax
 
     def test_init(self):
@@ -243,8 +246,10 @@ def test_scan_spole(self):
         tod = scs.scan(beam, return_tod=True, **chunk)
         self.assertEqual(tod.size, chunk['end'] - chunk['start'])
 
+        print("gnampfx")
         pix, nside_out, pa, hwp_ang = scs.scan(beam, return_point=True,
                                            **chunk)
+        print("gnampfy")
         self.assertEqual(pix.size, tod.size)
         self.assertEqual(nside, nside_out)
         self.assertEqual(pa.size, tod.size)
@@ -253,8 +258,10 @@ def test_scan_spole(self):
         # Turn on HWP
         scs.set_hwp_mod(mode='continuous', freq=1., start_ang=0)
         scs.rotate_hwp(**chunk)
+        print("gnampf1")
         tod2, pix2, nside_out2, pa2, hwp_ang2 = scs.scan(beam,
                         return_tod=True, return_point=True, **chunk)
+        print("gnampf2")
         np.testing.assert_almost_equal(pix, pix2)
         np.testing.assert_almost_equal(pix, pix2)
         np.testing.assert_almost_equal(pa, pa2)
@@ -1218,11 +1225,12 @@ def test_scan_spole_hwp_mueller(self):
         mmax = 2
         ra0=-10
         dec0=-57.5
-        fwhm = 200
-        nside = 128
+        fwhm = 0
+        nside = 1024
         az_throw = 10
-        polang = 20.
+        polang = 0.
 
+        print("Start")
         ces_opts = dict(ra0=ra0, dec0=dec0, az_throw=az_throw,
                         scan_speed=2.)
 
@@ -1233,10 +1241,11 @@ def test_scan_spole_hwp_mueller(self):
                                lmax=self.lmax, fwhm=fwhm,
                                polang=polang)
         beam = scs.beams[0][0]
+        print(beam)
         hwp_mueller = np.asarray([[1, 0, 0, 0],
-                                  [0, 1, 0, 0],
-                                  [0, 0, -1, 0],
-                                  [0, 0, 0, -1]])
+                                  [0, 0, 0, 0],
+                                  [0, 0, 0, 0],
+                                  [0, 0, 0, 0]])
         beam.hwp_mueller = hwp_mueller
         scs.init_detpair(self.alm, beam, nside_spin=nside,
                                    max_spin=mmax)
@@ -1251,8 +1260,10 @@ def test_scan_spole_hwp_mueller(self):
         # Turn on HWP
         scs.set_hwp_mod(mode='continuous', freq=1., start_ang=0)
         scs.rotate_hwp(**chunk)
+        print("Knampf1")
         tod, pix, nside_out, pa, hwp_ang = scs.scan(beam,
-                        return_tod=True, return_point=True, **chunk)
+                        return_tod=True, return_point=True, interp=True,  **chunk)
+        print("Knampf2")
 
         # Construct TOD manually.
         polang = beam.polang

From b873b477e67391f5c8a695e67f4b48a741f225fc Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Mon, 29 Nov 2021 12:27:02 +0100
Subject: [PATCH 2/7] s0a0 replacement works!

---
 beamconv/instrument.py      | 49 +++++++++++--------------------------
 tests/test_scan_strategy.py | 25 ++++++-------------
 2 files changed, 21 insertions(+), 53 deletions(-)

diff --git a/beamconv/instrument.py b/beamconv/instrument.py
index e1f7840..47dd32a 100644
--- a/beamconv/instrument.py
+++ b/beamconv/instrument.py
@@ -3214,11 +3214,7 @@ def scan(self, beam, add_to_tod=False, interp=False,
             # Convert qpoint output to input healpy (in-place),
             # needed for interpolation later.
             tools.radec2colatlong(ra, dec)
-            import matplotlib.pyplot as plt
-            plt.plot(ra)
-            plt.show()
-            plt.plot(dec)
-            plt.show()
+
         else:
             # In no interpolation is required, we can go straight
             # from quaternion to pix and pa.
@@ -3353,10 +3349,6 @@ def scan(self, beam, add_to_tod=False, interp=False,
                                     
             # Normal unpolarized part.
             if spinmaps['s0a0']['interpolators'] is not None:
-                print("Interpolator found!")
-                print("nspinmaps: ",  len(spinmaps['s0a0']['maps']))
-                print(spinmaps['s0a0']['s_vals'])
-                pa[()]=0.
                 if interp:
                     ptg = np.zeros((pix[0].shape[0], 3))
                     ptg[:, 0] = pix[1]
@@ -3368,26 +3360,8 @@ def scan(self, beam, add_to_tod=False, interp=False,
                     tptg = hp.pix2ang(nside,pix,nest=False)
                     ptg[:, 0] = tptg[0]
                     ptg[:, 1] = tptg[1]
-                ptg[:, 2] = np.pi
-                print(pa)
-                print(np.min(pa), np.max(pa))
-                tmp = spinmaps['s0a0']['interpolators'].interpol(ptg)
-                tod += 2*tmp[0]
-                xtmp = tod.copy()*0
-                self._scan_modulate_pa(xtmp, pix, pa,
-                                   spinmaps['s0a0']['maps'],
-                                   spinmaps['s0a0']['s_vals'],
-                                   reality=True, interp=interp)
-                import matplotlib.pyplot as plt
-  #              xtmp -= np.mean(xtmp)
-  #              tmp[0] -= np.mean(tmp[0])
-  #              tmp[0] *= np.max(xtmp)/np.max(tmp[0])
-                plt.plot(xtmp)
-                plt.plot(2*tmp[0])
-  #              plt.plot(ptg[:,0])
-  #              plt.plot(ptg[:,1])
-                plt.show()
-                print("blaaah",xtmp-2*tmp[0])
+                ptg[:, 2] = pa  # FIXME: np.radians(pa)?
+                tod += spinmaps['s0a0']['interpolators'].interpol(ptg)[0]
             else:
                 self._scan_modulate_pa(tod, pix, pa,
                                    spinmaps['s0a0']['maps'],
@@ -3714,15 +3688,20 @@ def _init_spinmaps(alm, blm, max_spin, nside,
                                                     beam_v=beam_v)
                 spinmap_dict['s0a0']['maps'] += ScanStrategy._spinmaps_real(
                     alm[3], blm_s0a0_v, spin_values_unpol, nside)
-                t_slm = np.array([alm[0],alm[3]])
-                t_blm = np.array([blm_s0a0[:nblm],blm_s0a0_v[:nblm]])
+                t_slm = np.array([alm[0].copy(),alm[3].copy()])
+                t_blm = np.array([blm_s0a0[:nblm].copy(),blm_s0a0_v[:nblm].copy()])
                 del blm_s0a0_v
             else:
-                t_slm = np.array([alm[0]])
-                t_blm = np.array([blm_s0a0[:nblm]])
-
+                t_slm = np.array([alm[0].copy()])
+                t_blm = np.array([blm_s0a0[:nblm].copy()])
+
+            lfac = np.sqrt((1.+2*np.arange(lmax+1.))/(4*np.pi))
+            ofs=0
+            for m in range(xspin+1):
+                t_blm[:, ofs:ofs+lmax+1-m] *= lfac[m:].reshape((1,-1))
+                ofs += lmax+1-m
             spinmap_dict['s0a0']['interpolators'] = ducc0.totalconvolve.Interpolator(
-                t_slm.copy(), t_blm.copy(), False, lmax, xspin, 1e-6, 2., 1)
+                t_slm, t_blm, False, lmax, xspin, 1e-11, 2., 1)
 
             spinmap_dict['s0a0']['s_vals'] = spin_values_unpol
             del blm_s0a0
diff --git a/tests/test_scan_strategy.py b/tests/test_scan_strategy.py
index ec805b4..271edb1 100644
--- a/tests/test_scan_strategy.py
+++ b/tests/test_scan_strategy.py
@@ -25,12 +25,9 @@ def rand_alm(lmax):
             alm += 1j * np.random.randn(hp.Alm.getsize(lmax))
             # Make m=0 modes real.
             alm[:lmax+1] = np.real(alm[:lmax+1])
-#            alm[:] = 0
-            alm[0]= 0
             return alm
 
         cls.alm = tuple([rand_alm(lmax) for i in range(3)])
-        cls.alm = (cls.alm[0], cls.alm[1]*0, cls.alm[2]*0)
         cls.lmax = lmax
 
     def test_init(self):
@@ -246,10 +243,8 @@ def test_scan_spole(self):
         tod = scs.scan(beam, return_tod=True, **chunk)
         self.assertEqual(tod.size, chunk['end'] - chunk['start'])
 
-        print("gnampfx")
         pix, nside_out, pa, hwp_ang = scs.scan(beam, return_point=True,
                                            **chunk)
-        print("gnampfy")
         self.assertEqual(pix.size, tod.size)
         self.assertEqual(nside, nside_out)
         self.assertEqual(pa.size, tod.size)
@@ -258,10 +253,8 @@ def test_scan_spole(self):
         # Turn on HWP
         scs.set_hwp_mod(mode='continuous', freq=1., start_ang=0)
         scs.rotate_hwp(**chunk)
-        print("gnampf1")
         tod2, pix2, nside_out2, pa2, hwp_ang2 = scs.scan(beam,
                         return_tod=True, return_point=True, **chunk)
-        print("gnampf2")
         np.testing.assert_almost_equal(pix, pix2)
         np.testing.assert_almost_equal(pix, pix2)
         np.testing.assert_almost_equal(pa, pa2)
@@ -1225,12 +1218,11 @@ def test_scan_spole_hwp_mueller(self):
         mmax = 2
         ra0=-10
         dec0=-57.5
-        fwhm = 0
-        nside = 1024
+        fwhm = 200
+        nside = 128
         az_throw = 10
-        polang = 0.
+        polang = 20.
 
-        print("Start")
         ces_opts = dict(ra0=ra0, dec0=dec0, az_throw=az_throw,
                         scan_speed=2.)
 
@@ -1241,11 +1233,10 @@ def test_scan_spole_hwp_mueller(self):
                                lmax=self.lmax, fwhm=fwhm,
                                polang=polang)
         beam = scs.beams[0][0]
-        print(beam)
         hwp_mueller = np.asarray([[1, 0, 0, 0],
-                                  [0, 0, 0, 0],
-                                  [0, 0, 0, 0],
-                                  [0, 0, 0, 0]])
+                                  [0, 1, 0, 0],
+                                  [0, 0, -1, 0],
+                                  [0, 0, 0, -1]])
         beam.hwp_mueller = hwp_mueller
         scs.init_detpair(self.alm, beam, nside_spin=nside,
                                    max_spin=mmax)
@@ -1260,10 +1251,8 @@ def test_scan_spole_hwp_mueller(self):
         # Turn on HWP
         scs.set_hwp_mod(mode='continuous', freq=1., start_ang=0)
         scs.rotate_hwp(**chunk)
-        print("Knampf1")
         tod, pix, nside_out, pa, hwp_ang = scs.scan(beam,
-                        return_tod=True, return_point=True, interp=True,  **chunk)
-        print("Knampf2")
+                        return_tod=True, return_point=True, **chunk)
 
         # Construct TOD manually.
         polang = beam.polang

From 61a0a7b8917bb6bf28a5d628a7f351deed5e26a4 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Mon, 29 Nov 2021 15:55:13 +0100
Subject: [PATCH 3/7] more total convolution

---
 beamconv/instrument.py    | 81 ++++++++++++++++-----------------------
 beamconv/totalconvolve.py | 58 ++++++++++++++++++++++++++++
 2 files changed, 90 insertions(+), 49 deletions(-)
 create mode 100644 beamconv/totalconvolve.py

diff --git a/beamconv/instrument.py b/beamconv/instrument.py
index 47dd32a..0220b43 100644
--- a/beamconv/instrument.py
+++ b/beamconv/instrument.py
@@ -9,9 +9,9 @@
 import numpy as np
 import qpoint as qp
 import healpy as hp
-import ducc0.totalconvolve
 
 from . import tools
+from . import totalconvolve
 from .detector import Beam
 
 class MPIBase(object):
@@ -3267,10 +3267,16 @@ def scan(self, beam, add_to_tod=False, interp=False,
 
         else:
             # New behaviour.
-            self._scan_modulate_pa(tod_c, pix, pa,
-                                   spinmaps['s2a4']['maps'],
-                                   spinmaps['s2a4']['s_vals'],
-                                   reality=False, interp=interp)
+
+            if spinmaps['s0a0']['interpolators'] is not None:
+                if interp:
+                    theta, phi = pix[1], pix[0]
+                else:
+                    nside = hp.npix2nside(len(spinmaps['s0a0']['maps'][0]))
+                    theta, phi = hp.pix2ang(nside,pix,nest=False)
+                psi = pa  # FIXME: np.radians(pa)?
+
+            tod_c += spinmaps['s2a4']['interpolators'].interpol(theta, phi, psi)
 
             # Modulate by HWP angle and polarization angle.
             expm2 = np.exp(1j * (4 * hwp_ang + 2 * np.radians(polang)))
@@ -3279,20 +3285,22 @@ def scan(self, beam, add_to_tod=False, interp=False,
             tod = np.real(tod_c).copy()
 
             tod_c *= 0.
-            self._scan_modulate_pa(tod_c, pix, pa,
-                                   spinmaps['s2a2']['maps'],
-                                   spinmaps['s2a2']['s_vals'],
-                                   reality=False, interp=interp)
+            # self._scan_modulate_pa(tod_c, pix, pa,
+                                   # spinmaps['s2a2']['maps'],
+                                   # spinmaps['s2a2']['s_vals'],
+                                   # reality=False, interp=interp)
+            tod_c += spinmaps['s2a2']['interpolators'].interpol(theta, phi, psi)
 
             expm2 = np.exp(1j * (2 * hwp_ang))
             tod_c *= expm2
             tod += np.real(tod_c)
 
             tod_c *= 0
-            self._scan_modulate_pa(tod_c, pix, pa,
-                                   spinmaps['s2a0']['maps'],
-                                   spinmaps['s2a0']['s_vals'],
-                                   reality=False, interp=interp)
+            # self._scan_modulate_pa(tod_c, pix, pa,
+                                   # spinmaps['s2a0']['maps'],
+                                   # spinmaps['s2a0']['s_vals'],
+                                   # reality=False, interp=interp)
+            tod_c += spinmaps['s2a0']['interpolators'].interpol(theta, phi, psi)
 
             tod += np.real(tod_c)
             del tod_c
@@ -3347,26 +3355,7 @@ def scan(self, beam, add_to_tod=False, interp=False,
             tod += tod_tmp
             del tod_tmp
                                     
-            # Normal unpolarized part.
-            if spinmaps['s0a0']['interpolators'] is not None:
-                if interp:
-                    ptg = np.zeros((pix[0].shape[0], 3))
-                    ptg[:, 0] = pix[1]
-                    ptg[:, 1] = pix[0]
-                else:
-                    nside = hp.npix2nside(len(spinmaps['s0a0']['maps'][0]))
-                    print("nside=",nside)
-                    ptg = np.zeros((len(pix), 3))
-                    tptg = hp.pix2ang(nside,pix,nest=False)
-                    ptg[:, 0] = tptg[0]
-                    ptg[:, 1] = tptg[1]
-                ptg[:, 2] = pa  # FIXME: np.radians(pa)?
-                tod += spinmaps['s0a0']['interpolators'].interpol(ptg)[0]
-            else:
-                self._scan_modulate_pa(tod, pix, pa,
-                                   spinmaps['s0a0']['maps'],
-                                   spinmaps['s0a0']['s_vals'],
-                                   reality=True, interp=interp)
+            tod += spinmaps['s0a0']['interpolators'].interpol(theta, phi, psi)
 
         if add_to_tod and hasattr(self, 'tod'):
             self.tod += tod
@@ -3668,9 +3657,6 @@ def _init_spinmaps(alm, blm, max_spin, nside,
         almB = alm[2]
 
 
-        xspin = 0 if symmetric else max_spin
-        nblm = ((xspin+1)*(xspin+2))//2 + (xspin+1)*(lmax-xspin)
-
         if hwp_mueller is not None:
 
             hwp_spin = tools.mueller2spin(hwp_mueller)
@@ -3682,26 +3668,17 @@ def _init_spinmaps(alm, blm, max_spin, nside,
             spinmap_dict['s0a0']['maps'] = ScanStrategy._spinmaps_real(
                 alm[0], blm_s0a0, spin_values_unpol, nside)
 
-
             if input_v:
                 blm_s0a0_v = ScanStrategy.blmxhwp(blm, hwp_spin, 's0a0_v',
                                                     beam_v=beam_v)
+                spinmap_dict['s0a0']['interpolators'] = totalconvolve.Interpolator_real(
+                    [alm[0], alm[3]], [blm_s0a0, blm_s0a0_v], lmax, 0 if symmetric else max_spin)
                 spinmap_dict['s0a0']['maps'] += ScanStrategy._spinmaps_real(
                     alm[3], blm_s0a0_v, spin_values_unpol, nside)
-                t_slm = np.array([alm[0].copy(),alm[3].copy()])
-                t_blm = np.array([blm_s0a0[:nblm].copy(),blm_s0a0_v[:nblm].copy()])
                 del blm_s0a0_v
             else:
-                t_slm = np.array([alm[0].copy()])
-                t_blm = np.array([blm_s0a0[:nblm].copy()])
-
-            lfac = np.sqrt((1.+2*np.arange(lmax+1.))/(4*np.pi))
-            ofs=0
-            for m in range(xspin+1):
-                t_blm[:, ofs:ofs+lmax+1-m] *= lfac[m:].reshape((1,-1))
-                ofs += lmax+1-m
-            spinmap_dict['s0a0']['interpolators'] = ducc0.totalconvolve.Interpolator(
-                t_slm, t_blm, False, lmax, xspin, 1e-11, 2., 1)
+                spinmap_dict['s0a0']['interpolators'] = totalconvolve.Interpolator_real(
+                    [alm[0]], [blm_s0a0], lmax, 0 if symmetric else max_spin)
 
             spinmap_dict['s0a0']['s_vals'] = spin_values_unpol
             del blm_s0a0
@@ -3712,6 +3689,8 @@ def _init_spinmaps(alm, blm, max_spin, nside,
             blmE, blmB = tools.spin2eb(blmm2, blmp2)
             spinmap_dict['s2a4']['maps'] = ScanStrategy._spinmaps_complex(
                 almE, almB, blmE, blmB, spin_values_pol, nside)
+            spinmap_dict['s2a4']['interpolators'] = totalconvolve.Interpolator_complex(
+                almE, almB, blmE, blmB, lmax, 0 if symmetric else max_spin)
             spinmap_dict['s2a4']['s_vals'] = spin_values_pol
 
             # s0a2.
@@ -3736,6 +3715,8 @@ def _init_spinmaps(alm, blm, max_spin, nside,
             blmE, blmB = tools.spin2eb(blmm2, blmp2)
             spinmap_dict['s2a2']['maps'] = ScanStrategy._spinmaps_complex(
                 almE, almB, blmE, blmB, spin_values_pol, nside)
+            spinmap_dict['s2a2']['interpolators'] = totalconvolve.Interpolator_complex(
+                almE, almB, blmE, blmB, lmax, 0 if symmetric else max_spin)
             spinmap_dict['s2a2']['s_vals'] = spin_values_pol
 
             # s2a0.
@@ -3744,6 +3725,8 @@ def _init_spinmaps(alm, blm, max_spin, nside,
             blmE, blmB = tools.spin2eb(blmm2, blmp2)
             spinmap_dict['s2a0']['maps'] = ScanStrategy._spinmaps_complex(
                 almE, almB, blmE, blmB, spin_values_pol, nside)
+            spinmap_dict['s2a0']['interpolators'] = totalconvolve.Interpolator_complex(
+                almE, almB, blmE, blmB, lmax, 0 if symmetric else max_spin)
             spinmap_dict['s2a0']['s_vals'] = spin_values_pol
 
         else:
diff --git a/beamconv/totalconvolve.py b/beamconv/totalconvolve.py
new file mode 100644
index 0000000..5faaad4
--- /dev/null
+++ b/beamconv/totalconvolve.py
@@ -0,0 +1,58 @@
+import ducc0.totalconvolve
+import numpy as np
+
+
+def _convert_blm(blm_in, lmax, kmax):
+    nblm = ((kmax+1)*(kmax+2))//2 + (kmax+1)*(lmax-kmax)
+    blm = np.array([x[:nblm].copy() for x in blm_in])
+    lfac = np.sqrt((1.+2*np.arange(lmax+1.))/(4*np.pi))
+    ofs=0
+    for m in range(kmax+1):
+       blm[:, ofs:ofs+lmax+1-m] *= lfac[m:].reshape((1,-1))
+       ofs += lmax+1-m
+    return blm
+def _convert_blm2(blm_in, lmax, kmax):
+    nblm = ((kmax+1)*(kmax+2))//2 + (kmax+1)*(lmax-kmax)
+    blm = blm_in.copy()
+    lfac = np.sqrt((1.+2*np.arange(lmax+1.))/(4*np.pi))
+    ofs=0
+    for m in range(kmax+1):
+       blm[:, ofs:ofs+lmax+1-m] = (-1)**m * np.conj(blm[:, ofs:ofs+lmax+1-m])
+       ofs += lmax+1-m
+    blm= blm[-1::-1]
+    return blm
+
+
+class Interpolator_real:
+    def __init__(self, slm, blm, lmax, kmax, epsilon=1e-11, nthreads=1):
+        _slm = np.array([x.copy() for x in slm])
+        _blm = _convert_blm(blm, lmax, kmax)
+        self._inter = ducc0.totalconvolve.Interpolator(
+            _slm, _blm, False, lmax, kmax, epsilon, 2., nthreads)
+
+    def interpol(self, theta, phi, psi):
+        ptg = np.zeros((theta.shape[0], 3))
+        ptg[:, 0] = theta
+        ptg[:, 1] = phi
+        ptg[:, 2] = psi
+        return self._inter.interpol(ptg)[0]
+
+
+class Interpolator_complex:
+    def __init__(self, slmE, slmB, blmE, blmB, lmax, kmax, epsilon=1e-11, nthreads=1):
+        _slm = np.array([slmE, slmB])
+        _blm = _convert_blm([blmE, blmB], lmax, kmax)
+        self._inter = ducc0.totalconvolve.Interpolator(
+            _slm, _blm, False, lmax, kmax, epsilon, 2., nthreads)
+        _blm2 = _convert_blm2(_blm, lmax, kmax)
+        self._inter2 = ducc0.totalconvolve.Interpolator(
+            _slm, _blm2, False, lmax, kmax, epsilon, 2., nthreads)
+
+    def interpol(self, theta, phi, psi):
+        ptg = np.zeros((theta.shape[0], 3))
+        ptg[:, 0] = theta
+        ptg[:, 1] = phi
+        ptg[:, 2] = psi
+        res = self._inter.interpol(ptg)
+        res2 = self._inter2.interpol(ptg)
+        return res[0]+1j*res2[0]

From 3eb714cbb0c4c9ad1f75ead279e1792670a545d6 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Mon, 29 Nov 2021 16:02:46 +0100
Subject: [PATCH 4/7] remove obsolescent healpy arguments

---
 beamconv/instrument.py      |  8 ++++----
 beamconv/test.py            | 20 ++++++++++----------
 tests/test_scan_strategy.py | 18 ++++++++----------
 3 files changed, 22 insertions(+), 24 deletions(-)

diff --git a/beamconv/instrument.py b/beamconv/instrument.py
index 0220b43..d3f981c 100644
--- a/beamconv/instrument.py
+++ b/beamconv/instrument.py
@@ -3274,7 +3274,7 @@ def scan(self, beam, add_to_tod=False, interp=False,
                 else:
                     nside = hp.npix2nside(len(spinmaps['s0a0']['maps'][0]))
                     theta, phi = hp.pix2ang(nside,pix,nest=False)
-                psi = pa  # FIXME: np.radians(pa)?
+                psi = pa
 
             tod_c += spinmaps['s2a4']['interpolators'].interpol(theta, phi, psi)
 
@@ -3899,7 +3899,7 @@ def _spinmaps_real(alm, blm, spin_values, nside):
             if s == 0: # Scalar transform.
 
                 flms = hp.almxfl(alm, bell, inplace=False)
-                func[sidx,:] = hp.alm2map(flms, nside, verbose=False)
+                func[sidx,:] = hp.alm2map(flms, nside)
 
             else: # Spin transforms.
 
@@ -3982,8 +3982,8 @@ def _spinmaps_complex(almE, almB, blmE, blmB, spin_values, nside):
 
             if s == 0:
                 # The (-1) factor for spin 0 is explained in HEALPix doc.
-                spinmaps = [hp.alm2map(-ps_flm_p, nside, verbose=False),
-                            hp.alm2map(ms_flm_m, nside, verbose=False)]
+                spinmaps = [hp.alm2map(-ps_flm_p, nside),
+                            hp.alm2map(ms_flm_m, nside)]
                 func_c[sidx,:] = spinmaps[0] + 1j * spinmaps[1]
 
             if s > 0:
diff --git a/beamconv/test.py b/beamconv/test.py
index a3cdafd..57dc5c6 100644
--- a/beamconv/test.py
+++ b/beamconv/test.py
@@ -122,8 +122,8 @@ def scan_bicep(lmax=700, mmax=5, fwhm=43, ra0=-10, dec0=-57.5,
 
         # plot smoothed input maps
         nside = hp.get_nside(maps[0])
-        hp.smoothalm(alm, fwhm=np.radians(fwhm/60.), verbose=False)
-        maps_raw = hp.alm2map(alm, nside, verbose=False)
+        hp.smoothalm(alm, fwhm=np.radians(fwhm/60.))
+        maps_raw = hp.alm2map(alm, nside)
 
         plot_iqu(maps_raw, img_out_path, 'raw_bicep',
                  sym_limits=[250, 5, 5],
@@ -265,8 +265,8 @@ def scan_atacama(lmax=700, mmax=5, fwhm=40,
 
         # plot smoothed input maps
         nside = hp.get_nside(maps[0])
-        hp.smoothalm(alm, fwhm=np.radians(fwhm/60.), verbose=False)
-        maps_raw = hp.alm2map(alm, nside, verbose=False)
+        hp.smoothalm(alm, fwhm=np.radians(fwhm/60.))
+        maps_raw = hp.alm2map(alm, nside)
 
         plot_iqu(maps_raw, img_out_path, 'raw_atacama',
                  sym_limits=[250, 5, 5],
@@ -925,8 +925,8 @@ def test_ghosts(lmax=700, mmax=5, fwhm=43, ra0=-10, dec0=-57.5,
 
         # plot smoothed input maps
         nside = hp.get_nside(maps[0])
-        hp.smoothalm(alm, fwhm=np.radians(fwhm/60.), verbose=False)
-        maps_raw = hp.alm2map(alm, nside, verbose=False)
+        hp.smoothalm(alm, fwhm=np.radians(fwhm/60.))
+        maps_raw = hp.alm2map(alm, nside)
 
         plot_iqu(maps_raw, '../scratch/img/', 'raw_ghost',
                  sym_limits=[250, 5, 5],
@@ -1162,8 +1162,8 @@ def idea_jon():
 
     # plot smoothed input maps
     nside = hp.get_nside(maps_g[0])
-    hp.smoothalm(alm, fwhm=np.radians(fwhm/60.), verbose=False)
-    maps_raw = hp.alm2map(alm, nside, verbose=False)
+    hp.smoothalm(alm, fwhm=np.radians(fwhm/60.))
+    maps_raw = hp.alm2map(alm, nside)
 
     plot_iqu(maps_raw, '../scratch/img/', 'raw_delta',
              sym_limits=[1, 1, 1],
@@ -1334,8 +1334,8 @@ def test_satellite_scan(lmax=700, mmax=2, fwhm=43,
 
         # plot smoothed input maps
         nside = hp.get_nside(maps[0])
-        hp.smoothalm(alm, fwhm=np.radians(fwhm/60.), verbose=False)
-        maps_raw = hp.alm2map(alm, nside, verbose=False)
+        hp.smoothalm(alm, fwhm=np.radians(fwhm/60.))
+        maps_raw = hp.alm2map(alm, nside)
 
         plot_iqu(maps_raw, '../scratch/img/', 'raw_satellite',
                  sym_limits=[250, 5, 5],
diff --git a/tests/test_scan_strategy.py b/tests/test_scan_strategy.py
index 271edb1..b203463 100644
--- a/tests/test_scan_strategy.py
+++ b/tests/test_scan_strategy.py
@@ -156,7 +156,7 @@ def test_init_detpair(self):
         # Since we have a infinitely narrow Gaussian the convolved
         # maps should just match the input (up to healpix quadrature
         # wonkyness).
-        input_map = hp.alm2map(self.alm, nside, verbose=False) # I, Q, U
+        input_map = hp.alm2map(self.alm, nside) # I, Q, U
         zero_map = np.zeros_like(input_map[0])
         np.testing.assert_array_almost_equal(input_map[0],
                                              func[0], decimal=6)
@@ -264,7 +264,7 @@ def test_scan_spole(self):
 
         # Construct TOD manually.
         polang = beam.polang
-        maps_sm = np.asarray(hp.alm2map(self.alm, nside, verbose=False,
+        maps_sm = np.asarray(hp.alm2map(self.alm, nside,
                                         fwhm=np.radians(beam.fwhm / 60.)))
 
         np.testing.assert_almost_equal(maps_sm[0],
@@ -346,7 +346,7 @@ def test_scan_spole_pol(self):
 
         # Construct TOD manually.
         polang = beam.polang
-        maps_sm = np.asarray(hp.alm2map(alm, nside, verbose=False,
+        maps_sm = np.asarray(hp.alm2map(alm, nside,
                                         fwhm=np.radians(beam.fwhm / 60.)))
 
         np.testing.assert_almost_equal(maps_sm[0],
@@ -407,9 +407,8 @@ def test_scan_spole_bin(self):
         # Solve for the maps.
         maps, cond = scs.solve_for_map(fill=np.nan)
 
-        alm = hp.smoothalm(self.alm, fwhm=np.radians(fwhm/60.),
-                     verbose=False)
-        maps_raw = np.asarray(hp.alm2map(self.alm, nside, verbose=False))
+        alm = hp.smoothalm(self.alm, fwhm=np.radians(fwhm/60.))
+        maps_raw = np.asarray(hp.alm2map(self.alm, nside))
 
         cond[~np.isfinite(cond)] = 10
 
@@ -1256,7 +1255,7 @@ def test_scan_spole_hwp_mueller(self):
 
         # Construct TOD manually.
         polang = beam.polang
-        maps_sm = np.asarray(hp.alm2map(self.alm, nside, verbose=False,
+        maps_sm = np.asarray(hp.alm2map(self.alm, nside,
                                         fwhm=np.radians(beam.fwhm / 60.)))
 
         np.testing.assert_almost_equal(maps_sm[0],
@@ -1326,9 +1325,8 @@ def test_scan_spole_bin_hwp_mueller(self):
         # Solve for the maps.
         maps, cond = scs.solve_for_map(fill=np.nan)
 
-        alm = hp.smoothalm(self.alm, fwhm=np.radians(fwhm/60.),
-                     verbose=False)
-        maps_raw = np.asarray(hp.alm2map(self.alm, nside, verbose=False))
+        alm = hp.smoothalm(self.alm, fwhm=np.radians(fwhm/60.))
+        maps_raw = np.asarray(hp.alm2map(self.alm, nside))
 
         cond[~np.isfinite(cond)] = 10
 

From 4562956a4eafa3acacc3be8c5838c299c3d3e93c Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Mon, 29 Nov 2021 16:13:15 +0100
Subject: [PATCH 5/7] cleanup

---
 beamconv/instrument.py | 33 +++++++++++++--------------------
 1 file changed, 13 insertions(+), 20 deletions(-)

diff --git a/beamconv/instrument.py b/beamconv/instrument.py
index d3f981c..1e7151d 100644
--- a/beamconv/instrument.py
+++ b/beamconv/instrument.py
@@ -3268,15 +3268,14 @@ def scan(self, beam, add_to_tod=False, interp=False,
         else:
             # New behaviour.
 
-            if spinmaps['s0a0']['interpolators'] is not None:
-                if interp:
-                    theta, phi = pix[1], pix[0]
-                else:
-                    nside = hp.npix2nside(len(spinmaps['s0a0']['maps'][0]))
-                    theta, phi = hp.pix2ang(nside,pix,nest=False)
-                psi = pa
+            if interp:
+                theta, phi = pix[1], pix[0]
+            else:
+                nside = hp.npix2nside(len(spinmaps['s0a0']['maps'][0]))
+                theta, phi = hp.pix2ang(nside,pix,nest=False)
+            psi = pa
 
-            tod_c += spinmaps['s2a4']['interpolators'].interpol(theta, phi, psi)
+            tod_c = spinmaps['s2a4']['interpolators'].interpol(theta, phi, psi)
 
             # Modulate by HWP angle and polarization angle.
             expm2 = np.exp(1j * (4 * hwp_ang + 2 * np.radians(polang)))
@@ -3284,23 +3283,13 @@ def scan(self, beam, add_to_tod=False, interp=False,
 
             tod = np.real(tod_c).copy()
 
-            tod_c *= 0.
-            # self._scan_modulate_pa(tod_c, pix, pa,
-                                   # spinmaps['s2a2']['maps'],
-                                   # spinmaps['s2a2']['s_vals'],
-                                   # reality=False, interp=interp)
-            tod_c += spinmaps['s2a2']['interpolators'].interpol(theta, phi, psi)
+            tod_c = spinmaps['s2a2']['interpolators'].interpol(theta, phi, psi)
 
             expm2 = np.exp(1j * (2 * hwp_ang))
             tod_c *= expm2
             tod += np.real(tod_c)
 
-            tod_c *= 0
-            # self._scan_modulate_pa(tod_c, pix, pa,
-                                   # spinmaps['s2a0']['maps'],
-                                   # spinmaps['s2a0']['s_vals'],
-                                   # reality=False, interp=interp)
-            tod_c += spinmaps['s2a0']['interpolators'].interpol(theta, phi, psi)
+            tod_c = spinmaps['s2a0']['interpolators'].interpol(theta, phi, psi)
 
             tod += np.real(tod_c)
             del tod_c
@@ -3700,12 +3689,16 @@ def _init_spinmaps(alm, blm, max_spin, nside,
             # Minus sign is because eb2spin applied to alm_I results in (-1) alm_I.
             spinmap_dict['s0a2']['maps'] = ScanStrategy._spinmaps_complex(
                 -alm[0], alm[0] * 0, blmE, blmB, spin_values_s0a2, nside)
+            spinmap_dict['s0a2']['interpolators0'] = totalconvolve.Interpolator_complex(
+                -alm[0], alm[0]*0, blmE, blmB, lmax, 0 if symmetric else max_spin)
 
             if input_v:
                 blmm2, blmp2 = ScanStrategy.blmxhwp(blm, hwp_spin, 's0a2_v')
                 blmE, blmB = tools.spin2eb(blmp2, blmm2)
                 spinmap_dict['s0a2']['maps'] += ScanStrategy._spinmaps_complex(
                     -alm[3], alm[3] * 0, blmE, blmB, spin_values_s0a2, nside)
+                spinmap_dict['s0a2']['interpolators3'] = totalconvolve.Interpolator_complex(
+                    -alm[3], alm[3]*0, blmE, blmB, lmax, 0 if symmetric else max_spin)
             spinmap_dict['s0a2']['s_vals'] = spin_values_s0a2
 
             # s2a2.

From 7fc9771f6d86128c95446f03075e0e4c8e9fd833 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Mon, 29 Nov 2021 16:13:56 +0100
Subject: [PATCH 6/7] cleanup

---
 beamconv/instrument.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/beamconv/instrument.py b/beamconv/instrument.py
index 1e7151d..08d59cf 100644
--- a/beamconv/instrument.py
+++ b/beamconv/instrument.py
@@ -3268,6 +3268,7 @@ def scan(self, beam, add_to_tod=False, interp=False,
         else:
             # New behaviour.
 
+            # obtain theta, phi, psi pointings
             if interp:
                 theta, phi = pix[1], pix[0]
             else:

From 4584a26c0bd4cb091427c775a29a398e751964c4 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Tue, 30 Nov 2021 09:48:51 +0100
Subject: [PATCH 7/7] start migrating old behaviour

---
 beamconv/instrument.py    | 49 +++++++++++++++++++++++----------------
 beamconv/totalconvolve.py | 16 ++++++++-----
 2 files changed, 39 insertions(+), 26 deletions(-)

diff --git a/beamconv/instrument.py b/beamconv/instrument.py
index 08d59cf..cb6b09f 100644
--- a/beamconv/instrument.py
+++ b/beamconv/instrument.py
@@ -3245,9 +3245,19 @@ def scan(self, beam, add_to_tod=False, interp=False,
         else:
             hwp_ang = self.hwp_ang
 
+        # obtain theta, phi, psi pointings
+        if interp:
+            theta, phi = pix[1], pix[0]
+        else:
+            nside = hp.npix2nside(len(spinmaps['s0a0']['maps'][0]))
+            theta, phi = hp.pix2ang(nside,pix,nest=False)
+        psi = pa
+
         # Find out if old or new HWP behaviour is desired.
         if set(spinmaps.keys()) == set(['s0a0', 's2a4']):
             # Old behaviour.
+# For some reason, this causes a single test to fail...
+#            tod_c += spinmaps['s2a4']['interpolators'].interpol(theta, phi, psi)
             self._scan_modulate_pa(tod_c, pix, pa,
                                    spinmaps['s2a4']['maps'],
                                    spinmaps['s2a4']['s_vals'],
@@ -3260,22 +3270,14 @@ def scan(self, beam, add_to_tod=False, interp=False,
             tod = np.real(tod_c) # Note, shares memory with tod_c.
 
             # Add unpolarized tod.
-            self._scan_modulate_pa(tod, pix, pa,
-                                   spinmaps['s0a0']['maps'],
-                                   spinmaps['s0a0']['s_vals'],
-                                   reality=True, interp=interp)
+            tod += spinmaps['s0a0']['interpolators'].interpol(theta, phi, psi)
+            # self._scan_modulate_pa(tod, pix, pa,
+                                   # spinmaps['s0a0']['maps'],
+                                   # spinmaps['s0a0']['s_vals'],
+                                   # reality=True, interp=interp)
 
         else:
             # New behaviour.
-
-            # obtain theta, phi, psi pointings
-            if interp:
-                theta, phi = pix[1], pix[0]
-            else:
-                nside = hp.npix2nside(len(spinmaps['s0a0']['maps'][0]))
-                theta, phi = hp.pix2ang(nside,pix,nest=False)
-            psi = pa
-
             tod_c = spinmaps['s2a4']['interpolators'].interpol(theta, phi, psi)
 
             # Modulate by HWP angle and polarization angle.
@@ -3662,13 +3664,13 @@ def _init_spinmaps(alm, blm, max_spin, nside,
                 blm_s0a0_v = ScanStrategy.blmxhwp(blm, hwp_spin, 's0a0_v',
                                                     beam_v=beam_v)
                 spinmap_dict['s0a0']['interpolators'] = totalconvolve.Interpolator_real(
-                    [alm[0], alm[3]], [blm_s0a0, blm_s0a0_v], lmax, 0 if symmetric else max_spin)
+                    [alm[0], alm[3]], [blm_s0a0, blm_s0a0_v], lmax, spin_values_unpol)
                 spinmap_dict['s0a0']['maps'] += ScanStrategy._spinmaps_real(
                     alm[3], blm_s0a0_v, spin_values_unpol, nside)
                 del blm_s0a0_v
             else:
                 spinmap_dict['s0a0']['interpolators'] = totalconvolve.Interpolator_real(
-                    [alm[0]], [blm_s0a0], lmax, 0 if symmetric else max_spin)
+                    [alm[0]], [blm_s0a0], lmax, spin_values_unpol)
 
             spinmap_dict['s0a0']['s_vals'] = spin_values_unpol
             del blm_s0a0
@@ -3680,7 +3682,7 @@ def _init_spinmaps(alm, blm, max_spin, nside,
             spinmap_dict['s2a4']['maps'] = ScanStrategy._spinmaps_complex(
                 almE, almB, blmE, blmB, spin_values_pol, nside)
             spinmap_dict['s2a4']['interpolators'] = totalconvolve.Interpolator_complex(
-                almE, almB, blmE, blmB, lmax, 0 if symmetric else max_spin)
+                almE, almB, blmE, blmB, lmax, spin_values_pol)
             spinmap_dict['s2a4']['s_vals'] = spin_values_pol
 
             # s0a2.
@@ -3691,7 +3693,7 @@ def _init_spinmaps(alm, blm, max_spin, nside,
             spinmap_dict['s0a2']['maps'] = ScanStrategy._spinmaps_complex(
                 -alm[0], alm[0] * 0, blmE, blmB, spin_values_s0a2, nside)
             spinmap_dict['s0a2']['interpolators0'] = totalconvolve.Interpolator_complex(
-                -alm[0], alm[0]*0, blmE, blmB, lmax, 0 if symmetric else max_spin)
+                -alm[0], alm[0]*0, blmE, blmB, lmax, spin_values_s0a2)
 
             if input_v:
                 blmm2, blmp2 = ScanStrategy.blmxhwp(blm, hwp_spin, 's0a2_v')
@@ -3699,7 +3701,7 @@ def _init_spinmaps(alm, blm, max_spin, nside,
                 spinmap_dict['s0a2']['maps'] += ScanStrategy._spinmaps_complex(
                     -alm[3], alm[3] * 0, blmE, blmB, spin_values_s0a2, nside)
                 spinmap_dict['s0a2']['interpolators3'] = totalconvolve.Interpolator_complex(
-                    -alm[3], alm[3]*0, blmE, blmB, lmax, 0 if symmetric else max_spin)
+                    -alm[3], alm[3]*0, blmE, blmB, lmax, spin_values_s0a2)
             spinmap_dict['s0a2']['s_vals'] = spin_values_s0a2
 
             # s2a2.
@@ -3710,7 +3712,7 @@ def _init_spinmaps(alm, blm, max_spin, nside,
             spinmap_dict['s2a2']['maps'] = ScanStrategy._spinmaps_complex(
                 almE, almB, blmE, blmB, spin_values_pol, nside)
             spinmap_dict['s2a2']['interpolators'] = totalconvolve.Interpolator_complex(
-                almE, almB, blmE, blmB, lmax, 0 if symmetric else max_spin)
+                almE, almB, blmE, blmB, lmax, spin_values_pol)
             spinmap_dict['s2a2']['s_vals'] = spin_values_pol
 
             # s2a0.
@@ -3720,7 +3722,7 @@ def _init_spinmaps(alm, blm, max_spin, nside,
             spinmap_dict['s2a0']['maps'] = ScanStrategy._spinmaps_complex(
                 almE, almB, blmE, blmB, spin_values_pol, nside)
             spinmap_dict['s2a0']['interpolators'] = totalconvolve.Interpolator_complex(
-                almE, almB, blmE, blmB, lmax, 0 if symmetric else max_spin)
+                almE, almB, blmE, blmB, lmax, spin_values_pol)
             spinmap_dict['s2a0']['s_vals'] = spin_values_pol
 
         else:
@@ -3733,8 +3735,13 @@ def _init_spinmaps(alm, blm, max_spin, nside,
                 alm[0], blm[0], spin_values_unpol, nside)
 
             if beam_v and input_v:
+                spinmap_dict['s0a0']['interpolators'] = totalconvolve.Interpolator_real(
+                    [alm[0], alm[3]], [blm[0], blm[3]], lmax, spin_values_unpol)
                 spinmap_dict['s0a0']['maps'] += ScanStrategy._spinmaps_real(
                     alm[3], blm[3], spin_values_unpol, nside)
+            else:
+                spinmap_dict['s0a0']['interpolators'] = totalconvolve.Interpolator_real(
+                    [alm[0]], [blm[0]], lmax, spin_values_unpol)
 
             spinmap_dict['s0a0']['s_vals'] = spin_values_unpol
 
@@ -3743,6 +3750,8 @@ def _init_spinmaps(alm, blm, max_spin, nside,
             spinmap_dict['s2a4']['maps'] = ScanStrategy._spinmaps_complex(
                 almE, almB, blmE, blmB, spin_values_pol, nside)
             spinmap_dict['s2a4']['s_vals'] = spin_values_pol
+            spinmap_dict['s2a4']['interpolators'] = totalconvolve.Interpolator_complex(
+                almE, almB, blmE, blmB, lmax, spin_values_pol)
 
         return spinmap_dict
 
diff --git a/beamconv/totalconvolve.py b/beamconv/totalconvolve.py
index 5faaad4..527d83f 100644
--- a/beamconv/totalconvolve.py
+++ b/beamconv/totalconvolve.py
@@ -2,13 +2,16 @@
 import numpy as np
 
 
-def _convert_blm(blm_in, lmax, kmax):
+def _convert_blm(blm_in, lmax, spin_values):
+    kmax = max(spin_values)
     nblm = ((kmax+1)*(kmax+2))//2 + (kmax+1)*(lmax-kmax)
     blm = np.array([x[:nblm].copy() for x in blm_in])
     lfac = np.sqrt((1.+2*np.arange(lmax+1.))/(4*np.pi))
     ofs=0
     for m in range(kmax+1):
        blm[:, ofs:ofs+lmax+1-m] *= lfac[m:].reshape((1,-1))
+       if (not m in spin_values):
+           blm[:, ofs:ofs+lmax+1-m] = 0
        ofs += lmax+1-m
     return blm
 def _convert_blm2(blm_in, lmax, kmax):
@@ -24,11 +27,11 @@ def _convert_blm2(blm_in, lmax, kmax):
 
 
 class Interpolator_real:
-    def __init__(self, slm, blm, lmax, kmax, epsilon=1e-11, nthreads=1):
+    def __init__(self, slm, blm, lmax, spin_values, epsilon=1e-11, nthreads=1):
         _slm = np.array([x.copy() for x in slm])
-        _blm = _convert_blm(blm, lmax, kmax)
+        _blm = _convert_blm(blm, lmax, spin_values)
         self._inter = ducc0.totalconvolve.Interpolator(
-            _slm, _blm, False, lmax, kmax, epsilon, 2., nthreads)
+            _slm, _blm, False, lmax, max(spin_values), epsilon, 2., nthreads)
 
     def interpol(self, theta, phi, psi):
         ptg = np.zeros((theta.shape[0], 3))
@@ -39,9 +42,10 @@ def interpol(self, theta, phi, psi):
 
 
 class Interpolator_complex:
-    def __init__(self, slmE, slmB, blmE, blmB, lmax, kmax, epsilon=1e-11, nthreads=1):
+    def __init__(self, slmE, slmB, blmE, blmB, lmax, spin_values, epsilon=1e-11, nthreads=1):
         _slm = np.array([slmE, slmB])
-        _blm = _convert_blm([blmE, blmB], lmax, kmax)
+        _blm = _convert_blm([blmE, blmB], lmax, spin_values)
+        kmax = max(spin_values)
         self._inter = ducc0.totalconvolve.Interpolator(
             _slm, _blm, False, lmax, kmax, epsilon, 2., nthreads)
         _blm2 = _convert_blm2(_blm, lmax, kmax)