Commit 4251ce8f authored by Martin Reinecke's avatar Martin Reinecke
Browse files

make DH grid available within pysharp

parent c8409114
......@@ -14,9 +14,6 @@ lmax = 2047
# maximum m. For SHTOOLS this is alway equal to lmax, if I understand correctly.
mmax = lmax
# Number of iso-latitude rings required for Gauss-Legendre grid
nlat = lmax+1
# Number of pixels per ring. Must be >=2*lmax+1, but I'm choosing a larger
# number for which the FFT is faster.
nlon = 4096
......@@ -43,6 +40,12 @@ alm[0:lmax+1].imag = 0.
# describe the a_lm array to the job
job.set_triangular_alm_info(lmax, mmax)
print("testing Gauss-Legendre grid")
# Number of iso-latitude rings required for Gauss-Legendre grid
nlat = lmax+1
# describe the Gauss-Legendre geometry to the job
job.set_Gauss_geometry(nlat, nlon)
......@@ -63,3 +66,30 @@ print("time for map analysis: {}s".format(time()-t0))
# make sure input was recovered accurately
assert_allclose(alm, alm2)
print("testing Driscoll-Healy grid")
# Number of iso-latitude rings required for Gauss-Legendre grid
nlat = 2*lmax
# describe the Gauss-Legendre geometry to the job
job.set_DH_geometry(nlat, nlon)
# go from a_lm to map
t0=time()
map = job.alm2map(alm)
print("time for map synthesis: {}s".format(time()-t0))
# map is a 1D real-valued array with (nlat*nlon) entries. It can be reshaped
# to (nlat, nlon) for plotting.
# Libsharp woks on "1D" maps because it apso supports pixelizations that varying
# number of pixels on each iso-latitude ring, which cannot be represented by 2D
# arrays (e.g. Healpix)
t0=time()
alm2 = job.map2alm(map)
print("time for map analysis: {}s".format(time()-t0))
# make sure input was recovered accurately
assert_allclose(alm, alm2)
......@@ -84,6 +84,13 @@ template<typename T> class py_sharpjob
npix_=nrings*nphi;
ginfo = sharp_make_ecp_geom_info (nrings, nphi, 0., 1, nphi);
}
void set_DH_geometry(int64_t nrings, int64_t nphi)
{
MR_assert(nrings>1,"bad nrings value");
MR_assert(nphi>0,"bad nphi value");
npix_=nrings*nphi;
ginfo = sharp_make_dh_geom_info (nrings, nphi, 0., 1, nphi);
}
void set_triangular_alm_info (int64_t lmax, int64_t mmax)
{
MR_assert(mmax>=0,"negative mmax");
......@@ -172,6 +179,8 @@ PYBIND11_MODULE(pysharp, m)
"nside"_a)
.def("set_ECP_geometry", &py_sharpjob<double>::set_ECP_geometry,
"nrings"_a, "nphi"_a)
.def("set_DH_geometry", &py_sharpjob<double>::set_ECP_geometry,
"nrings"_a, "nphi"_a)
.def("set_triangular_alm_info",
&py_sharpjob<double>::set_triangular_alm_info, "lmax"_a, "mmax"_a)
.def("n_alm", &py_sharpjob<double>::n_alm)
......
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