Open mreineck opened 1 year ago
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.
I'm telling a lie ... hwp_spin[2,0]
is indeed accessed in blmxhwp
, but hwp_spin[1,0]
isn't.
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)))
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
.
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
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
instrument.py
to blmm2 *= hwp_spin[2,0] * np.sqrt(2)
instrument.py
to blmm2_v *= hwp_spin[3,2] * np.sqrt(2)
psi
angle to np.pi-psi
for MuellerConvolver
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()
I finally found some time to look at this.
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.
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()
See #21
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?
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()
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?
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.
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
If I see correctly, you also need to flip the sign of psi
and alpha
.
In the particular case of this test, the psi
flip is equivalent to the pi-psi
tweak that I need to make in the MuellerConvolver
comparison (your beam is invariant to rotation by pi
). But the alpha
flip and the changes to the Mueller matrix I can't explain ... so either I made a mistake (or bad assumption) when implementing your equations for MuellerConvolver
, or there is some inconsistency between them and your analytical expression.
One more data point: there seems to b a difference between what healpy
and beamconv
consider to be a Gaussian beam with FWHM=0: the E and B components differ by a factor of -2 in both implementations. Perhaps this explains some of the remaining discrepancies. (I don't know what the effective blm
of Marta's analytical formula are.)
Here is a code snippet which demonstrates the discrepancies. If I did anything wrong during the import of the healpy
beam to beamconv
, please let me know!
import numpy as np
import healpy as hp
import beamconv
lmax = 10
# generate Gaussian beam directly in beamconv
beam1 = beamconv.Beam(btype='Gaussian', fwhm=0, lmax=lmax)
blm1 = np.array(beam1.blm)
# generate Gaussian beam with healpy and import it to beamconv
tmp = hp.blm_gauss(0., lmax, pol=True)
# prepare PO beam file
tmp2 = np.zeros((3, hp.Alm.getsize(lmax=lmax)), dtype=np.complex128)
tmp2[:,:tmp.shape[1]] = tmp
# convert E/B to spin
blmm, blmp = beamconv.tools.eb2spin(tmp2[1],tmp2[2])
tmp2[1] = blmm
tmp2[2] = blmp
np.save("temp_beam.npy", tmp2)
beam2 = beamconv.Beam(btype='PO', lmax=lmax, mmax=lmax, deconv_q=True, normalize=False, po_file="temp_beam.npy")
blm2 = np.array(beam2.blm)
# NaNs are OK, but any actual number different from 1 indicates a problem
print(blm1/blm2)
Hi:) I believe it's the first time that I write here. I'm Marco Bortolami, PhD student in Ferrara, now visiting Okayama working with Yuya Nagano (@okayamayu8) on beam systematics, among other things. We compared an independent convolution code that he wrote in Julia (Condor) with ducc0/MuellerConvolver, see here (if you cannot access, tell me and I will copy the results here).
This could let us conclude that Mueller Convolver is correct and there are some small bugs on the implementation of the equations on beamconv, as Yuya implemented exactly the equations on Adri's notes.
However, we would like to discuss about equations (72) and (76) of Adri's notes, as we do not understand completely the m±4
components.
Thanks to the hard work by @martamonelli, I now have a very compact code that allows direct comparison between the
beamconv
andMuellerConvolver
HWP simulators.My current setup is:
lmax
)I can get very good agreement between the output of
beamconv
andMuellerConvolver
, ifpsi
angle forMuellerConvolver
(This seems to be a convention issue only.)blm_E
andblm_B
by-np.sqrt(2)
forMuellerConvolver
(I have no idea where this factor comes from)mueller[0,2] = -mueller[2,0]
The last point seems to indicate some asymmetry within
beamconv
. If I runbeamconv
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
beamconv
, I think.test.tar.gz