sht_upsample_demo.py 2.54 KB
 Martin Reinecke committed Jun 03, 2020 1 ``````import ducc_0_1.sht as sht `````` Martin Reinecke committed Feb 19, 2020 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ``````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 alm = np.random.uniform(-1., 1., nalm) + 1j*np.random.uniform(-1., 1., nalm) # make a_lm with m==0 real-valued alm[0:lmax+1].imag = 0. `````` Martin Reinecke committed Jun 03, 2020 20 ``````job = sht.sharpjob_d() `````` Martin Reinecke committed Feb 19, 2020 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 ``````# 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 03, 2020 36 ``````map2 = sht.upsample_to_cc(map.reshape((nlat,nlon)), nlat2, False, False) `````` Martin Reinecke committed Feb 19, 2020 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 ``````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 03, 2020 58 ``````map2 = sht.upsample_to_cc(map.reshape((nlat,nlon)), nlat2, True, True) `````` Martin Reinecke committed Feb 19, 2020 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 ``````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 03, 2020 81 ``````map2 = sht.upsample_to_cc(map.reshape((nlat,nlon)), nlat2, False, True) `````` Martin Reinecke committed Feb 19, 2020 82 83 84 85 86 87 88 89 ``````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))``````