Imported latest snapshot of cosmotoolbox

external/cosmotool/sample/test_healpix_calls.cpp

@@ -0,0 +1,25 @@
+#include <iostream>
+#include "fourier/healpix.hpp"
+
+using namespace CosmoTool;
+using namespace std;
+
+int main()
+{
+  HealpixFourierTransform<double> dft(128,2*128,2*128, 40);
+  HealpixUtilityFunction<double> utils;
+  long Npix = dft.realSpace().size();
+
+  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()*Npix/(4*M_PI) << endl;
+  cout << "Spectrum analysis" << endl;
+  HealpixUtilityFunction<double>::Spectrum_ptr s = utils.estimateSpectrumFromMap(dft.fourierSpace());
+
+  dft.synthesis();
+  cout << "Resynthesis dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl;
+
+  return 0;
+
+}

external/cosmotool/src/cosmopower.cpp

@@ -0,0 +1,306 @@
+#include <cassert>
+#include <vector>
+#include <algorithm>
+#include <iostream>
+#include <cmath>
+#include <fstream>
+#include <gsl/gsl_integration.h>
+#include "cosmopower.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+#define USE_GSL
+#define TOLERANCE 1e-6
+#define NUM_ITERATION 8000
+
+
+CosmoPower::CosmoPower()
+{
+  eval = &CosmoPower::powerEfstathiou;
+  
+  n = 1.0;
+  K0 = 1;
+  V_LG_CMB = 627;
+
+  CMB_VECTOR[0] = 56.759;
+  CMB_VECTOR[1] = -540.02;
+  CMB_VECTOR[2] = 313.50;
+
+  h = 0.719;
+  SIGMA8 = 0.77;
+  OMEGA_B = 0.043969;
+  OMEGA_C = 0.21259;
+
+  Theta_27 = 2.728/2.7;
+
+  updateCosmology();
+}
+
+/*
+ * This is \hat{tophat}
+ */
+static double tophatFilter(double u)
+{
+  if (u != 0)    
+    return 3 / (u*u*u) * (sin(u) - u * cos(u));
+  else
+    return 1;
+}
+
+static double powC(double q, double alpha_c)
+{
+  return 14.2 / alpha_c + 386 / (1 + 69.9 * pow(q, 1.08));
+}
+
+static double T_tilde_0(double q, double alpha_c, double beta_c)
+{
+  double a = log(M_E + 1.8 * beta_c * q);
+  return a / ( a  + powC(q, alpha_c) * q * q);
+}
+
+static double j_0(double x)
+{
+  if (x == 0)
+    return 1.0;
+  return sin(x)/x;
+}
+
+static double powG(double y)
+{
+  double a = sqrt(1 + y);
+
+  return y * (-6 * a + (2 + 3 * y) *log((a + 1)/(a - 1)));
+}
+ 
+double CosmoPower::powerEfstathiou(double k)
+{
+  double a = 6.4/Gamma0;
+  double b = 3/Gamma0;
+  double c = 1.7/Gamma0;
+  double nu = 1.13;
+  
+  double f = (a*k) + pow(b*k,1.5) + pow(c*k,2);
+
+  // EFSTATHIOU  ET AL. (1992)
+  return normPower * pow(k,n) * pow(1+pow(f,nu),(-2/nu));
+}
+
+double CosmoPower::powerHuWiggles(double k)
+{
+  // EISENSTEIN ET HU (1998)
+  // FULL POWER SPECTRUM WITH BARYONS AND WIGGLES
+
+  double k_silk = 1.6 * pow(OMEGA_B * h * h, 0.52) * pow(OmegaEff, 0.73) * (1 + pow(10.4 * OmegaEff, -0.95));
+  double z_eq = 2.50e4 * OmegaEff * pow(Theta_27, -4);
+  double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
+  double f = 1 / (1 + pow(k * s / 5.4, 4));
+  double k_eq = 7.46e-2 * OmegaEff * pow(Theta_27, -2);
+  double a1 = pow(46.9 * OmegaEff, 0.670) * (1 + pow(32.1 * OmegaEff, -0.532));
+  double a2 = pow(12.0 * OmegaEff, 0.424) * (1 + pow(45.0 * OmegaEff, -0.582));
+  double alpha_c = pow(a1, -OMEGA_B/ OMEGA_0) * pow(a2, -pow(OMEGA_B / OMEGA_0, 3));  
+
+  double q = k / (13.41 * k_eq);
+  double b1_betac = 0.944 * 1/(1 + pow(458 * OmegaEff, -0.708));
+  double b2_betac = pow(0.395 * OmegaEff, -0.0266);
+  double beta_c = 1/ ( 1 + b1_betac * (pow(OMEGA_C / OMEGA_0, b2_betac) - 1)   );
+  double T_c = f * T_tilde_0(q, 1, beta_c) + (1 - f) * T_tilde_0(q, alpha_c, beta_c);
+
+  double b1_zd = 0.313 * pow(OmegaEff, -0.419) * (1 + 0.607 * pow(OmegaEff, 0.674));
+  double b2_zd = 0.238 * pow(OmegaEff, 0.223);
+  double z_d = 1291 * pow(OmegaEff, 0.251) / (1 + 0.659 * pow(OmegaEff, 0.828)) * (1 + b1_zd * pow(OmegaEff, b2_zd));
+  double R_d = 31.5 * OMEGA_B * h * h * pow(Theta_27, -4) * 1e3 / z_d;
+
+  double alpha_b = 2.07 * k_eq * s * pow(1 + R_d, -0.75) * powG((1 + z_eq)/(1 + z_d));
+  double beta_b = 0.5 + OMEGA_B / OMEGA_0 + (3 - 2 * OMEGA_B / OMEGA_0) * sqrt(pow(17.2 * OmegaEff, 2) + 1);
+  double beta_node = 8.41 * pow(OmegaEff, 0.435);
+  double s_tilde = s * pow(1 + pow(beta_node / (k * s), 3), -1./3);
+
+  double T_b = (T_tilde_0(q, 1, 1) / (1 + pow(k * s / 5.2, 2)) + alpha_b / (1 + pow(beta_b / (k * s), 3)) * exp(-pow(k/k_silk, 1.4))) * j_0(k * s_tilde);  
+
+  double T_k = OMEGA_B/OMEGA_0 * T_b + OMEGA_C/OMEGA_0 * T_c;  
+
+  return normPower * pow(k,n) * T_k * T_k;
+}
+
+double CosmoPower::powerHuBaryons(double k)
+{
+  double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
+  double alpha_Gamma = 1 - 0.328 * log(431 * OmegaEff) * OMEGA_B / OMEGA_0 + 0.38 * log(22.3 * OmegaEff) * pow(OMEGA_B / OMEGA_0, 2);
+  double GammaEff = OMEGA_0 * h * (alpha_Gamma + (1 - alpha_Gamma)/(1 + pow(0.43 * k * s, 4)));
+  double q = k/(h*GammaEff) * pow(Theta_27, 2);
+  double L_0 = log(2 * M_E + 1.8 * q);
+  double C_0 = 14.2 + 731 / (1 + 62.5 * q);
+  double T0 = L_0 / (L_0 + C_0 * q * q);
+  
+  return normPower * pow(k,n) * T0 * T0;
+}
+
+double CosmoPower::powerOld(double k)
+{
+  static const double l = 1/(Omega * h*h);
+  static const double alpha = 1.7 * l, beta = 9.0 * pow(l, 1.5), gamma = l*l;
+  return normPower * pow(k,n) * pow(1 + alpha * k +  beta * pow(k,1.5) + gamma *k*k,-2);
+}
+
+double CosmoPower::powerSugiyama(double k)
+{
+  double q = k * Theta_27*Theta_27 / (OmegaEff * exp(-OMEGA_B - sqrt(h/0.5)*OMEGA_B/OMEGA_0));
+  double L0 = log(2*M_E + 1.8 * q);
+  double C0 = 14.2 + 731 / (1 + 62.5 * q);
+  double T_k = L0 / (L0 + C0 * q*q);
+
+  return normPower * pow(k,n) * T_k * T_k;
+}
+
+double CosmoPower::powerBardeen(double k)
+{
+  double q = k / (OmegaEff);
+
+  double poly = 1 + 3.89 * q + pow(16.1*q,2) + pow(5.46*q,3) + pow(6.71*q,4);
+  double T_k = log(1+2.34*q)/(2.34*q) * pow(poly,-0.25);
+
+  return normPower * pow(k,n) * T_k * T_k;
+}
+
+double CosmoPower::powerBDM(double k)
+{
+  k /= h*h;
+  double alpha1 = 190;
+  double Gmu = 4.6;
+  double alpha2 = 11.5;
+  double alpha3 = 11;
+  double alpha4 = 12.55;
+  double alpha5 = 0.0004;
+  return normPower*k*alpha1*alpha1*Gmu*Gmu/(1+(alpha2*k)+pow(alpha3*k,2)+pow(alpha4*k,3))*pow(1+pow(alpha5/k,2), -2);
+}
+
+double CosmoPower::powerTest(double k)
+{
+  return 1/(1+k*k);
+}
+
+/*
+ * This function computes the normalization of the power spectrum. It requests
+ * a sigma8 (density fluctuations within 8 Mpc/h)
+ */
+static double gslPowSpecNorm(double k, void *params)
+{
+  CosmoPower *c = (CosmoPower *)params;
+
+  return c->integrandNormalize(k);
+}
+
+double CosmoPower::integrandNormalize(double x)
+{
+  double k = (1-x)/x;
+  double f = tophatFilter(k*8.0/h);
+  return power(k)*k*k*f*f/(x*x);
+}
+
+void CosmoPower::normalize()
+{
+  double normVal = 0;
+  double abserr;
+  gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION);
+  gsl_function f;
+
+  f.function = gslPowSpecNorm;
+  f.params = this;
+
+  normPower = 1;
+
+  ofstream ff("PP_k.txt");
+  for (int i = 0; i < 100; i++)
+    {
+       double k = pow(10.0, 4.0*i/100.-2);
+       ff << k << " " << power(k) << endl;
+    }
+
+ // gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &normVal, &abserr);
+  gsl_integration_qag(&f, 0, 1, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &normVal, &abserr);
+  gsl_integration_workspace_free(w);
+
+  normVal /= (2*M_PI*M_PI);
+
+  normPower =  SIGMA8*SIGMA8/normVal;
+}
+
+void CosmoPower::updateCosmology()
+{
+  OMEGA_0 = OMEGA_B+OMEGA_C;
+  Omega = OMEGA_0;
+  beta = pow(OMEGA_0, 5./9);
+  OmegaEff = OMEGA_0 * h * h;
+  Gamma0 = OMEGA_0 * h * h;
+  omega_B = OMEGA_B * h * h;
+  omega_C = OMEGA_C * h * h;
+}
+
+void CosmoPower::updatePhysicalCosmology()
+{
+  OMEGA_B = omega_B / (h*h);
+  OMEGA_C = omega_C / (h*h);
+  OMEGA_0 = Gamma0 / (h*h);
+  beta = pow(OMEGA_0, 5./9);
+}
+
+double CosmoPower::eval_theta_theta(double k)
+{
+  // Jennings (2012) fit
+  double P_deltadelta = power(k);
+
+  static const double alpha0 = -12480.5, alpha1 = 1.824, alpha2 = 2165.87, alpha3=1.796;
+  if (k > 0.3)
+    return 0;
+  double r =(alpha0*sqrt(P_deltadelta) + alpha1*P_deltadelta*P_deltadelta)/(alpha2 + alpha3*P_deltadelta);
+  assert(P_deltadelta > 0);
+  
+  if (r < 0)
+    return 0;
+  return r;
+}
+
+double CosmoPower::power(double k)
+{
+  return (this->*eval)(k);
+}
+
+
+void CosmoPower::setFunction(CosmoFunction f)
+{
+  switch (f)
+    {
+    case POWER_EFSTATHIOU:
+      eval = &CosmoPower::powerEfstathiou;
+      break;
+    case HU_WIGGLES:
+      eval = &CosmoPower::powerHuWiggles;
+      break;
+    case HU_BARYON:
+      eval = &CosmoPower::powerHuBaryons;
+      break;
+    case OLD_POWERSPECTRUM:
+      eval = &CosmoPower::powerOld;
+      break;
+    case POWER_BARDEEN:
+      eval = &CosmoPower::powerBardeen;
+      break;
+    case POWER_SUGIYAMA:
+      eval = &CosmoPower::powerSugiyama;
+      break;
+    case POWER_BDM:
+      eval = &CosmoPower::powerBDM;
+      break;
+    case POWER_TEST:
+      eval = &CosmoPower::powerTest;
+      break;
+    default:
+      abort();
+    }
+}
+
+void CosmoPower::setNormalization(double A_K)
+{
+  normPower = A_K/power(0.002);
+}

external/cosmotool/src/cosmopower.hpp

@@ -0,0 +1,74 @@
+#ifndef _COSMOPOWER_HPP
+#define _COSMOPOWER_HPP
+
+namespace CosmoTool {
+
+  class CosmoPower
+  {
+  public:
+    // PRIMARY VARIABLES
+    double n;
+    double K0;
+    double V_LG_CMB;
+
+    double CMB_VECTOR[3];
+
+    double h;
+    double SIGMA8;
+    double OMEGA_B;
+    double OMEGA_C;
+    double omega_B;
+    double omega_C;
+    double Theta_27;
+    
+    // DERIVED VARIABLES    
+    double OMEGA_0;
+    double Omega;
+    double beta;
+    double OmegaEff;
+    double Gamma0;
+    double normPower;
+
+    enum CosmoFunction
+      {
+	POWER_EFSTATHIOU,
+	HU_WIGGLES,
+	HU_BARYON,
+	OLD_POWERSPECTRUM,
+	POWER_BARDEEN,
+	POWER_SUGIYAMA,
+	POWER_BDM,
+	POWER_TEST
+      };
+
+    CosmoPower();
+
+    void setFunction(CosmoFunction f);
+
+    void updateCosmology();
+    void updatePhysicalCosmology();
+    void normalize();
+    void setNormalization(double A_K);
+
+ 
+    double eval_theta_theta(double k);
+    double power(double k);
+
+    double integrandNormalize(double k);
+  private:
+    double (CosmoPower::*eval)(double);
+
+    double powerEfstathiou(double k);
+    double powerHuWiggles(double k);
+    double powerHuBaryons(double k);
+    double powerOld(double k);
+    double powerBardeen(double k);
+    double powerSugiyama(double k);
+    double powerBDM(double k);
+    double powerTest(double k);
+
+  };
+
+};
+
+#endif

external/cosmotool/src/fourier/details/euclidian_maps.hpp

@@ -0,0 +1,232 @@
+#ifndef __DETAILS_EUCLIDIAN_MAPS
+#define __DETAILS_EUCLIDIAN_MAPS
+
+
+namespace CosmoTool
+{
+  
+  template<typename T>
+  class EuclidianFourierMapBase: public FourierMap<T>
+  {
+  public:
+    typedef std::vector<int> DimArray;
+  private:
+    boost::shared_ptr<T> m_data;
+    DimArray m_dims;
+    long m_size;
+  public:
+    
+    EuclidianFourierMapBase(boost::shared_ptr<T> indata, const DimArray& indims)
+    {
+      m_data = indata;
+      m_dims = indims;
+      m_size = 1;
+      for (int i = 0; i < m_dims.size(); i++)
+	m_size *= m_dims[i];
+    }
+
+    virtual ~EuclidianFourierMapBase()
+    {
+    }
+
+    const DimArray& getDims() const { return m_dims; }
+
+    virtual const T *data() const { return m_data.get(); }
+    virtual T *data() { return m_data.get(); }
+    virtual long size() const  { return m_size; } 
+    
+    virtual FourierMap<T> *copy() const
+    {
+      FourierMap<T> *m = this->mimick();
+      m->eigen() = this->eigen();
+      return m;
+    }
+
+  };
+
+  template<typename T>
+  class EuclidianFourierMapReal: public EuclidianFourierMapBase<T>
+  {
+  public:
+    typedef typename EuclidianFourierMapBase<T>::DimArray DimArray;
+
+    EuclidianFourierMapReal(boost::shared_ptr<T> indata, const DimArray& indims)
+      : EuclidianFourierMapBase<T>(indata, indims) 
+    {}
+
+    virtual FourierMap<T> *mimick() const
+    {
+      return new EuclidianFourierMapReal<T>(
+		      boost::shared_ptr<T>((T *)fftw_malloc(sizeof(T)*this->size()), 
+					   std::ptr_fun(fftw_free)),
+		      this->getDims());
+    }
+
+    virtual T dot_product(const FourierMap<T>& other) const
+      throw(std::bad_cast)
+    {
+      const EuclidianFourierMapReal<T>& m2 = dynamic_cast<const EuclidianFourierMapReal<T>&>(other);
+      if (this->size() != m2.size())
+        throw std::bad_cast();
+
+      return (this->eigen()*m2.eigen()).sum();
+    }
+  };
+
+  template<typename T>
+  class EuclidianFourierMapComplex: public EuclidianFourierMapBase<std::complex<T> >
+  {
+  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), 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];
+	  alleven = alleven && ((indims[q]%2)==0);
+	}
+    }
+
+    virtual FourierMap<std::complex<T> > *mimick() const
+    {
+      return 
+	new EuclidianFourierMapComplex<T>(
+		   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);
+    }
+
+    const std::vector<double>& get_delta_k() const
+    {
+      return this->delta_k;
+    }
+
+    template<typename Array, typename Array2>
+    void get_Kvec(const Array& ik, Array2& kvec) const
+    {
+      const DimArray& dims = this->getDims();
+      assert(ik.size() == dims.size());
+      assert(kvec.size() == dims.size());
+
+      kvec[0] = ik[0] * delta_k[0];
+      for (int q = 1; q < ik.size(); q++)
+	{
+	  int dk = ik[q];
+	  if (dk > dims[q]/2)
+	    dk = dk - dims[q];
+
+          kvec[q] = dk*delta_k[q];
+        }
+    }
+
+    template<typename Array2>
+    void get_Kvec(long p, Array2& kvec) const
+    {
+      const DimArray& dims = this->getDims();
+      DimArray d(delta_k.size());
+      get_IKvec(p, d);
+      get_Kvec(d, kvec);
+    }
+
+    void get_IKvec(long p, DimArray& ikvec) const
+    {
+      const DimArray& dims = this->getDims();
+      assert(dims.size()==ikvec.size());
+      for (int q = 0; q < ikvec.size(); q++)
+       {
+         ikvec[q] = p%dims[q];
+         p = (p-ikvec[q])/dims[q];
+       }
+    }
+
+
+    template<typename Array>
+    double get_K(const Array& ik) const
+    {
+      const DimArray& dims = this->getDims();
+      assert(ik.size() == dims.size());
+      double k2 = 0;
+      k2 += CosmoTool::square(ik[0]*delta_k[0]);
+
+      for (int q = 1; q < ik.size(); q++)
+	{
+	  int dk = ik[q];
+
+	  if (dk > dims[q]/2)
+	    dk = dk - dims[q];
+	  
+	  k2 += CosmoTool::square(delta_k[q]*dk);
+	}
+      return std::sqrt(k2);
+    }
+
+    double get_K(long p) const
+    {
+      const DimArray& dims = this->getDims();
+      DimArray d(delta_k.size());
+      for (int q = 0; q < d.size(); q++)
+       {
+         d[q] = p%dims[q];
+         p = (p-d[q])/dims[q];
+       }
+      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)
+    {
+      const EuclidianFourierMapComplex<T>& m2 = dynamic_cast<const EuclidianFourierMapComplex<T>&>(other);
+      if (this->size() != m2.size())
+        throw std::bad_cast();
+      
+      const std::complex<T> *d1 = this->data();
+      const std::complex<T> *d2 = m2.data();
+      const DimArray& dims = this->getDims();
+      int N0 = dims[0] + (even0 ? 0 : 1);
+      std::complex<T> result = 0;
+
+      for (long q0 = 1; q0 < N0-1; q0++)
+	{
+	  for (long p = 0; p < plane_size; p++)
+	    {
+	      long idx = q0+dims[0]*p;
+	      assert(idx < this->size());
+	      result += 2*(conj(d1[idx]) * d2[idx]).real();
+	    }
+	}
+      if (even0)
+	{
+	  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;
+    }
+
+  };
+};
+
+#endif

external/cosmotool/src/fourier/details/euclidian_spectrum_1d.hpp

@@ -0,0 +1,180 @@
+#ifndef __DETAILS_EUCLIDIAN_SPECTRUM_1D
+#define __DETAILS_EUCLIDIAN_SPECTRUM_1D
+
+#include  <iostream>
+#include <boost/function.hpp>
+
+
+namespace CosmoTool
+{
+  template<typename T>
+  class EuclidianOperator
+  {
+  public:
+    typedef boost::function1<T, T> Function;
+    
+    Function base, op;
+    T operator()(T k) {
+      return op(base(k));
+    }
+  };
+
+  template<typename T>
+  class EuclidianSpectrum_1D: public SpectrumFunction<T>
+  {
+  public:
+    typedef boost::function1<T, T> Function;
+  protected:
+    Function f;
+
+    static T msqrt(T a) { return std::sqrt(a); }
+  public:
+    typedef typename SpectrumFunction<T>::FourierMapType FourierMapType;
+    typedef typename SpectrumFunction<T>::SpectrumFunctionPtr SpectrumFunctionPtr;
+    typedef boost::shared_ptr<FourierMapType> ptr_map;
+
+    EuclidianSpectrum_1D(Function P)
+      : f(P)
+    {
+    }
+
+    void newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const;
+
+    SpectrumFunctionPtr copy() const {
+      return SpectrumFunctionPtr(new EuclidianSpectrum_1D(f));
+    }
+
+    void sqrt() {
+      EuclidianOperator<T> o;
+      o.base = f;
+      o.op = &EuclidianSpectrum_1D<T>::msqrt;
+      f = (Function(o));
+    }
+
+    void mul(FourierMapType& m) const;
+    void mul_sqrt(FourierMapType& m) const;
+    void mul_inv(FourierMapType& m) const;
+    void mul_inv_sqrt(FourierMapType& m) const;
+  };
+
+
+  template<typename T>
+  void EuclidianSpectrum_1D<T>::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const
+  {
+    typedef EuclidianFourierMapComplex<T> MapT;
+    typedef typename EuclidianSpectrum_1D<T>::ptr_map ptr_map;
+    typedef typename MapT::DimArray DimArray;
+
+    MapT& rand_map = dynamic_cast<MapT&>(out_map);
+
+    std::complex<T> *d = rand_map.data();
+    long idx;
+    const DimArray& dims = rand_map.getDims();
+    const std::vector<double>& delta_k = rand_map.get_delta_k();
+    long plane_size;
+    bool alleven = rand_map.allDimensionsEven();
+    double V = 1;
+
+    for (int p = 0; p < delta_k.size(); p++)
+      V *= (2*M_PI/delta_k[p]);
+
+    for (long p = 1; p < rand_map.size(); p++)
+      {
+	double A_k = std::sqrt(0.5*V*f(rand_map.get_K(p)));
+	d[p] = std::complex<T>(gsl_ran_gaussian(rng, A_k),
+			       gsl_ran_gaussian(rng, A_k));
+      }
+    // Generate the mean value
+    d[0] = std::complex<T>(gsl_ran_gaussian(rng, std::sqrt(V*f(0))), 0);
+
+    if (!rand_map.firstDimensionEven())
+      return;
+
+    // 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)
+      return;
+    
+    plane_size = 1;
+    for (int q = 1; q < dims.size(); q++)
+      {
+	plane_size *= dims[q];
+      }
+
+    for (long p = 1; p < plane_size/2; p++)
+      {
+	long q = (p+1)*dims[0]-1;
+	long q2 = (plane_size-p+1)*dims[0]-1;
+	assert(q < plane_size*dims[0]);
+	assert(q2 < plane_size*dims[0]);
+	d[q] = conj(d[q2]);
+      }
+
+    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);
+      }
+  }
+
+  template<typename T>
+  void EuclidianSpectrum_1D<T>::mul(FourierMapType& m) const
+  {
+    EuclidianFourierMapComplex<T>& m_c = dynamic_cast<EuclidianFourierMapComplex<T>&>(m);
+    std::complex<T> *d = m.data();
+
+    for (long p = 0; p < m_c.size(); p++)
+      d[p] *= f(m_c.get_K(p));
+  }
+
+  template<typename T>
+  void EuclidianSpectrum_1D<T>::mul_sqrt(FourierMapType& m) const
+  {
+    EuclidianFourierMapComplex<T>& m_c = dynamic_cast<EuclidianFourierMapComplex<T>&>(m);
+    std::complex<T> *d = m.data();
+
+    for (long p = 0; p < m_c.size(); p++)
+      d[p] *= std::sqrt(f(m_c.get_K(p)));
+  }
+
+  template<typename T>
+  void EuclidianSpectrum_1D<T>::mul_inv(FourierMapType& m) const
+  {
+    EuclidianFourierMapComplex<T>& m_c = dynamic_cast<EuclidianFourierMapComplex<T>&>(m);
+    std::complex<T> *d = m.data();
+
+    for (long p = 0; p < m_c.size(); p++)
+     {
+        T A = f(m_c.get_K(p));
+        if (A==0)
+          d[p] = 0;
+        else
+          d[p] /= A;
+     }
+  }
+
+  template<typename T>
+  void EuclidianSpectrum_1D<T>::mul_inv_sqrt(FourierMapType& m) const
+  {
+    EuclidianFourierMapComplex<T>& m_c = dynamic_cast<EuclidianFourierMapComplex<T>&>(m);
+    std::complex<T> *d = m.data();
+
+    for (long p = 0; p < m_c.size(); p++)
+      {
+        T A = std::sqrt(f(m_c.get_K(p)));
+        if (A == 0)
+          d[p] = 0;
+        else
+          d[p] /= A;
+      }
+  }
+
+};
+
+
+#endif

external/cosmotool/src/fourier/details/euclidian_spectrum_1d_bin.hpp

@@ -0,0 +1,65 @@
+#ifndef __DETAILS_EUCLIDIAN_SPECTRUM_1D_BIN
+#define __DETAILS_EUCLIDIAN_SPECTRUM_1D_BIN
+
+#include <boost/bind.hpp>
+#include <cmath>
+
+namespace CosmoTool
+{
+  template<typename T>
+  class EuclidianSpectrum_1D_Binned: public EuclidianSpectrum_1D<T>
+  {
+  protected:
+    T *m_data;
+    long m_size;
+    T m_kmin, m_kmax, m_logdeltak;
+    std::vector<T> m_dof;
+  public:
+    typedef typename SpectrumFunction<T>::FourierMapType FourierMapType;
+    typedef typename SpectrumFunction<T>::SpectrumFunctionPtr SpectrumFunctionPtr;
+    typedef boost::shared_ptr<FourierMapType> ptr_map;
+
+    T interpolate_spectrum(T k)
+    {
+      T q = std::log(k/m_kmin)/m_logdeltak;
+      long ik = std::floor(q);
+
+      if (ik >= m_size-1)
+	return m_data[m_size-1];
+      else if (ik < 0)
+	return k/m_kmin*m_data[0];
+
+      return std::exp((q-ik)*m_data[ik+1] + (1-(q-ik))*m_data[ik]);
+    }
+
+    EuclidianSpectrum_1D_Binned(int numBin, T kmin, T kmax)
+      : EuclidianSpectrum_1D<T>(boost::bind(&EuclidianSpectrum_1D_Binned::interpolate_spectrum, this, _1))
+    {
+      m_data = new T[numBin];
+      m_kmin = kmin;
+      m_kmax = kmax;
+      m_size = numBin;
+      m_logdeltak = std::log(m_kmax/m_kmin);
+    }
+
+    SpectrumFunctionPtr copy() const
+    {
+      EuclidianSpectrum_1D_Binned *s = new EuclidianSpectrum_1D_Binned(m_size, m_kmin, m_kmax);
+      std::copy(m_data, m_data+m_size, s->m_data);
+      return SpectrumFunctionPtr(s);
+    }
+
+    void set_dof(std::vector<T>& dof_array) 
+    {
+      assert(m_size == dof_array.size());
+      m_dof = dof_array;
+    }
+
+    const T *data() const { return m_data; }
+    long size() const { return m_size; }
+    const T *dof() const { return &m_dof[0]; }
+  };
+
+};
+
+#endif

external/cosmotool/src/fourier/details/euclidian_transform.hpp

@@ -0,0 +1,181 @@
+#ifndef __DETAILS_EUCLIDIAN_TRANSFORM
+#define __DETAILS_EUCLIDIAN_TRANSFORM
+
+namespace CosmoTool
+{
+  
+  template<typename T>
+  class EuclidianFourierTransform: public FourierTransform<T>
+  {
+  public:
+    typedef typename EuclidianFourierMapBase<T>::DimArray DimArray;
+  private:
+    typedef FFTW_Calls<T> calls;
+    EuclidianFourierMapReal<T> *realMap;
+    EuclidianFourierMapComplex<T> *fourierMap;
+    typename calls::plan_type m_analysis, m_synthesis;
+    double volume;
+    long N, Nc;
+    DimArray m_dims, m_dims_hc;
+    std::vector<double> m_L;
+  public:
+    EuclidianFourierTransform(const DimArray& dims, const std::vector<double>& L)
+    {
+      assert(L.size() == dims.size());
+      std::vector<double> dk(L.size());
+      std::vector<int> swapped_dims(dims.size());
+      
+      m_dims = dims;
+      m_dims_hc = dims;
+      m_dims_hc[0] = dims[0]/2+1;
+      m_L = L;
+
+      N = 1;
+      Nc = 1;
+      volume = 1;
+      for (int i = 0; i < dims.size(); i++)
+	{
+	  N *= dims[i];
+          Nc *= m_dims_hc[i];
+	  volume *= L[i];
+	  dk[i] = 2*M_PI/L[i];
+          swapped_dims[dims.size()-1-i] = dims[i];
+	}
+
+      realMap = new EuclidianFourierMapReal<T>(
+		   boost::shared_ptr<T>(calls::alloc_real(N),
+					std::ptr_fun(calls::free)),
+		   m_dims);
+      fourierMap = new EuclidianFourierMapComplex<T>(
+		   boost::shared_ptr<std::complex<T> >((std::complex<T>*)calls::alloc_complex(Nc),
+						       std::ptr_fun(calls::free)), 
+		   dims[0], m_dims_hc, dk);
+      m_analysis = calls::plan_dft_r2c(dims.size(), &swapped_dims[0],
+				       realMap->data(), (typename calls::complex_type *)fourierMap->data(),
+				       FFTW_DESTROY_INPUT|FFTW_MEASURE);
+      m_synthesis = calls::plan_dft_c2r(dims.size(), &swapped_dims[0],
+					(typename calls::complex_type *)fourierMap->data(), realMap->data(),
+					FFTW_DESTROY_INPUT|FFTW_MEASURE);
+    }
+    
+    virtual ~EuclidianFourierTransform()
+    {
+      delete realMap;
+      delete fourierMap;
+      calls::destroy_plan(m_synthesis);
+      calls::destroy_plan(m_analysis);
+    }
+
+    void synthesis()
+    {
+      calls::execute(m_synthesis);
+      realMap->scale(1/volume);
+    }
+    
+    void analysis()
+    {
+      calls::execute(m_analysis);
+      fourierMap->scale(volume/N);
+    }
+    
+    void synthesis_conjugate()
+    {
+      calls::execute(m_analysis);
+      fourierMap->scale(1/volume);
+    }
+    
+    void analysis_conjugate()
+    {
+      calls::execute(m_synthesis);
+      realMap->scale(volume/N);
+    }
+
+    const FourierMap<std::complex<T> >& fourierSpace() const
+    {
+      return *fourierMap;
+    }
+
+    FourierMap<std::complex<T> >& fourierSpace()
+    {
+      return *fourierMap;
+    }
+
+    const FourierMap<T>& realSpace() const
+    {
+      return *realMap;
+    }
+
+    FourierMap<T>& realSpace()
+    {
+      return *realMap;
+    }
+
+    FourierTransform<T> *mimick() const
+    {
+      return new EuclidianFourierTransform(m_dims, m_L);
+    }
+  };
+
+  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>
+  {
+  private:
+    template<typename T2>
+    static std::vector<T2> make_2d_vector(T2 a, T2 b)
+    {
+      T2 arr[2] = { a, b};
+      return std::vector<T2>(&arr[0], &arr[2]);
+    }
+  public:
+    EuclidianFourierTransform_2d(int Nx, int Ny, double Lx, double Ly)
+      : EuclidianFourierTransform<T>(make_2d_vector<int>(Nx, Ny), make_2d_vector<double>(Lx, Ly))
+    {
+    }
+
+    virtual ~EuclidianFourierTransform_2d() {}
+
+  };
+
+  template<typename T>
+  class EuclidianFourierTransform_3d: public EuclidianFourierTransform<T>
+  {
+  private:
+    template<typename T2>
+    static std::vector<T2> make_3d_vector(T2 a, T2 b, T2 c)
+    {
+      T2 arr[3] = { a, b, c};
+      return std::vector<T2>(&arr[0], &arr[3]);
+    }
+
+  public:
+    EuclidianFourierTransform_3d(int Nx, int Ny, int Nz, double Lx, double Ly, double Lz)
+      : EuclidianFourierTransform<T>(make_3d_vector<int>(Nx, Ny, Nz), make_3d_vector<double>(Lx, Ly, Lz))
+    {
+    }
+
+    virtual ~EuclidianFourierTransform_3d() {}
+  };
+
+};
+
+#endif

external/cosmotool/src/fourier/details/healpix_alms.hpp

@@ -0,0 +1,76 @@
+#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_ALM_HPP
+#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_ALM_HPP
+
+namespace CosmoTool
+{
+  template<typename T>
+  class HealpixFourierALM: public FourierMap<std::complex<T> >
+  {
+  private:
+    std::complex<T> *alms;
+    long m_size;
+    long Lmax_, Mmax_, TVal_;
+    Eigen::aligned_allocator<std::complex<T> > alloc;
+  public:
+    typedef unsigned long LType;
+
+    LType Lmax() const { return Lmax_; }
+    LType Mmax() const { return Mmax_; }
+
+    LType Num_Alms() const
+    {
+      return ((Mmax_+1)*(Mmax_+2))/2 + (Mmax_+1)*(Lmax_-Mmax_);
+    }
+
+    LType index_l0(LType m) const
+    {
+      return ((m*(TVal_-m))/2);
+    }
+
+    LType index(LType l, LType m) const
+    {
+      return index_l0(m) + l;
+    }
+    
+    HealpixFourierALM(LType lmax, LType mmax)
+      : Lmax_(lmax), Mmax_(mmax), TVal_(2*lmax+1)
+    {
+      m_size = Num_Alms();
+      alms = alloc.allocate(m_size);
+    }
+
+    virtual ~HealpixFourierALM()
+    {
+      alloc.deallocate(alms, m_size);
+    }
+
+    virtual const std::complex<T>* data() const { return alms; }
+    virtual std::complex<T> * data() { return alms;} 
+    virtual long size() const { return m_size; }
+
+    virtual FourierMap<std::complex<T> > *mimick() const
+    {
+      return new HealpixFourierALM<T>(Lmax_, Mmax_);
+    }
+
+    virtual std::complex<T> dot_product(const FourierMap<std::complex<T> >& other) const
+      throw(std::bad_cast)
+    {
+      const HealpixFourierALM<T>& mfm = dynamic_cast<const HealpixFourierALM<T>&>(other);
+      typedef typename FourierMap<std::complex<T> >::MapType MapType;
+      std::complex<T> S;
+
+      if (m_size != mfm.m_size)
+        throw std::bad_cast();
+
+      MapType m1(alms, m_size);
+      MapType m2(mfm.alms, mfm.m_size);
+      
+      S = (m1.block(0,0,1,Lmax_+1).conjugate() * m2.block(0,0,1,Lmax_+1)).sum();
+      S += std::complex<T>(2,0)*(m1.block(0,1+Lmax_,1,m_size-1-Lmax_).conjugate() * m2.block(0,1+Lmax_,1,m_size-1-Lmax_)).sum();
+      return S;
+    }
+  };
+};
+
+#endif

external/cosmotool/src/fourier/details/healpix_map.hpp

@@ -0,0 +1,54 @@
+#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_MAP_HPP
+#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_MAP_HPP
+
+namespace CosmoTool
+{
+
+  template<typename T>
+  class HealpixFourierMap: public FourierMap<T>
+  {
+  private:
+    T *m_data;
+    long Npix, m_Nside;
+    Eigen::aligned_allocator<T> alloc;
+  public:
+    HealpixFourierMap(long nSide)
+      : Npix(12*nSide*nSide), m_Nside(nSide)
+    {
+      m_data = alloc.allocate(Npix);
+    }
+
+    virtual ~HealpixFourierMap()
+    {
+      alloc.deallocate(m_data, Npix);
+    }
+
+    long Nside() const { return m_Nside; }
+    virtual const T* data() const { return m_data; }
+    virtual T *data() { return m_data; }
+    virtual long size() const { return Npix; }
+
+    virtual T dot_product(const FourierMap<T>& other) const
+      throw(std::bad_cast)
+    {
+      typedef typename FourierMap<T>::MapType MapType;
+
+      const HealpixFourierMap<T>& mfm = dynamic_cast<const HealpixFourierMap<T>&>(other);
+      if (Npix != mfm.size())
+        throw std::bad_cast();
+
+      MapType m1(m_data, Npix);
+      MapType m2(mfm.m_data, mfm.Npix);
+      
+      return (m1*m2).sum();
+    }
+    
+    virtual FourierMap<T> *mimick() const
+    {
+      return new HealpixFourierMap<T>(m_Nside);
+    }
+  };
+
+};
+
+#endif

external/cosmotool/src/fourier/details/healpix_spectrum.hpp

@@ -0,0 +1,144 @@
+#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HPP
+#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HPP
+
+namespace CosmoTool
+{
+  template<typename T>
+  class HealpixSpectrum: public SpectrumFunction<T>
+  {
+  protected:
+    std::vector<T> cls;
+    int *m_dof;
+  public:
+    typedef typename SpectrumFunction<T>::FourierMapType FourierMapType;
+    typedef boost::shared_ptr<FourierMapType> ptr_map;
+    typedef typename SpectrumFunction<T>::SpectrumFunctionPtr SpectrumFunctionPtr;
+    typedef unsigned long LType;
+
+    HealpixSpectrum(LType Lmax)
+      : cls(Lmax+1), m_dof(new int[Lmax+1]) 
+    {
+      for (int l = 0; l <= Lmax; l++)
+	m_dof[l] = l + 1;
+    }
+
+    T *data() { return &cls[0]; }
+    const T *data() const { return &cls[0]; }
+    const int *dof() const { return m_dof; } 
+    
+    void set_dof(LType l, int dof) { m_dof[l] = dof; }
+
+    LType Lmax() const { return cls.size()-1; }
+    long size() const { return cls.size(); }
+
+    void newRandomFourier(gsl_rng *rng, FourierMapType& like_map) const; 
+
+    SpectrumFunctionPtr copy() const {
+      HealpixSpectrum<T> *s = new HealpixSpectrum<T>(Lmax());
+      s->cls = cls;
+      return SpectrumFunctionPtr(s);
+    }
+
+    void sqrt() {
+      std::transform(cls.begin(), cls.end(), cls.begin(), std::ptr_fun<T,T>(std::sqrt));
+    }
+
+    void mul(FourierMapType& m) const;
+    void mul_sqrt(FourierMapType& m) const;
+    void mul_inv(FourierMapType& m) const;
+    void mul_inv_sqrt(FourierMapType& m) const;
+  };
+
+  template<typename T>
+  void HealpixSpectrum<T>::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const
+  { 
+    HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(out_map);
+    long lmaxGen = std::min(cls.size()-1, alms.Lmax());
+    std::complex<T> *new_data = alms.data();
+
+    for (LType l = 0; l <= lmaxGen; l++)
+      {
+        double Al = std::sqrt(cls[l]);
+
+        new_data[alms.index(l,0)] = gsl_ran_gaussian(rng, Al);
+        Al *= M_SQRT1_2;
+        for (LType m = 1; m <= std::min(l,alms.Mmax()); m++)
+          {
+            std::complex<T>& c = new_data[alms.index(l,m)];
+            c.real() = gsl_ran_gaussian(rng, Al); 
+            c.imag() = gsl_ran_gaussian(rng, Al); 
+          }
+      }
+  } 
+
+ template<typename T>
+ void HealpixSpectrum<T>::mul(FourierMapType& like_map) const
+  {
+    HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(like_map);
+    std::complex<T> *data = alms.data();
+
+    for (LType l = 0; l <= alms.Lmax(); l++)
+      {
+        double Al = cls[l];
+
+        for (LType m = 0; m <= std::min(l,alms.Mmax()); m++)
+          {
+            data[alms.index(l,m)] *= Al;
+          }
+      }
+  }
+
+ template<typename T>
+ void HealpixSpectrum<T>::mul_sqrt(FourierMapType& like_map) const
+  {
+    HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(like_map);
+    std::complex<T> *data = alms.data();
+
+    for (LType l = 0; l <= alms.Lmax(); l++)
+      {
+        double Al = std::sqrt(cls[l]);
+
+        for (LType m = 0; m <= std::min(l,alms.Mmax()); m++)
+          {
+            data[alms.index(l,m)] *= Al;
+          }
+      }
+  }
+
+ template<typename T>
+ void HealpixSpectrum<T>::mul_inv(FourierMapType& like_map) const
+  {
+    HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(like_map);
+    std::complex<T> *data = alms.data();
+
+    for (LType l = 0; l <= alms.Lmax(); l++)
+      {
+        double Al = (cls[l] <= 0) ? 0 : (1/cls[l]);
+
+        for (LType m = 0; m <= std::min(l,alms.Mmax()); m++)
+          {
+            data[alms.index(l,m)] *= Al;
+          }
+      }
+  }
+
+ template<typename T>
+ void HealpixSpectrum<T>::mul_inv_sqrt(FourierMapType& like_map) const
+  {
+    HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(like_map);
+    std::complex<T> *data = alms.data();
+
+    for (LType l = 0; l <= alms.Lmax(); l++)
+      {
+        double Al = (cls[l] <= 0) ? 0 : std::sqrt(1/cls[l]);
+
+        for (LType m = 0; m <= std::min(l,alms.Mmax()); m++)
+          {
+            data[alms.index(l,m)] *= Al;
+          }
+      }
+  }
+
+};
+
+#endif

external/cosmotool/src/fourier/details/healpix_transform.hpp

@@ -0,0 +1,97 @@
+#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_TRANSFORM_HPP
+#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_TRANSFORM_HPP
+
+namespace CosmoTool
+{
+
+  template<typename T> struct HealpixJobHelper__ {};
+
+  template<> struct HealpixJobHelper__<double>
+    { enum {val=1}; };
+
+  template<> struct HealpixJobHelper__<float>
+    { enum {val=0}; };
+
+
+  template<typename T>
+  class HealpixFourierTransform: public FourierTransform<T>
+  {
+  private:
+    sharp_alm_info *ainfo;
+    sharp_geom_info *ginfo;
+    HealpixFourierMap<T> realMap;
+    HealpixFourierALM<T> fourierMap;
+    int m_iterate;
+  public:
+    HealpixFourierTransform(long nSide, long Lmax, long Mmax, int iterate = 0)
+      : realMap(nSide), fourierMap(Lmax, Mmax), ainfo(0), ginfo(0), m_iterate(iterate)
+    {
+      sharp_make_healpix_geom_info (nSide, 1, &ginfo);
+      sharp_make_triangular_alm_info (Lmax, Mmax, 1, &ainfo);
+    }
+
+    virtual ~HealpixFourierTransform()
+    {
+      sharp_destroy_geom_info(ginfo);
+      sharp_destroy_alm_info(ainfo);
+    }
+
+    virtual const FourierMap<std::complex<T> >& fourierSpace() const { return fourierMap; }
+
+    virtual FourierMap<std::complex<T> >& fourierSpace() { return fourierMap; }
+
+    virtual const FourierMap<T>& realSpace() const { return realMap; }
+
+    virtual FourierMap<T>& realSpace() { return realMap; }
+    
+    virtual FourierTransform<T> *mimick() const 
+    {
+      return new HealpixFourierTransform<T>(realMap.Nside(), fourierMap.Lmax(), fourierMap.Mmax());
+    }
+
+    virtual void analysis()
+    {
+      void *aptr=reinterpret_cast<void *>(fourierMap.data()), *mptr=reinterpret_cast<void *>(realMap.data());
+
+      sharp_execute (SHARP_MAP2ALM, 0, 0, &aptr, &mptr, ginfo, ainfo, 1,
+        HealpixJobHelper__<T>::val,0,0,0);
+      for (int i = 0; i < m_iterate; i++)
+        {
+          HealpixFourierMap<T> tmp_map(realMap.Nside());          
+          void *tmp_ptr=reinterpret_cast<void *>(tmp_map.data());
+          typename HealpixFourierMap<T>::MapType m0 = tmp_map.eigen();
+          typename HealpixFourierMap<T>::MapType m1 = realMap.eigen();
+
+          sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &tmp_ptr, ginfo, ainfo, 1,
+            HealpixJobHelper__<T>::val,0,0,0);
+          m0 = m1 - m0;
+          sharp_execute (SHARP_MAP2ALM, 0, 1, &aptr, &tmp_ptr, ginfo, ainfo, 1,
+            HealpixJobHelper__<T>::val,0,0,0);
+        }
+    }
+
+    virtual void synthesis()
+    {
+      void *aptr=reinterpret_cast<void *>(fourierMap.data()), *mptr=reinterpret_cast<void *>(realMap.data());
+
+      sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &mptr, ginfo, ainfo, 1,
+        HealpixJobHelper__<T>::val,0,0,0);
+    }
+
+    virtual void analysis_conjugate()
+    {
+      synthesis();
+      realMap.scale(4*M_PI/realMap.size());
+    }
+
+    virtual void synthesis_conjugate()
+    {
+      analysis();
+      fourierMap.scale(realMap.size()/(4*M_PI));
+    }
+    
+  };
+
+};
+
+#endif

external/cosmotool/src/fourier/details/healpix_utility.hpp

@@ -0,0 +1,63 @@
+#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HP
+#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HP
+
+namespace CosmoTool
+{
+
+  template<typename T>
+  class HealpixUtilityFunction: public MapUtilityFunction<T>
+  {
+  public:
+    typedef typename MapUtilityFunction<T>::Spectrum Spectrum;
+    typedef typename MapUtilityFunction<T>::Spectrum_ptr Spectrum_ptr;
+    typedef typename MapUtilityFunction<T>::FMap FMap;
+    typedef typename HealpixSpectrum<T>::LType LType; 
+    
+    Spectrum_ptr estimateSpectrumFromMap(const FMap& m) const
+    {
+      const HealpixFourierALM<T>& alms = dynamic_cast<const HealpixFourierALM<T>&>(m);
+      LType Lmax = alms.Lmax(), Mmax = alms.Mmax();
+      HealpixSpectrum<T> *spectrum = new HealpixSpectrum<T>(Lmax);
+      T *data_cls = spectrum->data();
+      const std::complex<T> *data_alms = alms.data();
+
+      // make an estimate of the cls
+      std::fill(data_cls, data_cls+Lmax+1, 0);
+      LType q = 0;
+      for (LType m = 0; m <= Mmax; ++m )
+	{
+	  for (LType l = m; l <= Lmax; ++l)
+	    {
+	      // Triangular storage
+	      data_cls[l] += std::norm(data_alms[l+q]);
+	    }
+	  q += Lmax+1-m;
+	}
+      assert(q==alms.size());
+      
+      for (LType l = 0; l <= Lmax; ++l)
+	{
+	  int dof = 1 + std::min(l, Mmax);
+	  spectrum->set_dof(l, dof);
+	  data_cls[l] /= dof;
+	}
+
+      return Spectrum_ptr(spectrum);
+    }
+
+    Spectrum_ptr newSpectrumFromRaw(T *data, long size,
+				    const Spectrum& like_spec) const
+    {
+      const HealpixSpectrum<T>& in_spec = dynamic_cast<const HealpixSpectrum<T>&>(like_spec);
+      HealpixSpectrum<T> *new_spectrum = new HealpixSpectrum<T>(in_spec.Lmax());
+      T *out_d = new_spectrum->data();
+
+      std::copy(data, data + min(size,new_spectrum->size()), out_d);
+
+      return Spectrum_ptr(new_spectrum);
+    }
+  };
+
+};
+
+#endif