Comparison between beamconv and ducc0/MuellerConvolver

[mreineck]

Thanks to the hard work by @martamonelli, I now have a very compact code that allows direct comparison between the beamconv and MuellerConvolver HWP simulators.

My current setup is:

• polarised Gaussian beam (TEB)
• totally random sky a_lm (up to given lmax)
• totally random vectors of theta, phi, psi, alpha
• totally random Mueller matrix

I can get very good agreement between the output of beamconv and MuellerConvolver, if

• I flip the sign of the psi angle for MuellerConvolver (This seems to be a convention issue only.)
• I multiply blm_E and blm_B by -np.sqrt(2) for MuellerConvolver (I have no idea where this factor comes from)
• I make sure that mueller[0,2] = -mueller[2,0]

The last point seems to indicate some asymmetry within beamconv. If I run beamconv with a Mueller matrix that is zero except for entry [2,0], I always get a result that is exactly zero, independent of pointing,beam or sky. The same is not true for entry [0,2].

Once we understand the remaining differences, I'd like to extend the tests by

• going to non-Gaussian beams. For that I'd need a way to supply beam a_lm directly to beamconv, I think.
• adding the V component.

test.tar.gz

[mreineck]

Looking a bit deeper into the problem concerning Mueller matrix entry [2,0]: When converting the Mueller matrix to the "C" matrix from he notes, which is (I think) the hwp_spin matrix in the code, the only two entries affected by the original [2,0] entry are hwp_spin[1,0] and hwp_spin[2,0].

Grepping the sources seems to indicate that these two entries are never accessed by the code. If that is indeed the case, I think that's a bug.

[mreineck]

I'm telling a lie ... hwp_spin[2,0] is indeed accessed in blmxhwp, but hwp_spin[1,0] isn't.

[mreineck]

OK, let's focus on one thing first ...
The code below uses a Mueller matrix where only entry [2,0] is nonzero. For this, beamconv produces a signal that is always identically zero. Should this be the case?

import numpy as np
import healpy as hp

def nalm(lmax, mmax):
    return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)


def random_alm(lmax, mmax, spin, ncomp, rng):
    res = rng.uniform(-1., 1., (ncomp, nalm(lmax, mmax))) \
     + 1j*rng.uniform(-1., 1., (ncomp, nalm(lmax, mmax)))
    # make a_lm with m==0 real-valued
    res[:, 0:lmax+1].imag = 0.
    ofs=0
    for s in range(spin):
        res[:, ofs:ofs+spin-s] = 0.
        ofs += lmax+1-s
    return res

# code by Marta to get beamconv results for user-specified angles
def get_beamconv_values(lmax, fwhm_arcmin, slm, ptg, hwp_angles, mueller):
    from beamconv import Beam, ScanStrategy, tools
    import qpoint as qp

    # set up beam and HWP mueller matrix (identity, i.e. no HWP)
    beam = Beam(btype='Gaussian', fwhm=fwhm_arcmin, lmax=lmax)
    beam.hwp_mueller = mueller

    nsamp = ptg.shape[0]
    duration = 1       #scan duration in seconds
    sample_rate = nsamp  #samples per second

    # from (theta,phi) to (ra,dec) convention
    # also, all angles are converted in degrees
    ra = np.degrees(ptg[:,1])
    dec = 90. - np.degrees(ptg[:,0])
    psi = np.degrees(ptg[:,2])

    # calculate the quaternion
    q_bore_array = qp.QPoint().radecpa2quat(ra, dec, psi)

    def ctime_test(**kwargs):
        start = kwargs.pop('start')
        end = kwargs.pop('end')
        nsamp = end - start
        ctime0 = 0
        ctime = ctime0 + sample_rate*np.arange(nsamp)
        return ctime
    
    def q_bore_test(**kwargs):
        start = kwargs.pop('start')
        end = kwargs.pop('end')
        return q_bore_array[start:end]
   
    S = ScanStrategy(duration=duration, sample_rate=sample_rate, external_pointing=True, use_l2_scan=False)
    S.add_to_focal_plane(beam, combine=False)
    S.set_hwp_mod(mode='stepped', freq=sample_rate, angles=hwp_angles*180/np.pi)
    S.allocate_maps(nside=8)
    S.scan_instrument_mpi(slm, save_tod=True, ctime_func=ctime_test, q_bore_func=q_bore_test,
                      ctime_kwargs={'useless':0}, q_bore_kwargs={'useless':0},nside_spin=1024)

    return S.data(S.chunks[0], beam=beam, data_type='tod').copy()

# set up sky alm (monopole only)
rng = np.random.default_rng(np.random.SeedSequence(42))
lmax = 30

np.random.seed(10)
slm = random_alm(lmax, lmax, 2, 3, rng)

# Mueller matrix
mueller = np.zeros((4,4))
mueller[2,0]=1

fwhm_arcmin=32.
nptg=100
# completely random pointings
ptg = np.empty((nptg,3))
ptg[:,0]=np.random.uniform(0,np.pi,size=(nptg,))    # theta
ptg[:,1]=np.random.uniform(0,2*np.pi,size=(nptg,))  # phi
ptg[:,2]=np.random.uniform(0,2*np.pi,size=(nptg,))  # psi
hwp_angles = np.random.uniform(0,2*np.pi,size=(nptg,))  # alpha

# get the signal from beamconv

signal_beamconv = get_beamconv_values(lmax=lmax, fwhm_arcmin=fwhm_arcmin, slm=slm, ptg=ptg, hwp_angles=hwp_angles, mueller=mueller)

print("maximum absolute signal:", np.max(np.abs(signal_beamconv)))

[mreineck]

If I change line 4148 of instrument.py from
blmm2 *= hwp_spin[0,2] * np.sqrt(2)
to
blmm2 *= hwp_spin[2,0] * np.sqrt(2)

then I get perfect agreement between both codes.
This is probably not the right fix, but maybe it helps someone more knowledgeable than mself to understand what is actually going wrong. It's totally possible that the bug is in mueller_convolver.

[AdriJD]

Hi Martin, thank you for the detailed examples. It's not immediately clear to me what is going wrong here. I am trying your example to see if I can figure it out. I'll get back to you before the end of the week

[mreineck]

Here is my latest version, with more generic calls to beamconv and simpler description of the discrepancies.
With this version, I get complete agreement, provided

• I change line 4148 in instrument.py to blmm2 *= hwp_spin[2,0] * np.sqrt(2)
• I change line 4154 in instrument.py to blmm2_v *= hwp_spin[3,2] * np.sqrt(2)
• I change the psi angle to np.pi-psi for MuellerConvolver
• I use Gaussian beam a_lm for the E and B components of the beam

The point about psi feels very much as if I'm missing a complex conjugation somewhere ... I just can't find it :-(
The last point may actually be closely related to the third; apparently the results only agree at the moment when the beam E and B components .

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import healpy as hp
import mueller_convolver
import ducc0
from time import time

def nalm(lmax, mmax):
    return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)

def make_full_random_alm(lmax, mmax, rng):
    res = rng.uniform(-1., 1., (4, nalm(lmax, mmax))) \
     + 1j*rng.uniform(-1., 1., (4, nalm(lmax, mmax)))
    # make a_lm with m==0 real-valued
    res[:, 0:lmax+1].imag = 0.
    ofs=0
    # components 1 and 2 are spin-2, fix them accordingly
    spin=2
    for s in range(spin):
        res[1:3, ofs:ofs+spin-s] = 0.
        ofs += lmax+1-s
    return res

# code by Marta to get beamconv results for user-specified angles
def get_beamconv_values(lmax, kmax, slm, ptg, hwp_angles, mueller, beam_file):
    from beamconv import Beam, ScanStrategy, tools
    import qpoint as qp

    # set up beam and HWP mueller matrix (identity, i.e. no HWP)
    beam = Beam(btype='PO', lmax=lmax, mmax=lmax, deconv_q=True, normalize=False, po_file=beam_file, hwp_mueller=mueller)

    nsamp = ptg.shape[0]
    duration = 1       #scan duration in seconds
    sample_rate = nsamp  #samples per second

    # from (theta,phi) to (ra,dec) convention
    # also, all angles are converted in degrees
    ra = np.degrees(ptg[:,1])
    dec = 90. - np.degrees(ptg[:,0])
    psi = np.degrees(ptg[:,2])

    # calculate the quaternion
    q_bore_array = qp.QPoint().radecpa2quat(ra, dec, psi)

    def ctime_test(**kwargs):
        start = kwargs.pop('start')
        end = kwargs.pop('end')
        nsamp = end - start
        ctime0 = 0
        ctime = ctime0 + sample_rate*np.arange(nsamp)
        return ctime
    
    def q_bore_test(**kwargs):
        start = kwargs.pop('start')
        end = kwargs.pop('end')
        return q_bore_array[start:end]

    S = ScanStrategy(duration=duration, sample_rate=sample_rate, external_pointing=True, use_l2_scan=False)
    S.add_to_focal_plane(beam, combine=False)
    S.set_hwp_mod(mode='stepped', freq=sample_rate, angles=hwp_angles*180/np.pi)
    S.allocate_maps(nside=8)
    S.scan_instrument_mpi(slm.copy(), save_tod=True, ctime_func=ctime_test, q_bore_func=q_bore_test,
                      ctime_kwargs={'useless':0}, q_bore_kwargs={'useless':0},nside_spin=1024, interp=False, input_v=True, beam_v=True, max_spin=kmax+4)

    return S.data(S.chunks[0], beam=beam, data_type='tod').copy()

def save_blm_for_beamconv (blm, beam_file, lmax, kmax):
    import beamconv
    blm2 = np.zeros((blm.shape[0], hp.Alm.getsize(lmax=lmax)), dtype=np.complex128)
    blm2[:,:blm.shape[1]] = blm
    blmm, blmp = beamconv.tools.eb2spin(blm2[1],blm2[2])
    blm2[1] = blmm
    blm2[2] = blmp
    np.save(beam_file, blm2)

# set up sky alm (monopole only)
rng = np.random.default_rng(np.random.SeedSequence(42))
lmax = 30
kmax=2

np.random.seed(10)
slm =make_full_random_alm(lmax, lmax, rng)

# completely random Mueller matrix
mueller = np.random.uniform(-1,1,size=(4,4))


# completely random beam
blm = make_full_random_alm(lmax, kmax, rng)

# FIXME if I don't do this, I see discrepancies. Something is still odd
# with the spin +-2 components
blm[1:3] = hp.blm_gauss(0., lmax, pol=True)[1:3]

save_blm_for_beamconv(blm, "temp_beam.npy", lmax, kmax)

nptg=100
# completely random pointings
ptg = np.empty((nptg,3))
ptg[:,0]=np.random.uniform(0,np.pi,size=(nptg,))    # theta
ptg[:,1]=np.random.uniform(0,2*np.pi,size=(nptg,))  # phi
ptg[:,2]=np.random.uniform(0,2*np.pi,size=(nptg,))  # psi
hwp_angles = np.random.uniform(0,2*np.pi,size=(nptg,))  # alpha

# get the signal from beamconv

t0=time()
signal_beamconv = get_beamconv_values(lmax=lmax, kmax=kmax, slm=slm, ptg=ptg, hwp_angles=hwp_angles, mueller=mueller, beam_file="temp_beam.npy")
print("time for beamconv:", time()-t0)

# Now do the same thing with MuellerConvolver

# FIXME: I don't understand this angle transformation.
# It feels as if we are missing a conjugation somewhere, but I don't know where.
ptg[:,2] = np.pi-ptg[:,2]

t0=time()
fullconv = mueller_convolver.MuellerConvolver(
    lmax,
    kmax,
    slm,
    blm,
    mueller,
    single_precision=False,
    epsilon=1e-7,
    ofactor=1.8,
    nthreads=1,
)
signal_muellerconvolver = fullconv.signal(ptg, hwp_angles)
print("time for MuellerConvolver:", time()-t0)

# L2 error
print("L2 error between results:", ducc0.misc.l2error(signal_beamconv, signal_muellerconvolver))
plt.plot(signal_beamconv)
plt.plot(signal_muellerconvolver)
plt.show()

[AdriJD]

I finally found some time to look at this.

• I think the np.pi-psi thing is just due to qpoint following the COSMO convention. I would keep psi unchanged for the MuellerConvolver call and just add psi = 180 - np.degrees(ptg[:,2]) inside of the get_beamconv_values function.
• Updating line 4148 makes sense. The line (hwp_spin[0,2] * np.sqrt(2)) that was in the code is indeed weird. Perhaps a left-over from when I was assuming symmetric Mueller matrices. The s0a2 term describes the coupling of the Stokes I sky to the instrument and thus it should not contain a term like hwp_spin[0,2] (which describes the coupling of Stokes Q and U to the instrument). If I update the s0a2 and s0a2_v terms as below:

        elif mode == 's0a2':
            blmm2, blmp2 = tools.shift_blm(blm[1], blm[2], 2, eb=False)                                                                                                                                                         
            blmm2 *= hwp_spin[2,0] * np.sqrt(2)
            blmp2 *= hwp_spin[1,0] * np.sqrt(2)

            return blmm2, blmp2

        elif mode == 's0a2_v':
            blmm2_v, blmp2_v = tools.shift_blm(blm[1], blm[2], 2, eb=False)                                                                                                                                                   
            blmm2_v *= hwp_spin[1,3] * np.sqrt(2) 
            blmp2_v *= hwp_spin[2,3] * np.sqrt(2)

The two codes agree for all Mueller matrices and all beams as long as the [1:3,0] and [1:3,-1] elements are zero. This means that codes do agree on all P->P leakage and all P->T leakage due to the HWP. So that is good news already. There is just something wrong with the T->P leakage.

If I make the following change:

        elif mode == 's0a2':
            blmm2, blmp2 = tools.shift_blm(blm[1], blm[2], 2, eb=False)                                                                                                                                                         
            blmm2 *= hwp_spin[2,0] * np.sqrt(2) * -1
            blmp2 *= hwp_spin[1,0] * np.sqrt(2)

            return blmm2, blmp2

        elif mode == 's0a2_v':
            blmm2_v, blmp2_v = tools.shift_blm(blm[1], blm[2], 2, eb=False)                                                                                                                                                   
            blmm2_v *= hwp_spin[1,3] * np.sqrt(2) * -1
            blmp2_v *= hwp_spin[2,3] * np.sqrt(2)

and use a symmetric beam for E and B, e.g. using the call to hp.blm_gauss from your script, the codes agree for all Mueller matrices. This is as far as I got. The minus signs in the above code block I don't fully understand yet. The fact something goes wrong only if there is an asymmetric component to the E and B beam components and only for the s0a2 and s0a2_v terms does narrow the search down.

[mreineck]

Thank you very much, @AdriJD!

I wasn't aware of the COSMO convention, this helps a lot.

Concerning the last of the discrepancies, I think I found a way to "fix" them. Not sure if this is the right fix, or if some of the adjustments need to go to MuellerConvolver instead, but I'll open a PR (just for comparison purposes) soon, so we can discuss this further.

I'm really, really happy that we probably have convergence now!

Here is the latest version of the comparison code, with a few more cleanups:

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import healpy as hp
import mueller_convolver
import ducc0

def nalm(lmax, mmax):
    return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)

def make_full_random_alm(lmax, mmax, rng):
    res = rng.uniform(-1., 1., (4, nalm(lmax, mmax))) \
     + 1j*rng.uniform(-1., 1., (4, nalm(lmax, mmax)))
    # make a_lm with m==0 real-valued
    res[:, 0:lmax+1].imag = 0.
    ofs=0
    # components 1 and 2 are spin-2, fix them accordingly
    spin=2
    for s in range(spin):
        res[1:3, ofs:ofs+spin-s] = 0.
        ofs += lmax+1-s
    return res

# code by Marta to get beamconv results for user-specified angles
def get_beamconv_values(lmax, kmax, slm, blm, ptg, hwp_angles, mueller):
    import beamconv
    import qpoint as qp

    # prepare PO beam file
    blm2 = np.zeros((blm.shape[0], hp.Alm.getsize(lmax=lmax)), dtype=np.complex128)
    blm2[:,:blm.shape[1]] = blm
    blmm, blmp = beamconv.tools.eb2spin(blm2[1],blm2[2])
    blm2[1] = blmm
    blm2[2] = blmp
    np.save("temp_beam.npy", blm2)

    # set up beam and HWP mueller matrix (identity, i.e. no HWP)
    beam = beamconv.Beam(btype='PO', lmax=lmax, mmax=lmax, deconv_q=True, normalize=False, po_file="temp_beam.npy", hwp_mueller=mueller)

    nsamp = ptg.shape[0]

    # from (theta,phi) to (ra,dec) convention
    # also, all angles are converted in degrees
    ra = np.degrees(ptg[:,1])
    dec = 90. - np.degrees(ptg[:,0])
    # Adjustment for difference in convention between qpoint and MuellerConvolver?
    psi = 180. - np.degrees(ptg[:,2])

    # calculate the quaternion
    q_bore_array = qp.QPoint().radecpa2quat(ra, dec, psi)

    def ctime_test(**kwargs):
        return np.zeros(kwargs.pop('end')-kwargs.pop('start'))
    
    def q_bore_test(**kwargs):
        return q_bore_array[kwargs.pop('start'):kwargs.pop('end')]

    S = beamconv.ScanStrategy(duration=nsamp, sample_rate=1, external_pointing=True, use_l2_scan=False)
    S.add_to_focal_plane(beam, combine=False)
    S.set_hwp_mod(mode='stepped', freq=1, angles=hwp_angles*180/np.pi)

    # determine nside_spin necessary for good accuracy
    nside_spin = 1
    while nside_spin < 4*lmax:
        nside_spin *= 2

    S.scan_instrument_mpi(slm.copy(), save_tod=True, ctime_func=ctime_test, q_bore_func=q_bore_test,
                      ctime_kwargs={'useless':0}, q_bore_kwargs={'useless':0},nside_spin=nside_spin, interp=True, input_v=True, beam_v=True, max_spin=kmax+4, binning=False, verbose=0)

    return S.data(S.chunks[0], beam=beam, data_type='tod').copy()


np.random.seed(10)
rng = np.random.default_rng(np.random.SeedSequence(42))
lmax = 30
kmax = 18

# completely random sky
slm =make_full_random_alm(lmax, lmax, rng)

# completely random Mueller matrix
mueller = np.random.uniform(-1,1,size=(4,4))

# completely random beam
blm = make_full_random_alm(lmax, kmax, rng)

nptg=100
# completely random pointings
ptg = np.empty((nptg,3))
ptg[:,0]=np.random.uniform(0,np.pi,size=(nptg,))    # theta
ptg[:,1]=np.random.uniform(0,2*np.pi,size=(nptg,))  # phi
ptg[:,2]=np.random.uniform(0,2*np.pi,size=(nptg,))  # psi
hwp_angles = np.random.uniform(0,2*np.pi,size=(nptg,))  # alpha

# get the signal from beamconv
signal_beamconv = get_beamconv_values(lmax=lmax, kmax=kmax, slm=slm, blm=blm, ptg=ptg, hwp_angles=hwp_angles, mueller=mueller)

# Now do the same thing with MuellerConvolver
fullconv = mueller_convolver.MuellerConvolver(
    lmax,
    kmax,
    slm,
    blm,
    mueller,
    single_precision=False,
    epsilon=1e-7,
    ofactor=1.8,
    nthreads=1,
)
signal_muellerconvolver = fullconv.signal(ptg, hwp_angles)

# L2 error
print("L2 error between results:", ducc0.misc.l2error(signal_beamconv, signal_muellerconvolver))
plt.plot(signal_beamconv)
plt.plot(signal_muellerconvolver)
plt.show()

[mreineck]

See #21

[mreineck]

Ah, maybe there is a clue: line 4040 of instrument.py reads

            blmE, blmB = tools.spin2eb(blmp2, blmm2)

while in most other places the argument order is blmm2, blmp2.
(Analogously in line 4047.)
There is an exception in the s2a4 code, but you explicitly comment on this ... so perhaps the above is an oversight?

[martamonelli]

Another test that might help. The code below plots beamconv's TOD together with the analytical expectation. To make them match, I need to flip the sign of mueller[0,2] and mueller[2,0] in the formulae. (This is with instrument.py taken from Martin's pull request.)

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import healpy as hp

from beamconv import Beam, ScanStrategy, tools
import qpoint as qp

def ctime_test(**kwargs):
    start = kwargs.pop('start')
    end = kwargs.pop('end')
    nsamp = end - start
    ctime0 = 0
    ctime = ctime0 + sample_rate*np.arange(nsamp)
    return ctime
    
def q_bore_test(**kwargs):
    start = kwargs.pop('start')
    end = kwargs.pop('end')
    return q_bore_array[start:end]

def rot_mat(theta):
    c, s = np.cos(2*theta), np.sin(2*theta)
    return np.array(((1,0,0),(0,c,s),(0,-s,c)))

##################################################################
# DEFINITIONS
##################################################################

#fixing random seed
np.random.seed(50)

#setting up sky alm
nside = 256
lmax = 2*nside

cls = np.loadtxt('ancillary/wmap7_r0p03_lensed_uK_ext.txt', unpack=True) #Cls in uK^2
ells, cls = cls[0], cls[1:]

slm = hp.synalm(cls, lmax=lmax, new=True)
input_maps = hp.alm2map(slm,nside=nside)

#setting up constant pointings
nptg=1000

ptg = np.empty((nptg,3))
nsamp = ptg.shape[0]
duration = 5                  #scan duration in seconds
sample_rate = nsamp/duration  #samples per second

#setting constant theta, phi, psi
ptg[:,0]=np.ones(nptg)*np.random.rand(1)*np.pi    # theta
ptg[:,1]=np.ones(nptg)*np.random.rand(1)*2*np.pi  # phi
ptg[:,2]=np.ones(nptg)*np.random.rand(1)*2*np.pi  # psi

pixs = hp.ang2pix(nside, ptg[:,0], ptg[:,1])

#from (theta,phi) to (ra,dec) convention
#also, all angles are converted in degrees
ra = np.degrees(ptg[:,1])
dec = 90. - np.degrees(ptg[:,0])
psi = np.degrees(ptg[:,2])

#calculate the boresight quaternion
q_bore_array = qp.QPoint().radecpa2quat(ra, dec, psi)

#setting up ideal beam and fixing HWP rotation frequency
fwhm_arcmin=0
hwp_freq = 2.001

##################################################################
# RUNNING BEAMCONV
##################################################################

#setting up beam and HWP mueller matrix (random fluctuations on top of ideal HWP)
mueller = np.diag([1.,1.,-1.,-1.])
mueller[0,0] += np.random.rand(1)*1e-2
mueller[0,1] += np.random.rand(1)*1e-2
mueller[0,2] += np.random.rand(1)*1e-2
mueller[1,0] += np.random.rand(1)*1e-2
mueller[1,1] += np.random.rand(1)*1e-2
mueller[1,2] += np.random.rand(1)*1e-2
mueller[2,0] += np.random.rand(1)*1e-2
mueller[2,1] += np.random.rand(1)*1e-2
mueller[2,2] += np.random.rand(1)*1e-2

beam = Beam(btype='Gaussian', fwhm=fwhm_arcmin, lmax=lmax)
beam.hwp_mueller = mueller

S = ScanStrategy(duration=duration, sample_rate=sample_rate, external_pointing=True, use_l2_scan=False)
S.add_to_focal_plane(beam, combine=False)
S.set_hwp_mod(mode='continuous', freq=hwp_freq, start_ang=0)
S.allocate_maps(nside=nside)
S.scan_instrument_mpi(slm, save_tod=True, ctime_func=ctime_test, q_bore_func=q_bore_test,
                  ctime_kwargs={'useless':0}, q_bore_kwargs={'useless':0})
     
tod2 = S.data(S.chunks[0], beam=beam, data_type='tod').copy()

##################################################################
# ANALYTICAL COMPUTATIONS
##################################################################

d2 = np.empty_like(tod2)

#stokes vector at the observed pixel
stokes = input_maps[:,pixs[0]]

#mimicking HWP's continuous rotation
start_ang = 0
hwp_angles = np.linspace(start_ang,
                   start_ang + 2 * np.pi * nsamp / (sample_rate / float(hwp_freq)),
                   num=nsamp, endpoint=False, dtype=float) # in radians (w = 2 pi freq)

#setting detector angle
xi = np.radians(0)

#redefining mueller matrix elements in order for the TODs to overlap
mueller_again = mueller
mueller_again[2,0] *= -1
mueller_again[0,2] *= -1

for i in np.arange(nptg):
    R_psi = rot_mat(-ptg[i,2])       #MINUS SIGN
    R_phi = rot_mat(-hwp_angles[i])  #MINUS SIGN
    R_xi = rot_mat(xi)
    stokes2 = np.dot(R_xi,np.dot(R_phi.transpose(),np.dot(mueller_again[:3,:3],np.dot(R_phi,np.dot(R_psi,stokes)))))
    d2[i] = stokes2[0]+stokes2[1]    #MISSING 1/2
    
#plotting everything    
#solid lines: analytical prediction (with tod_man on top);
#dashed lines: beamconv's output TOD
plt.plot(d2[0:150], color='teal', label='analytic formula')
plt.plot(tod2[0:150], linestyle='--', color='powderblue', label='beamconv')
plt.legend(loc='lower right')
plt.show()

[mreineck]

Interesting: if you flip the sign of the [1,2] and [2,1] matrix entries as well, the difference between the results drops again by many orders of magnitude.

So perhaps this is related to the "we multiply everything with spin !=0 by -1" convention?

[mreineck]

So perhaps this is related to the "we multiply everything with spin !=0 by -1" convention?

Sorry, after thinking about it for a few seconds, I guess this sentence doesn't make any sense.

Still, I think that [1,2] and [2,1] must also be flipped to get the best agreement.

[martamonelli]

By including V too, it looks like that also [2,3] and [3,2] should be flipped.

#redefining mueller matrix elements in order for the TODs to overlap
mueller_again = mueller
mueller_again[2,:] *= -1
mueller_again[:,2] *= -1