From 72f658c7cf206b7dd045b058dcf2e947421bd7bf Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <guilhem.lavaux@gmail.com>
Date: Sat, 10 Nov 2012 10:56:47 -0500
Subject: [PATCH] Fixed fft_calls test and dimensions of fourier space

---
 sample/test_fft_calls.cpp |  5 +++++
 src/fourier/euclidian.hpp | 12 ++++++++----
 2 files changed, 13 insertions(+), 4 deletions(-)

diff --git a/sample/test_fft_calls.cpp b/sample/test_fft_calls.cpp
index 35eeb12..d633c14 100644
--- a/sample/test_fft_calls.cpp
+++ b/sample/test_fft_calls.cpp
@@ -1,12 +1,17 @@
+#include <iostream>
 #include "fourier/euclidian.hpp"
 
 using namespace CosmoTool;
+using namespace std;
 
 int main()
 {
   EuclidianFourierTransform_2d<double> dft(128,128,1.0,1.0);
+  double volume = 128*128;
 
   dft.realSpace().eigen().setRandom();
   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;
 }
diff --git a/src/fourier/euclidian.hpp b/src/fourier/euclidian.hpp
index a29c623..fac27f6 100644
--- a/src/fourier/euclidian.hpp
+++ b/src/fourier/euclidian.hpp
@@ -179,7 +179,9 @@ namespace CosmoTool
 	{
 	  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++)
@@ -205,7 +207,7 @@ namespace CosmoTool
     typename calls::plan_type m_analysis, m_synthesis;
     double volume;
     long N;
-    DimArray m_dims;
+    DimArray m_dims, m_dims_hc;
     std::vector<double> m_L;
   public:
     EuclidianFourierTransform(const DimArray& dims, const std::vector<double>& L)
@@ -214,6 +216,8 @@ namespace CosmoTool
       std::vector<double> dk(L.size());
       
       m_dims = dims;
+      m_dims_hc = dims;
+      m_dims_hc[0] = dims[0]/2+1;
       m_L = L;
 
       N = 1;
@@ -228,11 +232,11 @@ namespace CosmoTool
       realMap = new EuclidianFourierMapReal<T>(
 		   boost::shared_ptr<T>(calls::alloc_real(N),
 					std::ptr_fun(calls::free)),
-		   dims);
+		   m_dims);
       fourierMap = new EuclidianFourierMapComplex<T>(
 		   boost::shared_ptr<std::complex<T> >((std::complex<T>*)calls::alloc_complex(N),
 						       std::ptr_fun(calls::free)), 
-		   dims, dk);
+		   m_dims_hc, dk);
       m_analysis = calls::plan_dft_r2c(dims.size(), &dims[0],
 				       realMap->data(), (typename calls::complex_type *)fourierMap->data(),
 				       FFTW_MEASURE);