diff --git a/src/fourier/base_types.hpp b/src/fourier/base_types.hpp index 669e382..9a02a7e 100644 --- a/src/fourier/base_types.hpp +++ b/src/fourier/base_types.hpp @@ -63,7 +63,14 @@ namespace CosmoTool eigen() += map2->eigen(); } - virtual FourierMap *copy() const = 0; + virtual FourierMap *copy() const + { + FourierMap *m = this->mimick(); + + m->eigen() = this->eigen(); + return m; + } + virtual FourierMap *mimick() const = 0; }; diff --git a/src/fourier/healpix.hpp b/src/fourier/healpix.hpp new file mode 100644 index 0000000..1759a78 --- /dev/null +++ b/src/fourier/healpix.hpp @@ -0,0 +1,135 @@ +#ifndef __COSMOTOOL_FOURIER_HEALPIX_HPP +#define __COSMOTOOL_FOURIER_HEALPIX_HPP + +#include +#include +#include +#include +#include + +namespace CosmoTool +{ + + template + class HealpixFourierMap: public FourierMap, public Healpix_Base + { + private: + T *m_data; + Eigen::aligned_allocator alloc; + public: + HealpixFourierMap(long nSide) + : Healpix_Base(RING, nSide, SET_NSIDE) + { + m_data = alloc.allocate(Npix); + } + + virtual ~HealpixFourierMap() + { + alloc.deallocate(m_data); + } + + virtual const T* data() const { return m_data; } + virtual T *data() { return m_data; } + virtual long size() const { return Npix(); } + + virtual FourierMap *mimick() const + { + return new HealpixFourierMap(Nside()); + } + }; + + + template + class HealpixFourierALM: public FourierMap >, public Alm_Base + { + private: + std::complex *alms; + long size; + Eigen::aligned_allocator > alloc; + public: + HealpixFourierALM(long Lmax, long Mmax) + : Alm_Base(Lmax, Mmax) + { + size = Num_Alms(Lmax, Mmax); + alms = alloc.allocate(size); + } + + virtual ~HealpixFourierALM() + { + alloc.deallocate(alms); + } + + virtual const T* data() const { return alms; } + virtual T * data() { return alms;} + virtual long size() const { return size; } + + virtual FourierMap *mimick() const + { + return new HealpixFourierALM(Lmax(), Mmax()); + } + }; + + template + class HealpixFourierTransform: public FourierTransform + { + private: + HealpixFourierMap realMap; + HealpixFourierALM fourierMap; + psht_joblist jobs; + public: + HealpixFourierTransform(long nSide, long Lmax, long Mmax) + : realMap(nSide), fourierMap(Lmax, Mmax) + { + jobs.set_Healpix_geometry(nSide); + jobs.set_triangular_alm_info(Lmax, Mmax); + } + + virtual ~HealpixFourierTransform() {} + + virtual const FourierMap >& fourierSpace() const { return fourierMap; } + + virtual FourierMap >& fourierSpace() { return fourierMap; } + + virtual const FourierMap& realSpace() const { return realMap; } + + virtual FourierMap& realSpace() { return realMap; } + + virtual FourierTransform *mimick() const + { + return new HealpixFourierTransform(realMap.Nside(), fourierMap.Lmax(), fourierMap.Mmax()); + } + + virtual void analysis() + { + jobs.add_map2alm(realMap.data(), + reinterpret_cast *>(fourierMap.data()), + false); + jobs.execute(); + jobs.clear_jobs(); + } + + virtual void synthesis() + { + jobs.add_alm2map(reinterpret_cast *>(fourierMap.data()), + realMap.data(), + false); + jobs.execute(); + jobs.clear_jobs(); + } + + virtual void analysis_conjugate() + { + synthesis(); + realMap.scale(4*M_PI/realMap.Npix()); + } + + virtual void synthesis_conjugate() + { + analysis(); + fourierMap.scale(realMap.Npix()/(4*M_PI)); + } + + }; +}; + +#endif