guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fixed fft_calls test and dimensions of fourier space

Browse Source
This commit is contained in:
Guilhem Lavaux 2012-11-10 10:56:47 -05:00
parent 414edbd28a
commit 72f658c7cf
2 changed files with 13 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
sample/test_fft_calls.cpp
View File

@ -1,12 +1,17 @@
#include <iostream>
#include "fourier/euclidian.hpp" #include "fourier/euclidian.hpp"
using namespace CosmoTool; using namespace CosmoTool;
using namespace std;
int main() int main()
{ {
EuclidianFourierTransform_2d<double> dft(128,128,1.0,1.0); EuclidianFourierTransform_2d<double> dft(128,128,1.0,1.0);
double volume = 128*128;
dft.realSpace().eigen().setRandom(); dft.realSpace().eigen().setRandom();
dft.analysis(); dft.analysis();
cout << "Map dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl;
cout << "Fourier dot-product = " << dft.fourierSpace().dot_product(dft.fourierSpace()).real()*volume << endl;
return 0; return 0;
} }

12
src/fourier/euclidian.hpp
View File

@ -179,7 +179,9 @@ namespace CosmoTool
{ {
for (long p = 0; p < plane_size; p++) for (long p = 0; p < plane_size; p++)
{ {
result += 2*(conj(d1[q0+N0*p]) * d2[q0+N0*p]).real(); long idx = q0+N0*p;
assert(idx < this->size());
result += 2*(conj(d1[idx]) * d2[idx]).real();
} }
} }
for (long p = 0; p < plane_size; p++) for (long p = 0; p < plane_size; p++)
@ -205,7 +207,7 @@ namespace CosmoTool
typename calls::plan_type m_analysis, m_synthesis; typename calls::plan_type m_analysis, m_synthesis;
double volume; double volume;
long N; long N;
DimArray m_dims; DimArray m_dims, m_dims_hc;
std::vector<double> m_L; std::vector<double> m_L;
public: public:
EuclidianFourierTransform(const DimArray& dims, const std::vector<double>& L) EuclidianFourierTransform(const DimArray& dims, const std::vector<double>& L)
@ -214,6 +216,8 @@ namespace CosmoTool
std::vector<double> dk(L.size()); std::vector<double> dk(L.size());
m_dims = dims; m_dims = dims;
m_dims_hc = dims;
m_dims_hc[0] = dims[0]/2+1;
m_L = L; m_L = L;
N = 1; N = 1;
@ -228,11 +232,11 @@ namespace CosmoTool
realMap = new EuclidianFourierMapReal<T>( realMap = new EuclidianFourierMapReal<T>(
boost::shared_ptr<T>(calls::alloc_real(N), boost::shared_ptr<T>(calls::alloc_real(N),
std::ptr_fun(calls::free)), std::ptr_fun(calls::free)),
dims); m_dims);
fourierMap = new EuclidianFourierMapComplex<T>( fourierMap = new EuclidianFourierMapComplex<T>(
boost::shared_ptr<std::complex<T> >((std::complex<T>*)calls::alloc_complex(N), boost::shared_ptr<std::complex<T> >((std::complex<T>*)calls::alloc_complex(N),
std::ptr_fun(calls::free)), std::ptr_fun(calls::free)),
dims, dk); m_dims_hc, dk);
m_analysis = calls::plan_dft_r2c(dims.size(), &dims[0], m_analysis = calls::plan_dft_r2c(dims.size(), &dims[0],
realMap->data(), (typename calls::complex_type *)fourierMap->data(), realMap->data(), (typename calls::complex_type *)fourierMap->data(),
FFTW_MEASURE); FFTW_MEASURE);