sht_upsample_demo.py 2.59 KB
 Martin Reinecke committed Jun 03, 2020 1 ``````import ducc_0_1.sht as sht `````` Martin Reinecke committed Jun 05, 2020 2 ``````import ducc_0_1.misc as misc `````` Martin Reinecke committed Feb 19, 2020 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ``````import numpy as np from time import time def _l2error(a, b): return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2)) lmax = 1023 mmax = lmax print("Generating spherical harmonic coefficients up to {}".format(lmax)) # number of required a_lm coefficients nalm = ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax) # get random a_lm `````` Martin Reinecke committed Jun 08, 2020 17 18 ``````rng = np.random.default_rng(42) alm = rng.uniform(-1., 1., nalm) + 1j*rng.uniform(-1., 1., nalm) `````` Martin Reinecke committed Feb 19, 2020 19 20 21 ``````# make a_lm with m==0 real-valued alm[0:lmax+1].imag = 0. `````` Martin Reinecke committed Jun 03, 2020 22 ``````job = sht.sharpjob_d() `````` Martin Reinecke committed Feb 19, 2020 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 ``````# describe the a_lm array to the job job.set_triangular_alm_info(lmax, mmax) nlon = 2*mmax+2 nlat = lmax+1 # lmax+1 iso-latitude rings, first ring at 0.5*pi/nrings print("Converting them to Fejer1 grid with {}x{} points".format(nlat, nlon)) job.set_fejer1_geometry(nlat, nlon) t0=time() map = job.alm2map(alm) print("time for map synthesis: {}s".format(time()-t0)) nlat2 = 2*lmax+3 t0=time() `````` Martin Reinecke committed Jun 05, 2020 38 ``````map2 = misc.upsample_to_cc(map.reshape((nlat,nlon)), nlat2, False, False) `````` Martin Reinecke committed Feb 19, 2020 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 ``````print("time for upsampling: {}s".format(time()-t0)) job.set_cc_geometry(nlat2, nlon) t0=time() alm2 = job.map2alm(map2.reshape((-1,))) print("time for map analysis: {}s".format(time()-t0)) # make sure input was recovered accurately print("L2 error: ", _l2error(alm,alm2)) nlon = 2*mmax+2 nlat = lmax+2 # lmax+2 iso-latitude rings, first ring at north pole print("Converting them to Clenshaw-Curtis grid with {}x{} points".format(nlat, nlon)) job.set_cc_geometry(nlat, nlon) t0=time() map = job.alm2map(alm) print("time for map synthesis: {}s".format(time()-t0)) nlat2 = 2*lmax+3 t0=time() `````` Martin Reinecke committed Jun 05, 2020 60 ``````map2 = misc.upsample_to_cc(map.reshape((nlat,nlon)), nlat2, True, True) `````` Martin Reinecke committed Feb 19, 2020 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 ``````print("time for upsampling: {}s".format(time()-t0)) job.set_cc_geometry(nlat2, nlon) t0=time() alm2 = job.map2alm(map2.reshape((-1,))) print("time for map analysis: {}s".format(time()-t0)) # make sure input was recovered accurately print("L2 error: ", _l2error(alm,alm2)) nlon = 2*mmax+2 nlat = lmax+1 # pixelisation according to https://arxiv.org/abs/1110.6298 print("Converting them to McEwen-Wiaux grid with {}x{} points".format(nlat, nlon)) job.set_mw_geometry(nlat, nlon) t0=time() map = job.alm2map(alm) print("time for map synthesis: {}s".format(time()-t0)) nlat2 = 2*lmax+3 t0=time() `````` Martin Reinecke committed Jun 05, 2020 83 ``````map2 = misc.upsample_to_cc(map.reshape((nlat,nlon)), nlat2, False, True) `````` Martin Reinecke committed Feb 19, 2020 84 85 86 87 88 89 90 91 ``````print("time for upsampling: {}s".format(time()-t0)) job.set_cc_geometry(nlat2, nlon) t0=time() alm2 = job.map2alm(map2.reshape((-1,))) print("time for map analysis: {}s".format(time()-t0)) # make sure input was recovered accurately print("L2 error: ", _l2error(alm,alm2))``````