Fixed K computation. Handle fully the Nyquist plane depending on the parity of dimensions.

sample/test_fft_calls.cpp

@@ -9,14 +9,16 @@ using namespace std;
 
 double spectrum_generator(double k)
 {
-  return 1/(0.1+pow(k, 3.0));
+  if (k==0)
+    return 0;
+  return 1/(0.01+pow(k, 3.0));
 }
 
 int main()
 {
-  EuclidianFourierTransform_2d<double> dft(128,128,1.0,1.0);
+  EuclidianFourierTransform_2d<double> dft(1024,1024,1.0,1.0);
   EuclidianSpectrum_1D<double> spectrum(spectrum_generator);
-  double volume = 128*128;
+  double volume = 1024*1024;
   gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
 
   dft.realSpace().eigen().setRandom();
@@ -33,7 +35,7 @@ int main()
   dft.fourierSpace() = *m.get();
   dft.synthesis();
 
-  uint32_t dims[2] = { 128, 128 };
+  uint32_t dims[2] = { 1024, 1024 };
   CosmoTool::saveArray("generated_map.nc", dft.realSpace().data(), dims, 2);
   
   return 0;

src/fourier/euclidian.hpp

@@ -109,19 +109,26 @@ namespace CosmoTool
   protected:
     typedef boost::shared_ptr<std::complex<T> > ptr_t;
     std::vector<double> delta_k;
+    int m_dim0;
+    bool even0, alleven;
     long plane_size;
   public:
     typedef typename EuclidianFourierMapBase<std::complex<T> >::DimArray DimArray;
 
     EuclidianFourierMapComplex(ptr_t indata,
+			       int dim0,
 			       const DimArray& indims, 
 			       const std::vector<double>& dk)
-      : EuclidianFourierMapBase<std::complex<T> >(indata, indims), delta_k(dk)
+      : EuclidianFourierMapBase<std::complex<T> >(indata, indims), delta_k(dk), m_dim0(dim0), even0((dim0 % 2)==0)
     {
       assert(dk.size() == indims.size()); 
       plane_size = 1;
+      alleven = true;
       for (int q = 1; q < indims.size(); q++)
-	plane_size *= indims[q];
+	{
+	  plane_size *= indims[q];
+	  alleven = alleven && ((indims[q]%2)==0);
+	}
     }
 
     virtual FourierMap<std::complex<T> > *mimick() const
@@ -131,6 +138,7 @@ namespace CosmoTool
 		   ptr_t((std::complex<T> *)
 			 fftw_malloc(sizeof(std::complex<T>)*this->size()), 
 			 std::ptr_fun(fftw_free)),
+		   m_dim0,
 		   this->getDims(),
 		   this->delta_k);
     }
@@ -141,7 +149,9 @@ namespace CosmoTool
       const DimArray& dims = this->getDims();
       assert(ik.size() == dims.size());
       double k2 = 0;
-      for (int q = 0; q < ik.size(); q++)
+      k2 += CosmoTool::square(ik[0]*delta_k[0]);
+
+      for (int q = 1; q < ik.size(); q++)
 	{
 	  int dk = ik[q];
 
@@ -165,6 +175,9 @@ namespace CosmoTool
       return get_K(d);
     }
 
+    bool allDimensionsEven() const { return alleven; }
+    bool firstDimensionEven() const { return even0; }
+
     virtual std::complex<T> dot_product(const FourierMap<std::complex<T> >& other) const
       throw(std::bad_cast)
     {
@@ -175,7 +188,7 @@ namespace CosmoTool
       const std::complex<T> *d1 = this->data();
       const std::complex<T> *d2 = m2.data();
       const DimArray& dims = this->getDims();
-      int N0 = dims[0];
+      int N0 = dims[0] + (even0 ? 0 : 1);
       std::complex<T> result = 0;
 
       for (long q0 = 1; q0 < N0-1; q0++)
@@ -187,11 +200,14 @@ namespace CosmoTool
 	      result += 2*(conj(d1[idx]) * d2[idx]).real();
 	    }
 	}
-      for (long p = 0; p < plane_size; p++)
+      if (!even0)
 	{
-	  long q0 = N0*p, q1 = (p+1)*N0-1;
-	  result += conj(d1[q0]) * d2[q0];
-	  result += conj(d1[q1]) * d2[q1];
+	  for (long p = 0; p < plane_size; p++)
+	    {
+	      long q0 = N0*p, q1 = (p+1)*N0-1;
+	      result += conj(d1[q0]) * d2[q0];
+	      result += conj(d1[q1]) * d2[q1];
+	    }
 	}
       return result;
     }
@@ -239,7 +255,7 @@ namespace CosmoTool
       fourierMap = new EuclidianFourierMapComplex<T>(
 		   boost::shared_ptr<std::complex<T> >((std::complex<T>*)calls::alloc_complex(N),
 						       std::ptr_fun(calls::free)), 
-		   m_dims_hc, dk);
+		   dims[0], m_dims_hc, dk);
       m_analysis = calls::plan_dft_r2c(dims.size(), &dims[0],
 				       realMap->data(), (typename calls::complex_type *)fourierMap->data(),
 				       FFTW_MEASURE);
@@ -306,6 +322,26 @@ namespace CosmoTool
     }
   };
 
+  template<typename T>
+  class EuclidianFourierTransform_1d: public EuclidianFourierTransform<T>
+  {
+  private:
+    template<typename T2>
+    static std::vector<T2> make_1d_vector(T2 a)
+    {
+      T2 arr[2] = { a};
+      return std::vector<T2>(&arr[0],&arr[1]);
+    }
+  public:
+    EuclidianFourierTransform_1d(int Nx, double Lx)
+      : EuclidianFourierTransform<T>(make_1d_vector<int>(Nx), make_1d_vector<double>(Lx))
+    {
+    }
+
+    virtual ~EuclidianFourierTransform_1d() {}
+
+  };
+
   template<typename T>
   class EuclidianFourierTransform_2d: public EuclidianFourierTransform<T>
   {
@@ -363,6 +399,7 @@ namespace CosmoTool
     long idx;
     const DimArray& dims = rand_map.getDims();
     long plane_size;
+    bool alleven = rand_map.allDimensionsEven();
 
     for (long p = 1; p < m_c.size(); p++)
       {
@@ -371,16 +408,18 @@ namespace CosmoTool
 			       gsl_ran_gaussian(rng, A_k));
       }
     // Generate the mean value
-    d[0] = std::complex<T>(std::sqrt(f(0)), 0);
+    d[0] = std::complex<T>(gsl_ran_gaussian(rng, std::sqrt(f(0))), 0);
+
+    if (!rand_map.firstDimensionEven())
+      return ret_map;
+
     // Correct the Nyquist plane
     idx = dims[0]-1; // Stick to the last element of the first dimension
+    d[idx] = std::complex<T>(d[idx].real() + d[idx].imag(), 0);
     // 1D is special case
     if (dims.size() == 1)
-      {
-	d[idx] = std::complex<T>(d[idx].real() + d[idx].imag(), 0);
-        return ret_map;
-      }
-
+      return ret_map;
+    
     plane_size = 1;
     for (int q = 1; q < dims.size(); q++)
       {
@@ -395,8 +434,16 @@ namespace CosmoTool
 	assert(q2 < plane_size*dims[0]);
 	d[q] = conj(d[q2]);
       }
-    long q = dims[0];
-    d[q] = std::complex<T>(d[q].real() + d[q].imag());
+
+    if (alleven)
+      {
+	long q = 0;
+	for (int i = dims.size()-1; i >= 1; i--)
+	  q = dims[i]*q + dims[i]/2;
+	q += dims[0]-1;
+	d[q] = std::complex<T>(d[q].real()+d[q].imag(),0);
+      }
+
     return ret_map;
   }