diff --git a/beamconv/instrument.py b/beamconv/instrument.py index cac106a..96ade1d 100644 --- a/beamconv/instrument.py +++ b/beamconv/instrument.py @@ -1811,11 +1811,15 @@ def scan_instrument_mpi(self, alm, ground_alm=None, verbose=1, binning=True, # nside_spin from init_spinmaps kwargs (if not None). max_spin = kwargs.pop('max_spin', None) nside_spin = kwargs.pop('nside_spin', None) + mu_con_hwp = kwargs.pop('mu_con_hwp', True) + mu_con_spin = kwargs.pop('mu_con_spin', True) spinmaps_opts = {} if max_spin: spinmaps_opts.update({'max_spin' : max_spin}) if nside_spin: spinmaps_opts.update({'nside_spin' : nside_spin}) + spinmaps_opts.update({'mu_con_hwp' : mu_con_hwp}) + spinmaps_opts.update({'mu_con_spin' : mu_con_spin}) if verbose and self.mpi_rank == 0: print('Scanning with {:d} detectors'.format( @@ -3799,7 +3803,7 @@ def _scan_modulate_pa(tod, pix, pa, maps, s_vals, def init_spinmaps(self, alm, beam, ground_alm=None, max_spin=5, nside_spin=256, - symmetric=False, input_v=False, beam_v=False): + symmetric=False, input_v=False, beam_v=False, mu_con_hwp=True, mu_con_spin=True): ''' Compute appropriate spinmaps for beam and all its ghosts. @@ -3827,6 +3831,12 @@ def init_spinmaps(self, alm, beam, ground_alm=None, max_spin=5, nside_spin=256, include the 4th alm component beam_v : bool include the 4th blm component + mu_con_hwp : bool + Switch between the Mueller convolver matching + HWP tweak and the old (pre-July 2023) behavior + mu_con_spin : bool + Switch between the Mueller convolver matching + spin +-2 tweak and the old (pre-July 2023) behavior Notes @@ -3852,7 +3862,7 @@ def init_spinmaps(self, alm, beam, ground_alm=None, max_spin=5, nside_spin=256, spinmap_dict = self._init_spinmaps(alm, blm, max_s, nside_spin, symmetric=beam.symmetric, hwp_mueller=beam.hwp_mueller, input_v=input_v, - beam_v=beam_v) + beam_v=beam_v, mu_con_hwp=mu_con_hwp, mu_con_spin=mu_con_spin) # Names: s0a0, s0a2, s2a0, s2a2, s2a4. # s refers to spin value under psi, a to spin value under HWP rot. @@ -3886,20 +3896,20 @@ def init_spinmaps(self, alm, beam, ground_alm=None, max_spin=5, nside_spin=256, spinmap_dict = self._init_spinmaps(alm, blm, max_s, nside_spin, symmetric=ghost.symmetric, hwp_mueller=ghost.hwp_mueller, - input_v=input_v, beam_v=beam_v) + input_v=input_v, beam_v=beam_v, mu_con_hwp=mu_con_hwp, mu_con_spin=mu_con_spin) self.spinmaps['ghosts'][u] = spinmap_dict if ground_alm is not None: self.spinmaps['ground'] = self._init_spinmaps(ground_alm, blm, max_s, nside_spin, symmetric=beam.symmetric, hwp_mueller=beam.hwp_mueller, input_v=input_v, - beam_v=beam_v) + beam_v=beam_v, mu_con_hwp=mu_con_hwp, mu_con_spin=mu_con_spin) @staticmethod def _init_spinmaps(alm, blm, max_spin, nside, symmetric=False, hwp_mueller=None, - input_v=False, beam_v=False): + input_v=False, beam_v=False, mu_con_hwp=True, mu_con_spin=True): ''' Compute convolution of map with different spin modes of the beam. @@ -3928,6 +3938,13 @@ def _init_spinmaps(alm, blm, max_spin, nside, beam_v : bool include the 4th blm component if it exists + mu_con_hwp : bool + Switch between the Mueller convolver matching + HWP tweak and the old (pre-July 2023) behavior + mu_con_spin : bool + Switch between the Mueller convolver matching + spin +-2 tweak and the old (pre-July 2023) behavior + returns ------- @@ -4036,15 +4053,21 @@ def _init_spinmaps(alm, blm, max_spin, nside, # s0a2. spinmap_dict['s0a2'] = {} - blmm2, blmp2 = ScanStrategy.blmxhwp(blm, hwp_spin, 's0a2') - blmE, blmB = tools.spin2eb(blmp2, blmm2) + blmm2, blmp2 = ScanStrategy.blmxhwp(blm, hwp_spin, 's0a2', mu_con_hwp=mu_con_hwp) + if mu_con_spin==False: + blmE, blmB = tools.spin2eb(blmp2, blmm2) + else: + blmE, blmB = tools.spin2eb(blmm2, blmp2) # 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) if input_v: - blmm2, blmp2 = ScanStrategy.blmxhwp(blm, hwp_spin, 's0a2_v') - blmE, blmB = tools.spin2eb(blmp2, blmm2) + blmm2, blmp2 = ScanStrategy.blmxhwp(blm, hwp_spin, 's0a2_v', mu_con_hwp=mu_con_hwp) + if mu_con_spin==False: + blmE, blmB = tools.spin2eb(blmp2, blmm2) + else: + blmE, blmB = tools.spin2eb(blmm2, blmp2) spinmap_dict['s0a2']['maps'] += ScanStrategy._spinmaps_complex( -alm[3], alm[3] * 0, blmE, blmB, spin_values_s0a2, nside) spinmap_dict['s0a2']['s_vals'] = spin_values_s0a2 @@ -4090,7 +4113,7 @@ def _init_spinmaps(alm, blm, max_spin, nside, return spinmap_dict @staticmethod - def blmxhwp(blm, hwp_spin, mode, beam_v=False): + def blmxhwp(blm, hwp_spin, mode, beam_v=False, mu_con_hwp=True): ''' Arguments --------- @@ -4103,6 +4126,9 @@ def blmxhwp(blm, hwp_spin, mode, beam_v=False): Pick between 's0a0', 's0a0_v', s2a4', 's0a2', 's0a2_v', 's2a2' or 's2a0'. beam_v : bool include the 4th blm component when constructing the spinmaps + mu_con_hwp : bool + Switch between the Mueller convolver matching + HWP tweak and the old (pre-July 2023) behavior Returns ------- @@ -4145,13 +4171,19 @@ def blmxhwp(blm, hwp_spin, mode, beam_v=False): elif mode == 's0a2': blmm2, blmp2 = tools.shift_blm(blm[1], blm[2], 2, eb=False) - blmm2 *= hwp_spin[0,2] * np.sqrt(2) + if mu_con_hwp == False: + blmm2 *= hwp_spin[0,2] * np.sqrt(2) + else: + blmm2 *= hwp_spin[1,0] * np.sqrt(2) blmp2 *= hwp_spin[2,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[3,2] * np.sqrt(2) + if mu_con_hwp == False: + blmm2_v *= hwp_spin[3,2] * np.sqrt(2) + else: + blmm2_v *= hwp_spin[1,3] * np.sqrt(2) blmp2_v *= hwp_spin[2,3] * np.sqrt(2) return blmm2_v, blmp2_v @@ -4235,7 +4267,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. @@ -4318,8 +4350,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/comparison_script.py b/comparison_script.py new file mode 100644 index 0000000..923acf2 --- /dev/null +++ b/comparison_script.py @@ -0,0 +1,180 @@ +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 + +def get_gauss_beam_from_beamconv(fwhm, lmax): + import beamconv + blmT, blmm2 = beamconv.tools.gauss_blm(fwhm*180*60/np.pi, lmax, pol=True) + res = np.zeros((4,blmT.shape[0]), dtype=np.complex128) + blmE, blmB = beamconv.tools.spin2eb(blmm2, blmm2*0, spin=2) + res[0] = blmT + res[1] = blmE + res[2] = blmB + res[3] = blmT # correct? + return res + +def blm_gauss_new(fwhm, lmax, pol=False): + fwhm = float(fwhm) + lmax = int(lmax) + pol = bool(pol) + mmax = 2 if pol else 0 + ncomp = 3 if pol else 1 + nval = hp.Alm.getsize(lmax, mmax) + + if mmax > lmax: + raise ValueError("lmax value too small") + + blm = np.zeros((ncomp, nval), dtype=np.complex128) + sigmasq = fwhm * fwhm / (8 * np.log(2.0)) + + for l in range(lmax+1): + blm[0, hp.Alm.getidx(lmax, l, 0)] = np.exp(-0.5*sigmasq*l*(l+1)) + + if pol: + for l in range(2, lmax+1): + blm[1, hp.Alm.getidx(lmax, l, 2)] = np.exp(-0.5 * sigmasq * (l*(l+1)-4)) + blm[2] = 1j * blm[1] + + return blm + +# blm_gauss_new times sqrt((2*l+1)/(4pi)) +def Blm_gauss_new(fwhm, lmax, pol=False): + blm = blm_gauss_new(fwhm, lmax, pol) + for l in range(lmax+1): + blm[0, hp.Alm.getidx(lmax, l, 0)] *= np.sqrt((2*l+1) / (4*np.pi)) + + if pol: + for l in range(2, lmax+1): + blm[1:3, hp.Alm.getidx(lmax, l, 2)] *= np.sqrt((2*l+1) / (4*np.pi)) + + return blm + +# code by Marta to get beamconv results for user-specified angles +def get_beamconv_values(lmax, kmax, slm, blm, ptg, hwp_angles, mueller, + mu_con_hwp, mu_con_spin): + 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 to 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, mu_con_hwp=mu_con_hwp, mu_con_spin=mu_con_spin) + + 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)) +#mueller[1:3,0]=mueller[1:3,-1] = 0 +#mueller[0,2]=mueller[2,0] = 0 + + +# completely random beam +blm = make_full_random_alm(lmax, kmax, rng) +# ... or use a Gauss beam +#blmtmp = blm_gauss_new(np.radians(10.), lmax, True) +#blm *= 0 +#blm [0:3, 0:blmtmp.shape[1]] = blmtmp + +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_FF = get_beamconv_values(lmax=lmax, kmax=kmax, slm=slm, blm=blm, ptg=ptg, hwp_angles=hwp_angles, mueller=mueller, mu_con_hwp=False, mu_con_spin=False) +signal_beamconv_FT = get_beamconv_values(lmax=lmax, kmax=kmax, slm=slm, blm=blm, ptg=ptg, hwp_angles=hwp_angles, mueller=mueller, mu_con_hwp=False, mu_con_spin=True) +signal_beamconv_TF = get_beamconv_values(lmax=lmax, kmax=kmax, slm=slm, blm=blm, ptg=ptg, hwp_angles=hwp_angles, mueller=mueller, mu_con_hwp=True, mu_con_spin=False) +signal_beamconv_TT = get_beamconv_values(lmax=lmax, kmax=kmax, slm=slm, blm=blm, ptg=ptg, hwp_angles=hwp_angles, mueller=mueller, mu_con_hwp=True, mu_con_spin=True) + +# Now do the same thing with MuellerConvolver +fullconv = mueller_convolver.MuellerConvolver( + lmax=lmax, + kmax=kmax, + slm=slm, + blm=blm, + mueller=mueller, + single_precision=False, + epsilon=1e-7, + nthreads=1, +) +signal_muellerconvolver = fullconv.signal(ptg=ptg, alpha=hwp_angles) + +# L2 error +print("L2 error to original beamconv:", ducc0.misc.l2error(signal_beamconv_FF, signal_muellerconvolver)) +print("L2 error to beamconv + HWP:", ducc0.misc.l2error(signal_beamconv_TF, signal_muellerconvolver)) +print("L2 error to beamconv + spin:", ducc0.misc.l2error(signal_beamconv_FT, signal_muellerconvolver)) +print("L2 error to beamconv + HWP + spin:", ducc0.misc.l2error(signal_beamconv_TT, signal_muellerconvolver)) +plt.plot(signal_beamconv_FF, label="beamconv, original") +plt.plot(signal_beamconv_TF, label="beamconv + HWP") +plt.plot(signal_beamconv_FT, label="beamconv + spin") +plt.plot(signal_beamconv_TT, label="beamconv + HWP + spin") +plt.plot(signal_muellerconvolver, label="MuellerConvolver") +plt.legend() +plt.show() diff --git a/mueller_convolver.py b/mueller_convolver.py new file mode 100644 index 0000000..24dd871 --- /dev/null +++ b/mueller_convolver.py @@ -0,0 +1,313 @@ +import ducc0 +import numpy as np + + +__all__ = ["MuellerConvolver"] + + +def nalm(lmax, mmax): + return ((mmax+1) * (mmax+2)) // 2 + (mmax+1) * (lmax-mmax) + + +# Adri 2020 A25/A35 +def mueller_to_C(mueller): + T = np.zeros((4, 4), dtype=np.complex128) + T[0, 0] = T[3, 3] = 1.0 + T[1, 1] = T[2, 1] = 1.0 / np.sqrt(2.0) + T[1, 2] = 1j / np.sqrt(2.0) + T[2, 2] = -1j / np.sqrt(2.0) + return T.dot(mueller.dot(np.conj(T.T))) + + +def truncate_blm(inp, lmax, kmax, epsilon=1e-10): + limit = epsilon * np.max(np.abs(inp)) + out = [] + for i in range(len(inp)): + maxk = -1 + idx = 0 + for k in range(kmax + 1): + if np.max(np.abs(inp[i, :, idx:idx+lmax+1-k])) > limit: + maxk = k + idx += lmax+1-k + # print("component",i,"maxk=",maxk) + if maxk == -1: + out.append(None) + else: + out.append((inp[i, :, : nalm(lmax, maxk)].copy(), maxk)) + return out + + +class MuellerConvolver: + """Class for computing convolutions between arbitrary beams and skies in the + presence of an optical element with arbitrary Mueller matrix in front of the + detector. + + Most of the expressions in this code are derived from + Duivenvoorden et al. 2021, MNRAS 502, 4526 + (https://arxiv.org/abs/2012.10437) + + Parameters + ---------- + lmax : int + maximum l moment of the provided sky and beam a_lm + kmax : int + maximum m moment of the provided beam a_lm + slm : numpy.ndarray((n_comp, n_slm), dtype=complex) + input sky a_lm + ncomp can be 1, 3, or 4, for T, TEB, TEBV components, respectively. + The components have the a_lm format used by healpy + blm : numpy.ndarray((n_comp, n_blm), dtype=complex) + input beam a_lm + ncomp can be 1, 3, or 4, for T, TEB, TEBV components, respectively. + The components have the a_lm format used by healpy + mueller : np.ndarray((4,4), dtype=np.float64) + Mueller matrix of the optical element in front of the detector + single_precision : bool + if True, store internal data in single precision, else double precision + epsilon : float + desired accuracy for the interpolation; a typical value is 1e-4 + npoints : int + total number of irregularly spaced points you want to use this object for + (only used for performance fine-tuning) + sigma_min, sigma_max: float + minimum and maximum allowed oversampling factors + 1.2 <= sigma_min < sigma_max <= 2.5 + nthreads : int + the number of threads to use for computation + """ + + # Very simple class to store a_lm that allow negative m values + class AlmPM: + def __init__(self, lmax, mmax): + if lmax < 0 or mmax < 0 or lmax < mmax: + raise ValueError("bad parameters") + self._lmax, self._mmax = lmax, mmax + self._data = np.zeros((2*mmax+1, lmax+1), dtype=np.complex128) + + def __getitem__(self, lm): + l, m = lm + + if isinstance(l, slice): + if l.step is not None or l.start < 0 or l.stop-1 > self._lmax: + print(l, m) + raise ValueError("out of bounds read access") + else: + if l < 0 or l > self._lmax: # or abs(m) > l: + print(l, m) + raise ValueError("out of bounds read access") + # if we are asked for elements outside our m range, return 0 + if m < -self._mmax or m > self._mmax: + return 0.0 + 0j + return self._data[m + self._mmax, l] + + def __setitem__(self, lm, val): + l, m = lm + if isinstance(l, slice): + if l.step is not None or l.start < 0 or l.stop-1 > self._lmax: + print(l, m) + raise ValueError("out of bounds write access") + else: + if l < 0 or l > self._lmax or abs(m) > l or abs(m) > self._mmax: + print(l, m) + raise ValueError("out of bounds write access") + self._data[m+self._mmax, l] = val + + def mueller_tc_prep(self, blm, mueller, lmax, mmax): + ncomp = blm.shape[0] + + # convert input blm to T/P/P*/V blm + blm2 = [self.AlmPM(lmax, mmax+4) for _ in range(4)] + idx = 0 + for m in range(mmax + 1): + sign = (-1)**m + lrange = slice(m, lmax+1) + idxrange = slice(idx, idx+lmax+1-m) + # T component + blm2[0][lrange, m] = blm[0, idxrange] + blm2[0][lrange, -m] = np.conj(blm[0, idxrange]) * sign + # V component + if ncomp > 3: + blm2[3][lrange, m] = blm[3, idxrange] + blm2[3][lrange, -m] = np.conj(blm[3, idxrange]) * sign + # E/B components + if ncomp > 2: + # Adri's notes [10], spin +2 + blm2[1][lrange, m] = -blm[1, idxrange] - 1j*blm[2, idxrange] + # Adri's notes [9], spin -2 + blm2[2][lrange, m] = -blm[1, idxrange] + 1j*blm[2, idxrange] + # negative m + # Adri's notes [2] + blm2[1][lrange, -m] = np.conj(blm2[2][lrange, m]) * sign + blm2[2][lrange, -m] = np.conj(blm2[1][lrange, m]) * sign + idx += lmax+1-m + + C = mueller_to_C(mueller) + + # compute the blm for the full beam+Mueller matrix system at angles + # n*pi/5 for n in [0; 5[ + sqrt2 = np.sqrt(2.0) + nbeam = 5 + inc = 4 + res = np.zeros((nbeam, ncomp, nalm(lmax, mmax+inc)), dtype=self._ctype) + blm_eff = [self.AlmPM(lmax, mmax+4) for _ in range(4)] + + for ibeam in range(nbeam): + alpha = ibeam * np.pi / nbeam + e2ia = np.exp(2j*alpha) + e2iac = np.exp(-2j*alpha) + e4ia = np.exp(4j*alpha) + e4iac = np.exp(-4j*alpha) + # FIXME: do I need to calculate anything for negative m? + for m in range(-mmax-4, mmax+4+1): + lrange = slice(abs(m), lmax+1) + # T component, Marta notes [4a] + blm_eff[0][lrange, m] = ( + C[0, 0] * blm2[0][lrange, m] + + C[3, 0] * blm2[3][lrange, m] + + 1.0/sqrt2 * ((C[1, 0]*e2ia ) * blm2[2][lrange, m+2] + +(C[2, 0]*e2iac) * blm2[1][lrange, m-2])) + # E/B components, Marta notes [4b,c] + blm_eff[1][lrange, m] = ( + (sqrt2*e2iac) * (C[0, 1] * blm2[0][lrange, m+2] + + C[3, 1] * blm2[3][lrange, m+2]) + + (C[2, 1]*e4iac) * blm2[2][lrange, m+4] + + C[1, 1] * blm2[1][lrange, m]) + blm_eff[2][lrange, m] = ( + (sqrt2*e2ia) * (C[0, 2] * blm2[0][lrange, m-2] + + C[3, 2] * blm2[3][lrange, m-2]) + + (C[1, 2]*e4ia) * blm2[1][lrange, m-4] + + C[2, 2] * blm2[2][lrange, m]) + # V component, Marta notes [4d] + blm_eff[3][lrange, m] = ( + C[0, 3] * blm2[0][lrange, m] + + C[3, 3] * blm2[3][lrange, m] + + 1.0/sqrt2 * ((C[1, 3]*e2ia ) * blm2[2][lrange, m+2] + +(C[2, 3]*e2iac) * blm2[1][lrange, m-2])) + + # TEMPORARY sanity check ... + def check(xp, xm, lrange, m): + sign = (-1)**m + diff = xp[lrange, m] - sign*np.conj(xm[lrange, -m]) + return np.max(np.abs(diff)) <= 1e-4 + + for m in range(0, mmax+4+1): + sign = (-1)**m + lrange = slice(abs(m), lmax+1) + if not check(blm_eff[0], blm_eff[0], lrange, m): + raise RuntimeError("error T") + if not check(blm_eff[1], blm_eff[2], lrange, m): + raise RuntimeError("error 12") + if not check(blm_eff[2], blm_eff[1], lrange, m): + raise RuntimeError("error 21") + if not check(blm_eff[3], blm_eff[3], lrange, m): + raise RuntimeError("error V") + # ... up to here + + # back to original TEBV b_lm format + idx = 0 + for m in range(mmax+inc+1): + lrange = slice(m, lmax+1) + idxrange = slice(idx, idx+lmax+1-m) + # T component + res[ibeam, 0, idxrange] = blm_eff[0][lrange, m] + # P/P* components + if ncomp > 2: + # Adri's notes [10] + res[ibeam, 1, idxrange] = -0.5 * ( + blm_eff[1][lrange, m] + blm_eff[2][lrange, m] + ) + res[ibeam, 2, idxrange] = 0.5j * ( + blm_eff[1][lrange, m] - blm_eff[2][lrange, m] + ) + # V component + if ncomp > 3: + res[ibeam, 3, idxrange] = blm_eff[3][lrange, m] + idx += lmax+1-m + + return res + + # "Fourier transform" the blm at different alpha to obtain + # blm(alpha) = out[0] + cos(2*alpha)*out[1] + sin(2*alpha)*out[2] + # + cos(4*alpha)*out[3] + sin(4*alpha)*out[4] + def pseudo_fft(self, inp): + out = np.zeros((5, inp.shape[1], inp.shape[2]), dtype=self._ctype) + out[0] = 0.2 * (inp[0]+inp[1]+inp[2]+inp[3]+inp[4]) + # FIXME: I'm not absolutely sure about the sign of the angles yet + c1, s1 = np.cos(2*np.pi/5), np.sin(2*np.pi/5) + c2, s2 = np.cos(4*np.pi/5), np.sin(4*np.pi/5) + out[1] = 0.4 * (inp[0] + c1*(inp[1]+inp[4]) + c2*(inp[2]+inp[3])) + out[2] = 0.4 * ( s1*(inp[1]-inp[4]) + s2*(inp[2]-inp[3])) + out[3] = 0.4 * (inp[0] + c2*(inp[1]+inp[4]) + c1*(inp[2]+inp[3])) + out[4] = 0.4 * ( s2*(inp[1]-inp[4]) - s1*(inp[2]-inp[3])) + # Alternative way via real FFT + # out2 = inp.copy() + # out2 = out2.view(np.float64) + # out2 = ducc0.fft.r2r_fftpack( + # out2, real2hermitian=True, forward=False, axes=(0,), out=out2 + # ) + # out2[0] *=0.2 + # out2[1:] *=0.4 + # out2 = out2.view(np.complex128) + # print(np.max(np.abs(out-out2))) + return out + + def __init__(self, *, lmax, kmax, slm, blm, mueller, single_precision=True, + epsilon=1e-4, npoints=1000000000,sigma_min=1.2, sigma_max=2.5, + nthreads=1): + self._ftype = np.float32 if single_precision else np.float64 + self._ctype = np.complex64 if single_precision else np.complex128 + self._slm = slm.astype(self._ctype) + self._lmax = lmax + self._kmax = kmax + tmp = self.mueller_tc_prep(blm, mueller, lmax, kmax) + tmp = self.pseudo_fft(tmp) + + # construct the five interpolators for the individual components + # All sets of blm are checked up to which kmax they contain significant + # coefficients, and the interpolator is chosen accordingly + tmp = truncate_blm(tmp, self._lmax, self._kmax + 4) + + self._inter = [] + intertype = (ducc0.totalconvolve.Interpolator_f + if self._ctype == np.complex64 + else ducc0.totalconvolve.Interpolator) + for i in range(5): + if tmp[i] is not None: # component is not zero + self._inter.append( + intertype(self._slm, tmp[i][0], False, self._lmax, + tmp[i][1], epsilon=epsilon, npoints=npoints, + sigma_min=sigma_min, sigma_max=sigma_max, + nthreads=nthreads)) + else: # we can ignore this component entirely + self._inter.append(None) + + def signal(self, *, ptg, alpha): + """Computes the convolved signal for a set of pointings and HWP angles. + + Parameters + ---------- + ptg : numpy.ndarray((nptg, 3), dtype=float) + the input pointings in radians, in (theta, phi, psi) order + alpha : numpy.ndarray((nptg,), dtype=np.float) + the HWP angles in radians + + Returns + ------- + signal : numpy.ndarray((nptg,), dtype=np.float) + the signal measured by the detector + """ + ptg = ptg.astype(self._ftype) + alpha = alpha.astype(self._ftype) + if self._inter[0] is not None: + res = self._inter[0].interpol(ptg)[0] + else: + res = np.zeros(ptg.shape[0], dtype=self._ftype) + if self._inter[1] is not None: + res += np.cos(2*alpha) * self._inter[1].interpol(ptg)[0] + if self._inter[2] is not None: + res += np.sin(2*alpha) * self._inter[2].interpol(ptg)[0] + if self._inter[3] is not None: + res += np.cos(4*alpha) * self._inter[3].interpol(ptg)[0] + if self._inter[4] is not None: + res += np.sin(4*alpha) * self._inter[4].interpol(ptg)[0] + return res diff --git a/munich4.py b/munich4.py new file mode 100644 index 0000000..02f73df --- /dev/null +++ b/munich4.py @@ -0,0 +1,407 @@ +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt + +print('Loading modules...') +import os, copy +import numpy as np +import healpy as hp +import qpoint as qp +from beamconv import ScanStrategy, Beam, tools, plot_tools +print('..done') + + +matplotlib.rcParams.update({'font.size':10}) +matplotlib.rcParams.update({'xtick.direction':'in'}) +matplotlib.rcParams.update({'ytick.direction':'in'}) +matplotlib.rcParams.update({'xtick.top':True}) +matplotlib.rcParams.update({'ytick.right':True}) +matplotlib.rcParams.update({'legend.fontsize':8}) +lfs = 14 +hdl = 2 + +opj = os.path.join +nside = 128 +dpi = 300 + +cw = [np.array([78., 121., 165.])/255., +np.array([241., 143., 59.])/255., +np.array([224., 88., 91.])/255., +np.array([119., 183., 178.])/255., +np.array([90., 161., 85.])/255., +np.array([237., 201., 88.])/255., +np.array([175., 122., 160.])/255., +np.array([254., 158., 168.])/255., +np.array([156., 117., 97.])/255., +np.array([186., 176., 172.])/255.] + +M = qp.QMap() +M.init_dest(nside=nside, pol=False, reset=True) + + +# def get_cls(no_B=False, no_T=False): + +# wmap_dir = '/Users/jon/git_repos/beamconv/ancillary/' +# fname='wmap7_r0p03_lensed_uK_ext.txt' + +# cls = np.loadtxt(opj(wmap_dir, fname), unpack=True) # Cl in uK^2 +# cls[3] *= 1.25 +# print(cls[2][:50]) +# if no_B: +# print('NO B') +# data = np.loadtxt('lensing.txt', unpack=True) +# lmax = len(cls[0]) +# cls[3][1:] = data[1][:lmax-1] +# cls[3][0] = 0.0 + +# if no_T: +# cls[1][:] = 0.0 + +# print(cls[2][:50]) +# print(np.shape(cls)) + +# return cls[0], cls[1:] + +def short_ctime(**kwargs): + + start = kwargs.pop('start') + end = kwargs.pop('end') + cidx = kwargs.pop('cidx') + + ctime = np.linspace(0, 10., 1000) + + return ctime + +def get_theta(nsamp=1000): + + # return np.linspace(0, 4 * np.pi / 6., nsamp) + return np.pi / 2.0 * np.ones(nsamp) + +def get_phi(nsamp=1000): + + # return np.linspace(0, np.pi/6., nsamp) + return np.zeros(nsamp) + +def short_scan(nsamp=1000, **kwargs): + + theta = get_theta(nsamp) + phi = get_phi(nsamp) + + # theta = np.zeros(1000) + # phi = np.zeros(1000) + + vec = hp.ang2vec(theta, phi) + # Ra and Dec for qpoint + ra, dec = hp.vec2ang(vec, lonlat=True) + + plt.plot(theta) + plt.savefig(opj('img/', 'theta.png')) + plt.close() + + # pa = np.linspace(0, 360, 1000) + pa = np.zeros_like(phi) + + q_bore = M.radecpa2quat(ra, dec, pa) + + return q_bore + +def scan2(use_dust=True, lmax=300, tag='', + no_T=False, no_E=False, no_B=False, cmap='RdBu_r'): + + nside_spin = 128 + mmax = 4 + mlen = 10 + sample_rate = 100.0 + nside_out = 128 + verbose = 1 + ctime0 = 1510000000 + + # fwhm = 19.03863 + fwhm = 60. + + # ell, cls = get_cls(no_B=False) + + np.random.seed(42) + + # alm = hp.synalm(cls, lmax=lmax, new=True, verbose=True) # uK + # alm = np.asarray(alm) + + alm = np.zeros((3, int(lmax*(lmax+1)/2 + lmax +1))) + 0j + + idx1 = hp.Alm.getidx(lmax, 1, -1) + idx2 = hp.Alm.getidx(lmax, 1, 0) + idx3 = hp.Alm.getidx(lmax, 1, 1) + idx4 = hp.Alm.getidx(lmax, 2, -2) + idx5 = hp.Alm.getidx(lmax, 2, -1) + idx6 = hp.Alm.getidx(lmax, 2, 0) + idx7 = hp.Alm.getidx(lmax, 2, 1) + idx8 = hp.Alm.getidx(lmax, 2, 2) + + print('Idx 1: {}'.format(idx1)) + print('Idx 2: {}'.format(idx2)) + print('Idx 3: {}'.format(idx3)) + print('Idx 4: {}'.format(idx4)) + print('Idx 5: {}'.format(idx5)) + print('Idx 6: {}'.format(idx6)) + print('Idx 7: {}'.format(idx7)) + print('Idx 8: {}'.format(idx8)) + + alm[0][:] = 0.0 + alm[0][0] = 1.0 + alm[1][2] = 1.0 + alm[2][2] = 0.0 + alm[1][3:] = 0.0 + alm[2][3:] = 0.0 + + # print(np.shape(alm)) + # print(np.shape(alm2)) + + # print(alm[0][:]) + # print(alm2[0][:]) + # print(alm[0][0]) + # print(alm2[0][0]) + + + maps_raw = hp.alm2map(alm, 256) + + sym_limits = [250, 0.5, 0.5] + limits = [[0, 0], [-0.4, 0.0], [-0.4, 0]] + + + cmap = copy.copy(plt.cm.get_cmap(cmap)) + cmap.set_over(cmap(1.0), 1.0) + cmap.set_under(cmap(0.0), 1.0) + + + theta = np.linspace(0, 4 * np.pi / 6., 1000) + phi = np.linspace(0, np.pi/6., 1000) + + # theta = np.zeros(1000) + # phi = np.zeros(1000) + + nsamp =1000 + theta = get_theta(nsamp) + phi = get_phi(nsamp) + + vec = hp.ang2vec(theta, phi) + ra, dec = hp.vec2ang(vec, lonlat=True) + + plt.plot(theta) + plt.savefig(opj('img/', 'theta.png')) + plt.close() + + plt.plot(phi) + plt.savefig(opj('img/', 'phi.png')) + plt.close() + + plt.plot(ra) + plt.savefig(opj('img/', 'ra.png')) + plt.close() + + plt.plot(dec) + plt.savefig(opj('img/', 'dec.png')) + plt.close() + + print('Maximum value in map: {:.3f} uK'.format(1e6 * np.max(maps_raw[1].flatten()))) + + # plot_tools.plot_iqu(maps_raw, 'img/', 'input_spher'+tag, + # plot_func=projview, + # sym_limits=[250, 3, 3], cbar=False, #no_limits=True, + # # limits=[[-50, 200], [-10, 40], [-7, 7]], cbar=False, + # projection='nsper', lon_0=0, lat_0=40, pmin=0.5, pmax=99.5, + # cmap=cmap, bgcolor='white', unit='Signal [uK]', cbar_extend='both', + # diverging=False, dpi=300, transparent=True) + + print('setting up') + sat = ScanStrategy(mlen, external_pointing=True, sample_rate=sample_rate, + location='space', ctime0=ctime0, nside_out=nside_out) + + scan_opts = dict(q_bore_func=short_scan, + ctime_func=short_ctime, + q_bore_kwargs=dict(jitter_amp=0.), + ctime_kwargs=dict(), + max_spin=mmax, + nside_spin=nside_spin, + verbose=verbose, + binning=False, + interp=True) + + ############### Gaussian scanning + sat.create_focal_plane(nrow=1, ncol=1, fov=0, + fwhm=fwhm, lmax=lmax, mmax=mmax, amplitude=1., + no_pairs=True) + + ############### Continuous ideal HWP + hfreq = 0.1 + varphi = 0.0 + sat.set_hwp_mod(mode='continuous', freq=hfreq, varphi=varphi) + sat.create_focal_plane(nrow=1, ncol=1, fov=0, + fwhm=fwhm, lmax=lmax, mmax=mmax, amplitude=1., + no_pairs=True, combine=False) + + print('scanning with Gaussian w/ideal HWP') + #print(sat.hwp_dict) + + mueller = np.diag([1., 1., -1., -1.]) +# mueller = np.random.uniform(-1,1,size=(4,4)) + + # mueller = np.diag([0, 0, 0, 0.]) + for beami in sat.beams: + print(beami[0]) + print(type(beami[0])) + beami[0].hwp_mueller = mueller + beami[1].hwp_mueller = mueller + + sat.scan_instrument_mpi(alm, mu_con_hwp=True, mu_con_spin=True, **scan_opts) + gaussian_tod_whwp = sat.tod.copy() + hwp_ang = sat.hwp_ang.copy() + hwp_ang4 = 4 * 180 * hwp_ang / np.pi + hwp_ang4 = hwp_ang4 % 360 + + hwp_ang2 = 2 * 180 * hwp_ang / np.pi + hwp_ang2 = hwp_ang2 % 360 + + hwp_ang = 180 * hwp_ang / np.pi + + from scipy.signal import argrelextrema + idx2 = argrelextrema(hwp_ang2, np.less)[0] + idx4 = argrelextrema(hwp_ang4, np.less)[0] + + plt.plot(hwp_ang) + plt.plot(hwp_ang2) + plt.plot(hwp_ang4) + plt.savefig(opj('img/', 'hwp_ang.png')) + plt.close() + + sat.set_hwp_mod(mode='continuous', freq=hfreq, varphi=varphi) + sat.scan_instrument_mpi(alm, mu_con_hwp=True, mu_con_spin=False, **scan_opts) + gaussian_tod_whwp_nospin = sat.tod.copy() + + sat.set_hwp_mod(mode='continuous', freq=hfreq, varphi=varphi) + sat.scan_instrument_mpi(alm, mu_con_hwp=False, mu_con_spin=False, **scan_opts) + gaussian_tod_whwp_orig = sat.tod.copy() + + + #sat.scan_instrument_mpi(alm_noB, **scan_opts) + #gaussian_tod_whwp_noB = sat.tod.copy() + + ############## 1BR HWP + hwp_model='1BR' + sat.create_focal_plane(nrow=1, ncol=1, fov=0, + fwhm=fwhm, lmax=lmax, mmax=mmax, amplitude=1., + no_pairs=True, combine=False) + + for beami in sat.beams: + print(beami[0]) + print(type(beami[0])) + beami[0].set_hwp_mueller(model_name=hwp_model) + beami[1].set_hwp_mueller(model_name=hwp_model) + + print('scanning with Gaussian') + #print(sat.hwp_dict) + sat.set_hwp_mod(mode='continuous', freq=hfreq, varphi=varphi) + sat.scan_instrument_mpi(alm, mu_con_hwp=True, mu_con_spin=True, **scan_opts) + gaussian_tod_w1BR = sat.tod.copy() + + sat.set_hwp_mod(mode='continuous', freq=hfreq, varphi=varphi) + sat.scan_instrument_mpi(alm, mu_con_hwp=True, mu_con_spin=False, **scan_opts) + gaussian_tod_w1BR_nospin = sat.tod.copy() + + sat.set_hwp_mod(mode='continuous', freq=hfreq, varphi=varphi) + sat.scan_instrument_mpi(alm, mu_con_hwp=False, mu_con_spin=False, **scan_opts) + gaussian_tod_w1BR_orig = sat.tod.copy() + + ############## 3BR HWP + hwp_model='3BR' + sat.create_focal_plane(nrow=1, ncol=1, fov=0, + fwhm=fwhm, lmax=lmax, mmax=mmax, amplitude=1., + no_pairs=True, combine=False) + + for beami in sat.beams: + print(beami[0]) + print(type(beami[0])) + beami[0].set_hwp_mueller(model_name=hwp_model) + beami[1].set_hwp_mueller(model_name=hwp_model) + + print('scanning with Gaussian') + sat.set_hwp_mod(mode='continuous', freq=hfreq, varphi=varphi) + sat.scan_instrument_mpi(alm, mu_con_hwp=True, mu_con_spin=True, **scan_opts) + gaussian_tod_w3BR = sat.tod.copy() + + sat.set_hwp_mod(mode='continuous', freq=hfreq, varphi=varphi) + sat.scan_instrument_mpi(alm, mu_con_hwp=True, mu_con_spin=False, **scan_opts) + gaussian_tod_w3BR_nospin = sat.tod.copy() + + sat.set_hwp_mod(mode='continuous', freq=hfreq, varphi=varphi) + sat.scan_instrument_mpi(alm, mu_con_hwp=False, mu_con_spin=False, **scan_opts) + gaussian_tod_w3BR_orig = sat.tod.copy() + + ############## Plotting + fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(8, 6)) + + print('plotting') + + ax1.plot(gaussian_tod_whwp, color=cw[0], label='Gaussian') + ax1.plot(gaussian_tod_whwp_nospin, color=cw[0], ls='--', label='Gaussian (no spin)') + ax1.plot(gaussian_tod_whwp_orig, color=cw[0], ls=':', label='Gaussian (orig)') + ax1.plot(gaussian_tod_w1BR, color=cw[1], label='1BR') + ax1.plot(gaussian_tod_w1BR_nospin, color=cw[1], ls='--', label='1BR (no spin)') + ax1.plot(gaussian_tod_w1BR_orig, color=cw[1], ls=':', label='1BR (orig)') + ax1.plot(gaussian_tod_w3BR, color=cw[2], label='3BR') + ax1.plot(gaussian_tod_w3BR_nospin, color=cw[2], ls='--', label='3BR (no spin)') + ax1.plot(gaussian_tod_w3BR_orig, color=cw[2], ls=':', label='3BR (orig)') + + # print('RMS amplitude (Gaussian) {:.6f}'.format(np.std(gaussian_tod))) + # print('RMS amplitude (Gaussian noB) {:.6f}'.format(np.std(gaussian_tod_noB))) + # print('RMS amplitude (EG) {:.6f}'.format(np.std(eg_tod))) + # print('RMS amplitude (PO) {:.6f}'.format(np.std(po_tod))) + + ax2.plot(gaussian_tod_whwp - gaussian_tod_w1BR, color=cw[0], + label='Gaussian w/HWP - Gaussian w/HWP 1BR') + ax2.plot(gaussian_tod_whwp - gaussian_tod_w3BR, color=cw[1], + label='Gaussian w/HWP - Gaussian w/HWP 3BR') + + ax3.plot(gaussian_tod_whwp - gaussian_tod_whwp_nospin, color=cw[0], + label='Ideal: Both - No spin') + ax3.plot(gaussian_tod_whwp - gaussian_tod_whwp_orig, color=cw[0], ls='--', + label='Ideal: Both - Orig') + + ax3.plot(gaussian_tod_w1BR - gaussian_tod_w1BR_nospin, color=cw[1], + label='1BR: Both - No spin') + ax3.plot(gaussian_tod_w1BR - gaussian_tod_w1BR_orig, color=cw[1], ls='--', + label='1BR: Both - Orig') + + ax3.plot(gaussian_tod_w3BR - gaussian_tod_w3BR_nospin, color=cw[2], + label='3BR: Both - No spin') + ax3.plot(gaussian_tod_w3BR - gaussian_tod_w3BR_orig, color=cw[2], ls='--', + label='3BR: Both - Orig') + + for idx in idx4: + + ax3.axvline(x=idx, color='black', alpha=0.3) + + ax1.legend(borderaxespad=0.1, ncol=3, loc=2, + facecolor='white', edgecolor='None', framealpha=0, + handlelength=3, fontsize=8) + ax2.legend(borderaxespad=0.1, ncol=3, loc=2, + facecolor='white', edgecolor='None', framealpha=0, + handlelength=3, fontsize=8) + ax3.legend(borderaxespad=0.1, ncol=3, loc=2, + facecolor='white', edgecolor='None', framealpha=0, + handlelength=3, fontsize=8) + + ax1.set_ylabel('Signal [uK]') + ax2.set_ylabel('Signal [uK]') + ax3.set_ylabel('Signal [uK]') + ax2.set_xlabel('Sample number') + plt.savefig(opj('img/','hwp_tod.png'), + bbox_inches='tight', dpi=300) + plt.close() + + +if __name__ == '__main__': + + scan2(tag='munich', cmap='RdBu_r') + + +