From e41c26a46c58630d06ffa7d33eebe8a9b9e74d51 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 19 Feb 2013 11:18:20 -0500 Subject: [PATCH] Imported latest snapshot of cosmotoolbox --- .../cosmotool/sample/test_healpix_calls.cpp | 25 ++ external/cosmotool/src/cosmopower.cpp | 306 ++++++++++++++++++ external/cosmotool/src/cosmopower.hpp | 74 +++++ .../src/fourier/details/euclidian_maps.hpp | 232 +++++++++++++ .../fourier/details/euclidian_spectrum_1d.hpp | 180 +++++++++++ .../details/euclidian_spectrum_1d_bin.hpp | 65 ++++ .../fourier/details/euclidian_transform.hpp | 181 +++++++++++ .../src/fourier/details/healpix_alms.hpp | 76 +++++ .../src/fourier/details/healpix_map.hpp | 54 ++++ .../src/fourier/details/healpix_spectrum.hpp | 144 +++++++++ .../src/fourier/details/healpix_transform.hpp | 97 ++++++ .../src/fourier/details/healpix_utility.hpp | 63 ++++ 12 files changed, 1497 insertions(+) create mode 100644 external/cosmotool/sample/test_healpix_calls.cpp create mode 100644 external/cosmotool/src/cosmopower.cpp create mode 100644 external/cosmotool/src/cosmopower.hpp create mode 100644 external/cosmotool/src/fourier/details/euclidian_maps.hpp create mode 100644 external/cosmotool/src/fourier/details/euclidian_spectrum_1d.hpp create mode 100644 external/cosmotool/src/fourier/details/euclidian_spectrum_1d_bin.hpp create mode 100644 external/cosmotool/src/fourier/details/euclidian_transform.hpp create mode 100644 external/cosmotool/src/fourier/details/healpix_alms.hpp create mode 100644 external/cosmotool/src/fourier/details/healpix_map.hpp create mode 100644 external/cosmotool/src/fourier/details/healpix_spectrum.hpp create mode 100644 external/cosmotool/src/fourier/details/healpix_transform.hpp create mode 100644 external/cosmotool/src/fourier/details/healpix_utility.hpp diff --git a/external/cosmotool/sample/test_healpix_calls.cpp b/external/cosmotool/sample/test_healpix_calls.cpp new file mode 100644 index 0000000..7cc4272 --- /dev/null +++ b/external/cosmotool/sample/test_healpix_calls.cpp @@ -0,0 +1,25 @@ +#include +#include "fourier/healpix.hpp" + +using namespace CosmoTool; +using namespace std; + +int main() +{ + HealpixFourierTransform dft(128,2*128,2*128, 40); + HealpixUtilityFunction 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::Spectrum_ptr s = utils.estimateSpectrumFromMap(dft.fourierSpace()); + + dft.synthesis(); + cout << "Resynthesis dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl; + + return 0; + +} diff --git a/external/cosmotool/src/cosmopower.cpp b/external/cosmotool/src/cosmopower.cpp new file mode 100644 index 0000000..19aec97 --- /dev/null +++ b/external/cosmotool/src/cosmopower.cpp @@ -0,0 +1,306 @@ +#include +#include +#include +#include +#include +#include +#include +#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); +} diff --git a/external/cosmotool/src/cosmopower.hpp b/external/cosmotool/src/cosmopower.hpp new file mode 100644 index 0000000..1edf5b6 --- /dev/null +++ b/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 diff --git a/external/cosmotool/src/fourier/details/euclidian_maps.hpp b/external/cosmotool/src/fourier/details/euclidian_maps.hpp new file mode 100644 index 0000000..a7c51b5 --- /dev/null +++ b/external/cosmotool/src/fourier/details/euclidian_maps.hpp @@ -0,0 +1,232 @@ +#ifndef __DETAILS_EUCLIDIAN_MAPS +#define __DETAILS_EUCLIDIAN_MAPS + + +namespace CosmoTool +{ + + template + class EuclidianFourierMapBase: public FourierMap + { + public: + typedef std::vector DimArray; + private: + boost::shared_ptr m_data; + DimArray m_dims; + long m_size; + public: + + EuclidianFourierMapBase(boost::shared_ptr 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 *copy() const + { + FourierMap *m = this->mimick(); + m->eigen() = this->eigen(); + return m; + } + + }; + + template + class EuclidianFourierMapReal: public EuclidianFourierMapBase + { + public: + typedef typename EuclidianFourierMapBase::DimArray DimArray; + + EuclidianFourierMapReal(boost::shared_ptr indata, const DimArray& indims) + : EuclidianFourierMapBase(indata, indims) + {} + + virtual FourierMap *mimick() const + { + return new EuclidianFourierMapReal( + boost::shared_ptr((T *)fftw_malloc(sizeof(T)*this->size()), + std::ptr_fun(fftw_free)), + this->getDims()); + } + + virtual T dot_product(const FourierMap& other) const + throw(std::bad_cast) + { + const EuclidianFourierMapReal& m2 = dynamic_cast&>(other); + if (this->size() != m2.size()) + throw std::bad_cast(); + + return (this->eigen()*m2.eigen()).sum(); + } + }; + + template + class EuclidianFourierMapComplex: public EuclidianFourierMapBase > + { + protected: + typedef boost::shared_ptr > ptr_t; + std::vector delta_k; + int m_dim0; + bool even0, alleven; + long plane_size; + public: + typedef typename EuclidianFourierMapBase >::DimArray DimArray; + + EuclidianFourierMapComplex(ptr_t indata, + int dim0, + const DimArray& indims, + const std::vector& dk) + : EuclidianFourierMapBase >(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 > *mimick() const + { + return + new EuclidianFourierMapComplex( + ptr_t((std::complex *) + fftw_malloc(sizeof(std::complex)*this->size()), + std::ptr_fun(fftw_free)), + m_dim0, + this->getDims(), + this->delta_k); + } + + const std::vector& get_delta_k() const + { + return this->delta_k; + } + + template + 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 + 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 + 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 dot_product(const FourierMap >& other) const + throw(std::bad_cast) + { + const EuclidianFourierMapComplex& m2 = dynamic_cast&>(other); + if (this->size() != m2.size()) + throw std::bad_cast(); + + const std::complex *d1 = this->data(); + const std::complex *d2 = m2.data(); + const DimArray& dims = this->getDims(); + int N0 = dims[0] + (even0 ? 0 : 1); + std::complex 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 diff --git a/external/cosmotool/src/fourier/details/euclidian_spectrum_1d.hpp b/external/cosmotool/src/fourier/details/euclidian_spectrum_1d.hpp new file mode 100644 index 0000000..fbc4f7f --- /dev/null +++ b/external/cosmotool/src/fourier/details/euclidian_spectrum_1d.hpp @@ -0,0 +1,180 @@ +#ifndef __DETAILS_EUCLIDIAN_SPECTRUM_1D +#define __DETAILS_EUCLIDIAN_SPECTRUM_1D + +#include +#include + + +namespace CosmoTool +{ + template + class EuclidianOperator + { + public: + typedef boost::function1 Function; + + Function base, op; + T operator()(T k) { + return op(base(k)); + } + }; + + template + class EuclidianSpectrum_1D: public SpectrumFunction + { + public: + typedef boost::function1 Function; + protected: + Function f; + + static T msqrt(T a) { return std::sqrt(a); } + public: + typedef typename SpectrumFunction::FourierMapType FourierMapType; + typedef typename SpectrumFunction::SpectrumFunctionPtr SpectrumFunctionPtr; + typedef boost::shared_ptr 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 o; + o.base = f; + o.op = &EuclidianSpectrum_1D::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 + void EuclidianSpectrum_1D::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const + { + typedef EuclidianFourierMapComplex MapT; + typedef typename EuclidianSpectrum_1D::ptr_map ptr_map; + typedef typename MapT::DimArray DimArray; + + MapT& rand_map = dynamic_cast(out_map); + + std::complex *d = rand_map.data(); + long idx; + const DimArray& dims = rand_map.getDims(); + const std::vector& 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(gsl_ran_gaussian(rng, A_k), + gsl_ran_gaussian(rng, A_k)); + } + // Generate the mean value + d[0] = std::complex(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(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(d[q].real()+d[q].imag(),0); + } + } + + template + void EuclidianSpectrum_1D::mul(FourierMapType& m) const + { + EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); + std::complex *d = m.data(); + + for (long p = 0; p < m_c.size(); p++) + d[p] *= f(m_c.get_K(p)); + } + + template + void EuclidianSpectrum_1D::mul_sqrt(FourierMapType& m) const + { + EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); + std::complex *d = m.data(); + + for (long p = 0; p < m_c.size(); p++) + d[p] *= std::sqrt(f(m_c.get_K(p))); + } + + template + void EuclidianSpectrum_1D::mul_inv(FourierMapType& m) const + { + EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); + std::complex *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 + void EuclidianSpectrum_1D::mul_inv_sqrt(FourierMapType& m) const + { + EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); + std::complex *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 diff --git a/external/cosmotool/src/fourier/details/euclidian_spectrum_1d_bin.hpp b/external/cosmotool/src/fourier/details/euclidian_spectrum_1d_bin.hpp new file mode 100644 index 0000000..ff1c98d --- /dev/null +++ b/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 +#include + +namespace CosmoTool +{ + template + class EuclidianSpectrum_1D_Binned: public EuclidianSpectrum_1D + { + protected: + T *m_data; + long m_size; + T m_kmin, m_kmax, m_logdeltak; + std::vector m_dof; + public: + typedef typename SpectrumFunction::FourierMapType FourierMapType; + typedef typename SpectrumFunction::SpectrumFunctionPtr SpectrumFunctionPtr; + typedef boost::shared_ptr 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(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& 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 diff --git a/external/cosmotool/src/fourier/details/euclidian_transform.hpp b/external/cosmotool/src/fourier/details/euclidian_transform.hpp new file mode 100644 index 0000000..8faba32 --- /dev/null +++ b/external/cosmotool/src/fourier/details/euclidian_transform.hpp @@ -0,0 +1,181 @@ +#ifndef __DETAILS_EUCLIDIAN_TRANSFORM +#define __DETAILS_EUCLIDIAN_TRANSFORM + +namespace CosmoTool +{ + + template + class EuclidianFourierTransform: public FourierTransform + { + public: + typedef typename EuclidianFourierMapBase::DimArray DimArray; + private: + typedef FFTW_Calls calls; + EuclidianFourierMapReal *realMap; + EuclidianFourierMapComplex *fourierMap; + typename calls::plan_type m_analysis, m_synthesis; + double volume; + long N, Nc; + DimArray m_dims, m_dims_hc; + std::vector m_L; + public: + EuclidianFourierTransform(const DimArray& dims, const std::vector& L) + { + assert(L.size() == dims.size()); + std::vector dk(L.size()); + std::vector 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( + boost::shared_ptr(calls::alloc_real(N), + std::ptr_fun(calls::free)), + m_dims); + fourierMap = new EuclidianFourierMapComplex( + boost::shared_ptr >((std::complex*)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 >& fourierSpace() const + { + return *fourierMap; + } + + FourierMap >& fourierSpace() + { + return *fourierMap; + } + + const FourierMap& realSpace() const + { + return *realMap; + } + + FourierMap& realSpace() + { + return *realMap; + } + + FourierTransform *mimick() const + { + return new EuclidianFourierTransform(m_dims, m_L); + } + }; + + template + class EuclidianFourierTransform_1d: public EuclidianFourierTransform + { + private: + template + static std::vector make_1d_vector(T2 a) + { + T2 arr[2] = { a}; + return std::vector(&arr[0],&arr[1]); + } + public: + EuclidianFourierTransform_1d(int Nx, double Lx) + : EuclidianFourierTransform(make_1d_vector(Nx), make_1d_vector(Lx)) + { + } + + virtual ~EuclidianFourierTransform_1d() {} + + }; + + template + class EuclidianFourierTransform_2d: public EuclidianFourierTransform + { + private: + template + static std::vector make_2d_vector(T2 a, T2 b) + { + T2 arr[2] = { a, b}; + return std::vector(&arr[0], &arr[2]); + } + public: + EuclidianFourierTransform_2d(int Nx, int Ny, double Lx, double Ly) + : EuclidianFourierTransform(make_2d_vector(Nx, Ny), make_2d_vector(Lx, Ly)) + { + } + + virtual ~EuclidianFourierTransform_2d() {} + + }; + + template + class EuclidianFourierTransform_3d: public EuclidianFourierTransform + { + private: + template + static std::vector make_3d_vector(T2 a, T2 b, T2 c) + { + T2 arr[3] = { a, b, c}; + return std::vector(&arr[0], &arr[3]); + } + + public: + EuclidianFourierTransform_3d(int Nx, int Ny, int Nz, double Lx, double Ly, double Lz) + : EuclidianFourierTransform(make_3d_vector(Nx, Ny, Nz), make_3d_vector(Lx, Ly, Lz)) + { + } + + virtual ~EuclidianFourierTransform_3d() {} + }; + +}; + +#endif diff --git a/external/cosmotool/src/fourier/details/healpix_alms.hpp b/external/cosmotool/src/fourier/details/healpix_alms.hpp new file mode 100644 index 0000000..849781f --- /dev/null +++ b/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 + class HealpixFourierALM: public FourierMap > + { + private: + std::complex *alms; + long m_size; + long Lmax_, Mmax_, TVal_; + Eigen::aligned_allocator > 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* data() const { return alms; } + virtual std::complex * data() { return alms;} + virtual long size() const { return m_size; } + + virtual FourierMap > *mimick() const + { + return new HealpixFourierALM(Lmax_, Mmax_); + } + + virtual std::complex dot_product(const FourierMap >& other) const + throw(std::bad_cast) + { + const HealpixFourierALM& mfm = dynamic_cast&>(other); + typedef typename FourierMap >::MapType MapType; + std::complex 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(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 diff --git a/external/cosmotool/src/fourier/details/healpix_map.hpp b/external/cosmotool/src/fourier/details/healpix_map.hpp new file mode 100644 index 0000000..6b96974 --- /dev/null +++ b/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 + class HealpixFourierMap: public FourierMap + { + private: + T *m_data; + long Npix, m_Nside; + Eigen::aligned_allocator 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& other) const + throw(std::bad_cast) + { + typedef typename FourierMap::MapType MapType; + + const HealpixFourierMap& mfm = dynamic_cast&>(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 *mimick() const + { + return new HealpixFourierMap(m_Nside); + } + }; + +}; + +#endif diff --git a/external/cosmotool/src/fourier/details/healpix_spectrum.hpp b/external/cosmotool/src/fourier/details/healpix_spectrum.hpp new file mode 100644 index 0000000..cfeebc9 --- /dev/null +++ b/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 + class HealpixSpectrum: public SpectrumFunction + { + protected: + std::vector cls; + int *m_dof; + public: + typedef typename SpectrumFunction::FourierMapType FourierMapType; + typedef boost::shared_ptr ptr_map; + typedef typename SpectrumFunction::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 *s = new HealpixSpectrum(Lmax()); + s->cls = cls; + return SpectrumFunctionPtr(s); + } + + void sqrt() { + std::transform(cls.begin(), cls.end(), cls.begin(), std::ptr_fun(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 + void HealpixSpectrum::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const + { + HealpixFourierALM& alms = dynamic_cast&>(out_map); + long lmaxGen = std::min(cls.size()-1, alms.Lmax()); + std::complex *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& c = new_data[alms.index(l,m)]; + c.real() = gsl_ran_gaussian(rng, Al); + c.imag() = gsl_ran_gaussian(rng, Al); + } + } + } + + template + void HealpixSpectrum::mul(FourierMapType& like_map) const + { + HealpixFourierALM& alms = dynamic_cast&>(like_map); + std::complex *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 + void HealpixSpectrum::mul_sqrt(FourierMapType& like_map) const + { + HealpixFourierALM& alms = dynamic_cast&>(like_map); + std::complex *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 + void HealpixSpectrum::mul_inv(FourierMapType& like_map) const + { + HealpixFourierALM& alms = dynamic_cast&>(like_map); + std::complex *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 + void HealpixSpectrum::mul_inv_sqrt(FourierMapType& like_map) const + { + HealpixFourierALM& alms = dynamic_cast&>(like_map); + std::complex *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 diff --git a/external/cosmotool/src/fourier/details/healpix_transform.hpp b/external/cosmotool/src/fourier/details/healpix_transform.hpp new file mode 100644 index 0000000..bc542e7 --- /dev/null +++ b/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 struct HealpixJobHelper__ {}; + + template<> struct HealpixJobHelper__ + { enum {val=1}; }; + + template<> struct HealpixJobHelper__ + { enum {val=0}; }; + + + template + class HealpixFourierTransform: public FourierTransform + { + private: + sharp_alm_info *ainfo; + sharp_geom_info *ginfo; + HealpixFourierMap realMap; + HealpixFourierALM 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 >& 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() + { + void *aptr=reinterpret_cast(fourierMap.data()), *mptr=reinterpret_cast(realMap.data()); + + sharp_execute (SHARP_MAP2ALM, 0, 0, &aptr, &mptr, ginfo, ainfo, 1, + HealpixJobHelper__::val,0,0,0); + for (int i = 0; i < m_iterate; i++) + { + HealpixFourierMap tmp_map(realMap.Nside()); + void *tmp_ptr=reinterpret_cast(tmp_map.data()); + typename HealpixFourierMap::MapType m0 = tmp_map.eigen(); + typename HealpixFourierMap::MapType m1 = realMap.eigen(); + + sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &tmp_ptr, ginfo, ainfo, 1, + HealpixJobHelper__::val,0,0,0); + m0 = m1 - m0; + sharp_execute (SHARP_MAP2ALM, 0, 1, &aptr, &tmp_ptr, ginfo, ainfo, 1, + HealpixJobHelper__::val,0,0,0); + } + } + + virtual void synthesis() + { + void *aptr=reinterpret_cast(fourierMap.data()), *mptr=reinterpret_cast(realMap.data()); + + sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &mptr, ginfo, ainfo, 1, + HealpixJobHelper__::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 diff --git a/external/cosmotool/src/fourier/details/healpix_utility.hpp b/external/cosmotool/src/fourier/details/healpix_utility.hpp new file mode 100644 index 0000000..d989dd8 --- /dev/null +++ b/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 + class HealpixUtilityFunction: public MapUtilityFunction + { + public: + typedef typename MapUtilityFunction::Spectrum Spectrum; + typedef typename MapUtilityFunction::Spectrum_ptr Spectrum_ptr; + typedef typename MapUtilityFunction::FMap FMap; + typedef typename HealpixSpectrum::LType LType; + + Spectrum_ptr estimateSpectrumFromMap(const FMap& m) const + { + const HealpixFourierALM& alms = dynamic_cast&>(m); + LType Lmax = alms.Lmax(), Mmax = alms.Mmax(); + HealpixSpectrum *spectrum = new HealpixSpectrum(Lmax); + T *data_cls = spectrum->data(); + const std::complex *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& in_spec = dynamic_cast&>(like_spec); + HealpixSpectrum *new_spectrum = new HealpixSpectrum(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