usage.py 2.1 KB
 Martin Reinecke committed Apr 16, 2020 1 2 ``````# Short usage demo for the interpol_ng module `````` Martin Reinecke committed Apr 14, 2020 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 ``````import pyinterpol_ng import numpy as np # establish a random number generator rng = np.random.default_rng(np.random.SeedSequence(42)) def nalm(lmax, mmax): return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax) def random_alm(lmax, mmax, ncomp): res = rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \ + 1j*rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) # make a_lm with m==0 real-valued res[0:lmax+1,:].imag = 0. return res # simulation parameters lmax=1024 # highest l and m moment for sky, highest l moment for beam kmax=13 # highest m moment for beam ncomp=3 # T, E and B # get random sky a_lm # the a_lm arrays follow the same conventions as those in healpy slm = random_alm(lmax, lmax, ncomp) # build beam a_lm blm = random_alm(lmax, kmax, ncomp) # build random pointings npnt = 1000000 # play with this to measure performance ptg = rng.uniform(0., 1., (npnt,3)) ptg[:,0]*=np.pi # theta ptg[:,1]*=2*np.pi # phi ptg[:,2]*=2*np.pi # psi # build the interpolator # For a "classic" CMB experiment we want to combine the T, E and B signals, # so we set `separate` to `False` print("classic interpolator setup...") inter_classic = pyinterpol_ng.PyInterpolator( slm,blm,separate=False,lmax=lmax, kmax=kmax, epsilon=1e-4, nthreads=2) print("...done") # get interpolated values print("interpolating...") res = inter_classic.interpol(ptg) print("...done") # res is an array of shape (nptg, 1). # If we had specified `separate=True`, it would be of shape(nptg, 3). print(res.shape) `````` Martin Reinecke committed Apr 23, 2020 60 61 ``````# Since the interpolator object holds large data structures, it should be # deleted once it is no longer needed `````` Martin Reinecke committed Apr 14, 2020 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 ``````del inter_classic # Now the same thing for an experiment with HWP. In this case we need the # interpolated T, E and B signals separate. separate = True print("HWP interpolator setup...") inter_hwp = pyinterpol_ng.PyInterpolator( slm,blm,separate=True,lmax=lmax, kmax=kmax, epsilon=1e-4, nthreads=2) print("...done") # get interpolated values print("interpolating...") res = inter_hwp.interpol(ptg) print("...done") # now res has shape(nptg,3) print(res.shape)``````