Commit bd30feab authored by Martin Reinecke's avatar Martin Reinecke

cleanup

parent ef3bd6c6
import ducc_0_1.totalconvolve as totalconvolve
import numpy as np
import ducc_0_1.sht as sht
import ducc_0_1.misc as misc
import time
import matplotlib.pyplot as plt
np.random.seed(48)
def nalm(lmax, mmax):
return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
def random_alm(lmax, mmax, ncomp):
res = np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
+ 1j*np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp))
# make a_lm with m==0 real-valued
res[0:lmax+1,:].imag = 0.
return res
def compress_alm(alm,lmax):
res = np.empty(2*len(alm)-lmax-1, dtype=np.float64)
res[0:lmax+1] = alm[0:lmax+1].real
res[lmax+1::2] = np.sqrt(2)*alm[lmax+1:].real
res[lmax+2::2] = np.sqrt(2)*alm[lmax+1:].imag
return res
def myalmdot(a1,a2,lmax,mmax,spin):
return np.vdot(compress_alm(a1,lmax),compress_alm(np.conj(a2),lmax))
def convolve(alm1, alm2, lmax):
job = sht.sharpjob_d()
job.set_triangular_alm_info(lmax, lmax)
job.set_gauss_geometry(lmax+1, 2*lmax+1)
map = job.alm2map(alm1)*job.alm2map(alm2)
job.set_triangular_alm_info(0,0)
return job.map2alm(map)[0]*np.sqrt(4*np.pi)
lmax=60
kmax=13
ncomp=1
separate=True
ncomp2 = ncomp if separate else 1
# 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)
t0=time.time()
# build interpolator object for slm and blm
foo = totalconvolve.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=1e-4, nthreads=2)
print("setup time: ",time.time()-t0)
nth = lmax+1
nph = 2*lmax+1
# compute a convolved map at a fixed psi and compare it to a map convolved
# "by hand"
ptg = np.zeros((nth,nph,3))
ptg[:,:,0] = (np.pi*(0.5+np.arange(nth))/nth).reshape((-1,1))
ptg[:,:,1] = (2*np.pi*(0.5+np.arange(nph))/nph).reshape((1,-1))
ptg[:,:,2] = np.pi*0.2
t0=time.time()
# do the actual interpolation
bar=foo.interpol(ptg.reshape((-1,3))).reshape((nth,nph,ncomp2))
print("interpolation time: ", time.time()-t0)
plt.subplot(2,2,1)
plt.imshow(bar[:,:,0])
bar2 = np.zeros((nth,nph))
blmfull = np.zeros(slm.shape)+0j
blmfull[0:blm.shape[0],:] = blm
for ith in range(nth):
rbeamth=misc.rotate_alm(blmfull[:,0], lmax, ptg[ith,0,2],ptg[ith,0,0],0)
for iph in range(nph):
rbeam=misc.rotate_alm(rbeamth, lmax, 0, 0, ptg[ith,iph,1])
bar2[ith,iph] = convolve(slm[:,0], rbeam, lmax).real
plt.subplot(2,2,2)
plt.imshow(bar2)
plt.subplot(2,2,3)
plt.imshow(bar2-bar[:,:,0])
plt.show()
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment