WIP transition to monofonic infrastructure

src/TransferFunction.hh

@@ -1,410 +0,0 @@
-/*
- This file is part of MUSIC -
- a tool to generate initial conditions for cosmological simulations
- 
- Copyright (C) 2008-12  Oliver Hahn, ojha@gmx.de
- 
- This program is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
- 
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- GNU General Public License for more details.
- 
- You should have received a copy of the GNU General Public License
- along with this program.  If not, see <http://www.gnu.org/licenses/>.
-*/
-
-#ifndef __TRANSFERFUNCTION_HH
-#define __TRANSFERFUNCTION_HH
-
-#include <vector>
-#include <sstream>
-#include <fstream>
-#include <iostream>
-#include <cmath>
-#include <stdexcept>
-
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_spline.h>
-#include <gsl/gsl_sf_gamma.h>
-
-#include "Numerics.hh"
-#include "general.hh"
-
-#include <complex>
-
-#define NZERO_Q
-
-typedef std::complex<double> complex;
-
-//! Abstract base class for transfer functions
-/*!
-    This class implements a purely virtual interface that can be
-    used to derive instances implementing various transfer functions.
-*/ 
-class TransferFunction{
-public:
-	Cosmology m_Cosmology;
-	
-public:
-  
-	TransferFunction( Cosmology acosm ) : m_Cosmology( acosm ) { };
-	virtual double compute( double k ) = 0;
-	virtual ~TransferFunction(){ };
-	virtual double get_kmax( void ) = 0;
-	virtual double get_kmin( void ) = 0;
-};
-
-class TransferFunction_real
-{
-	
-public:
-	gsl_interp_accel *accp, *accn;
-	gsl_spline *splinep, *splinen;
-	double Tr0_, Tmin_, Tmax_, Tscale_;
-	double rneg_, rneg2_;
-	static TransferFunction *ptf_;
-	static double nspec_;
-	
-protected:
-	
-	double krgood( double mu, double q, double dlnr, double kr )
-	{
-		double krnew = kr;
-		complex cdgamma, zm, zp;
-		double arg, iarg, xm, xp, y;
-		gsl_sf_result g_a, g_p;
-		
-		xp = 0.5*(mu+1.0+q);
-		xm = 0.5*(mu+1.0-q);
-		y = M_PI/(2.0*dlnr);
-		zp=complex(xp,y);
-		zm=complex(xm,y);
-		
-		gsl_sf_lngamma_complex_e (zp.real(), zp.imag(), &g_a, &g_p);
-		zp=std::polar(exp(g_a.val),g_p.val);
-		double zpa = g_p.val;
-
-		gsl_sf_lngamma_complex_e (zm.real(), zm.imag(), &g_a, &g_p);
-		zm=std::polar(exp(g_a.val),g_p.val);
-		double zma = g_p.val;
-		
-		arg=log(2.0/kr)/dlnr+(zpa+zma)/M_PI;
-		iarg=(double)((int)(arg + 0.5));
-		
-		if( arg!=iarg )
-			krnew=kr*exp((arg-iarg)*dlnr);
-		
-		return krnew;
-	}
-	
-	void transform( double pnorm, double dplus, unsigned N, double q, std::vector<double>& rr, std::vector<double>& TT )
-	{
-		const double mu = 0.5;
-		double qmin = 1.0e-6, qmax = 1.0e+6;
-		
-		q = 0.0;
-		
-		N = 16384;
-		
-#ifdef NZERO_Q
-		//q = 0.4;
-		q = 0.2;
-#endif
-		
-		double kmin = qmin, kmax=qmax;
-		double rmin = qmin, rmax = qmax;
-		double k0 = exp(0.5*(log(kmax)+log(kmin)));
-		double r0 = exp(0.5*(log(rmax)+log(rmin)));
-		double L = log(rmax)-log(rmin);
-		double k0r0 = k0*r0;
-		double dlnk = L/N, dlnr = L/N;
-		
-		double sqrtpnorm = sqrt(pnorm);
-		
-		double dir = 1.0;
-		
-		double fftnorm = 1.0/N;
-		
-		complex_t in[N], out[N];
-		fftw_plan p,ip;
-		
-		//... perform anti-ringing correction from Hamilton (2000)
-		k0r0 = krgood( mu, q, dlnr, k0r0 );
-
-		std::ofstream ofsk("transfer_k.txt");
-		double sum_in = 0.0;
-		for( unsigned i=0; i<N; ++i )
-		{
-			
-			double k = k0*exp(((int)i - (int)N/2+1) * dlnk);
-			//double k = k0*exp(((int)i - (int)N/2) * dlnk);
-			//double k = k0*exp(ii * dlnk);
-			
-			//... some constants missing ...//
-			in[i].re = dplus*sqrtpnorm*ptf_->compute( k )*pow(k,0.5*nspec_)*pow(k,1.5-q);
-			in[i].im = 0.0;
-			
-			sum_in += in[i].re;
-			ofsk << std::setw(16) << k <<std::setw(16) << in[i].re << std::endl;
-		}
-		ofsk.close();
-		
-		
-		p = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE);
-		ip = fftw_create_plan(N, FFTW_BACKWARD, FFTW_ESTIMATE);
-		
-		//fftw_one(p, in, out);
-		fftw_one(p, in, out);
-		
-		//... compute the Hankel transform by convolution with the Bessel function
-		for( unsigned i=0; i<N; ++i )
-		{
-			int ii=i;
-			if( ii > (int)N/2 )
-				ii -= N;
-			
-#ifndef NZERO_Q
-			double y=ii*M_PI/L;
-			complex zp((mu+1.0)*0.5,y);
-			gsl_sf_result g_a, g_p;
-			gsl_sf_lngamma_complex_e(zp.real(), zp.imag(), &g_a, &g_p);
-			
-			double arg = 2.0*(log(2.0/k0r0)*y+g_p.val);
-			complex cu = complex(out[i].re,out[i].im)*std::polar(1.0,arg);
-			out[i].re = cu.real()*fftnorm;
-			out[i].im = cu.imag()*fftnorm;
-			
-#else		
-			//complex x(dir*q, (double)ii*2.0*M_PI/L);
-			complex x(dir*q, (double)ii*2.0*M_PI/L);
-			gsl_sf_result g_a, g_p;
-			
-			complex g1, g2, garg, U, phase;						
-			complex twotox = pow(complex(2.0,0.0),x);
-			
-			/////////////////////////////////////////////////////////
-			//.. evaluate complex Gamma functions
-			
-			garg = 0.5*(mu+1.0+x);
-			gsl_sf_lngamma_complex_e (garg.real(), garg.imag(), &g_a, &g_p);
-			g1 = std::polar(exp(g_a.val),g_p.val);
-
-			
-			garg = 0.5*(mu+1.0-x);
-			gsl_sf_lngamma_complex_e (garg.real(), garg.imag(), &g_a, &g_p);
-			g2 = std::polar(exp(g_a.val),g_p.val);
-
-			/////////////////////////////////////////////////////////
-			//.. compute U
-			
-			if( (fabs(g2.real()) < 1e-19 && fabs(g2.imag()) < 1e-19) )
-			{
-				//std::cerr << "Warning : encountered possible singularity in TransferFunction_real::transform!\n";
-				g1 = 1.0; g2 = 1.0;
-			}
-			
-			
-			U = twotox * g1 / g2;
-			phase = pow(complex(k0r0,0.0),complex(0.0,2.0*M_PI*(double)ii/L));
-			
-			complex cu = complex(out[i].re,out[i].im)*U*phase*fftnorm;
-			
-			out[i].re = cu.real();
-			out[i].im = cu.imag();
-			
-			if( (out[i].re != out[i].re)||(out[i].im != out[i].im) )
-			{	std::cerr << "NaN @ i=" << i << ", U= " << U << ", phase = " << phase << ", g1 = " << g1 << ", g2 = " << g2 << std::endl;
-				std::cerr << "mu+1+q = " << mu+1.0+q << std::endl;
-				//break;
-			}
-			
-#endif
-
-		}
-			
-		/*out[N/2].im = 0.0;
-		out[N/2+1].im = 0.0;
-		out[N/2+1].re = out[N/2].re;
-		out[N/2].im = 0.0;*/
-		
-		fftw_one(ip, out, in);
-		
-		rr.assign(N,0.0);
-		TT.assign(N,0.0);
-		
-		r0 = k0r0/k0;
-		
-		for( unsigned i=0; i<N; ++i )
-		{
-			int ii = i;
-			ii -= N/2-1;
-			//ii -= N/2;
-			//if( ii>N/2)
-			//	ii-=N;
-			
-			
-			
-			double r = r0*exp(-ii*dlnr);
-			rr[N-i-1] = r;
-			TT[N-i-1] = 4.0*M_PI* sqrt(M_PI/2.0) *  in[i].re*pow(r,-(1.5+q));
-			
-			//TT[N-i-1] = 4.0*M_PI* sqrt(M_PI/2.0) *  in[i].re*exp( -dir*(q+1.5)*ii*dlnr +q*log(k0r0))/r0;
-			
-			//rr[i] = r;
-			//TT[i] = 4.0*M_PI* sqrt(M_PI/2.0) *  in[i].re*pow(r,-(1.5+q));
-			
-		}
-		
-		
-		{
-			std::ofstream ofs("transfer_real_new.txt");
-			for( unsigned i=0; i<N; ++i )
-			{
-				int ii = i;
-				ii -= N/2-1;
-				
-				double r = r0*exp(-ii*dlnr);//r0*exp(ii*dlnr);
-				double T = 4.0*M_PI* sqrt(M_PI/2.0) *  in[i].re*pow(r,-(1.5+q));
-				ofs << r << "\t\t" << T << "\t\t" << in[i].im << std::endl;
-			}
-		}
-		
-
-		fftw_destroy_plan(p);
-		fftw_destroy_plan(ip);
-	}
-	
-public:
-	TransferFunction_real( TransferFunction *tf, double nspec, double pnorm, double dplus, double rmin, double rmax, double knymax, unsigned nr )
-	{
-				
-		ptf_ = tf;
-		nspec_ = nspec;
-	
-		double q = 0.8;
-		
-		std::vector<double> r,T,xp,yp,xn,yn;
-		
-		transform( pnorm, dplus, nr, q, r, T );
-		
-		//... determine r=0 zero component by integrating up to the Nyquist frequency
-		gsl_integration_workspace * wp; 
-		gsl_function F;
-		wp = gsl_integration_workspace_alloc(20000);
-		F.function = &call_wrapper;
-		double par[2]; par[0] = dplus*sqrt(pnorm); //par[1] = M_PI/kny;
-		F.params = (void*)par;
-		double error;
-		
-		//#warning factor of sqrt(1.5) needs to be adjusted for non-equilateral boxes
-		//.. need a factor sqrt( 2*kny^2_x + 2*kny^2_y + 2*kny^2_z )/2 = sqrt(3/2)kny (in equilateral case)
-		gsl_integration_qag (&F, 0.0, sqrt(1.5)*knymax, 0, 1e-8, 20000, GSL_INTEG_GAUSS21, wp, &Tr0_, &error); 
-		//Tr0_ = 0.0;
-		gsl_integration_workspace_free(wp);
-				
-		
-		for( unsigned i=0; i<r.size(); ++i )
-		{
-			// spline positive and negative part separately
-			/*if( T[i] > 0.0 )
-			{
-				xp.push_back( 2.0*log10(r[i]) );
-				yp.push_back( log10(T[i]) );
-				rneg_ = r[i];
-				rneg2_ = rneg_*rneg_;
-			}else {
-				xn.push_back( 2.0*log10(r[i]) );
-				yn.push_back( log10(-T[i]) );
-			}*/
-			
-			
-			if( r[i] > rmin && r[i] < rmax )
-			{
-				xp.push_back( 2.0*log10(r[i]) );
-				yp.push_back( log10(fabs(T[i])) );
-				xn.push_back( 2.0*log10(r[i]) );
-				if( T[i] >= 0.0 ) 
-					yn.push_back( 1.0 );
-				else
-					yn.push_back( -1.0 );
-				
-				
-				//ofs << std::setw(16) << xp.back() << std::setw(16) << yp.back() << std::endl;
-			}
-			
-		}
-		
-
-		
-		
-		
-		accp = gsl_interp_accel_alloc ();
-		accn = gsl_interp_accel_alloc ();
-		
-		//... spline interpolation is only marginally slower here
-		splinep = gsl_spline_alloc (gsl_interp_cspline, xp.size() );
-		splinen = gsl_spline_alloc (gsl_interp_cspline, xn.size() );
-
-		//... set up everything for spline interpolation
-		gsl_spline_init (splinep, &xp[0], &yp[0], xp.size() );
-		gsl_spline_init (splinen, &xn[0], &yn[0], xn.size() );		
-		
-
-		
-		
-		{
-			double dlogr = (log10(rmax)-log10(rmin))/100;
-			std::ofstream ofs("transfer_splinep.txt");			
-			
-			for( int i=0; i< 100; ++i ) 
-			{
-				double r = rmin*pow(10.0,i*dlogr);
-				ofs << std::setw(16) << r << std::setw(16) << compute_real(r*r) << std::endl;
-			}
-		}
-		
-	}
-	
-	static double call_wrapper( double k, void *arg )
-	{
-		double *a = (double*)arg;
-		return 4.0*M_PI*a[0]*ptf_->compute( k )*pow(k,0.5*nspec_)*k*k;
-	}
-	
-	~TransferFunction_real()
-	{
-		gsl_spline_free (splinen);
-		gsl_interp_accel_free (accn);
-		gsl_spline_free (splinep);
-		gsl_interp_accel_free (accp);
-
-	}
-	
-	inline double compute_real( double r2 ) const
-	{
-		const double EPS = 1e-8;
-		const double Reps2 = EPS*EPS;
-		
-		if( r2 <Reps2 )
-			return Tr0_;
-		double q;
-		/*if( r2 < rneg2_ )
-			q = pow(10.0,gsl_spline_eval (splinep, log10(r2), accp));
-		else
-			q = -pow(10.0,gsl_spline_eval(splinen, log10(r2), accn));*/
-		
-		double logr2 = log10(r2);
-		q = pow(10.0,gsl_spline_eval(splinep, logr2, accp));
-		double sign = 1.0;
-		if( gsl_spline_eval(splinen, logr2, accn) < 0.0 )
-			sign = -1.0;
-		return q*sign;
-	}
-};
-
-#endif

src/constraints.cc

@@ -31,7 +31,7 @@ double dsigma2_tophat( double k, void *pparams )
 	
 	double w = 3.0*(sin(x)-x*cos(x))/(x*x*x);
 	
-	double tfk = ptf->compute(k,total);
+	double tfk = ptf->compute(k,delta_matter);
 	return k*k * w*w * pow(k,nspect) * tfk*tfk;
 }
 
@@ -48,7 +48,7 @@ double dsigma2_gauss( double k, void *pparams )
 	
 	double w = exp(-x*x*0.5);
 	
-	double tfk = ptf->compute(k,total);
+	double tfk = ptf->compute(k,delta_matter);
 	return k*k * w*w * pow(k,nspect) * tfk*tfk;
 }
 
@@ -76,15 +76,15 @@ void compute_sigma_tophat( config_file& cf, transfer_function *ptf, double R, st
 	z.clear();
 	sigma.clear();
 	
-	cosmology cosm( cf );
+	cosmology::parameters cosm( cf );
-	CosmoCalc ccalc( cosm, ptf );
+	cosmology::calculator ccalc( cf );
 	
 	double zmin = 0.0, zmax = 200.0;
 	int nz = 100;
 	for( int i=0; i <nz; ++i )
 		z.push_back( zmax - i*(zmax-zmin)/(nz-1.0) );
 	
-	double D0 = ccalc.CalcGrowthFactor(1.0);
+	double D0 = ccalc.get_growth_factor(1.0);
 
 	double sigma8 = cf.get_value<double>("cosmology","sigma_8"); 
 	double nspec = cf.get_value<double>("cosmology","nspec");
@@ -109,7 +109,7 @@ void compute_sigma_tophat( config_file& cf, transfer_function *ptf, double R, st
 		params[2] = reinterpret_cast<char*> (&nspec);
 		
 		double sig = sqrt(4.0*M_PI*integrate( &dsigma2_tophat, 1e-4, 1e4, reinterpret_cast<void*>(params) ));
-		double Dz  = ccalc.CalcGrowthFactor(1./(1.+z[i]));
+		double Dz  = ccalc.get_growth_factor(1./(1.+z[i]));
 		sigma.push_back( sig*sigma8/sigma0*Dz/D0 );
 	}
 }
@@ -119,15 +119,15 @@ void compute_sigma_gauss( config_file& cf, transfer_function *ptf, double R, std
 	z.clear();
 	sigma.clear();
 	
-	cosmology cosm( cf );
+	cosmology::parameters cosm( cf );
-	CosmoCalc ccalc( cosm, ptf );
+	cosmology::calculator ccalc( cf );
 	
 	double zmin = 0.0, zmax = 200.0;
 	int nz = 100;
 	for( int i=0; i <nz; ++i )
 		z.push_back( zmax - i*(zmax-zmin)/(nz-1.0) );
 	
-	double D0 = ccalc.CalcGrowthFactor(1.0);
+	double D0 = ccalc.get_growth_factor(1.0);
 	
 	double sigma8 = cf.get_value<double>("cosmology","sigma_8"); 
 	double nspec = cf.get_value<double>("cosmology","nspec");
@@ -152,7 +152,7 @@ void compute_sigma_gauss( config_file& cf, transfer_function *ptf, double R, std
 		params[2] = reinterpret_cast<char*> (&nspec);
 		
 		double sig = sqrt(4.0*M_PI*integrate( &dsigma2_gauss, 1e-4, 1e4, reinterpret_cast<void*>(params) ));
-		double Dz  = ccalc.CalcGrowthFactor(1./(1.+z[i]));
+		double Dz  = ccalc.get_growth_factor(1./(1.+z[i]));
 		
 		//std::cerr << z[i] << "    " << sig << std::endl;
 		sigma.push_back( sig*sigma8/sigma0*Dz/D0 );
@@ -163,8 +163,8 @@ void compute_sigma_gauss( config_file& cf, transfer_function *ptf, double R, std
 constraint_set::constraint_set( config_file& cf, transfer_function *ptf )
 : pcf_( &cf ), ptf_( ptf )
 {
-	pcosmo_ = new Cosmology( cf );
+	pcosmo_ = new cosmology::parameters( cf );
-	pccalc_ = new CosmoCalc( *pcosmo_, ptf_ );
+	pccalc_ = new cosmology::calculator( cf );
 	dplus0_ = 1.0;//pccalc_->CalcGrowthFactor( 1.0 );
 	
 	
@@ -218,9 +218,9 @@ constraint_set::constraint_set( config_file& cf, transfer_function *ptf )
 				double zcoll = cf.get_value<double>( "constraints", temp1 );
 				new_c.Rg = pow((mass/pow(2.*M_PI,1.5)/rhom),1./3.);
 				
-				new_c.sigma = 1.686/(pccalc_->CalcGrowthFactor(1./(1.+zcoll))/pccalc_->CalcGrowthFactor(1.0));
+				new_c.sigma = 1.686/(pccalc_->get_growth_factor(1./(1.+zcoll))/pccalc_->get_growth_factor(1.0));
                 music::ilog.Print("sigma of constraint : %g", new_c.sigma );
-                new_c.sigma *=pccalc_->CalcGrowthFactor(astart)/pccalc_->CalcGrowthFactor(1.0);
+                new_c.sigma *=pccalc_->get_growth_factor(astart)/pccalc_->get_growth_factor(1.0);
 				music::ilog.Print("Constraint %d : halo with %g h-1 M_o",i,pow(2.*M_PI,1.5)*rhom*pow(new_c.Rg,3));
 			}
 			else if( new_c.type == peak )
@@ -324,7 +324,7 @@ void constraint_set::wnoise_constr_corr( double dx, size_t nx, size_t ny, size_t
 					
 					double k = sqrt(iix*iix+iiy*iiy+iiz*iiz)*(double)nx/lsub;
 					
-					double T = ptf_->compute(k,total);
+					double T = ptf_->compute(k,delta_matter);
 					double Pk = pnorm*T*T*pow(k,nspec)*d3k;
 					
 					size_t q = ((size_t)ix*ny+(size_t)iy)*nzp+(size_t)iz;
@@ -407,7 +407,7 @@ void constraint_set::wnoise_constr_corr( double dx, complex_t* cw, size_t nx, si
 					iiz *= 2.0*M_PI/nx;
 					
 					double k = sqrt(iix*iix+iiy*iiy+iiz*iiz)*(double)nx/lsub;
-					double T = ptf_->compute(k,total);
+					double T = ptf_->compute(k,delta_matter);
 					
 					std::complex<double> v(std::conj(eval_constr(i,iix,iiy,iiz)));
 					
@@ -472,7 +472,7 @@ void constraint_set::icov_constr( double dx, size_t nx, size_t ny, size_t nz, ma
 						iiz *= 2.0*M_PI/nx;
 						
 						double k = sqrt(iix*iix+iiy*iiy+iiz*iiz)*(double)nx/lsub;
-						double T = ptf_->compute(k,total);
+						double T = ptf_->compute(k,delta_matter);
 						std::complex<double> v(std::conj(eval_constr(i,iix,iiy,iiz)));
 						v *= eval_constr(j,iix,iiy,iiz);
 						v *= pnorm * pow(k,nspec) * T * T * d3k;

src/constraints.hh

@@ -17,7 +17,7 @@
 #include <general.hh>
 #include <config_file.hh>
 #include <transfer_function.hh>
-#include <cosmology.hh>
+#include <cosmology_calculator.hh>
 
 //! matrix class serving as a gsl wrapper
 class matrix
@@ -119,8 +119,8 @@ protected:
 	config_file *pcf_;
 	std::vector<constraint> cset_;
 	transfer_function *ptf_;
-	CosmoCalc *pccalc_;
+	const cosmology::calculator *pccalc_;
-	Cosmology *pcosmo_;
+	const cosmology::parameters *pcosmo_;
 	double dplus0_;
 	unsigned constr_level_;
 

src/convolution_kernel.cc

@@ -157,7 +157,7 @@ class kernel_k : public kernel
 {
 protected:
 	/**/
-	double boxlength_, patchlength_, nspec_, pnorm_, volfac_, kfac_, kmax_;
+	double boxlength_, patchlength_, nspec_, pnorm_, kfac_, kmax_;
 	TransferFunction_k *tfk_;
 
 public:
@@ -165,9 +165,8 @@ public:
 		: kernel(cf, ptf, refh, type)
 	{
 		boxlength_ = pcf_->get_value<double>("setup", "boxlength");
-		nspec_ = pcf_->get_value<double>("cosmology", "nspec");
+		nspec_ = ptf->cosmo_params_["n_s"];
-		pnorm_ = pcf_->get_value<double>("cosmology", "pnorm");
+		pnorm_ = ptf->cosmo_params_["pnorm"];
-		volfac_ = 1.0; //pow(boxlength,3)/pow(2.0*M_PI,3);
 		kfac_ = 2.0 * M_PI / boxlength_;
 		kmax_ = kfac_ / 2;
 		tfk_ = new TransferFunction_k(type_, ptf_, nspec_, pnorm_);
@@ -231,7 +230,7 @@ public:
 		for (size_t i = 0; i < len; ++i)
 		{
 			double kk = kfac_ * in_k[i];
-			out_Tk[i] = volfac_ * tfk_->compute(kk);
+			out_Tk[i] = tfk_->compute(kk);
 		}
 	}
 

src/cosmology.hh

@@ -1,224 +0,0 @@
-/*
- 
- cosmology.hh - This file is part of MUSIC -
- a code to generate multi-scale initial conditions 
- for cosmological simulations 
- 
- Copyright (C) 2010  Oliver Hahn
- 
- */
-
-#ifndef _COSMOLOGY_HH
-#define _COSMOLOGY_HH
-
-
-#include "transfer_function.hh"
-#include "mesh.hh"
-#include "general.hh"
-
-/*!
- * @class CosmoCalc
- * @brief provides functions to compute cosmological quantities
- *
- * This class provides member functions to compute cosmological quantities
- * related to the Friedmann equations and linear perturbation theory
- */
-class CosmoCalc
-{
-public:
-	//! data structure to store cosmological parameters
-	Cosmology m_Cosmology;
-	
-	//! pointer to an instance of a transfer function plugin
-	transfer_function_plugin *m_pTransferFunction;
-	
-	
-	//! constructor for a cosmology calculator object
-	/*!
-	 * @param acosmo a cosmological parameters structure
-	 * @param pTransferFunction pointer to an instance of a transfer function object
-	 */
-	 
-	CosmoCalc( const Cosmology acosmo, transfer_function_plugin *pTransferFunction )
-	{
-		m_Cosmology = acosmo;
-		m_pTransferFunction = pTransferFunction;
-	}
-	
-	//! returns the amplitude of amplitude of the power spectrum
-	/*!
-	 * @param k the wave number in h/Mpc
-	 * @param a the expansion factor of the universe
-	 * @returns power spectrum amplitude for wave number k at time a
-	 */
-	inline real_t Power( real_t k, real_t a ){
-		real_t m_Dplus    = CalcGrowthFactor( a );
-		real_t m_DplusOne = CalcGrowthFactor( 1.0 );
-		real_t m_pNorm = ComputePNorm( 1e4 );
-		m_Dplus    /= m_DplusOne;
-		m_DplusOne = 1.0;
-		real_t scale = m_Dplus/m_DplusOne;
-		return m_pNorm*scale*scale*TransferSq(k)*pow((double)k,(double)m_Cosmology.nspect);
-	}
-
-  inline static double H_of_a( double a, void *Params )
-  {
-    Cosmology *cosm = (Cosmology*)Params;
-    double a2 = a*a;
-    double Ha = sqrt(cosm->Omega_m/(a2*a) + cosm->Omega_k/a2
-		     + cosm->Omega_DE * pow(a,-3.*(1.+cosm->w_0+cosm->w_a)) * exp(-3.*(1.0-a)*cosm->w_a) );
-    return Ha;
-  }
-
-  inline static double Hprime_of_a( double a, void *Params )
-  {
-    Cosmology *cosm = (Cosmology*)Params;
-    double a2 = a*a;
-    double H  = H_of_a( a, Params );
-    double Hprime = 1/(a*H) * ( -1.5 * cosm->Omega_m / (a2*a) - cosm->Omega_k / a2
-      - 1.5 * cosm->Omega_DE * pow( a, -3.*(1.+cosm->w_0+cosm->w_a) ) * exp( -3.*(1.0-a)*cosm->w_a )
-      * ( 1. + cosm->w_0 + (1.-a) * cosm->w_a ) );
-    return Hprime;
-  }
-
-
-  //! Integrand used by function CalcGrowthFactor to determine the linear growth factor D+
-  inline static double GrowthIntegrand( double a, void *Params )
-  {
-    double Ha = a * H_of_a( a, Params );
-    return 2.5/( Ha * Ha * Ha );
-  }
-  
-    //! Computes the linear theory growth factor D+
-    /*! Function integrates over member function GrowthIntegrand and computes
-    *                      /a
-    *   D+(a) = 5/2 H(a) * |  [a'^3 * H(a')^3]^(-1) da'
-    *                      /0
-    */
-    real_t CalcGrowthFactor( real_t a )
-    {
-        real_t integral = integrate( &GrowthIntegrand, 0.0, a, (void*)&m_Cosmology );
-        return H_of_a( a, (void*)&m_Cosmology ) * integral;
-    }
-
-    //! Compute the factor relating particle displacement and velocity
-    /*! Function computes
-    *
-    *  vfac = a^2 * H(a) * dlogD+ / d log a = a^2 * H'(a) + 5/2 * [ a * D+(a) * H(a) ]^(-1)
-    *
-    */
-        real_t CalcVFact( real_t a )
-        {
-        real_t Dp = CalcGrowthFactor( a );
-        real_t H  = H_of_a( a, (void*)&m_Cosmology );
-        real_t Hp = Hprime_of_a( a, (void*)&m_Cosmology );
-        real_t a2 = a*a;
-
-        return ( a2 * Hp + 2.5 / ( a * Dp * H ) ) * 100.0;
-    }
-		
-	
-    //! Integrand for the sigma_8 normalization of the power spectrum
-    /*! Returns the value of the primordial power spectrum multiplied with 
-     the transfer function and the window function of 8 Mpc/h at wave number k */
-    static double dSigma8( double k, void *Params )
-    {
-        if( k<=0.0 )
-            return 0.0f;
-        
-        transfer_function *ptf = (transfer_function *)Params;
-        
-        double x = k*8.0;
-        double w = 3.0*(sin(x)-x*cos(x))/(x*x*x);
-        static double nspect = (double)ptf->cosmo_.nspect;
-        
-        double tf = ptf->compute(k, total);
-        
-        //... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1
-        return k*k * w*w * pow((double)k,(double)nspect) * tf*tf;
-        
-    }
-    
-    //! Integrand for the sigma_8 normalization of the power spectrum
-	/*! Returns the value of the primordial power spectrum multiplied with 
-	 the transfer function and the window function of 8 Mpc/h at wave number k */
-	static double dSigma8_0( double k, void *Params )
-	{
-		if( k<=0.0 )
-			return 0.0f;
-		
-		transfer_function *ptf = (transfer_function *)Params;
-		
-		double x = k*8.0;
-		double w = 3.0*(sin(x)-x*cos(x))/(x*x*x);
-		static double nspect = (double)ptf->cosmo_.nspect;
-		
-		double tf = ptf->compute(k, total0);
-		
-		//... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1
-		return k*k * w*w * pow((double)k,(double)nspect) * tf*tf;
-		
-	}
-	
-	
-	//! Computes the square of the transfer function
-	/*! Function evaluates the supplied transfer function m_pTransferFunction
-	 * and returns the square of its value at wave number k
-	 * @param k wave number at which to evaluate the transfer function
-	 */
-	inline real_t TransferSq( real_t k ){
-		//.. parameter supplied transfer function
-		real_t tf1 = m_pTransferFunction->compute(k, total);
-		return tf1*tf1;
-	}
-	
-	
-	//! Computes the normalization for the power spectrum
-	/*!
-	 * integrates the power spectrum to fix the normalization to that given
-	 * by the sigma_8 parameter
-	 */
-	real_t ComputePNorm( real_t kmax )
-	{
-		real_t sigma0, kmin;
-		kmax = m_pTransferFunction->get_kmax();//m_Cosmology.H0/8.0;
-		kmin = m_pTransferFunction->get_kmin();//0.0;
-        
-        if( !m_pTransferFunction->tf_has_total0() )
-            sigma0 = 4.0 * M_PI * integrate( &dSigma8, (double)kmin, (double)kmax, (void*)m_pTransferFunction );
-		else
-            sigma0 = 4.0 * M_PI * integrate( &dSigma8_0, (double)kmin, (double)kmax, (void*)m_pTransferFunction );
-		
-        return m_Cosmology.sigma8*m_Cosmology.sigma8/sigma0;
-	}
-	
-};
-
-
-//! compute the jeans sound speed
-/*! given a density in g/cm^-3 and a mass in g it gives back the sound
- *  speed in cm/s for which the input mass is equal to the jeans mass
- *  @param rho density 
- *  @param mass mass scale
- *  @returns jeans sound speed
- */
-inline double jeans_sound_speed( double rho, double mass )
-{
-	const double G = 6.67e-8;
-	return pow( 6.0*mass/M_PI*sqrt(rho)*pow(G,1.5), 1.0/3.0 );
-}
-
-//! computes the density from the potential using the Laplacian
-void compute_Lu_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order=4 );
-
-//! computes the 2nd order density perturbations using also off-diagonal terms in the potential Hessian 
-void compute_LLA_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order=4 );
-
-//! computes the source term for the 2nd order perturbations in the displacements
-void compute_2LPT_source( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order=4 );
-
-void compute_2LPT_source_FFT( config_file& cf_, const grid_hierarchy& u, grid_hierarchy& fnew );
-
-
-#endif // _COSMOLOGY_HH
-

src/cosmology_calculator.hh

@@ -17,11 +17,11 @@
 #pragma once
 
 #include <array>
-#include <vec.hh>
+#include "math/vec.hh"
 
 #include <cosmology_parameters.hh>
 #include <physical_constants.hh>
-#include <transfer_function_plugin.hh>
+#include <transfer_function.hh>
 #include <math/ode_integrate.hh>
 #include <logger.hh>
 
@@ -30,7 +30,7 @@
 #include <gsl/gsl_integration.h>
 #include <gsl/gsl_errno.h>
 
-namespace cosmology_new
+namespace cosmology
 {
 
 /*!
@@ -48,7 +48,7 @@ public:
     cosmology::parameters cosmo_param_;
 
     //! pointer to an instance of a transfer function plugin
-    std::unique_ptr<TransferFunction_plugin> transfer_function_;
+    std::unique_ptr<transfer_function_plugin> transfer_function_;
 
 private:
     static constexpr double REL_PRECISION = 1e-10;
@@ -90,7 +90,7 @@ private:
      */
     void compute_growth( std::vector<double>& tab_a, std::vector<double>& tab_D, std::vector<double>& tab_f )
     {
-        using v_t = vec_t<3, double>;
+        using v_t = vec_t<3,double>;
 
         // set ICs, very deep in radiation domination
         const double a0 = 1e-10;
@@ -177,9 +177,10 @@ public:
         Dplus_target_ = D_of_a_( atarget_ ) / Dnow_;
 
         music::ilog << "Linear growth factors: D+_target = " << Dplus_target_ << ", D+_start = " << Dplus_start_ << std::endl;
+        music::ilog << "-------------------------------------------------------------------------------" << std::endl;
 
         // set up transfer functions and compute normalisation
-        transfer_function_ = select_TransferFunction_plugin(cf, cosmo_param_);
+        transfer_function_ = select_transfer_function_plugin(cf, cosmo_param_);
         transfer_function_->intialise();
         if( !transfer_function_->tf_isnormalised_ ){
             cosmo_param_.set("pnorm", this->compute_pnorm_from_sigma8() );
@@ -190,9 +191,14 @@ public:
         }
         cosmo_param_.set("sqrtpnorm", std::sqrt(cosmo_param_["pnorm"]));
 
-        music::ilog << std::setw(32) << std::left << "TF supports distinct CDM+baryons"
+        // if (!transfer_function_->tf_is_distinct())
+        //     music::wlog << " - WARNING: The selected transfer function does not support" << std::endl
+        //                 << "            distinct amplitudes for baryon and DM fields!" << std::endl
+        //                 << "            Perturbation amplitudes will be identical!" << std::endl;
+
+        music::ilog << std::setw(32) << std::left << " . TF supports distinct CDM+baryons"
                     << " : " << (transfer_function_->tf_is_distinct() ? "yes" : "no") << std::endl;
-        music::ilog << std::setw(32) << std::left << "TF maximum wave number"
+        music::ilog << std::setw(32) << std::left << " . TF maximum wave number"
                     << " : " << transfer_function_->get_kmax() << " h/Mpc" << std::endl;
 
         m_n_s_ = cosmo_param_["n_s"];

src/cosmology_parameters.cc

@@ -21,7 +21,7 @@
  * @brief namespace encapsulating all things cosmology
  * 
  */
-namespace cosmology_new{
+namespace cosmology{
 
 //! we store here the preset cosmological paramters
 parameters::defaultmmap_t parameters::default_pmaps_

src/cosmology_parameters.hh

@@ -23,7 +23,7 @@
 #include <config_file.hh>
 #include <general.hh>
 
-namespace cosmology_new
+namespace cosmology
 {
     //! structure for cosmological parameters
     class parameters
@@ -50,6 +50,19 @@ namespace cosmology_new
             return it->second;
         }
 
+        //! get routine for cosmological parameter key-value pairs
+        double& get(const std::string &key)
+        {
+            auto it = pmap_.find(key);
+            if (it == pmap_.end())
+            {
+                auto errmsg = std::string("Cosmological parameter \'") + key + std::string("\' does not exist in internal list.");
+                music::elog << errmsg << std::endl;
+                throw std::runtime_error(errmsg.c_str());
+            }
+            return it->second;
+        }
+
         //! set routine for cosmological parameter key-value pairs
         void set(const std::string &key, const double value)
         {
@@ -69,6 +82,9 @@ namespace cosmology_new
         //! shortcut get routine for cosmological parameter key-value pairs through bracket operator
         inline double operator[](const std::string &key) const { return this->get(key); }
 
+        //! shortcut get routine for cosmological parameter key-value pairs through bracket operator
+        inline double& operator[](const std::string &key) { return this->get(key); }
+
         //! default constructor does nothing
         parameters() {}
 

src/densities.cc

@@ -253,9 +253,10 @@ void fft_interpolate(m1 &V, m2 &v, bool from_basegrid = false)
 /*******************************************************************************************/
 /*******************************************************************************************/
 
-void GenerateDensityUnigrid(config_file &cf, transfer_function *ptf, tf_type type,
+void GenerateDensityUnigrid(config_file &cf, const cosmology::calculator* cc, tf_type type,
 							refinement_hierarchy &refh, noise_generator &rand, grid_hierarchy &delta, bool smooth, bool shift)
 {
+	auto ptf = cc->transfer_function_.get();
 	unsigned levelmin, levelmax, levelminPoisson;
 
 	levelminPoisson = cf.get_value<unsigned>("setup", "levelmin");
@@ -308,10 +309,11 @@ void GenerateDensityUnigrid(config_file &cf, transfer_function *ptf, tf_type typ
 /*******************************************************************************************/
 /*******************************************************************************************/
 
-void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type type,
+void GenerateDensityHierarchy(config_file &cf, const cosmology::calculator* cc, tf_type type,
 							  refinement_hierarchy &refh, noise_generator &rand,
 							  grid_hierarchy &delta, bool smooth, bool shift)
 {
+	auto ptf = cc->transfer_function_.get();
 	unsigned levelmin, levelmax, levelminPoisson;
 	std::vector<long> rngseeds;
 	std::vector<std::string> rngfnames;

src/densities.hh

@@ -15,16 +15,14 @@
 
 #include "general.hh"
 #include "config_file.hh"
-// #include "density_grid.hh"
 #include "random.hh"
-#include "cosmology.hh"
 #include "transfer_function.hh"
 #include "general.hh"
 
-void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type type,
+void GenerateDensityHierarchy(config_file &cf, const cosmology::calculator* cc, tf_type type,
 															refinement_hierarchy &refh, noise_generator &rand, grid_hierarchy &delta, bool smooth, bool shift);
 
-void GenerateDensityUnigrid(config_file &cf, transfer_function *ptf, tf_type type,
+void GenerateDensityUnigrid(config_file &cf, const cosmology::calculator*, tf_type type,
 														refinement_hierarchy &refh, noise_generator &rand, grid_hierarchy &delta, bool smooth, bool shift);
 
 void normalize_density(grid_hierarchy &delta);

src/general.hh

@@ -12,9 +12,10 @@
 
 #include <logger.hh>
 #include <config_file.hh>
-
+#include <memory>
 #include <cassert>
 #include <omp.h>
+#include <complex>
 
 #include <fftw3.h>
 
@@ -35,6 +36,8 @@
   #define FFTW_PREFIX fftwl
 #endif
 
+using ccomplex_t = std::complex<real_t>;
+
 #define FFTW_GEN_NAME_PRIM(a, b) a##_##b
 #define FFTW_GEN_NAME(a, b) FFTW_GEN_NAME_PRIM(a, b)
 #define FFTW_API(x) FFTW_GEN_NAME(FFTW_PREFIX, x)
@@ -46,14 +49,13 @@ using fftw_plan_t = FFTW_GEN_NAME(FFTW_PREFIX, plan);
 
 #include <vector>
 #include <array>
-using vec3_t = std::array<real_t,3>;
 
 namespace CONFIG
 {
 // extern int MPI_thread_support;
-// extern int MPI_task_rank;
+extern int MPI_task_rank;
-// extern int MPI_task_size;
+extern int MPI_task_size;
-// extern bool MPI_ok;
+extern bool MPI_ok;
 // extern bool MPI_threads_ok;
 extern bool FFTW_threads_ok;
 extern int num_threads;
@@ -80,7 +82,8 @@ inline T POW4( T a ){
 
 
 //! structure for cosmological parameters
-typedef struct cosmology{
+/*
+typedef struct cosmology_old{
   double 
     Omega_m,		//!< baryon+dark matter density
     Omega_b,		//!< baryon matter density
@@ -133,7 +136,7 @@ typedef struct cosmology{
   {
     
   }
-}Cosmology;
+}Cosmology;*/
 
 //! basic box/grid/refinement structure parameters
 typedef struct {

src/main.cc

@@ -44,7 +44,9 @@ extern "C"
 #include <densities.hh>
 
 #include <convolution_kernel.hh>
-#include <cosmology.hh>
+#include <perturbation_theory.hh>
+#include <cosmology_parameters.hh>
+#include <cosmology_calculator.hh>
 #include <transfer_function.hh>
 
 #define THE_CODE_NAME "music!"
@@ -53,9 +55,9 @@ extern "C"
 // initialise with "default" values
 namespace CONFIG{
 // int  MPI_thread_support = -1;
-// int  MPI_task_rank = 0;
+int  MPI_task_rank = 0;
-// int  MPI_task_size = 1;
+int  MPI_task_size = 1;
-// bool MPI_ok = false;
+bool MPI_ok = false;
 // bool MPI_threads_ok = false;
 bool FFTW_threads_ok = false;
 int  num_threads = 1;
@@ -67,7 +69,7 @@ namespace music
 
 	struct framework
 	{
-		transfer_function *the_transfer_function;
+		// transfer_function *the_transfer_function;
 		// poisson_solver *the_poisson_solver;
 		config_file *the_config_file;
 		refinement_hierarchy *the_refinement_hierarchy;
@@ -76,13 +78,15 @@ namespace music
 }
 
 //... declare static class members here
-transfer_function *TransferFunction_real::ptf_ = NULL;
+// transfer_function *TransferFunction_real::ptf_ = NULL;
 transfer_function *TransferFunction_k::ptf_ = NULL;
 tf_type TransferFunction_k::type_;
-tf_type TransferFunction_real::type_;
+// tf_type TransferFunction_real::type_;
-real_t TransferFunction_real::nspec_ = -1.0;
+// real_t TransferFunction_real::nspec_ = -1.0;
 real_t TransferFunction_k::nspec_ = -1.0;
 
+std::unique_ptr<cosmology::calculator>  the_cosmo_calc;
+
 //... prototypes for routines used in main driver routine
 void splash(void);
 void modify_grid_for_TF(const refinement_hierarchy &rh_full, refinement_hierarchy &rh_TF, config_file &cf);
@@ -95,17 +99,24 @@ void splash(void)
 {
 
 	music::ilog << std::endl
-			<< "    __    __     __  __     ______     __     ______      " << std::endl
+			<< "           __    __     __  __     ______     __     ______      " << std::endl
-			<< "   /\\ \"-./  \\   /\\ \\/\\ \\   /\\  ___\\   /\\ \\   /\\  ___\\  " << std::endl
+			<< "          /\\ \"-./  \\   /\\ \\/\\ \\   /\\  ___\\   /\\ \\   /\\  ___\\  " << std::endl
-			<< "   \\ \\ \\-./\\ \\  \\ \\ \\_\\ \\  \\ \\___  \\  \\ \\ \\  \\ \\ \\____ " << std::endl
+			<< "          \\ \\ \\-./\\ \\  \\ \\ \\_\\ \\  \\ \\___  \\  \\ \\ \\  \\ \\ \\____ " << std::endl
-			<< "    \\ \\_\\ \\ \\_\\  \\ \\_____\\  \\/\\_____\\  \\ \\_\\  \\ \\_____\\ " << std::endl
+			<< "           \\ \\_\\ \\ \\_\\  \\ \\_____\\  \\/\\_____\\  \\ \\_\\  \\ \\_____\\ " << std::endl
-			<< "     \\/_/  \\/_/   \\/_____/   \\/_____/   \\/_/   \\/_____/ " << std::endl << std::endl
+			<< "            \\/_/  \\/_/   \\/_____/   \\/_____/   \\/_/   \\/_____/ " << std::endl << std::endl
-			<< "                            this is " << THE_CODE_NAME << " version " << THE_CODE_VERSION << std::endl << std::endl;
+			<< "                                   this is " << THE_CODE_NAME << " version " << THE_CODE_VERSION << std::endl << std::endl;
 
-#if defined(CMAKE_BUILD)
+// git and versioning info:
-	music::ilog.Print("Version built from git rev.: %s, tag: %s, branch: %s", GIT_REV, GIT_TAG, GIT_BRANCH);
+    music::ilog << "Version: git rev.: " << GIT_REV << ", tag: " << GIT_TAG << ", branch: " << GIT_BRANCH << std::endl;
+    
+    // Compilation CMake configuration, time etc info:
+    music::ilog << "This " << CMAKE_BUILDTYPE_STR << " build was compiled at " << __TIME__ << " on " <<  __DATE__ << std::endl;
+
+#ifdef __GNUC__
+    music::ilog << "Compiled with GNU C++ version " << __VERSION__ <<std::endl;
+#else
+    music::ilog << "Compiled with " << __VERSION__ << std::endl;
 #endif
-	music::ilog << std::endl << std::endl;
 }
 
 void modify_grid_for_TF(const refinement_hierarchy &rh_full, refinement_hierarchy &rh_TF, config_file &cf)
@@ -378,6 +389,7 @@ int main(int argc, const char *argv[])
 		splash();
 		std::cout << " This version is compiled with the following plug-ins:\n";
 
+		cosmology::print_ParameterSets();
 		print_region_generator_plugins();
 		print_transfer_function_plugins();
 		print_RNG_plugins();
@@ -407,8 +419,6 @@ int main(int argc, const char *argv[])
 	config_file cf(argv[1]);
 	std::string tfname, randfname, temp;
 	bool force_shift(false);
-	double boxlength;
-
 
 	//------------------------------------------------------------------------------
 	//... init multi-threading
@@ -424,7 +434,6 @@ int main(int argc, const char *argv[])
 	//... initialize some parameters about grid set-up
 	//------------------------------------------------------------------------------
 
-	boxlength = cf.get_value<double>("setup", "boxlength");
 	lbase = cf.get_value<unsigned>("setup", "levelmin");
 	lmax = cf.get_value<unsigned>("setup", "levelmax");
 	lbaseTF = cf.get_value_safe<unsigned>("setup", "levelmin_TF", lbase);
@@ -463,39 +472,42 @@ int main(int argc, const char *argv[])
 			do_LLA = cf.get_value_safe<bool>("setup", "use_LLA", false),
 			do_counter_mode = cf.get_value_safe<bool>("setup", "zero_zoom_velocity", false);
 
-	transfer_function_plugin *the_transfer_function_plugin = select_transfer_function_plugin(cf);
+	the_cosmo_calc              = std::make_unique<cosmology::calculator>(cf);
 
-	cosmology cosmo(cf);
+	bool tf_has_velocities = the_cosmo_calc.get()->transfer_function_.get()->tf_has_velocities();
+	//--------------------------------------------------------------------------------------------------------
+	//! starting redshift
+	const real_t zstart = cf.get_value<double>("setup", "zstart");
+	const real_t astart = 1.0/(1.0+zstart);
+	
+	music::ilog << "- starting at a=" << 1.0/(1.0+zstart) << std::endl;
 
-	music::ilog << "- starting at a=" << cosmo.astart << std::endl;
+	double cosmo_dplus = the_cosmo_calc->get_growth_factor(astart) / the_cosmo_calc->get_growth_factor(1.0);
-
+	double cosmo_vfact = the_cosmo_calc->get_vfact(astart);
-	CosmoCalc ccalc(cosmo, the_transfer_function_plugin);
-	cosmo.pnorm = ccalc.ComputePNorm(2.0 * M_PI / boxlength);
-	cosmo.dplus = ccalc.CalcGrowthFactor(cosmo.astart) / ccalc.CalcGrowthFactor(1.0);
-	cosmo.vfact = ccalc.CalcVFact(cosmo.astart);
-
-	if (!the_transfer_function_plugin->tf_has_total0())
-		cosmo.pnorm *= cosmo.dplus * cosmo.dplus;
 
+	if (!the_cosmo_calc.get()->transfer_function_.get()->tf_has_total0()){
+		the_cosmo_calc->cosmo_param_["pnorm"] *= cosmo_dplus * cosmo_dplus;
+	}
+	double cosmo_pnorm = the_cosmo_calc->cosmo_param_["pnorm"];
 	//... directly use the normalisation via a parameter rather than the calculated one
-	cosmo.pnorm = cf.get_value_safe<double>("setup", "force_pnorm", cosmo.pnorm);
+	cosmo_pnorm = cf.get_value_safe<double>("setup", "force_pnorm", cosmo_pnorm);
 
 	double vfac2lpt = 1.0;
 
-	if (the_transfer_function_plugin->tf_velocity_units() && do_baryons)
+	if (the_cosmo_calc->transfer_function_->tf_velocity_units() && do_baryons)
 	{
-		vfac2lpt = cosmo.vfact; // if the velocities are in velocity units, we need to divide by vfact for the 2lPT term
+		vfac2lpt = cosmo_vfact; // if the velocities are in velocity units, we need to divide by vfact for the 2lPT term
-		cosmo.vfact = 1.0;
+		cosmo_vfact = 1.0;
 	}
 
-	//
+	
 	{
 		char tmpstr[128];
-		snprintf(tmpstr, 128, "%.12g", cosmo.pnorm);
+		snprintf(tmpstr, 128, "%.12g", cosmo_pnorm);
 		cf.insert_value("cosmology", "pnorm", tmpstr);
-		snprintf(tmpstr, 128, "%.12g", cosmo.dplus);
+		snprintf(tmpstr, 128, "%.12g", cosmo_dplus);
 		cf.insert_value("cosmology", "dplus", tmpstr);
-		snprintf(tmpstr, 128, "%.12g", cosmo.vfact);
+		snprintf(tmpstr, 128, "%.12g", cosmo_vfact);
 		cf.insert_value("cosmology", "vfact", tmpstr);
 	}
 
@@ -505,16 +517,13 @@ int main(int argc, const char *argv[])
 	//... determine run parameters
 	//------------------------------------------------------------------------------
 
-	if (!the_transfer_function_plugin->tf_is_distinct() && do_baryons)
+	
-		music::wlog << " - WARNING: The selected transfer function does not support" << std::endl
-							<< "            distinct amplitudes for baryon and DM fields!" << std::endl
-							<< "            Perturbation amplitudes will be identical!" << std::endl;
 
 
 	//------------------------------------------------------------------------------
 	//... start up the random number generator plugin
 	//... see if we need to set some grid building constraints
-	noise_generator rand( cf, the_transfer_function_plugin );
+	noise_generator rand( cf );
 
 	//------------------------------------------------------------------------------
 	//... determine the refinement hierarchy
@@ -549,13 +558,11 @@ int main(int argc, const char *argv[])
 	music::ilog << "   GENERATING WHITE NOISE\n";
 	music::ilog << "-------------------------------------------------------------------------------" << std::endl;
 	music::ilog << "Computing white noise..." << std::endl;
-	// rand_gen rand(cf, rh_TF, the_transfer_function_plugin);
 	rand.initialize_for_grid_structure( rh_TF );
 
 	//------------------------------------------------------------------------------
 	//... initialize the Poisson solver
 	//------------------------------------------------------------------------------
-	// bool bdefd	= cf.get_value_safe<bool> ( "poisson" , "fft_fine", true );
 	bool bdefd = true; // we set this by default and don't allow it to be changed outside any more
 	bool bglass = cf.get_value_safe<bool>("output", "glass", false);
 	bool bsph = cf.get_value_safe<bool>("setup", "do_SPH", false) && do_baryons;
@@ -610,11 +617,11 @@ int main(int argc, const char *argv[])
 			music::ulog.Print("Computing dark matter displacements...");
 
 			grid_hierarchy f(nbnd); //, u(nbnd);
-			tf_type my_tf_type = cdm;
+			tf_type my_tf_type = delta_cdm;
-			if (!do_baryons || !the_transfer_function_plugin->tf_is_distinct())
+			if (!do_baryons)
-				my_tf_type = total;
+				my_tf_type = delta_matter;
 
-			GenerateDensityHierarchy(cf, the_transfer_function_plugin, my_tf_type, rh_TF, rand, f, false, false);
+			GenerateDensityHierarchy(cf, the_cosmo_calc.get(), my_tf_type, rh_TF, rand, f, false, false);
 			coarsen_density(rh_Poisson, f, bspectral_sampling);
 			f.add_refinement_mask(rh_Poisson.get_coord_shift());
 
@@ -678,7 +685,7 @@ int main(int argc, const char *argv[])
 				music::ilog << "   COMPUTING BARYON DENSITY\n";
 				music::ilog << "-------------------------------------------------------------------------------" << std::endl;
 				music::ulog.Print("Computing baryon density...");
-				GenerateDensityHierarchy(cf, the_transfer_function_plugin, baryon, rh_TF, rand, f, false, bbshift);
+				GenerateDensityHierarchy(cf, the_cosmo_calc.get(), delta_baryon, rh_TF, rand, f, false, bbshift);
 				coarsen_density(rh_Poisson, f, bspectral_sampling);
 				f.add_refinement_mask(rh_Poisson.get_coord_shift());
 				normalize_density(f);
@@ -741,17 +748,17 @@ int main(int argc, const char *argv[])
 			//------------------------------------------------------------------------------
 			//... velocities
 			//------------------------------------------------------------------------------
-			if ((!the_transfer_function_plugin->tf_has_velocities() || !do_baryons) && !bsph)
+			if ((!tf_has_velocities || !do_baryons) && !bsph)
 			{
 				music::ilog << "===============================================================================" << std::endl;
 				music::ilog << "   COMPUTING VELOCITIES\n";
 				music::ilog << "-------------------------------------------------------------------------------" << std::endl;
 				music::ulog.Print("Computing velocitites...");
 
-				if (do_baryons || the_transfer_function_plugin->tf_has_velocities())
+				if (do_baryons || tf_has_velocities)
 				{
 					music::ulog.Print("Generating velocity perturbations...");
-					GenerateDensityHierarchy(cf, the_transfer_function_plugin, vtotal, rh_TF, rand, f, false, false);
+					GenerateDensityHierarchy(cf, the_cosmo_calc.get(), theta_cdm, rh_TF, rand, f, false, false);
 					coarsen_density(rh_Poisson, f, bspectral_sampling);
 					f.add_refinement_mask(rh_Poisson.get_coord_shift());
 					normalize_density(f);
@@ -779,7 +786,7 @@ int main(int argc, const char *argv[])
 						the_poisson_solver->gradient(icoord, u, data_forIO);
 
 					//... multiply to get velocity
-					data_forIO *= cosmo.vfact;
+					data_forIO *= cosmo_vfact;
 
 					//... velocity kick to keep refined region centered?
 
@@ -796,7 +803,7 @@ int main(int argc, const char *argv[])
 					coarsen_density(rh_Poisson, data_forIO, false);
 
 					// add counter velocity-mode
-					if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact );
+					if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo_vfact );
 
 					music::ulog.Print("Writing CDM velocities");
 					the_output_plugin->write_dm_velocity(icoord, data_forIO);
@@ -821,7 +828,7 @@ int main(int argc, const char *argv[])
 
 				//... we do baryons and have velocity transfer functions, or we do SPH and not to shift
 				//... do DM first
-				GenerateDensityHierarchy(cf, the_transfer_function_plugin, vcdm, rh_TF, rand, f, false, false);
+				GenerateDensityHierarchy(cf, the_cosmo_calc.get(), theta_cdm, rh_TF, rand, f, false, false);
 				coarsen_density(rh_Poisson, f, bspectral_sampling);
 				f.add_refinement_mask(rh_Poisson.get_coord_shift());
 				normalize_density(f);
@@ -851,7 +858,7 @@ int main(int argc, const char *argv[])
 						the_poisson_solver->gradient(icoord, u, data_forIO);
 
 					//... multiply to get velocity
-					data_forIO *= cosmo.vfact;
+					data_forIO *= cosmo_vfact;
 
 					double sigv = compute_finest_sigma(data_forIO);
 					music::ilog.Print("sigma of %c-velocity of high-res DM is %f", 'x' + icoord, sigv);
@@ -866,7 +873,7 @@ int main(int argc, const char *argv[])
 					coarsen_density(rh_Poisson, data_forIO, false);
 
 					// add counter velocity mode
-					if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact );
+					if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo_vfact );
 
 					music::ulog.Print("Writing CDM velocities");
 					the_output_plugin->write_dm_velocity(icoord, data_forIO);
@@ -880,7 +887,7 @@ int main(int argc, const char *argv[])
 				music::ilog << "-------------------------------------------------------------------------------" << std::endl;
 				music::ulog.Print("Computing baryon velocitites...");
 				//... do baryons
-				GenerateDensityHierarchy(cf, the_transfer_function_plugin, vbaryon, rh_TF, rand, f, false, bbshift);
+				GenerateDensityHierarchy(cf, the_cosmo_calc.get(), theta_baryon, rh_TF, rand, f, false, bbshift);
 				coarsen_density(rh_Poisson, f, bspectral_sampling);
 				f.add_refinement_mask(rh_Poisson.get_coord_shift());
 				normalize_density(f);
@@ -910,7 +917,7 @@ int main(int argc, const char *argv[])
 						the_poisson_solver->gradient(icoord, u, data_forIO);
 
 					//... multiply to get velocity
-					data_forIO *= cosmo.vfact;
+					data_forIO *= cosmo_vfact;
 
 					double sigv = compute_finest_sigma(data_forIO);
 					music::ilog.Print("sigma of %c-velocity of high-res baryons is %f", 'x' + icoord, sigv);
@@ -925,7 +932,7 @@ int main(int argc, const char *argv[])
 					coarsen_density(rh_Poisson, data_forIO, false);
 
 					// add counter velocity mode
-					if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact );
+					if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo_vfact );
 
 					music::ulog.Print("Writing baryon velocities");
 					the_output_plugin->write_gas_velocity(icoord, data_forIO);
@@ -946,13 +953,13 @@ int main(int argc, const char *argv[])
 
 			grid_hierarchy f(nbnd), u1(nbnd), u2LPT(nbnd), f2LPT(nbnd);
 
-			tf_type my_tf_type = vcdm;
+			tf_type my_tf_type = theta_cdm;
 			bool dm_only = !do_baryons;
-			if (!do_baryons || !the_transfer_function_plugin->tf_has_velocities())
+			if (!do_baryons || !tf_has_velocities)
-				my_tf_type = total;
+				my_tf_type = theta_matter;
 
 			music::ilog << "===============================================================================" << std::endl;
-			if (my_tf_type == total)
+			if (my_tf_type == theta_matter)
 			{
 				music::ilog << "   COMPUTING VELOCITIES" << std::endl;
 			}
@@ -962,7 +969,7 @@ int main(int argc, const char *argv[])
 			}
 			music::ilog << "-------------------------------------------------------------------------------" << std::endl;
 
-			GenerateDensityHierarchy(cf, the_transfer_function_plugin, my_tf_type, rh_TF, rand, f, false, false);
+			GenerateDensityHierarchy(cf, the_cosmo_calc.get(), my_tf_type, rh_TF, rand, f, false, false);
 			coarsen_density(rh_Poisson, f, bspectral_sampling);
 			f.add_refinement_mask(rh_Poisson.get_coord_shift());
 			normalize_density(f);
@@ -1028,7 +1035,7 @@ int main(int argc, const char *argv[])
 				else
 					the_poisson_solver->gradient(icoord, u1, data_forIO);
 
-				data_forIO *= cosmo.vfact;
+				data_forIO *= cosmo_vfact;
 
 				double sigv = compute_finest_sigma(data_forIO);
 
@@ -1051,7 +1058,7 @@ int main(int argc, const char *argv[])
 				music::ulog.Print("Writing CDM velocities");
 				the_output_plugin->write_dm_velocity(icoord, data_forIO);
 
-				if (do_baryons && !the_transfer_function_plugin->tf_has_velocities() && !bsph)
+				if (do_baryons && !tf_has_velocities && !bsph)
 				{
 					music::ulog.Print("Writing baryon velocities");
 					the_output_plugin->write_gas_velocity(icoord, data_forIO);
@@ -1061,14 +1068,14 @@ int main(int argc, const char *argv[])
 			if (!dm_only)
 				u1.deallocate();
 
-			if (do_baryons && (the_transfer_function_plugin->tf_has_velocities() || bsph))
+			if (do_baryons && (tf_has_velocities || bsph))
 			{
 				music::ilog << "===============================================================================" << std::endl;
 				music::ilog << "   COMPUTING BARYON VELOCITIES" << std::endl;
 				music::ilog << "-------------------------------------------------------------------------------" << std::endl;
 				music::ulog.Print("Computing baryon displacements...");
 
-				GenerateDensityHierarchy(cf, the_transfer_function_plugin, vbaryon, rh_TF, rand, f, false, bbshift);
+				GenerateDensityHierarchy(cf, the_cosmo_calc.get(), theta_baryon, rh_TF, rand, f, false, bbshift);
 				coarsen_density(rh_Poisson, f, bspectral_sampling);
 				f.add_refinement_mask(rh_Poisson.get_coord_shift());
 				normalize_density(f);
@@ -1126,7 +1133,7 @@ int main(int argc, const char *argv[])
 					else
 						the_poisson_solver->gradient(icoord, u1, data_forIO);
 
-					data_forIO *= cosmo.vfact;
+					data_forIO *= cosmo_vfact;
 
 					double sigv = compute_finest_sigma(data_forIO);
 
@@ -1162,11 +1169,11 @@ int main(int argc, const char *argv[])
 			if (!dm_only)
 			{
 				// my_tf_type is cdm if do_baryons==true, total otherwise
-				my_tf_type = cdm;
+				my_tf_type = delta_cdm;
-				if (!do_baryons || !the_transfer_function_plugin->tf_is_distinct())
+				if (!do_baryons || !the_cosmo_calc->transfer_function_->tf_is_distinct())
-					my_tf_type = total;
+					my_tf_type = delta_matter;
 
-				GenerateDensityHierarchy(cf, the_transfer_function_plugin, my_tf_type, rh_TF, rand, f, false, false);
+				GenerateDensityHierarchy(cf, the_cosmo_calc.get(), my_tf_type, rh_TF, rand, f, false, false);
 				coarsen_density(rh_Poisson, f, bspectral_sampling);
 				f.add_refinement_mask(rh_Poisson.get_coord_shift());
 				normalize_density(f);
@@ -1248,7 +1255,7 @@ int main(int argc, const char *argv[])
 				coarsen_density(rh_Poisson, data_forIO, false);
 
 				// add counter mode
-				if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo.vfact );
+				if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo_vfact );
 
 				music::ulog.Print("Writing CDM displacements");
 				the_output_plugin->write_dm_position(icoord, data_forIO);
@@ -1264,7 +1271,7 @@ int main(int argc, const char *argv[])
 				music::ilog << "-------------------------------------------------------------------------------" << std::endl;
 				music::ulog.Print("Computing baryon density...");
 
-				GenerateDensityHierarchy(cf, the_transfer_function_plugin, baryon, rh_TF, rand, f, true, false);
+				GenerateDensityHierarchy(cf, the_cosmo_calc.get(), delta_baryon, rh_TF, rand, f, true, false);
 				coarsen_density(rh_Poisson, f, bspectral_sampling);
 				f.add_refinement_mask(rh_Poisson.get_coord_shift());
 				normalize_density(f);
@@ -1307,7 +1314,7 @@ int main(int argc, const char *argv[])
 				music::ilog << "-------------------------------------------------------------------------------" << std::endl;
 				music::ulog.Print("Computing baryon displacements...");
 
-				GenerateDensityHierarchy(cf, the_transfer_function_plugin, baryon, rh_TF, rand, f, false, bbshift);
+				GenerateDensityHierarchy(cf, the_cosmo_calc.get(), delta_baryon, rh_TF, rand, f, false, bbshift);
 				coarsen_density(rh_Poisson, f, bspectral_sampling);
 				f.add_refinement_mask(rh_Poisson.get_coord_shift());
 				normalize_density(f);
@@ -1365,7 +1372,7 @@ int main(int argc, const char *argv[])
 					coarsen_density(rh_Poisson, data_forIO, false);
 
 					// add counter mode
-					if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo.vfact );
+					if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo_vfact );
 
 
 					music::ulog.Print("Writing baryon displacements");
@@ -1397,10 +1404,12 @@ int main(int argc, const char *argv[])
 		music::ulog.Print("Wrote output file \'%s\'.", outfname.c_str());
 	}
 
+
+
 	//------------------------------------------------------------------------------
 	//... clean up
 	//------------------------------------------------------------------------------
-	delete the_transfer_function_plugin;
+	// delete the_transfer_function_plugin;
 	delete the_poisson_solver;
 
 	if( CONFIG::FFTW_threads_ok )

src/math/interpolate.hh

@@ -0,0 +1,85 @@
+// This file is part of monofonIC (MUSIC2)
+// A software package to generate ICs for cosmological simulations
+// Copyright (C) 2020 by Oliver Hahn
+// 
+// monofonIC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+// 
+// monofonIC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+// 
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#pragma once
+
+#include <vector>
+#include <cassert>
+#include <gsl/gsl_spline.h>
+#include <gsl/gsl_errno.h>
+
+template <bool logx, bool logy, bool periodic>
+class interpolated_function_1d
+{
+
+private:
+  bool isinit_;
+  std::vector<double> data_x_, data_y_;
+  gsl_interp_accel *gsl_ia_;
+  gsl_spline *gsl_sp_;
+
+  void deallocate()
+  {
+    gsl_spline_free(gsl_sp_);
+    gsl_interp_accel_free(gsl_ia_);
+  }
+
+public:
+  interpolated_function_1d(const interpolated_function_1d &) = delete;
+
+  interpolated_function_1d() : isinit_(false){}
+
+  interpolated_function_1d(const std::vector<double> &data_x, const std::vector<double> &data_y)
+  : isinit_(false)
+  {
+    static_assert(!(logx & periodic),"Class \'interpolated_function_1d\' cannot both be periodic and logarithmic in x!");
+    this->set_data( data_x, data_y );
+  }
+
+  ~interpolated_function_1d()
+  {
+    if (isinit_) this->deallocate();
+  }
+
+  void set_data(const std::vector<double> &data_x, const std::vector<double> &data_y)
+  {
+    data_x_ = data_x;
+    data_y_ = data_y;
+    
+    assert(data_x_.size() == data_y_.size());
+    assert(data_x_.size() > 5);
+
+    if (logx) for (auto &d : data_x_) d = std::log(d);
+    if (logy) for (auto &d : data_y_) d = std::log(d);
+
+    if (isinit_) this->deallocate();
+
+    gsl_ia_ = gsl_interp_accel_alloc();
+    gsl_sp_ = gsl_spline_alloc(periodic ? gsl_interp_cspline_periodic : gsl_interp_cspline, data_x_.size());
+    gsl_spline_init(gsl_sp_, &data_x_[0], &data_y_[0], data_x_.size());
+
+    isinit_ = true;
+  }
+
+  double operator()(double x) const noexcept
+  {
+    assert( isinit_ && !(logx&&x<=0.0) );
+    const double xa = logx ? std::log(x) : x;
+    const double y(gsl_spline_eval(gsl_sp_, xa, gsl_ia_));
+    return logy ? std::exp(y) : y;
+  }
+};

src/math/mat3.hh

@@ -0,0 +1,175 @@
+// This file is part of monofonIC (MUSIC2)
+// A software package to generate ICs for cosmological simulations
+// Copyright (C) 2020 by Oliver Hahn
+// 
+// monofonIC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+// 
+// monofonIC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+// 
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#pragma once
+
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_eigen.h>
+
+#include <math/vec3.hh>
+
+//! class for 3x3 matrix calculations
+template<typename T>
+class mat3_t{
+protected:
+    std::array<T,9> data_;
+    std::array<double,9> data_double_;
+    gsl_matrix_view m_;
+    gsl_vector *eval_;
+    gsl_matrix *evec_;
+	gsl_eigen_symmv_workspace * wsp_;
+    bool bdid_alloc_gsl_;
+						
+    void init_gsl(){
+        // allocate memory for GSL operations if we haven't done so yet
+        if( !bdid_alloc_gsl_ )
+        {
+            if( typeid(T)!=typeid(double) ){
+                m_ = gsl_matrix_view_array ( &data_double_[0], 3, 3);
+            }else{
+                m_ = gsl_matrix_view_array ( (double*)(&data_[0]), 3, 3); // this should only ever be called for T==double so cast is to avoid constexpr ifs from C++17
+            }
+            eval_ = gsl_vector_alloc (3);
+            evec_ = gsl_matrix_alloc (3, 3);
+            wsp_ = gsl_eigen_symmv_alloc (3);
+            bdid_alloc_gsl_ = true;
+        }
+
+        if( typeid(T)!=typeid(double) ){
+            for( int i=0; i<9; ++i ) data_double_[i] = double(data_[i]);
+        }
+    }
+
+    void free_gsl(){
+        // free memory for GSL operations if it was allocated
+        if( bdid_alloc_gsl_ )
+        {
+            gsl_eigen_symmv_free (wsp_);
+            gsl_vector_free (eval_);
+            gsl_matrix_free (evec_);
+        }
+    }
+
+public:
+
+    mat3_t()
+    : bdid_alloc_gsl_(false) 
+    {}
+
+    //! copy constructor
+    mat3_t( const mat3_t<T> &m)
+    : data_(m.data_), bdid_alloc_gsl_(false) 
+    {}
+    
+    //! move constructor
+    mat3_t( mat3_t<T> &&m)
+    : data_(std::move(m.data_)), bdid_alloc_gsl_(false) 
+    {}
+
+    //! construct mat3_t from initializer list
+    template<typename ...E>
+    mat3_t(E&&...e) 
+    : data_{{std::forward<E>(e)...}}, bdid_alloc_gsl_(false)
+    {}
+
+    mat3_t<T>& operator=(const mat3_t<T>& m) noexcept{
+        data_ = m.data_;
+        return *this;
+    }
+
+    mat3_t<T>& operator=(const mat3_t<T>&& m) noexcept{
+        data_ = std::move(m.data_);
+        return *this;
+    }
+
+    //! destructor
+    ~mat3_t(){
+        this->free_gsl();
+    }
+    
+    //! bracket index access to vector components
+    T &operator[](size_t i) noexcept { return data_[i];}
+    
+    //! const bracket index access to vector components
+    const T &operator[](size_t i) const noexcept { return data_[i]; }
+
+    //! matrix 2d index access
+    T &operator()(size_t i, size_t j) noexcept { return data_[3*i+j]; }
+
+    //! const matrix 2d index access
+    const T &operator()(size_t i, size_t j) const noexcept { return data_[3*i+j]; }
+
+    //! in-place addition
+    mat3_t<T>& operator+=( const mat3_t<T>& rhs ) noexcept{
+        for (size_t i = 0; i < 9; ++i) {
+           (*this)[i] += rhs[i];
+        }
+        return *this;
+    }
+
+    //! in-place subtraction
+    mat3_t<T>& operator-=( const mat3_t<T>& rhs ) noexcept{
+        for (size_t i = 0; i < 9; ++i) {
+           (*this)[i] -= rhs[i];
+        }
+        return *this;
+    }
+
+    void zero() noexcept{
+        for (size_t i = 0; i < 9; ++i) data_[i]=0;
+    }
+
+    void eigen( vec3_t<T>& evals, vec3_t<T>& evec1, vec3_t<T>& evec2, vec3_t<T>& evec3_t )
+    {
+        this->init_gsl();
+
+        gsl_eigen_symmv (&m_.matrix, eval_, evec_, wsp_);
+        gsl_eigen_symmv_sort (eval_, evec_, GSL_EIGEN_SORT_VAL_ASC);
+
+        for( int i=0; i<3; ++i ){
+            evals[i] = gsl_vector_get( eval_, i );
+            evec1[i] = gsl_matrix_get( evec_, i, 0 );
+            evec2[i] = gsl_matrix_get( evec_, i, 1 );
+            evec3_t[i] = gsl_matrix_get( evec_, i, 2 );
+        }
+    }
+};
+
+template<typename T>
+constexpr const mat3_t<T> operator+(const mat3_t<T> &lhs, const mat3_t<T> &rhs) noexcept
+{
+    mat3_t<T> result;
+    for (size_t i = 0; i < 9; ++i) {
+        result[i] = lhs[i] + rhs[i];
+    }
+    return result;
+}
+
+// matrix - vector multiplication
+template<typename T>
+inline vec3_t<T> operator*( const mat3_t<T> &A, const vec3_t<T> &v ) noexcept
+{
+    vec3_t<T> result;
+    for( int mu=0; mu<3; ++mu ){
+        result[mu] = 0.0;
+        for( int nu=0; nu<3; ++nu ){
+            result[mu] += A(mu,nu)*v[nu];
+        }
+    }
+    return result;
+}
+

src/math/ode_integrate.hh

@@ -0,0 +1,114 @@
+// This file is part of monofonIC (MUSIC2)
+// A software package to generate ICs for cosmological simulations
+// Copyright (C) 2020 by Oliver Hahn
+// 
+// monofonIC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+// 
+// monofonIC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+// 
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#pragma once
+
+//! ODE integration namespace
+namespace ode_integrate
+{
+
+// simple Runge-Kutta 4th order step without error estimate
+template <typename vector_t, typename function_t>
+inline void rk4_step(double h, double &t, vector_t &y, function_t f)
+{
+    vector_t k1(h * f(t, y));
+    vector_t k2(h * f(t + h / 2, y + k1 / 2));
+    vector_t k3(h * f(t + h / 2, y + k2 / 2));
+    vector_t k4(h * f(t + h, y + k3));
+    y += (k1 + 2 * k2 + 2 * k3 + k4) / 6;
+    t += h;
+}
+
+// Cash-Karp modified Runge-Kutta scheme, 5th order with 4th order error estimate
+// see Press & Teukolsky (1992): "Adaptive Stepsize Runge-Kutta Integration"
+// in Computers in Physics 6, 188 (1992); doi: 10.1063/1.4823060
+template <typename vector_t, typename function_t>
+inline vector_t ckrk5_step(double h, double &t, vector_t &y, function_t f)
+{
+  static constexpr double
+      a2 = 0.20,
+      a3 = 0.30, a4 = 0.60, a5 = 1.0, a6 = 0.8750,
+      b21 = 0.20,
+      b31 = 3.0 / 40.0, b32 = 9.0 / 40.0,
+      b41 = 0.30, b42 = -0.90, b43 = 1.20,
+      b51 = -11.0 / 54.0, b52 = 2.50, b53 = -70.0 / 27.0, b54 = 35.0 / 27.0,
+      b61 = 1631.0 / 55296.0, b62 = 175.0 / 512.0, b63 = 575.0 / 13824.0, b64 = 44275.0 / 110592.0, b65 = 253.0 / 4096.0,
+      c1 = 37.0 / 378.0, c3 = 250.0 / 621.0, c4 = 125.0 / 594.0, c6 = 512.0 / 1771.0,
+      dc1 = c1 - 2825.0 / 27648.0, dc3 = c3 - 18575.0 / 48384.0,
+      dc4 = c4 - 13525.0 / 55296.0, dc5 = -277.0 / 14336.0, dc6 = c6 - 0.250;
+
+  vector_t k1(h * f(t, y));
+  vector_t k2(h * f(t + a2 * h, y + b21 * k1));
+  vector_t k3(h * f(t + a3 * h, y + b31 * k1 + b32 * k2));
+  vector_t k4(h * f(t + a4 * h, y + b41 * k1 + b42 * k2 + b43 * k3));
+  vector_t k5(h * f(t + a5 * h, y + b51 * k1 + b52 * k2 + b53 * k3 + b54 * k4));
+  vector_t k6(h * f(t + a6 * h, y + b61 * k1 + b62 * k2 + b63 * k3 + b64 * k4 + b65 * k5));
+
+  y += c1 * k1 + c3 * k3 + c4 * k4 + c6 * k6;
+
+  return dc1 * k1 + dc3 * k3 + dc4 * k4 + dc5 * k5 + dc6 * k6;
+}
+
+// Adaptive step-size quality-controlled routine for ckrk5_step, see
+// Press & Teukolsky (1992): "Adaptive Stepsize Runge-Kutta Integration"
+// in Computers in Physics 6, 188 (1992); doi: 10.1063/1.4823060
+template <typename vector_t, typename function_t>
+inline void rk_step_qs(double htry, double &t, vector_t &y, vector_t &yscale, function_t f, double eps, double &hdid, double &hnext)
+{
+  static constexpr double SAFETY{0.9};
+  static constexpr double PSHRNK{-0.25};
+  static constexpr double PGROW{-0.2};
+  static constexpr double ERRCON{1.89e-4};
+
+  auto h(htry);
+  vector_t ytemp(y);
+  vector_t yerr;
+  double errmax;
+
+do_ckrk5trialstep:
+  yerr = ckrk5_step(h, t, ytemp, f);
+  errmax = 0.0;
+  for (size_t i = 0; i < yerr.size(); ++i)
+  {
+    errmax = std::max(errmax, std::fabs(yerr[i] / yscale[i]));
+  }
+  errmax = errmax / eps;
+  if (errmax > 1.0)
+  {
+    h *= std::max(0.1, SAFETY*std::pow(errmax, PSHRNK));
+    if (t + h == t)
+    {
+      std::cerr << "stepsize underflow in rkqs" << std::endl;
+      abort();
+    }
+    goto do_ckrk5trialstep;
+  }
+  else
+  {
+    if( errmax > ERRCON ){
+      hnext = h * SAFETY * std::pow(errmax, PGROW);
+    }else{
+      hnext = 5*h;
+    }
+    hdid = h;
+    t += h;
+    y = ytemp;
+  }
+}
+
+
+} // namespace ode_integrate

src/math/vec.hh

@@ -0,0 +1,153 @@
+// This file is part of monofonIC (MUSIC2)
+// A software package to generate ICs for cosmological simulations
+// Copyright (C) 2020 by Oliver Hahn
+// 
+// monofonIC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+// 
+// monofonIC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+// 
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#pragma once
+
+#include <array>
+
+//! implements general N-dim vectors of arbitrary primtive type with some arithmetic ops
+template <int N, typename T = double>
+struct vec_t
+{
+  std::array<T, N> data_;
+
+  vec_t() {}
+
+  vec_t(const vec_t<N, T> &v)
+      : data_(v.data_) {}
+
+  vec_t(vec_t<N, T> &&v)
+      : data_(std::move(v.data_)) {}
+
+  template <typename... E>
+  vec_t(E... e)
+      : data_{{std::forward<E>(e)...}}
+  {
+    static_assert(sizeof...(E) == N, "Brace-enclosed initialiser list doesn't match vec_t length!");
+  }
+
+  //! bracket index access to vector components
+  T &operator[](size_t i) noexcept { return data_[i]; }
+
+  //! const bracket index access to vector components
+  const T &operator[](size_t i) const noexcept { return data_[i]; }
+
+  // assignment operator
+  vec_t<N, T> &operator=(const vec_t<N, T> &v) noexcept
+  {
+    data_ = v.data_;
+    return *this;
+  }
+
+  //! implementation of summation of vec_t
+  vec_t<N, T> operator+(const vec_t<N, T> &v) const noexcept
+  {
+    vec_t<N, T> res;
+    for (int i = 0; i < N; ++i)
+      res[i] = data_[i] + v[i];
+    return res;
+  }
+
+  //! implementation of difference of vec_t
+  vec_t<N, T> operator-(const vec_t<N, T> &v) const noexcept
+  {
+    vec_t<N, T> res;
+    for (int i = 0; i < N; ++i)
+      res[i] = data_[i] - v[i];
+    return res;
+  }
+
+  //! implementation of unary negative
+  vec_t<N, T> operator-() const noexcept
+  {
+    vec_t<N, T> res;
+    for (int i = 0; i < N; ++i)
+      res[i] = -data_[i];
+    return res;
+  }
+
+  //! implementation of scalar multiplication
+  template <typename T2>
+  vec_t<N, T> operator*(T2 s) const noexcept
+  {
+    vec_t<N, T> res;
+    for (int i = 0; i < N; ++i)
+      res[i] = data_[i] * s;
+    return res;
+  }
+
+  //! implementation of scalar division
+  vec_t<N, T> operator/(T s) const noexcept
+  {
+    vec_t<N, T> res;
+    for (int i = 0; i < N; ++i)
+      res[i] = data_[i] / s;
+    return res;
+  }
+
+  //! takes the absolute value of each element
+  vec_t<N, T> abs(void) const noexcept
+  {
+    vec_t<N, T> res;
+    for (int i = 0; i < N; ++i)
+      res[i] = std::fabs(data_[i]);
+    return res;
+  }
+
+  //! implementation of implicit summation of vec_t
+  vec_t<N, T> &operator+=(const vec_t<N, T> &v) noexcept
+  {
+    for (int i = 0; i < N; ++i)
+      data_[i] += v[i];
+    return *this;
+  }
+
+  //! implementation of implicit subtraction of vec_t
+  vec_t<N, T> &operator-=(const vec_t<N, T> &v) noexcept
+  {
+    for (int i = 0; i < N; ++i)
+      data_[i] -= v[i];
+    return *this;
+  }
+
+  //! implementation of implicit scalar multiplication of vec_t
+  vec_t<N, T> &operator*=(T s) noexcept
+  {
+    for (int i = 0; i < N; ++i)
+      data_[i] *= s;
+    return *this;
+  }
+
+  //! implementation of implicit scalar division of vec_t
+  vec_t<N, T> &operator/=(T s) noexcept
+  {
+    for (int i = 0; i < N; ++i)
+      data_[i] /= s;
+    return *this;
+  }
+
+  size_t size(void) const noexcept { return N; }
+};
+
+//! multiplication with scalar
+template <typename T2, int N, typename T = double>
+inline vec_t<N, T> operator*(T2 s, const vec_t<N, T> &v)
+{
+  vec_t<N, T> res;
+  for (int i = 0; i < N; ++i)
+    res[i] = v[i] * s;
+  return res;
+}

src/math/vec3.hh

@@ -0,0 +1,140 @@
+// This file is part of monofonIC (MUSIC2)
+// A software package to generate ICs for cosmological simulations
+// Copyright (C) 2020 by Oliver Hahn
+// 
+// monofonIC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+// 
+// monofonIC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+// 
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#pragma once
+
+//! implements a simple class of 3-vectors of arbitrary scalar type
+template< typename T >
+class vec3_t{
+private:
+    //! holds the data
+    std::array<T,3> data_;
+    
+public: 
+    //! expose access to elements via references
+    T &x,&y,&z;
+
+    //! empty constructor
+    vec3_t()
+    : data_{{T(0),T(0),T(0)}},x(data_[0]),y(data_[1]),z(data_[2]){}
+
+    //! copy constructor
+    vec3_t( const vec3_t<T> &v)
+    : data_(v.data_), x(data_[0]),y(data_[1]),z(data_[2]){}
+
+    //! copy constructor for non-const reference, needed to avoid variadic template being called for non-const reference
+    vec3_t( vec3_t<T>& v)
+    : data_(v.data_), x(data_[0]),y(data_[1]),z(data_[2]){}
+
+    //! move constructor
+    vec3_t( vec3_t<T> &&v)
+    : data_(std::move(v.data_)), x(data_[0]), y(data_[1]), z(data_[2]){}
+
+    //! construct vec3_t from initializer list
+//    template<typename ...E>
+//    vec3_t(E&&...e) 
+//    : data_{{std::forward<E>(e)...}}, x{data_[0]}, y{data_[1]}, z{data_[2]}
+//    {}
+    vec3_t( T a, T b, T c )
+        : data_{{a,b,c}}, x(data_[0]), y(data_[1]), z(data_[2]){}
+  
+    explicit vec3_t( T a)
+        : data_{{a,a,a}}, x(data_[0]), y(data_[1]), z(data_[2]){}
+    
+    //! bracket index access to vector components
+    T &operator[](size_t i) noexcept{ return data_[i];}
+    
+    //! const bracket index access to vector components
+    const T &operator[](size_t i) const noexcept { return data_[i]; }
+
+    // assignment operator
+    vec3_t<T>& operator=( const vec3_t<T>& v ) noexcept { data_=v.data_; return *this; }
+
+    //! implementation of summation of vec3_t
+    vec3_t<T> operator+( const vec3_t<T>& v ) const noexcept{ return vec3_t<T>({x+v.x,y+v.y,z+v.z}); }
+
+    //! implementation of difference of vec3_t
+    vec3_t<T> operator-( const vec3_t<T>& v ) const noexcept{ return vec3_t<T>({x-v.x,y-v.y,z-v.z}); }
+
+    //! implementation of unary negative
+    vec3_t<T> operator-() const noexcept{ return vec3_t<T>({-x,-y,-z}); }
+
+    //! implementation of scalar multiplication
+    vec3_t<T> operator*( T s ) const noexcept{ return vec3_t<T>({x*s,y*s,z*s}); }
+
+    //! implementation of scalar division
+    vec3_t<T> operator/( T s ) const noexcept{ return vec3_t<T>({x/s,y/s,z/s}); }
+
+    //! implementation of += operator
+    vec3_t<T>& operator+=( const vec3_t<T>& v ) noexcept{ x+=v.x; y+=v.y; z+=v.z; return *this; }
+
+    //! implementation of -= operator
+    vec3_t<T>& operator-=( const vec3_t<T>& v ) noexcept{ x-=v.x; y-=v.y; z-=v.z; return *this; }
+
+    //! multiply with scalar
+    vec3_t<T>& operator*=( T s ) noexcept{ x*=s; y*=s; z*=s; return *this; }
+    
+    //! divide by scalar
+    vec3_t<T>& operator/=( T s ) noexcept{ x/=s; y/=s; z/=s; return *this; }
+
+    //! compute dot product with another vector
+    T dot(const vec3_t<T> &a) const noexcept
+    {
+        return data_[0] * a.data_[0] + data_[1] * a.data_[1] + data_[2] * a.data_[2];
+    }
+    
+    //! returns 2-norm squared of vector
+    T norm_squared(void) const noexcept { return this->dot(*this); }
+
+    //! returns 2-norm of vector
+    T norm(void) const noexcept { return std::sqrt( this->norm_squared() ); }
+
+    //! wrap absolute vector to box of size p
+    vec3_t<T>& wrap_abs( T p = 1.0 ) noexcept{
+        for( auto& x : data_ ) x = std::fmod( 2*p + x, p );
+        return *this;
+    }
+
+    //! wrap relative vector to box of size p
+    vec3_t<T>& wrap_rel( T p = 1.0 ) noexcept{
+        for( auto& x : data_ ) x = (x<-p/2)? x+p : (x>=p/2)? x-p : x;
+        return *this;
+    }
+
+    //! takes the absolute value of each element
+    vec3_t<T> abs(void) const noexcept
+    {
+        vec3_t<T> res;
+        for (int i = 0; i < 3; ++i)
+            res[i] = std::fabs(data_[i]);
+        return res;
+    }
+
+    //! ordering, allows 3d sorting of vec3_ts
+    bool operator<( const vec3_t<T>& o ) const noexcept{
+        if( x!=o.x ) return x<o.x?true:false;
+        if( y!=o.y ) return y<o.y?true:false;
+        if( z!=o.z ) return z<o.z?true:false;
+        return false;
+    }
+};
+
+//! multiplication with scalar
+template<typename T>
+vec3_t<T> operator*( T s, const vec3_t<T>& v ){
+    return vec3_t<T>({v.x*s,v.y*s,v.z*s});
+}

src/perturbation_theory.cc

@@ -8,10 +8,10 @@
 
  */
 
-#include "cosmology.hh"
+#include <perturbation_theory.hh>
-#include "mesh.hh"
+#include <mesh.hh>
-#include "mg_operators.hh"
+#include <mg_operators.hh>
-#include "general.hh"
+#include <general.hh>
 
 #define ACC(i, j, k) ((*u.get_grid((ilevel)))((i), (j), (k)))
 #define SQR(x) ((x) * (x))

src/perturbation_theory.hh

@@ -0,0 +1,40 @@
+/*
+ 
+ cosmology.hh - This file is part of MUSIC -
+ a code to generate multi-scale initial conditions 
+ for cosmological simulations 
+ 
+ Copyright (C) 2010  Oliver Hahn
+ 
+ */
+
+#pragma once
+
+#include "mesh.hh"
+#include "general.hh"
+
+
+//! compute the jeans sound speed
+/*! given a density in g/cm^-3 and a mass in g it gives back the sound
+ *  speed in cm/s for which the input mass is equal to the jeans mass
+ *  @param rho density 
+ *  @param mass mass scale
+ *  @returns jeans sound speed
+ */
+inline double jeans_sound_speed( double rho, double mass )
+{
+	const double G = 6.67e-8;
+	return pow( 6.0*mass/M_PI*sqrt(rho)*pow(G,1.5), 1.0/3.0 );
+}
+
+//! computes the density from the potential using the Laplacian
+void compute_Lu_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order=4 );
+
+//! computes the 2nd order density perturbations using also off-diagonal terms in the potential Hessian 
+void compute_LLA_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order=4 );
+
+//! computes the source term for the 2nd order perturbations in the displacements
+void compute_2LPT_source( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order=4 );
+
+void compute_2LPT_source_FFT( config_file& cf_, const grid_hierarchy& u, grid_hierarchy& fnew );
+

src/plugins/random_panphasia.cc

@@ -483,7 +483,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
         size_t idx = ((size_t)i * nxremap[1] + (size_t)j) * (nxremap[2] / 2 + 1) + (size_t)k;
 
         double fx(1.0), fy(1.0), fz(1.0), arg = 0.;
-        complex gx(0., 0.), gy(0., 0.), gz(0., 0.);
+        ccomplex_t gx(0., 0.), gy(0., 0.), gz(0., 0.);
 
         int ii(i), jj(j), kk(k);
         if (i > nxremap[0] / 2)
@@ -497,7 +497,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
         {
           arg = M_PI * (double)ii / (double)nxremap[0];
           fx = sin(arg) / arg;
-          gx = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
+          gx = ccomplex_t(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
         }
         else
         {
@@ -509,7 +509,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
         {
           arg = M_PI * (double)jj / (double)nxremap[1];
           fy = sin(arg) / arg;
-          gy = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
+          gy = ccomplex_t(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
         }
         else
         {
@@ -521,7 +521,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
         {
           arg = M_PI * (double)kk / (double)nxremap[2];
           fz = sin(arg) / arg;
-          gz = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
+          gz = ccomplex_t(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
         }
         else
         {
@@ -529,13 +529,13 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
           gz = 0.0;
         }
 
-        complex temp_comp = (fx + sqrt(3.0) * gx) * (fy + sqrt(3.0) * gy) * (fz + sqrt(3.0) * gz);
+        ccomplex_t temp_comp = (fx + sqrt(3.0) * gx) * (fy + sqrt(3.0) * gy) * (fz + sqrt(3.0) * gz);
         double magnitude = sqrt(1.0 - std::abs(temp_comp * temp_comp));
 
         if (abs(ii) != nxremap[0] / 2 && abs(jj) != nxremap[1] / 2 &&
             abs(kk) != nxremap[2] / 2)
         { // kkmax != nxremap[2]/2 ){
-          complex x, y0(RE(pc0[idx]), IM(pc0[idx])), y1(RE(pc1[idx]), IM(pc1[idx])), y2(RE(pc2[idx]), IM(pc2[idx])),
+          ccomplex_t x, y0(RE(pc0[idx]), IM(pc0[idx])), y1(RE(pc1[idx]), IM(pc1[idx])), y2(RE(pc2[idx]), IM(pc2[idx])),
               y3(RE(pc3[idx]), IM(pc3[idx])), y4(RE(pc4[idx]), IM(pc4[idx]));
 
           x = y0 * fx * fy * fz + sqrt(3.0) * (y1 * gx * fy * fz + y2 * fx * gy * fz + y3 * fx * fy * gz) +
@@ -670,7 +670,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
         size_t idx = ((size_t)i * nxremap[1] + (size_t)j) * (nxremap[2] / 2 + 1) + (size_t)k;
 
         double fx(1.0), fy(1.0), fz(1.0), arg = 0.;
-        complex gx(0., 0.), gy(0., 0.), gz(0., 0.);
+        ccomplex_t gx(0., 0.), gy(0., 0.), gz(0., 0.);
 
         int ii(i), jj(j), kk(k);
         if (i > nxremap[0] / 2)
@@ -684,7 +684,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
         {
           arg = M_PI * (double)ii / (double)nxremap[0];
           fx = sin(arg) / arg;
-          gx = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
+          gx = ccomplex_t(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
         }
         else
         {
@@ -696,7 +696,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
         {
           arg = M_PI * (double)jj / (double)nxremap[1];
           fy = sin(arg) / arg;
-          gy = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
+          gy = ccomplex_t(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
         }
         else
         {
@@ -708,7 +708,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
         {
           arg = M_PI * (double)kk / (double)nxremap[2];
           fz = sin(arg) / arg;
-          gz = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
+          gz = ccomplex_t(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
         }
         else
         {
@@ -719,7 +719,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
         if (abs(ii) != nxremap[0] / 2 && abs(jj) != nxremap[1] / 2 &&
             abs(kk) != nxremap[2] / 2)
         { // kkmax != nxremap[2]/2 ){
-          complex x, y1(RE(pc1[idx]), IM(pc1[idx])), y2(RE(pc2[idx]), IM(pc2[idx])), y3(RE(pc3[idx]), IM(pc3[idx])),
+          ccomplex_t x, y1(RE(pc1[idx]), IM(pc1[idx])), y2(RE(pc2[idx]), IM(pc2[idx])), y3(RE(pc3[idx]), IM(pc3[idx])),
               y4(RE(pc4[idx]), IM(pc4[idx]));
 
           x = 3.0 * (y1 * gx * gy * fz + y2 * fx * gy * gz + y3 * gx * fy * gz) + sqrt(27.0) * y4 * gx * gy * gz;
@@ -767,10 +767,10 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
 
             size_t idx2 = ((size_t)ir * nx_m[1] + (size_t)jr) * ((size_t)nx_m[2] / 2 + 1) + (size_t)kr;
 
-            complex x(RE(pc0[idx]), IM(pc0[idx]));
+            ccomplex_t x(RE(pc0[idx]), IM(pc0[idx]));
             double total_phase_shift;
             total_phase_shift = inter_grid_phase_adjustment_ * (double)(ii + jj + kk);
-            x = x * exp(complex(0.0, total_phase_shift));
+            x = x * exp(ccomplex_t(0.0, total_phase_shift));
             RE(pc1[idx2]) = x.real();
             IM(pc1[idx2]) = x.imag();
           }

src/plugins/transfer_bbks.cc

@@ -26,10 +26,12 @@ public:
 	 \param aCosm Structure of type Cosmology carrying the cosmological parameters
 	 \param bSugiyama flag whether the Sugiyama (1995) correction shall be applied (default=true)
 	 */
-	transfer_bbks_plugin( config_file& cf )
+	transfer_bbks_plugin( config_file& cf, const cosmology::parameters& cp )
-    : transfer_function_plugin( cf )
+    : transfer_function_plugin( cf, cp )
 	{  
-		double Omega0 = cosmo_.Omega_m;
+		const double Omega0 = cp["Omega_m"];
+		const double Omegab = cp["Omega_b"];
+		const double H0 = cp["H0"];
 		double FreeGamma = -1.0;
 		
 		bool bSugiyama(true);
@@ -43,9 +45,9 @@ public:
 		FreeGamma = pcf_->get_value_safe<double>( "cosmology", "gamma", FreeGamma );
 		
 		if( FreeGamma <= 0.0 ){
-			m_Gamma = Omega0*0.01*cosmo_.H0;
+			m_Gamma = Omega0*0.01*H0;
 			if( bSugiyama )
-				m_Gamma *= exp(-cosmo_.Omega_b*(1.0+sqrt(2.0*0.01*cosmo_.H0)/Omega0));
+				m_Gamma *= exp(-Omegab*(1.0+sqrt(2.0*0.01*H0)/Omega0));
 		}else
 			m_Gamma = FreeGamma;
 		

src/plugins/transfer_camb.cc

@@ -18,9 +18,9 @@ private:
   std::vector<double> m_tab_k, m_tab_Tk_tot, m_tab_Tk_cdm, m_tab_Tk_baryon;
   std::vector<double> m_tab_Tvk_tot, m_tab_Tvk_cdm, m_tab_Tvk_baryon;
   gsl_interp_accel *acc_tot, *acc_cdm, *acc_baryon;
-  gsl_interp_accel *acc_vtot, *acc_vcdm, *acc_vbaryon;
+  gsl_interp_accel *acc_vtot, *acc_theta_cdm, *acc_theta_baryon;
   gsl_spline *spline_tot, *spline_cdm, *spline_baryon;
-  gsl_spline *spline_vtot, *spline_vcdm, *spline_vbaryon;
+  gsl_spline *spline_vtot, *spline_theta_cdm, *spline_theta_baryon;
 
   double m_kmin, m_kmax, m_Omega_b, m_Omega_m, m_zstart;
   unsigned m_nlines;
@@ -162,12 +162,13 @@ private:
   }
 
 public:
-  transfer_CAMB_plugin(config_file &cf)
+  transfer_CAMB_plugin(config_file &cf, const cosmology::parameters& cp)
-  : transfer_function_plugin(cf)
+  : transfer_function_plugin(cf,cp)
   {
     m_filename_Tk = pcf_->get_value<std::string>("cosmology", "transfer_file");
-    m_Omega_m=cf.get_value<double>("cosmology","Omega_m"); //MvD
+		
-    m_Omega_b=cf.get_value<double>("cosmology","Omega_b"); //MvD
+    m_Omega_m=cp["Omega_m"]; //MvD
+    m_Omega_b=cp["Omega_b"]; //MvD
     m_zstart =cf.get_value<double>("setup","zstart"); //MvD
 
     read_table();
@@ -176,16 +177,16 @@ public:
     acc_cdm = gsl_interp_accel_alloc();
     acc_baryon = gsl_interp_accel_alloc();
     acc_vtot = gsl_interp_accel_alloc();    //>[150609SH: add]
-    acc_vcdm = gsl_interp_accel_alloc();    //>[150609SH: add]
+    acc_theta_cdm = gsl_interp_accel_alloc();    //>[150609SH: add]
-    acc_vbaryon = gsl_interp_accel_alloc(); //>[150609SH: add]
+    acc_theta_baryon = gsl_interp_accel_alloc(); //>[150609SH: add]
 
     spline_tot = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
     spline_cdm = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
     spline_baryon = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
     spline_vtot = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
-    spline_vcdm =
+    spline_theta_cdm =
         gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); //>[150609SH: add]
-    spline_vbaryon =
+    spline_theta_baryon =
         gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); //>[150609SH: add]
 
     gsl_spline_init(spline_tot, &m_tab_k[0], &m_tab_Tk_tot[0], m_tab_k.size());
@@ -194,9 +195,9 @@ public:
                     m_tab_k.size());
     gsl_spline_init(spline_vtot, &m_tab_k[0], &m_tab_Tvk_tot[0],
                     m_tab_k.size()); //>[150609SH: add]
-    gsl_spline_init(spline_vcdm, &m_tab_k[0], &m_tab_Tvk_cdm[0],
+    gsl_spline_init(spline_theta_cdm, &m_tab_k[0], &m_tab_Tvk_cdm[0],
                     m_tab_k.size()); //>[150609SH: add]
-    gsl_spline_init(spline_vbaryon, &m_tab_k[0], &m_tab_Tvk_baryon[0],
+    gsl_spline_init(spline_theta_baryon, &m_tab_k[0], &m_tab_Tvk_baryon[0],
                     m_tab_k.size()); //>[150609SH: add]
 
     tf_distinct_ = true; // [150612SH: different density between CDM v.s. Baryon]
@@ -208,15 +209,15 @@ public:
     gsl_spline_free(spline_cdm);
     gsl_spline_free(spline_baryon);
     gsl_spline_free(spline_vtot);
-    gsl_spline_free(spline_vcdm);    //>[150609SH: add]
+    gsl_spline_free(spline_theta_cdm);    //>[150609SH: add]
-    gsl_spline_free(spline_vbaryon); //>[150609SH: add]
+    gsl_spline_free(spline_theta_baryon); //>[150609SH: add]
 
     gsl_interp_accel_free(acc_tot);
     gsl_interp_accel_free(acc_cdm);
     gsl_interp_accel_free(acc_baryon);
     gsl_interp_accel_free(acc_vtot);
-    gsl_interp_accel_free(acc_vcdm);    //>[150609SH: add]
+    gsl_interp_accel_free(acc_theta_cdm);    //>[150609SH: add]
-    gsl_interp_accel_free(acc_vbaryon); //>[150609SH: add]
+    gsl_interp_accel_free(acc_theta_baryon); //>[150609SH: add]
   }
 
   // linear interpolation in log-log
@@ -230,31 +231,31 @@ public:
     double delk = lk - m_tab_k[n];
 
     switch (type) {
-    case cdm:
+    case delta_cdm:
       v1 = m_tab_Tk_cdm[n1];
       v2 = m_tab_Tk_cdm[n];
       return pow(10.0, (v2 - v1) / dk * (delk) + v2);
-    case baryon:
+    case delta_baryon:
       v1 = m_tab_Tk_baryon[n1];
       v2 = m_tab_Tk_baryon[n];
       if (m_linbaryoninterp)
         return std::max((v2 - v1) / dk * (delk) + v2, tiny);
       return pow(10.0, (v2 - v1) / dk * (delk) + v2);
-    case vtotal: //>[150609SH: add]
+    case theta_matter: //>[150609SH: add]
       v1 = m_tab_Tvk_tot[n1];
       v2 = m_tab_Tvk_tot[n];
       return pow(10.0, (v2 - v1) / dk * (delk) + v2);
-    case vcdm: //>[150609SH: add]
+    case theta_cdm: //>[150609SH: add]
       v1 = m_tab_Tvk_cdm[n1];
       v2 = m_tab_Tvk_cdm[n];
       return pow(10.0, (v2 - v1) / dk * (delk) + v2);
-    case vbaryon: //>[150609SH: add]
+    case theta_baryon: //>[150609SH: add]
       v1 = m_tab_Tvk_baryon[n1];
       v2 = m_tab_Tvk_baryon[n];
       if (m_linbaryoninterp)
         return std::max((v2 - v1) / dk * (delk) + v2, tiny);
       return pow(10.0, (v2 - v1) / dk * (delk) + v2);
-    case total:
+    case delta_matter:
       v1 = m_tab_Tk_tot[n1];
       v2 = m_tab_Tk_tot[n];
       return pow(10.0, (v2 - v1) / dk * (delk) + v2);
@@ -270,21 +271,21 @@ public:
     // use constant interpolation on the left side of the tabulated values
     if (k < m_kmin) {
       switch (type) {
-      case cdm:
+      case delta_cdm:
         return pow(10.0, m_tab_Tk_cdm[0]);
-      case baryon:
+      case delta_baryon:
         if (m_linbaryoninterp)
           return m_tab_Tk_baryon[0];
         return pow(10.0, m_tab_Tk_baryon[0]);
-      case vtotal:
+      case theta_matter:
         return pow(10.0, m_tab_Tvk_tot[0]);
-      case vcdm:
+      case theta_cdm:
         return pow(10.0, m_tab_Tvk_cdm[0]);
-      case vbaryon:
+      case theta_baryon:
         if (m_linbaryoninterp)
           return m_tab_Tvk_baryon[0];
         return pow(10.0, m_tab_Tvk_baryon[0]);
-      case total:
+      case delta_matter:
         return pow(10.0, m_tab_Tk_tot[0]);
       default:
         throw std::runtime_error(
@@ -297,21 +298,21 @@ public:
 
     double lk = log10(k);
     switch (type) {
-    case cdm:
+    case delta_cdm:
       return pow(10.0, gsl_spline_eval(spline_cdm, lk, acc_cdm));
-    case baryon:
+    case delta_baryon:
       if (m_linbaryoninterp)
         return gsl_spline_eval(spline_baryon, lk, acc_baryon);
       return pow(10.0, gsl_spline_eval(spline_baryon, lk, acc_baryon));
-    case vtotal:
+    case theta_matter:
       return pow(10.0, gsl_spline_eval(spline_vtot, lk, acc_vtot)); //MvD
-    case vcdm:
+    case theta_cdm:
-      return pow(10.0, gsl_spline_eval(spline_vcdm, lk, acc_vcdm));
+      return pow(10.0, gsl_spline_eval(spline_theta_cdm, lk, acc_theta_cdm));
-    case vbaryon:
+    case theta_baryon:
       if (m_linbaryoninterp)
-        return gsl_spline_eval(spline_vbaryon, lk, acc_vbaryon);
+        return gsl_spline_eval(spline_theta_baryon, lk, acc_theta_baryon);
-      return pow(10.0, gsl_spline_eval(spline_vbaryon, lk, acc_vbaryon));
+      return pow(10.0, gsl_spline_eval(spline_theta_baryon, lk, acc_theta_baryon));
-    case total:
+    case delta_matter:
       return pow(10.0, gsl_spline_eval(spline_tot, lk, acc_tot));
     default:
       throw std::runtime_error(

src/plugins/transfer_eisenstein.cc

@@ -1,9 +1,9 @@
 /*
- 
+
  transfer_eisenstein.cc - This file is part of MUSIC -
- a code to generate multi-scale initial conditions 
+ a code to generate multi-scale initial conditions
- for cosmological simulations 
+ for cosmological simulations
- 
+
  */
 
 #include "transfer_function.hh"
@@ -11,30 +11,29 @@
 // forward declaration of WDM class
 class transfer_eisenstein_wdm_plugin;
 
-
 struct eisenstein_transfer
 {
-  //Cosmology m_Cosmology;
+  // Cosmology m_Cosmology;
-  double  m_h0;
+  double m_h0;
-  double	omhh,		/* Omega_matter*h^2 */
+  double omhh,           /* Omega_matter*h^2 */
-  obhh,		/* Omega_baryon*h^2 */
+      obhh,              /* Omega_baryon*h^2 */
-  theta_cmb,	/* Tcmb in units of 2.7 K */
+      theta_cmb,         /* Tcmb in units of 2.7 K */
-  z_equality,	/* Redshift of matter-radiation equality, really 1+z */
+      z_equality,        /* Redshift of matter-radiation equality, really 1+z */
-  k_equality,	/* Scale of equality, in Mpc^-1 */
+      k_equality,        /* Scale of equality, in Mpc^-1 */
-  z_drag,		/* Redshift of drag epoch */
+      z_drag,            /* Redshift of drag epoch */
-  R_drag,		/* Photon-baryon ratio at drag epoch */
+      R_drag,            /* Photon-baryon ratio at drag epoch */
-  R_equality,	/* Photon-baryon ratio at equality epoch */
+      R_equality,        /* Photon-baryon ratio at equality epoch */
-  sound_horizon,	/* Sound horizon at drag epoch, in Mpc */
+      sound_horizon,     /* Sound horizon at drag epoch, in Mpc */
-  k_silk,		/* Silk damping scale, in Mpc^-1 */
+      k_silk,            /* Silk damping scale, in Mpc^-1 */
-  alpha_c,	/* CDM suppression */
+      alpha_c,           /* CDM suppression */
-  beta_c,		/* CDM log shift */
+      beta_c,            /* CDM log shift */
-  alpha_b,	/* Baryon suppression */
+      alpha_b,           /* Baryon suppression */
-  beta_b,		/* Baryon envelope shift */
+      beta_b,            /* Baryon envelope shift */
-  beta_node,	/* Sound horizon shift */
+      beta_node,         /* Sound horizon shift */
-  k_peak,		/* Fit to wavenumber of first peak, in Mpc^-1 */
+      k_peak,            /* Fit to wavenumber of first peak, in Mpc^-1 */
-  sound_horizon_fit,	/* Fit to sound horizon, in Mpc */
+      sound_horizon_fit, /* Fit to sound horizon, in Mpc */
-  alpha_gamma;	/* Gamma suppression in approximate TF */
+      alpha_gamma;       /* Gamma suppression in approximate TF */
-  
+
   //! private member function: sets internal quantities for Eisenstein & Hu fitting
   void TFset_parameters(double omega0hh, double f_baryon, double Tcmb)
   /* Set all the scalars quantities for Eisenstein & Hu 1997 fitting formula */
@@ -50,57 +49,58 @@ struct eisenstein_transfer
   {
     double z_drag_b1, z_drag_b2;
     double alpha_c_a1, alpha_c_a2, beta_c_b1, beta_c_b2, alpha_b_G, y;
-    
+
-    if (f_baryon<=0.0 || omega0hh<=0.0) {
+    if (f_baryon <= 0.0 || omega0hh <= 0.0)
+    {
       fprintf(stderr, "TFset_parameters(): Illegal input.\n");
       exit(1);
     }
     omhh = omega0hh;
-    obhh = omhh*f_baryon;
+    obhh = omhh * f_baryon;
-    if (Tcmb<=0.0) Tcmb=2.728;	/* COBE FIRAS */
+    if (Tcmb <= 0.0)
-    theta_cmb = Tcmb/2.7;
+      Tcmb = 2.728; /* COBE FIRAS */
-    
+    theta_cmb = Tcmb / 2.7;
-    z_equality = 2.50e4*omhh/POW4(theta_cmb);  /* Really 1+z */
+
-    k_equality = 0.0746*omhh/SQR(theta_cmb);
+    z_equality = 2.50e4 * omhh / POW4(theta_cmb); /* Really 1+z */
-    
+    k_equality = 0.0746 * omhh / SQR(theta_cmb);
-    z_drag_b1 = 0.313*pow((double)omhh,-0.419)*(1+0.607*pow((double)omhh,0.674));
+
-    z_drag_b2 = 0.238*pow((double)omhh,0.223);
+    z_drag_b1 = 0.313 * pow((double)omhh, -0.419) * (1 + 0.607 * pow((double)omhh, 0.674));
-    z_drag = 1291*pow(omhh,0.251)/(1+0.659*pow((double)omhh,0.828))*
+    z_drag_b2 = 0.238 * pow((double)omhh, 0.223);
-    (1+z_drag_b1*pow((double)obhh,(double)z_drag_b2));
+    z_drag = 1291 * pow(omhh, 0.251) / (1 + 0.659 * pow((double)omhh, 0.828)) *
-    
+             (1 + z_drag_b1 * pow((double)obhh, (double)z_drag_b2));
-    R_drag = 31.5*obhh/POW4(theta_cmb)*(1000/(1+z_drag));
+
-    R_equality = 31.5*obhh/POW4(theta_cmb)*(1000/z_equality);
+    R_drag = 31.5 * obhh / POW4(theta_cmb) * (1000 / (1 + z_drag));
-    
+    R_equality = 31.5 * obhh / POW4(theta_cmb) * (1000 / z_equality);
-    sound_horizon = 2./3./k_equality*sqrt(6./R_equality)*
+
-    log((sqrt(1+R_drag)+sqrt(R_drag+R_equality))/(1+sqrt(R_equality)));
+    sound_horizon = 2. / 3. / k_equality * sqrt(6. / R_equality) *
-    
+                    log((sqrt(1 + R_drag) + sqrt(R_drag + R_equality)) / (1 + sqrt(R_equality)));
-    k_silk = 1.6*pow((double)obhh,0.52)*pow((double)omhh,0.73)*(1+pow((double)10.4*omhh,-0.95));
+
-    
+    k_silk = 1.6 * pow((double)obhh, 0.52) * pow((double)omhh, 0.73) * (1 + pow((double)10.4 * omhh, -0.95));
-    alpha_c_a1 = pow((double)46.9*omhh,0.670)*(1+pow(32.1*omhh,-0.532));
+
-    alpha_c_a2 = pow((double)12.0*omhh,0.424)*(1+pow(45.0*omhh,-0.582));
+    alpha_c_a1 = pow((double)46.9 * omhh, 0.670) * (1 + pow(32.1 * omhh, -0.532));
-    alpha_c = pow(alpha_c_a1,-f_baryon)*
+    alpha_c_a2 = pow((double)12.0 * omhh, 0.424) * (1 + pow(45.0 * omhh, -0.582));
-    pow(alpha_c_a2,-CUBE(f_baryon));
+    alpha_c = pow(alpha_c_a1, -f_baryon) *
-    
+              pow(alpha_c_a2, -CUBE(f_baryon));
-    beta_c_b1 = 0.944/(1+pow(458*omhh,-0.708));
+
-    beta_c_b2 = pow(0.395*omhh, -0.0266);
+    beta_c_b1 = 0.944 / (1 + pow(458 * omhh, -0.708));
-    beta_c = 1.0/(1+beta_c_b1*(pow(1-f_baryon, beta_c_b2)-1));
+    beta_c_b2 = pow(0.395 * omhh, -0.0266);
-    
+    beta_c = 1.0 / (1 + beta_c_b1 * (pow(1 - f_baryon, beta_c_b2) - 1));
-    y = z_equality/(1+z_drag);
+
-    alpha_b_G = y*(-6.*sqrt(1+y)+(2.+3.*y)*log((sqrt(1+y)+1)/(sqrt(1+y)-1)));
+    y = z_equality / (1 + z_drag);
-    alpha_b = 2.07*k_equality*sound_horizon*pow(1+R_drag,-0.75)*alpha_b_G;
+    alpha_b_G = y * (-6. * sqrt(1 + y) + (2. + 3. * y) * log((sqrt(1 + y) + 1) / (sqrt(1 + y) - 1)));
-    
+    alpha_b = 2.07 * k_equality * sound_horizon * pow(1 + R_drag, -0.75) * alpha_b_G;
-    beta_node = 8.41*pow(omhh, 0.435);
+
-    beta_b = 0.5+f_baryon+(3.-2.*f_baryon)*sqrt(pow(17.2*omhh,2.0)+1);
+    beta_node = 8.41 * pow(omhh, 0.435);
-    
+    beta_b = 0.5 + f_baryon + (3. - 2. * f_baryon) * sqrt(pow(17.2 * omhh, 2.0) + 1);
-    k_peak = 2.5*3.14159*(1+0.217*omhh)/sound_horizon;
+
-    sound_horizon_fit = 44.5*log(9.83/omhh)/sqrt(1+10.0*pow(obhh,0.75));
+    k_peak = 2.5 * 3.14159 * (1 + 0.217 * omhh) / sound_horizon;
-    
+    sound_horizon_fit = 44.5 * log(9.83 / omhh) / sqrt(1 + 10.0 * pow(obhh, 0.75));
-    alpha_gamma = 1-0.328*log(431.0*omhh)*f_baryon + 0.38*log(22.3*omhh)*
+
-    SQR(f_baryon);
+    alpha_gamma = 1 - 0.328 * log(431.0 * omhh) * f_baryon + 0.38 * log(22.3 * omhh) * SQR(f_baryon);
-    
+
     return;
   }
-  
+
   //! private member function: computes transfer function for mode k (k in Mpc)
   inline double TFfit_onek(double k, double *tf_baryon, double *tf_cdm)
   /* Input: k -- Wavenumber at which to calculate transfer function, in Mpc^-1.
@@ -113,71 +113,75 @@ struct eisenstein_transfer
   /* Notes: Units are Mpc, not h^-1 Mpc. */
   {
     double T_c_ln_beta, T_c_ln_nobeta, T_c_C_alpha, T_c_C_noalpha;
-    double q, xx, xx_tilde;//, q_eff;
+    double q, xx, xx_tilde; //, q_eff;
     double T_c_f, T_c, s_tilde, T_b_T0, T_b, f_baryon, T_full;
-    //double T_0_L0, T_0_C0, T_0, gamma_eff;
+    // double T_0_L0, T_0_C0, T_0, gamma_eff;
-    //double T_nowiggles_L0, T_nowiggles_C0, T_nowiggles;
+    // double T_nowiggles_L0, T_nowiggles_C0, T_nowiggles;
-    
+
-    k = fabs(k);	/* Just define negative k as positive */
+    k = fabs(k); /* Just define negative k as positive */
-    if (k==0.0) {
+    if (k == 0.0)
-      if (tf_baryon!=NULL) *tf_baryon = 1.0;
+    {
-      if (tf_cdm!=NULL) *tf_cdm = 1.0;
+      if (tf_baryon != NULL)
+        *tf_baryon = 1.0;
+      if (tf_cdm != NULL)
+        *tf_cdm = 1.0;
       return 1.0;
     }
-    
+
-    q = k/13.41/k_equality;
+    q = k / 13.41 / k_equality;
-    xx = k*sound_horizon;
+    xx = k * sound_horizon;
-    
+
-    T_c_ln_beta = log(2.718282+1.8*beta_c*q);
+    T_c_ln_beta = log(2.718282 + 1.8 * beta_c * q);
-    T_c_ln_nobeta = log(2.718282+1.8*q);
+    T_c_ln_nobeta = log(2.718282 + 1.8 * q);
-    T_c_C_alpha = 14.2/alpha_c + 386.0/(1+69.9*pow(q,1.08));
+    T_c_C_alpha = 14.2 / alpha_c + 386.0 / (1 + 69.9 * pow(q, 1.08));
-    T_c_C_noalpha = 14.2 + 386.0/(1+69.9*pow(q,1.08));
+    T_c_C_noalpha = 14.2 + 386.0 / (1 + 69.9 * pow(q, 1.08));
-    
+
-    T_c_f = 1.0/(1.0+POW4(xx/5.4));
+    T_c_f = 1.0 / (1.0 + POW4(xx / 5.4));
-    T_c = T_c_f*T_c_ln_beta/(T_c_ln_beta+T_c_C_noalpha*SQR(q)) +
+    T_c = T_c_f * T_c_ln_beta / (T_c_ln_beta + T_c_C_noalpha * SQR(q)) +
-    (1-T_c_f)*T_c_ln_beta/(T_c_ln_beta+T_c_C_alpha*SQR(q));
+          (1 - T_c_f) * T_c_ln_beta / (T_c_ln_beta + T_c_C_alpha * SQR(q));
-    
+
-    s_tilde = sound_horizon*pow(1.+CUBE(beta_node/xx),-1./3.);
+    s_tilde = sound_horizon * pow(1. + CUBE(beta_node / xx), -1. / 3.);
-    xx_tilde = k*s_tilde;
+    xx_tilde = k * s_tilde;
-    
+
-    T_b_T0 = T_c_ln_nobeta/(T_c_ln_nobeta+T_c_C_noalpha*SQR(q));
+    T_b_T0 = T_c_ln_nobeta / (T_c_ln_nobeta + T_c_C_noalpha * SQR(q));
-    T_b = sin(xx_tilde)/(xx_tilde)*(T_b_T0/(1.+SQR(xx/5.2))+
+    T_b = sin(xx_tilde) / (xx_tilde) * (T_b_T0 / (1. + SQR(xx / 5.2)) + alpha_b / (1. + CUBE(beta_b / xx)) * exp(-pow(k / k_silk, 1.4)));
-                                    alpha_b/(1.+CUBE(beta_b/xx))*exp(-pow(k/k_silk,1.4)));
+
-    
+    f_baryon = obhh / omhh;
-    f_baryon = obhh/omhh;
+    T_full = f_baryon * T_b + (1 - f_baryon) * T_c;
-    T_full = f_baryon*T_b + (1-f_baryon)*T_c;
+
-    
     /* Now to store these transfer functions */
-    if (tf_baryon!=NULL) *tf_baryon = T_b;
+    if (tf_baryon != NULL)
-    if (tf_cdm!=NULL) *tf_cdm = T_c;
+      *tf_baryon = T_b;
+    if (tf_cdm != NULL)
+      *tf_cdm = T_c;
     return T_full;
   }
-  
+
   double fb_, fc_;
-  
+
   eisenstein_transfer()
-  {  }
-  
-  void set_parameters( const cosmology& cosmo, double Tcmb )
   {
-    m_h0 = cosmo.H0*0.01;
-    TFset_parameters( (cosmo.Omega_m)*cosmo.H0*cosmo.H0*(0.01*0.01),
-                     cosmo.Omega_b/cosmo.Omega_m, Tcmb);
-    
-    fb_ = cosmo.Omega_b/(cosmo.Omega_m);
-    fc_ = (cosmo.Omega_m-cosmo.Omega_b)/(cosmo.Omega_m) ;
   }
-  
+
-  inline double at_k( double k )
+  void set_parameters(const cosmology::parameters &cp, double Tcmb)
+  {
+    m_h0 = cp["H0"] * 0.01;
+    TFset_parameters((cp["Omega_m"]) * cp["H0"] * cp["H0"] * (0.01 * 0.01),
+                     cp["Omega_b"] / cp["Omega_m"], Tcmb);
+
+    fb_ = cp["Omega_b"]  / cp["Omega_m"] ;
+    fc_ = (cp["Omega_m"] - cp["Omega_b"]) / cp["Omega_m"] ;
+  }
+
+  inline double at_k(double k)
   {
     double tfb, tfcdm;
-    TFfit_onek( k*m_h0, &tfb, &tfcdm );
+    TFfit_onek(k * m_h0, &tfb, &tfcdm);
-    return fb_*tfb+fc_*tfcdm;
+    return fb_ * tfb + fc_ * tfcdm;
   }
 };
 
-
+//! Implementation of abstract base class TransferFunction for the Eisenstein & Hu transfer function
-//! Implementation of abstract base class TransferFunction for the Eisenstein & Hu transfer function 
 /*!
  This class implements the analytical fit to the matter transfer
  function by Eisenstein & Hu (1999). In fact it is their code.
@@ -185,196 +189,193 @@ struct eisenstein_transfer
 class transfer_eisenstein_plugin : public transfer_function_plugin
 {
 protected:
-	using transfer_function_plugin::cosmo_;
+  using transfer_function_plugin::cosmo_params_;
-   eisenstein_transfer etf_;
+  eisenstein_transfer etf_;
-	
-public:
-	//! Constructor for Eisenstein & Hu fitting for transfer function
-	/*!
-	 \param aCosm structure of type Cosmology carrying the cosmological parameters
-	 \param Tcmb mean temperature of the CMB fluctuations (defaults to
-	 Tcmb = 2.726 if not specified)
-	 */
-	transfer_eisenstein_plugin( config_file &cf )//Cosmology aCosm, double Tcmb = 2.726 )
-    :  transfer_function_plugin(cf)
-	{
-		double Tcmb = pcf_->get_value_safe("cosmology","Tcmb",2.726);
-		
-        etf_.set_parameters( cosmo_, Tcmb );
-		
-		tf_distinct_ = false;
-		tf_withvel_  = false;
-	}
-	
-	//! Computes the transfer function for k in Mpc/h by calling TFfit_onek
-	inline double compute( double k, tf_type type ){
-        return etf_.at_k( k );
-	}
-	
-	inline double get_kmin( void ){
-		return 1e-5;
-	}
-	
-	inline double get_kmax( void ){
-		return 1.e5;
-	}
-	
-};
 
+public:
+  //! Constructor for Eisenstein & Hu fitting for transfer function
+  /*!
+   \param aCosm structure of type Cosmology carrying the cosmological parameters
+   \param Tcmb mean temperature of the CMB fluctuations (defaults to
+   Tcmb = 2.726 if not specified)
+   */
+  transfer_eisenstein_plugin(config_file &cf, const cosmology::parameters& cp)
+      : transfer_function_plugin(cf, cp)
+  {
+    double Tcmb = pcf_->get_value_safe("cosmology", "Tcmb", 2.726);
+
+    etf_.set_parameters(cp, Tcmb);
+
+    tf_distinct_ = false;
+    tf_withvel_ = false;
+  }
+
+  //! Computes the transfer function for k in Mpc/h by calling TFfit_onek
+  inline double compute(double k, tf_type type)
+  {
+    return etf_.at_k(k);
+  }
+
+  inline double get_kmin(void)
+  {
+    return 1e-5;
+  }
+
+  inline double get_kmax(void)
+  {
+    return 1.e5;
+  }
+};
 
 #include <map>
 class transfer_eisenstein_wdm_plugin : public transfer_function_plugin
 {
 protected:
-	real_t m_WDMalpha, m_h0;
+  real_t m_WDMalpha, m_h0;
-	double omegam_, wdmm_, wdmgx_, wdmnu_, H0_, omegab_;
+  double omegam_, wdmm_, wdmgx_, wdmnu_, H0_, omegab_;
-    std::string type_;
+  std::string type_;
-    std::map< std::string, int > typemap_;
+  std::map<std::string, int> typemap_;
-  
+
-    eisenstein_transfer etf_;
+  eisenstein_transfer etf_;
-  
+
-    enum wdmtyp { wdm_bode, wdm_viel, wdm_bode_wrong=99};
+  enum wdmtyp
-    
+  {
+    wdm_bode,
+    wdm_viel,
+    wdm_bode_wrong = 99
+  };
+
 public:
-	transfer_eisenstein_wdm_plugin( config_file &cf )
+  transfer_eisenstein_wdm_plugin(config_file &cf, const cosmology::parameters& cp)
-	: transfer_function_plugin(cf), m_h0( cosmo_.H0*0.01 )
+      : transfer_function_plugin(cf,cp), m_h0(cp["H0"] * 0.01)
-	{
+  {
-        double Tcmb = pcf_->get_value_safe("cosmology","Tcmb",2.726);
+    double Tcmb = pcf_->get_value_safe("cosmology", "Tcmb", 2.726);
-        etf_.set_parameters( cosmo_, Tcmb );
+    etf_.set_parameters(cp, Tcmb);
-      
+
-        typemap_.insert( std::pair<std::string,int>( "BODE", wdm_bode ) );
+    typemap_.insert(std::pair<std::string, int>("BODE", wdm_bode));
-        typemap_.insert( std::pair<std::string,int>( "VIEL", wdm_viel ) ); // add the other types
+    typemap_.insert(std::pair<std::string, int>("VIEL", wdm_viel));             // add the other types
-        typemap_.insert( std::pair<std::string,int>( "BODE_WRONG", wdm_bode_wrong ) ); // add the other types
+    typemap_.insert(std::pair<std::string, int>("BODE_WRONG", wdm_bode_wrong)); // add the other types
-        
+
-		omegam_ = cf.get_value<double>("cosmology","Omega_m");
+    omegam_ = cp["Omega_m"];//cf.get_value<double>("cosmology", "Omega_m");
-		omegab_ = cf.get_value<double>("cosmology","Omega_b");
+    omegab_ = cp["Omega_b"];//cf.get_value<double>("cosmology", "Omega_b");
-		wdmm_   = cf.get_value<double>("cosmology","WDMmass");
+    wdmm_ = cf.get_value<double>("cosmology", "WDMmass");
-      
+
-		
+    H0_ = cp["H0"];//cf.get_value<double>("cosmology", "H0");
-        H0_     = cf.get_value<double>("cosmology","H0");
+    type_ = cf.get_value_safe<std::string>("cosmology", "WDMtftype", "BODE");
-        type_   = cf.get_value_safe<std::string>("cosmology","WDMtftype","BODE");
+
-        
+    // type_ = std::string( toupper( type_.c_str() ) );
-        //type_ = std::string( toupper( type_.c_str() ) );
+
-        
+    if (typemap_.find(type_) == typemap_.end())
-        if( typemap_.find( type_ ) == typemap_.end() )
+      throw std::runtime_error("unknown transfer function fit for WDM");
-            throw std::runtime_error("unknown transfer function fit for WDM");
+
-        
+    m_WDMalpha = 1.0;
-        m_WDMalpha = 1.0;
+
-        
+    switch (typemap_[type_])
-        switch( typemap_[type_] )
+    {
-        {
+    //... parameterisation from Bode et al. (2001), ApJ, 556, 93
-            //... parameterisation from Bode et al. (2001), ApJ, 556, 93
+    case wdm_bode:
-            case wdm_bode:
+      wdmnu_ = cf.get_value_safe<double>("cosmology", "WDMnu", 1.0);
-                    wdmnu_  = cf.get_value_safe<double>("cosmology","WDMnu",1.0);
+      wdmgx_ = cf.get_value_safe<double>("cosmology", "WDMg_x", 1.5);
-                    wdmgx_  = cf.get_value_safe<double>("cosmology","WDMg_x",1.5);
+      m_WDMalpha = 0.05 * pow(omegam_ / 0.4, 0.15) * pow(H0_ * 0.01 / 0.65, 1.3) * pow(wdmm_, -1.15) * pow(1.5 / wdmgx_, 0.29);
-                    m_WDMalpha = 0.05 * pow( omegam_/0.4,0.15)
+
-                        *pow(H0_*0.01/0.65,1.3)*pow(wdmm_,-1.15)
+      break;
-                        *pow(1.5/wdmgx_,0.29);
+
-                
+    //... parameterisation from Viel et al. (2005), Phys Rev D, 71
-                break;
+    case wdm_viel:
-            
+      wdmnu_ = cf.get_value_safe<double>("cosmology", "WDMnu", 1.12);
-            //... parameterisation from Viel et al. (2005), Phys Rev D, 71
+      m_WDMalpha = 0.049 * pow(omegam_ / 0.25, 0.11) * pow(H0_ * 0.01 / 0.7, 1.22) * pow(wdmm_, -1.11);
-            case wdm_viel:
+      break;
-                    wdmnu_  = cf.get_value_safe<double>("cosmology","WDMnu",1.12);
+
-                    m_WDMalpha = 0.049 * pow( omegam_/0.25,0.11)
+    //.... below is for historical reasons due to the buggy parameterisation
-                        *pow(H0_*0.01/0.7,1.22)*pow(wdmm_,-1.11);
+    //.... in early versions of MUSIC, but apart from H instead of h, Bode et al.
-                break;
+    case wdm_bode_wrong:
-            
+      wdmnu_ = cf.get_value_safe<double>("cosmology", "WDMnu", 1.0);
-            
+      wdmgx_ = cf.get_value_safe<double>("cosmology", "WDMg_x", 1.5);
-            //.... below is for historical reasons due to the buggy parameterisation
+      m_WDMalpha = 0.05 * pow(omegam_ / 0.4, 0.15) * pow(H0_ / 0.65, 1.3) * pow(wdmm_, -1.15) * pow(1.5 / wdmgx_, 0.29);
-            //.... in early versions of MUSIC, but apart from H instead of h, Bode et al.
+      break;
-            case wdm_bode_wrong:
+
-                    wdmnu_  = cf.get_value_safe<double>("cosmology","WDMnu",1.0);
+    default:
-                    wdmgx_  = cf.get_value_safe<double>("cosmology","WDMg_x",1.5);
+      wdmnu_ = cf.get_value_safe<double>("cosmology", "WDMnu", 1.0);
-                    m_WDMalpha = 0.05 * pow( omegam_/0.4,0.15)
+      wdmgx_ = cf.get_value_safe<double>("cosmology", "WDMg_x", 1.5);
-                        *pow(H0_/0.65,1.3)*pow(wdmm_,-1.15)
+      m_WDMalpha = 0.05 * pow(omegam_ / 0.4, 0.15) * pow(H0_ * 0.01 / 0.65, 1.3) * pow(wdmm_, -1.15) * pow(1.5 / wdmgx_, 0.29);
-                        *pow(1.5/wdmgx_,0.29);
+      break;
-                break;
-                
-            default:
-                    wdmnu_  = cf.get_value_safe<double>("cosmology","WDMnu",1.0);
-                    wdmgx_  = cf.get_value_safe<double>("cosmology","WDMg_x",1.5);
-                    m_WDMalpha = 0.05 * pow( omegam_/0.4,0.15)
-                        *pow(H0_*0.01/0.65,1.3)*pow(wdmm_,-1.15)
-                        *pow(1.5/wdmgx_,0.29);
-                break;
-        }
-        std::cerr << "WDM alpha = " << m_WDMalpha << std::endl;
-	}
-	
-	inline double compute( double k, tf_type type )
-	{
-		return etf_.at_k( k )*pow(1.0+pow(m_WDMalpha*k,2.0*wdmnu_),-5.0/wdmnu_);
-	}
-  
-    inline double get_kmin( void ){
-      return 1e-4;
     }
-    
+    std::cerr << "WDM alpha = " << m_WDMalpha << std::endl;
-    inline double get_kmax( void ){
+  }
-      return 1.e4;
+
-    }
+  inline double compute(double k, tf_type type)
-	
+  {
+    return etf_.at_k(k) * pow(1.0 + pow(m_WDMalpha * k, 2.0 * wdmnu_), -5.0 / wdmnu_);
+  }
+
+  inline double get_kmin(void)
+  {
+    return 1e-4;
+  }
+
+  inline double get_kmax(void)
+  {
+    return 1.e4;
+  }
 };
 
 // CDM Bino type WIMP small-scale damped spectrum from Green, Hofmann & Schwarz (2004)
 class transfer_eisenstein_cdmbino_plugin : public transfer_function_plugin
 {
 protected:
-	real_t m_h0;
+  real_t m_h0;
-	double omegam_, H0_, omegab_, mcdm_, Tkd_, kfs_, kd_;
+  double omegam_, H0_, omegab_, mcdm_, Tkd_, kfs_, kd_;
-    eisenstein_transfer etf_;
+  eisenstein_transfer etf_;
-    
+
 public:
-	transfer_eisenstein_cdmbino_plugin( config_file &cf )
+  transfer_eisenstein_cdmbino_plugin(config_file &cf, const cosmology::parameters &cp)
-	: transfer_function_plugin(cf), m_h0( cosmo_.H0*0.01 )
+      : transfer_function_plugin(cf, cp)
-	{
+  {
-        double Tcmb = pcf_->get_value_safe("cosmology","Tcmb",2.726);
+    double Tcmb = pcf_->get_value_safe("cosmology", "Tcmb", 2.726);
-        etf_.set_parameters( cosmo_, Tcmb );
+    etf_.set_parameters(cp, Tcmb);
-        
+
-		omegam_ = cf.get_value<double>("cosmology","Omega_m");
+    omegam_ = cp["Omega_m"];
-		omegab_ = cf.get_value<double>("cosmology","Omega_b");
+    omegab_ = cp["Omega_b"];
-		H0_     = cf.get_value<double>("cosmology","H0");
+    H0_ = cf.get_value<double>("cosmology", "H0");
-        
+    m_h0 = H0_ * 0.01;
-        mcdm_   = cf.get_value_safe<double>("cosmology","CDM_mass", 100.0); // bino particle mass in GeV
-        Tkd_    = cf.get_value_safe<double>("cosmology","CDM_Tkd", 33.0); // temperature at which CDM particle kinetically decouples (in MeV)
-        
-        kfs_    = 1.7e6 / m_h0 * sqrt( mcdm_ / 100. * Tkd_ / 30. ) / (1.0 + log( Tkd_ / 30. ) / 19.2 );
-        kd_     = 3.8e7 / m_h0 * sqrt( mcdm_ / 100. * Tkd_ / 30. );
-        
-        music::ilog.Print(" bino CDM: k_fs = %g, k_d = %g", kfs_, kd_ );
-		
-	}
-	
-	inline double compute( double k, tf_type type )
-	{
-        double kkfs = k/kfs_;
-        double kkfs2 = kkfs*kkfs;
-        double kkd2 = (k/kd_)*(k/kd_);
-        
-        // in principle the Green et al. (2004) works only up to k/k_fs < 1
-        // the fit crosses zero at (k/k_fs)**2 = 3/2, we just zero it there...
-        if( kkfs2 < 1.5 )
-            return etf_.at_k( k ) * (1.0-2.0/3.0*kkfs2) * exp( -kkfs2 - kkd2 );
-        else
-            return 0.0;
-	}
     
-    inline double get_kmin( void ){
+    mcdm_ = cf.get_value_safe<double>("cosmology", "CDM_mass", 100.0); // bino particle mass in GeV
-        return 1e-4;
+    Tkd_ = cf.get_value_safe<double>("cosmology", "CDM_Tkd", 33.0);    // temperature at which CDM particle kinetically decouples (in MeV)
-    }
+
-    
+    kfs_ = 1.7e6 / m_h0 * sqrt(mcdm_ / 100. * Tkd_ / 30.) / (1.0 + log(Tkd_ / 30.) / 19.2);
-    inline double get_kmax( void ){
+    kd_ = 3.8e7 / m_h0 * sqrt(mcdm_ / 100. * Tkd_ / 30.);
-        return 1.e8;
+
-    }
+    music::ilog.Print(" bino CDM: k_fs = %g, k_d = %g", kfs_, kd_);
-	
+  }
+
+  inline double compute(double k, tf_type type)
+  {
+    double kkfs = k / kfs_;
+    double kkfs2 = kkfs * kkfs;
+    double kkd2 = (k / kd_) * (k / kd_);
+
+    // in principle the Green et al. (2004) works only up to k/k_fs < 1
+    // the fit crosses zero at (k/k_fs)**2 = 3/2, we just zero it there...
+    if (kkfs2 < 1.5)
+      return etf_.at_k(k) * (1.0 - 2.0 / 3.0 * kkfs2) * exp(-kkfs2 - kkd2);
+    else
+      return 0.0;
+  }
+
+  inline double get_kmin(void)
+  {
+    return 1e-4;
+  }
+
+  inline double get_kmax(void)
+  {
+    return 1.e8;
+  }
 };
 
-
+namespace
-
+{
-namespace{
+  transfer_function_plugin_creator_concrete<transfer_eisenstein_plugin> creator("eisenstein");
-	transfer_function_plugin_creator_concrete< transfer_eisenstein_plugin > creator("eisenstein");
+  transfer_function_plugin_creator_concrete<transfer_eisenstein_wdm_plugin> creator2("eisenstein_wdm");
-	transfer_function_plugin_creator_concrete< transfer_eisenstein_wdm_plugin > creator2("eisenstein_wdm");
+  transfer_function_plugin_creator_concrete<transfer_eisenstein_cdmbino_plugin> creator3("eisenstein_cdmbino");
-    transfer_function_plugin_creator_concrete< transfer_eisenstein_cdmbino_plugin > creator3("eisenstein_cdmbino");
 }
-

src/plugins/transfer_eisenstein_suppressLSS.cc

@@ -1,243 +0,0 @@
-/*
- 
- transfer_eisenstein.cc - This file is part of MUSIC -
- a code to generate multi-scale initial conditions 
- for cosmological simulations 
- 
- */
-
-#include "transfer_function.hh"
-
-// forward declaration of WDM class
-class transfer_eisenstein_wdm_plugin;
-
-
-struct eisenstein_transfer
-{
-  //Cosmology m_Cosmology;
-  double  m_h0;
-  double	omhh,		/* Omega_matter*h^2 */
-  obhh,		/* Omega_baryon*h^2 */
-  theta_cmb,	/* Tcmb in units of 2.7 K */
-  z_equality,	/* Redshift of matter-radiation equality, really 1+z */
-  k_equality,	/* Scale of equality, in Mpc^-1 */
-  z_drag,		/* Redshift of drag epoch */
-  R_drag,		/* Photon-baryon ratio at drag epoch */
-  R_equality,	/* Photon-baryon ratio at equality epoch */
-  sound_horizon,	/* Sound horizon at drag epoch, in Mpc */
-  k_silk,		/* Silk damping scale, in Mpc^-1 */
-  alpha_c,	/* CDM suppression */
-  beta_c,		/* CDM log shift */
-  alpha_b,	/* Baryon suppression */
-  beta_b,		/* Baryon envelope shift */
-  beta_node,	/* Sound horizon shift */
-  k_peak,		/* Fit to wavenumber of first peak, in Mpc^-1 */
-  sound_horizon_fit,	/* Fit to sound horizon, in Mpc */
-  alpha_gamma;	/* Gamma suppression in approximate TF */
-  
-  //! private member function: sets internal quantities for Eisenstein & Hu fitting
-  void TFset_parameters(double omega0hh, double f_baryon, double Tcmb)
-  /* Set all the scalars quantities for Eisenstein & Hu 1997 fitting formula */
-  /* Input: omega0hh -- The density of CDM and baryons, in units of critical dens,
-   multiplied by the square of the Hubble constant, in units
-   of 100 km/s/Mpc */
-  /* 	  f_baryon -- The fraction of baryons to CDM */
-  /*        Tcmb -- The temperature of the CMB in Kelvin.  Tcmb<=0 forces use
-   of the COBE value of  2.728 K. */
-  /* Output: Nothing, but set many global variables used in TFfit_onek().
-   You can access them yourself, if you want. */
-  /* Note: Units are always Mpc, never h^-1 Mpc. */
-  {
-    double z_drag_b1, z_drag_b2;
-    double alpha_c_a1, alpha_c_a2, beta_c_b1, beta_c_b2, alpha_b_G, y;
-    
-    if (f_baryon<=0.0 || omega0hh<=0.0) {
-      fprintf(stderr, "TFset_parameters(): Illegal input.\n");
-      exit(1);
-    }
-    omhh = omega0hh;
-    obhh = omhh*f_baryon;
-    if (Tcmb<=0.0) Tcmb=2.728;	/* COBE FIRAS */
-    theta_cmb = Tcmb/2.7;
-    
-    z_equality = 2.50e4*omhh/POW4(theta_cmb);  /* Really 1+z */
-    k_equality = 0.0746*omhh/SQR(theta_cmb);
-    
-    z_drag_b1 = 0.313*pow((double)omhh,-0.419)*(1+0.607*pow((double)omhh,0.674));
-    z_drag_b2 = 0.238*pow((double)omhh,0.223);
-    z_drag = 1291*pow(omhh,0.251)/(1+0.659*pow((double)omhh,0.828))*
-    (1+z_drag_b1*pow((double)obhh,(double)z_drag_b2));
-    
-    R_drag = 31.5*obhh/POW4(theta_cmb)*(1000/(1+z_drag));
-    R_equality = 31.5*obhh/POW4(theta_cmb)*(1000/z_equality);
-    
-    sound_horizon = 2./3./k_equality*sqrt(6./R_equality)*
-    log((sqrt(1+R_drag)+sqrt(R_drag+R_equality))/(1+sqrt(R_equality)));
-    
-    k_silk = 1.6*pow((double)obhh,0.52)*pow((double)omhh,0.73)*(1+pow((double)10.4*omhh,-0.95));
-    
-    alpha_c_a1 = pow((double)46.9*omhh,0.670)*(1+pow(32.1*omhh,-0.532));
-    alpha_c_a2 = pow((double)12.0*omhh,0.424)*(1+pow(45.0*omhh,-0.582));
-    alpha_c = pow(alpha_c_a1,-f_baryon)*
-    pow(alpha_c_a2,-CUBE(f_baryon));
-    
-    beta_c_b1 = 0.944/(1+pow(458*omhh,-0.708));
-    beta_c_b2 = pow(0.395*omhh, -0.0266);
-    beta_c = 1.0/(1+beta_c_b1*(pow(1-f_baryon, beta_c_b2)-1));
-    
-    y = z_equality/(1+z_drag);
-    alpha_b_G = y*(-6.*sqrt(1+y)+(2.+3.*y)*log((sqrt(1+y)+1)/(sqrt(1+y)-1)));
-    alpha_b = 2.07*k_equality*sound_horizon*pow(1+R_drag,-0.75)*alpha_b_G;
-    
-    beta_node = 8.41*pow(omhh, 0.435);
-    beta_b = 0.5+f_baryon+(3.-2.*f_baryon)*sqrt(pow(17.2*omhh,2.0)+1);
-    
-    k_peak = 2.5*3.14159*(1+0.217*omhh)/sound_horizon;
-    sound_horizon_fit = 44.5*log(9.83/omhh)/sqrt(1+10.0*pow(obhh,0.75));
-    
-    alpha_gamma = 1-0.328*log(431.0*omhh)*f_baryon + 0.38*log(22.3*omhh)*
-    SQR(f_baryon);
-    
-    return;
-  }
-  
-  //! private member function: computes transfer function for mode k (k in Mpc)
-  inline double TFfit_onek(double k, double *tf_baryon, double *tf_cdm)
-  /* Input: k -- Wavenumber at which to calculate transfer function, in Mpc^-1.
-   *tf_baryon, *tf_cdm -- Input value not used; replaced on output if
-   the input was not NULL. */
-  /* Output: Returns the value of the full transfer function fitting formula.
-   This is the form given in Section 3 of Eisenstein & Hu (1997).
-   *tf_baryon -- The baryonic contribution to the full fit.
-   *tf_cdm -- The CDM contribution to the full fit. */
-  /* Notes: Units are Mpc, not h^-1 Mpc. */
-  {
-    double T_c_ln_beta, T_c_ln_nobeta, T_c_C_alpha, T_c_C_noalpha;
-    double q, xx, xx_tilde;//, q_eff;
-    double T_c_f, T_c, s_tilde, T_b_T0, T_b, f_baryon, T_full;
-    //double T_0_L0, T_0_C0, T_0, gamma_eff;
-    //double T_nowiggles_L0, T_nowiggles_C0, T_nowiggles;
-    
-    k = fabs(k);	/* Just define negative k as positive */
-    if (k==0.0) {
-      if (tf_baryon!=NULL) *tf_baryon = 1.0;
-      if (tf_cdm!=NULL) *tf_cdm = 1.0;
-      return 1.0;
-    }
-    
-    q = k/13.41/k_equality;
-    xx = k*sound_horizon;
-    
-    T_c_ln_beta = log(2.718282+1.8*beta_c*q);
-    T_c_ln_nobeta = log(2.718282+1.8*q);
-    T_c_C_alpha = 14.2/alpha_c + 386.0/(1+69.9*pow(q,1.08));
-    T_c_C_noalpha = 14.2 + 386.0/(1+69.9*pow(q,1.08));
-    
-    T_c_f = 1.0/(1.0+POW4(xx/5.4));
-    T_c = T_c_f*T_c_ln_beta/(T_c_ln_beta+T_c_C_noalpha*SQR(q)) +
-    (1-T_c_f)*T_c_ln_beta/(T_c_ln_beta+T_c_C_alpha*SQR(q));
-    
-    s_tilde = sound_horizon*pow(1.+CUBE(beta_node/xx),-1./3.);
-    xx_tilde = k*s_tilde;
-    
-    T_b_T0 = T_c_ln_nobeta/(T_c_ln_nobeta+T_c_C_noalpha*SQR(q));
-    T_b = sin(xx_tilde)/(xx_tilde)*(T_b_T0/(1.+SQR(xx/5.2))+
-                                    alpha_b/(1.+CUBE(beta_b/xx))*exp(-pow(k/k_silk,1.4)));
-    
-    f_baryon = obhh/omhh;
-    T_full = f_baryon*T_b + (1-f_baryon)*T_c;
-    
-    /* Now to store these transfer functions */
-    if (tf_baryon!=NULL) *tf_baryon = T_b;
-    if (tf_cdm!=NULL) *tf_cdm = T_c;
-    return T_full;
-  }
-  
-  double fb_, fc_;
-  
-  eisenstein_transfer()
-  {  }
-  
-  void set_parameters( const cosmology& cosmo, double Tcmb )
-  {
-    m_h0 = cosmo.H0*0.01;
-    TFset_parameters( (cosmo.Omega_m)*cosmo.H0*cosmo.H0*(0.01*0.01),
-                     cosmo.Omega_b/(cosmo.Omega_m-cosmo.Omega_b), Tcmb);
-    
-    fb_ = cosmo.Omega_b/(cosmo.Omega_m);
-    fc_ = (cosmo.Omega_m-cosmo.Omega_b)/(cosmo.Omega_m) ;
-  }
-  
-  inline double at_k( double k )
-  {
-    double tfb, tfcdm;
-    TFfit_onek( k*m_h0, &tfb, &tfcdm );
-    return fb_*tfb+fc_*tfcdm;
-  }
-};
-
-#include "cosmology.hh"
-
-//! Implementation of abstract base class TransferFunction for the Eisenstein & Hu transfer function with an additional suppression of large-scale power
-/*!
- This class implements the analytical fit to the matter transfer
- function by Eisenstein & Hu (1999). In fact it is their code.
- */
-class transfer_eisensteinS_plugin : public transfer_function_plugin
-{
-protected:
-	using transfer_function_plugin::cosmo_;
-   eisenstein_transfer etf_;
-    double ktrunc_, normfac_;
-    double dplus_;
-	
-public:
-	//! Constructor for Eisenstein & Hu fitting for transfer function
-	/*!
-	 \param aCosm structure of type Cosmology carrying the cosmological parameters
-	 \param Tcmb mean temperature of the CMB fluctuations (defaults to
-	 Tcmb = 2.726 if not specified)
-	 */
-	transfer_eisensteinS_plugin( config_file &cf )//Cosmology aCosm, double Tcmb = 2.726 )
-    :  transfer_function_plugin(cf)
-	{
-		double Tcmb = pcf_->get_value_safe<double>("cosmology","Tcmb",2.726);
-        //double boxlength = pcf_->get_value<double>("setup","boxlength");
-        ktrunc_ = pcf_->get_value<double>("cosmology","ktrunc");
-        normfac_ = 2.0/M_PI;
-		
-        etf_.set_parameters( cosmo_, Tcmb );
-		
-		tf_distinct_ = false;
-		tf_withvel_  = false;
-        tf_withtotal0_ = true;
-        
-        cosmology cosmo( cf );
-        
-        CosmoCalc ccalc(cosmo, this);
-        dplus_ = ccalc.CalcGrowthFactor( cosmo.astart )/ccalc.CalcGrowthFactor( 1.0 );
-	}
-	
-	//! Computes the transfer function for k in Mpc/h by calling TFfit_onek
-	inline double compute( double k, tf_type type ){
-        if( type == total0 )
-            return etf_.at_k( k )/dplus_;
-        return etf_.at_k( k ) * atan(k/ktrunc_)*normfac_;
-	}
-	
-	inline double get_kmin( void ){
-		return 1e-4;
-	}
-	
-	inline double get_kmax( void ){
-		return 1.e4;
-	}
-	
-};
-
-
-
-namespace{
-	transfer_function_plugin_creator_concrete< transfer_eisensteinS_plugin > creator("eisenstein_suppress");
-}
-

src/plugins/transfer_inflation.cc

@@ -1,64 +0,0 @@
-/*
- 
- tranfer_inflation.cc - This file is part of MUSIC -
- a code to generate multi-scale initial conditions 
- for cosmological simulations 
- 
- Copyright (C) 2010  Oliver Hahn
- 
- This program is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
- 
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- GNU General Public License for more details.
- 
- You should have received a copy of the GNU General Public License
- along with this program.  If not, see <http://www.gnu.org/licenses/>.
- 
- */
-
-#include "transfer_function.hh"
-
-class transfer_inflation_plugin : public transfer_function_plugin
-{
-protected:
-	
-	double ns2_;
-	
-public:
-	
-	transfer_inflation_plugin( config_file& cf ) 
-	: transfer_function_plugin( cf )
-	{ 
-		ns2_ = 0.5*cf.get_value<double>("cosmology","nspec");
-		tf_distinct_ = true;
-	}
-	
-	~transfer_inflation_plugin(){ };
-	
-	double compute( double k, tf_type type=baryon)
-	{
-		return pow(k,ns2_);
-	}
-	
-	double get_kmax( void )
-	{
-		return 1e10;
-	}
-	
-	double get_kmin( void )
-	{
-		return 1e-30;
-	}
-	
-};
-
-
-namespace{
-	transfer_function_plugin_creator_concrete< transfer_inflation_plugin > creator("inflation");
-}
-

src/plugins/transfer_linger++.cc

@@ -1,59 +1,63 @@
 /*
- 
+
  transfer_camb.cc - This file is part of MUSIC -
- a code to generate multi-scale initial conditions 
+ a code to generate multi-scale initial conditions
- for cosmological simulations 
+ for cosmological simulations
- 
+
  Copyright (C) 2010  Oliver Hahn
- 
+
  */
 
 #include "transfer_function.hh"
 
 class transfer_LINGERpp_plugin : public transfer_function_plugin
 {
-	
+
 private:
 	std::string m_filename_Pk, m_filename_Tk;
 	std::vector<double> m_tab_k, m_tab_Tk_tot, m_tab_Tk_cdm, m_tab_Tk_baryon, m_tab_Tvk_tot, m_tab_Tvk_cdm, m_tab_Tvk_baryon, m_tab_Tk_tot0;
-	gsl_interp_accel *acc_dtot, *acc_dcdm, *acc_dbaryon, *acc_vtot, *acc_vcdm, *acc_vbaryon, *acc_dtot0;
+	gsl_interp_accel *acc_dtot, *acc_dcdm, *acc_dbaryon, *acc_vtot, *acc_theta_cdm, *acc_theta_baryon, *acc_dtot0;
-	gsl_spline *spline_dtot, *spline_dcdm, *spline_dbaryon, *spline_vtot, *spline_vcdm, *spline_vbaryon, *spline_dtot0;
+	gsl_spline *spline_dtot, *spline_dcdm, *spline_dbaryon, *spline_vtot, *spline_theta_cdm, *spline_theta_baryon, *spline_dtot0;
-	
+
-	bool m_bnovrel;	
+	bool m_bnovrel;
 	bool m_bz0norm;
-	
+
-	void read_table( void ){
+	void read_table(void)
+	{
 #ifdef WITH_MPI
-		if( MPI::COMM_WORLD.Get_rank() == 0 ){
+		if (MPI::COMM_WORLD.Get_rank() == 0)
+		{
 #endif
-			std::cerr 
+			std::cerr
-			<< " - reading tabulated transfer function data from file \n"
+					<< " - reading tabulated transfer function data from file \n"
-			<< "    \'" << m_filename_Tk << "\'\n";
+					<< "    \'" << m_filename_Tk << "\'\n";
-			
+
 			std::string line;
-			std::ifstream ifs( m_filename_Tk.c_str() );
+			std::ifstream ifs(m_filename_Tk.c_str());
-			
+
-			if(! ifs.good() )
+			if (!ifs.good())
-				throw std::runtime_error("Could not find transfer function file \'"+m_filename_Tk+"\'");
+				throw std::runtime_error("Could not find transfer function file \'" + m_filename_Tk + "\'");
-			
+
 			m_tab_k.clear();
 			m_tab_Tk_tot.clear();
 			m_tab_Tk_cdm.clear();
 			m_tab_Tk_baryon.clear();
-            m_tab_Tvk_tot.clear();
+			m_tab_Tvk_tot.clear();
 			m_tab_Tvk_cdm.clear();
 			m_tab_Tvk_baryon.clear();
 			m_tab_Tk_tot0.clear();
-			
+
 			const double zero = 1e-10;
-			
+
-			while( !ifs.eof() ){
+			while (!ifs.eof())
-				getline(ifs,line);
+			{
-				
+				getline(ifs, line);
-				if(ifs.eof()) break;
+
-				
+				if (ifs.eof())
+					break;
+
 				std::stringstream ss(line);
-				
+
 				double k, Tkc, Tkb, Tktot, Tkvc, Tkvb, Tkvtot, Tktot0;
 				ss >> k;
 				ss >> Tktot;
@@ -61,290 +65,288 @@ private:
 				ss >> Tkb;
 				ss >> Tkvc;
 				ss >> Tkvb;
-                ss >> Tkvtot;
+				ss >> Tkvtot;
-                ss >> Tktot0;
+				ss >> Tktot0;
 
-		if( m_bnovrel )
+				if (m_bnovrel)
-		{
+				{
-			std::cerr << " - transfer_linger++ : disabling baryon-DM relative velocity\n";
+					std::cerr << " - transfer_linger++ : disabling baryon-DM relative velocity\n";
-			Tkvb = Tkvc;
+					Tkvb = Tkvc;
-		}		
+				}
-				Tktot = std::max(zero,Tktot);
+				Tktot = std::max(zero, Tktot);
-				Tkc   = std::max(zero,Tkc);
+				Tkc = std::max(zero, Tkc);
-				Tkb   = std::max(zero,Tkb);
+				Tkb = std::max(zero, Tkb);
-                Tkvtot= std::max(zero,Tkvtot);
+				Tkvtot = std::max(zero, Tkvtot);
-				Tkvc  = std::max(zero,Tkvc);
+				Tkvc = std::max(zero, Tkvc);
-				Tkvb  = std::max(zero,Tkvb);
+				Tkvb = std::max(zero, Tkvb);
-                Tktot0= std::max(zero,Tktot0);
+				Tktot0 = std::max(zero, Tktot0);
-								
+
-				m_tab_k.push_back( log10(k) );
+				m_tab_k.push_back(log10(k));
-				
+
-				m_tab_Tk_tot.push_back( log10(Tktot) );
+				m_tab_Tk_tot.push_back(log10(Tktot));
-				m_tab_Tk_baryon.push_back( log10(Tkb) );
+				m_tab_Tk_baryon.push_back(log10(Tkb));
-				m_tab_Tk_cdm.push_back( log10(Tkc) );
+				m_tab_Tk_cdm.push_back(log10(Tkc));
-                m_tab_Tvk_tot.push_back( log10(Tkvtot) );
+				m_tab_Tvk_tot.push_back(log10(Tkvtot));
-                m_tab_Tvk_cdm.push_back( log10(Tkvc) );
+				m_tab_Tvk_cdm.push_back(log10(Tkvc));
-				m_tab_Tvk_baryon.push_back( log10(Tkvb) );
+				m_tab_Tvk_baryon.push_back(log10(Tkvb));
-                m_tab_Tk_tot0.push_back( log10(Tktot0) );
+				m_tab_Tk_tot0.push_back(log10(Tktot0));
-				
 			}
-			
+
-			ifs.close();			
+			ifs.close();
 #ifdef WITH_MPI
 		}
-		
+
-		unsigned n=m_tab_k.size();
+		unsigned n = m_tab_k.size();
-		MPI::COMM_WORLD.Bcast( &n, 1, MPI_UNSIGNED, 0 );
+		MPI::COMM_WORLD.Bcast(&n, 1, MPI_UNSIGNED, 0);
-		
+
-		if( MPI::COMM_WORLD.Get_rank() > 0 ){
+		if (MPI::COMM_WORLD.Get_rank() > 0)
-			m_tab_k.assign(n,0);
+		{
-			m_tab_Tk_tot.assign(n,0);
+			m_tab_k.assign(n, 0);
-			m_tab_Tk_cdm.assign(n,0);
+			m_tab_Tk_tot.assign(n, 0);
-			m_tab_Tk_baryon.assign(n,0);
+			m_tab_Tk_cdm.assign(n, 0);
-			m_tab_Tvk_tot.assign(n,0);
+			m_tab_Tk_baryon.assign(n, 0);
-			m_tab_Tvk_cdm.assign(n,0);
+			m_tab_Tvk_tot.assign(n, 0);
-			m_tab_Tvk_baryon.assign(n,0);
+			m_tab_Tvk_cdm.assign(n, 0);
-            m_tab_Tk_tot0.assign(n,0);
+			m_tab_Tvk_baryon.assign(n, 0);
+			m_tab_Tk_tot0.assign(n, 0);
 		}
-		
+
-		MPI::COMM_WORLD.Bcast( &m_tab_k[0],  n, MPI_DOUBLE, 0 );
+		MPI::COMM_WORLD.Bcast(&m_tab_k[0], n, MPI_DOUBLE, 0);
-		MPI::COMM_WORLD.Bcast( &m_tab_Tk_tot[0], n, MPI_DOUBLE, 0 );
+		MPI::COMM_WORLD.Bcast(&m_tab_Tk_tot[0], n, MPI_DOUBLE, 0);
-		MPI::COMM_WORLD.Bcast( &m_tab_Tk_cdm[0], n, MPI_DOUBLE, 0 );
+		MPI::COMM_WORLD.Bcast(&m_tab_Tk_cdm[0], n, MPI_DOUBLE, 0);
-		MPI::COMM_WORLD.Bcast( &m_tab_Tk_baryon[0], n, MPI_DOUBLE, 0 );
+		MPI::COMM_WORLD.Bcast(&m_tab_Tk_baryon[0], n, MPI_DOUBLE, 0);
-        MPI::COMM_WORLD.Bcast( &m_tab_Tvk_tot[0], n, MPI_DOUBLE, 0 );
+		MPI::COMM_WORLD.Bcast(&m_tab_Tvk_tot[0], n, MPI_DOUBLE, 0);
-		MPI::COMM_WORLD.Bcast( &m_tab_Tvk_cdm[0], n, MPI_DOUBLE, 0 );
+		MPI::COMM_WORLD.Bcast(&m_tab_Tvk_cdm[0], n, MPI_DOUBLE, 0);
-		MPI::COMM_WORLD.Bcast( &m_tab_Tvk_baryon[0], n, MPI_DOUBLE, 0 );
+		MPI::COMM_WORLD.Bcast(&m_tab_Tvk_baryon[0], n, MPI_DOUBLE, 0);
-        MPI::COMM_WORLD.Bcast( &m_tab_Tk_tot0[0], n, MPI_DOUBLE, 0 );
+		MPI::COMM_WORLD.Bcast(&m_tab_Tk_tot0[0], n, MPI_DOUBLE, 0);
-		
+
 #endif
-		
 	}
-	
+
 public:
-	transfer_LINGERpp_plugin( config_file& cf )
+	transfer_LINGERpp_plugin(config_file &cf, const cosmology::parameters &cp)
-	: transfer_function_plugin( cf )
+			: transfer_function_plugin(cf, cp)
 	{
-		m_filename_Tk	= pcf_->get_value<std::string>("cosmology","transfer_file");
+		m_filename_Tk = pcf_->get_value<std::string>("cosmology", "transfer_file");
-		
+
 		//.. disable the baryon-CDM relative velocity (both follow the total matter potential)
-		m_bnovrel		= pcf_->get_value_safe<bool>("cosmology","no_vrel",false);
+		m_bnovrel = pcf_->get_value_safe<bool>("cosmology", "no_vrel", false);
-		
+
 		//.. normalize at z=0 rather than using the linearly scaled zini spectrum
 		//.. this can be different due to radiation still being non-negligible at
 		//.. high redshifts
-		m_bz0norm		= pcf_->get_value_safe<bool>("cosmology","z0norm",true);
+		m_bz0norm = pcf_->get_value_safe<bool>("cosmology", "z0norm", true);
-		
+
-		tf_distinct_   = true;
+		tf_distinct_ = true;
-		tf_withvel_    = true;
+		tf_withvel_ = true;
-		tf_velunits_   = true;
+		tf_velunits_ = true;
 
 		//.. normalize with z=0 spectrum rather than zini spectrum?
-		if( m_bz0norm )
+		if (m_bz0norm)
 			tf_withtotal0_ = true;
 		else
 			tf_withtotal0_ = false;
-		
+
-		
+		read_table();
-		read_table( );
+
-		
 		acc_dtot = gsl_interp_accel_alloc();
 		acc_dcdm = gsl_interp_accel_alloc();
 		acc_dbaryon = gsl_interp_accel_alloc();
-        acc_vtot = gsl_interp_accel_alloc();
+		acc_vtot = gsl_interp_accel_alloc();
-		acc_vcdm = gsl_interp_accel_alloc();
+		acc_theta_cdm = gsl_interp_accel_alloc();
-		acc_vbaryon = gsl_interp_accel_alloc();
+		acc_theta_baryon = gsl_interp_accel_alloc();
 		acc_dtot0 = gsl_interp_accel_alloc();
-		
+
-		spline_dtot = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
+		spline_dtot = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
-		spline_dcdm = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
+		spline_dcdm = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
-		spline_dbaryon = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
+		spline_dbaryon = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
-		spline_vtot = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
+		spline_vtot = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
-		spline_vcdm = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
+		spline_theta_cdm = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
-		spline_vbaryon = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
+		spline_theta_baryon = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
-		spline_dtot0 = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
+		spline_dtot0 = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
-		
+
-		gsl_spline_init (spline_dtot, &m_tab_k[0], &m_tab_Tk_tot[0], m_tab_k.size() );
+		gsl_spline_init(spline_dtot, &m_tab_k[0], &m_tab_Tk_tot[0], m_tab_k.size());
-		gsl_spline_init (spline_dcdm, &m_tab_k[0], &m_tab_Tk_cdm[0], m_tab_k.size() );
+		gsl_spline_init(spline_dcdm, &m_tab_k[0], &m_tab_Tk_cdm[0], m_tab_k.size());
-		gsl_spline_init (spline_dbaryon, &m_tab_k[0], &m_tab_Tk_baryon[0], m_tab_k.size() );
+		gsl_spline_init(spline_dbaryon, &m_tab_k[0], &m_tab_Tk_baryon[0], m_tab_k.size());
-		gsl_spline_init (spline_vtot, &m_tab_k[0], &m_tab_Tvk_tot[0], m_tab_k.size() );
+		gsl_spline_init(spline_vtot, &m_tab_k[0], &m_tab_Tvk_tot[0], m_tab_k.size());
-		gsl_spline_init (spline_vcdm, &m_tab_k[0], &m_tab_Tvk_cdm[0], m_tab_k.size() );
+		gsl_spline_init(spline_theta_cdm, &m_tab_k[0], &m_tab_Tvk_cdm[0], m_tab_k.size());
-		gsl_spline_init (spline_vbaryon, &m_tab_k[0], &m_tab_Tvk_baryon[0], m_tab_k.size() );
+		gsl_spline_init(spline_theta_baryon, &m_tab_k[0], &m_tab_Tvk_baryon[0], m_tab_k.size());
-		
+
-		if( tf_withtotal0_ )
+		if (tf_withtotal0_)
-			gsl_spline_init (spline_dtot0, &m_tab_k[0], &m_tab_Tk_tot0[0], m_tab_k.size() );
+			gsl_spline_init(spline_dtot0, &m_tab_k[0], &m_tab_Tk_tot0[0], m_tab_k.size());
 		else
-			gsl_spline_init (spline_dtot0, &m_tab_k[0], &m_tab_Tk_tot[0], m_tab_k.size() );
+			gsl_spline_init(spline_dtot0, &m_tab_k[0], &m_tab_Tk_tot[0], m_tab_k.size());
 	}
-	
+
 	~transfer_LINGERpp_plugin()
 	{
-		gsl_spline_free (spline_dtot);
+		gsl_spline_free(spline_dtot);
-		gsl_spline_free (spline_dcdm);
+		gsl_spline_free(spline_dcdm);
-		gsl_spline_free (spline_dbaryon);
+		gsl_spline_free(spline_dbaryon);
-        gsl_spline_free (spline_vtot);
+		gsl_spline_free(spline_vtot);
-		gsl_spline_free (spline_vcdm);
+		gsl_spline_free(spline_theta_cdm);
-		gsl_spline_free (spline_vbaryon);
+		gsl_spline_free(spline_theta_baryon);
-		gsl_spline_free (spline_dtot0);
+		gsl_spline_free(spline_dtot0);
-		
+
-		gsl_interp_accel_free (acc_dtot);
+		gsl_interp_accel_free(acc_dtot);
-		gsl_interp_accel_free (acc_dcdm);
+		gsl_interp_accel_free(acc_dcdm);
-		gsl_interp_accel_free (acc_dbaryon);
+		gsl_interp_accel_free(acc_dbaryon);
-		gsl_interp_accel_free (acc_vtot);
+		gsl_interp_accel_free(acc_vtot);
-		gsl_interp_accel_free (acc_vcdm);
+		gsl_interp_accel_free(acc_theta_cdm);
-		gsl_interp_accel_free (acc_vbaryon);
+		gsl_interp_accel_free(acc_theta_baryon);
-        gsl_interp_accel_free (acc_dtot0);
+		gsl_interp_accel_free(acc_dtot0);
 	}
-	
+
-  inline double extrap_left( double k, const tf_type& type ) 
+	inline double extrap_left(double k, const tf_type &type)
-  {
+	{
-    if( k<1e-8 )
+		if (k < 1e-8)
-      return 1.0;
+			return 1.0;
-    
+
-    double v1(1.0), v2(1.0);
+		double v1(1.0), v2(1.0);
-    switch( type )
+		switch (type)
-      {
+		{
-      case cdm:
+		case delta_cdm:
-	v1 = m_tab_Tk_cdm[0];
+			v1 = m_tab_Tk_cdm[0];
-	v2 = m_tab_Tk_cdm[1];
+			v2 = m_tab_Tk_cdm[1];
-	break;
+			break;
-      case baryon:
+		case delta_baryon:
-	v1 = m_tab_Tk_baryon[0];
+			v1 = m_tab_Tk_baryon[0];
-	v2 = m_tab_Tk_baryon[1];
+			v2 = m_tab_Tk_baryon[1];
-	break;
+			break;
-      case vtotal:
+		case theta_matter:
-	v1 = m_tab_Tvk_tot[0];
+			v1 = m_tab_Tvk_tot[0];
-	v2 = m_tab_Tvk_tot[1];
+			v2 = m_tab_Tvk_tot[1];
-	break;
+			break;
-      case vcdm:
+		case theta_cdm:
-	v1 = m_tab_Tvk_cdm[0];
+			v1 = m_tab_Tvk_cdm[0];
-	v2 = m_tab_Tvk_cdm[1];
+			v2 = m_tab_Tvk_cdm[1];
-	break;
+			break;
-      case vbaryon:
+		case theta_baryon:
-	v1 = m_tab_Tvk_baryon[0];
+			v1 = m_tab_Tvk_baryon[0];
-	v2 = m_tab_Tvk_baryon[1];
+			v2 = m_tab_Tvk_baryon[1];
-	break;
+			break;
-      case total: 
+		case delta_matter:
-	v1 = m_tab_Tk_tot[0];
+			v1 = m_tab_Tk_tot[0];
-	v2 = m_tab_Tk_tot[1];
+			v2 = m_tab_Tk_tot[1];
-	break;
+			break;
-      case total0:
+		case delta_matter0:
-	v1 = m_tab_Tk_tot0[0];
+			v1 = m_tab_Tk_tot0[0];
-	v2 = m_tab_Tk_tot0[1];
+			v2 = m_tab_Tk_tot0[1];
-	break;
+			break;
-	
+
-      default:
+		default:
-	throw std::runtime_error("Invalid type requested in transfer function evaluation");
+			throw std::runtime_error("Invalid type requested in transfer function evaluation");
-      }
+		}
-    
+
-    double lk = log10(k);
+		double lk = log10(k);
-    double dk = m_tab_k[1]-m_tab_k[0];
+		double dk = m_tab_k[1] - m_tab_k[0];
-    double delk = lk-m_tab_k[0];
+		double delk = lk - m_tab_k[0];
-    
+
-    //double xi = (v2-v1)/dk;
+		// double xi = (v2-v1)/dk;
-    return pow(10.0,(v2-v1)/dk*(delk)+v1);
+		return pow(10.0, (v2 - v1) / dk * (delk) + v1);
-  }
-  
-  inline double extrap_right( double k, const tf_type& type ) 
-  {
-    double v1(1.0), v2(1.0);
-    
-    int n=m_tab_k.size()-1, n1=n-1;
-    switch( type )
-      {
-      case cdm:
-	v1 = m_tab_Tk_cdm[n1];
-	v2 = m_tab_Tk_cdm[n];
-	break;
-      case baryon:
-	v1 = m_tab_Tk_baryon[n1];
-	v2 = m_tab_Tk_baryon[n];
-	break;
-      case vtotal:
-	v1 = m_tab_Tvk_tot[n1];
-	v2 = m_tab_Tvk_tot[n];
-	break;
-      case vcdm:
-	v1 = m_tab_Tvk_cdm[n1];
-	v2 = m_tab_Tvk_cdm[n];
-	break;
-      case vbaryon:
-	v1 = m_tab_Tvk_baryon[n1];
-	v2 = m_tab_Tvk_baryon[n];
-	break;
-      case total: 
-	v1 = m_tab_Tk_tot[n1];
-	v2 = m_tab_Tk_tot[n];
-	break;
-      case total0:
-	v1 = m_tab_Tk_tot0[n1];
-	v2 = m_tab_Tk_tot0[n];
-	break;
-	
-      default:
-	throw std::runtime_error("Invalid type requested in transfer function evaluation");
-      }
-    
-    double lk = log10(k);
-    double dk = m_tab_k[n]-m_tab_k[n1];
-    double delk = lk-m_tab_k[n];
-    
-    //double xi = (v2-v1)/dk;
-    return pow(10.0,(v2-v1)/dk*(delk)+v2);
-  }
-  
-  inline double compute( double k, tf_type type ){
-    
-    double lk = log10(k);
-    
-    //if( lk<m_tab_k[1])
-    //	return 1.0;
-    
-    //if( lk>m_tab_k[m_tab_k.size()-2] );
-    //	return m_tab_Tk_cdm[m_tab_k.size()-2]/k/k;
-    
-    if( k<get_kmin() )
-      return extrap_left(k, type );
-    
-    if( k>get_kmax() )
-      return extrap_right(k,type );
-    
-    
-    switch( type )
-      {
-      case cdm:
-	return pow(10.0, gsl_spline_eval (spline_dcdm, lk, acc_dcdm) );
-      case baryon:
-	return pow(10.0, gsl_spline_eval (spline_dbaryon, lk, acc_dbaryon) );
-      case vtotal:
-	return pow(10.0, gsl_spline_eval (spline_vtot, lk, acc_vtot) );
-      case vcdm:
-	return pow(10.0, gsl_spline_eval (spline_vcdm, lk, acc_vcdm) );
-      case vbaryon:
-	return pow(10.0, gsl_spline_eval (spline_vbaryon, lk, acc_vbaryon) );
-      case total: 
-	return pow(10.0, gsl_spline_eval (spline_dtot, lk, acc_dtot) );
-      case total0:
-	return pow(10.0, gsl_spline_eval (spline_dtot0, lk, acc_dtot0) );
-        
-      default:
-	throw std::runtime_error("Invalid type requested in transfer function evaluation");
-      }
-    
-    return 1.0;
-  }
-	
-	inline double get_kmin( void ){
-		return pow(10.0,m_tab_k[0]);
 	}
-	
+
-	inline double get_kmax( void ){
+	inline double extrap_right(double k, const tf_type &type)
-		return pow(10.0,m_tab_k[m_tab_k.size()-1]);
+	{
+		double v1(1.0), v2(1.0);
+
+		int n = m_tab_k.size() - 1, n1 = n - 1;
+		switch (type)
+		{
+		case delta_cdm:
+			v1 = m_tab_Tk_cdm[n1];
+			v2 = m_tab_Tk_cdm[n];
+			break;
+		case delta_baryon:
+			v1 = m_tab_Tk_baryon[n1];
+			v2 = m_tab_Tk_baryon[n];
+			break;
+		case theta_matter:
+			v1 = m_tab_Tvk_tot[n1];
+			v2 = m_tab_Tvk_tot[n];
+			break;
+		case theta_cdm:
+			v1 = m_tab_Tvk_cdm[n1];
+			v2 = m_tab_Tvk_cdm[n];
+			break;
+		case theta_baryon:
+			v1 = m_tab_Tvk_baryon[n1];
+			v2 = m_tab_Tvk_baryon[n];
+			break;
+		case delta_matter:
+			v1 = m_tab_Tk_tot[n1];
+			v2 = m_tab_Tk_tot[n];
+			break;
+		case delta_matter0:
+			v1 = m_tab_Tk_tot0[n1];
+			v2 = m_tab_Tk_tot0[n];
+			break;
+
+		default:
+			throw std::runtime_error("Invalid type requested in transfer function evaluation");
+		}
+
+		double lk = log10(k);
+		double dk = m_tab_k[n] - m_tab_k[n1];
+		double delk = lk - m_tab_k[n];
+
+		// double xi = (v2-v1)/dk;
+		return pow(10.0, (v2 - v1) / dk * (delk) + v2);
+	}
+
+	inline double compute(double k, tf_type type)
+	{
+
+		double lk = log10(k);
+
+		// if( lk<m_tab_k[1])
+		//	return 1.0;
+
+		// if( lk>m_tab_k[m_tab_k.size()-2] );
+		//	return m_tab_Tk_cdm[m_tab_k.size()-2]/k/k;
+
+		if (k < get_kmin())
+			return extrap_left(k, type);
+
+		if (k > get_kmax())
+			return extrap_right(k, type);
+
+		switch (type)
+		{
+		case delta_cdm:
+			return pow(10.0, gsl_spline_eval(spline_dcdm, lk, acc_dcdm));
+		case delta_baryon:
+			return pow(10.0, gsl_spline_eval(spline_dbaryon, lk, acc_dbaryon));
+		case theta_matter:
+			return pow(10.0, gsl_spline_eval(spline_vtot, lk, acc_vtot));
+		case theta_cdm:
+			return pow(10.0, gsl_spline_eval(spline_theta_cdm, lk, acc_theta_cdm));
+		case theta_baryon:
+			return pow(10.0, gsl_spline_eval(spline_theta_baryon, lk, acc_theta_baryon));
+		case delta_matter:
+			return pow(10.0, gsl_spline_eval(spline_dtot, lk, acc_dtot));
+		case delta_matter0:
+			return pow(10.0, gsl_spline_eval(spline_dtot0, lk, acc_dtot0));
+
+		default:
+			throw std::runtime_error("Invalid type requested in transfer function evaluation");
+		}
+
+		return 1.0;
+	}
+
+	inline double get_kmin(void)
+	{
+		return pow(10.0, m_tab_k[0]);
+	}
+
+	inline double get_kmax(void)
+	{
+		return pow(10.0, m_tab_k[m_tab_k.size() - 1]);
 	}
-	
 };
 
-namespace{
+namespace
-	transfer_function_plugin_creator_concrete< transfer_LINGERpp_plugin > creator("linger++");
+{
+	transfer_function_plugin_creator_concrete<transfer_LINGERpp_plugin> creator("linger++");
 }
-
-

src/plugins/transfer_music.cc

@@ -16,8 +16,8 @@ class transfer_MUSIC_plugin : public transfer_function_plugin
 private:
 	std::string m_filename_Pk, m_filename_Tk;
 	std::vector<double> m_tab_k, m_tab_Tk_tot, m_tab_Tk_cdm, m_tab_Tk_baryon, m_tab_Tvk_cdm, m_tab_Tvk_baryon;
-	gsl_interp_accel *acc_dtot, *acc_dcdm, *acc_dbaryon, *acc_vcdm, *acc_vbaryon;
+	gsl_interp_accel *acc_dtot, *acc_dcdm, *acc_dbaryon, *acc_theta_cdm, *acc_theta_baryon;
-	gsl_spline *spline_dtot, *spline_dcdm, *spline_dbaryon, *spline_vcdm, *spline_vbaryon;
+	gsl_spline *spline_dtot, *spline_dcdm, *spline_dbaryon, *spline_theta_cdm, *spline_theta_baryon;
 	
 	
 	
@@ -117,8 +117,8 @@ private:
 	}
 	
 public:
-	transfer_MUSIC_plugin( config_file& cf )
+	transfer_MUSIC_plugin( config_file& cf, const cosmology::parameters& cp )
-	: transfer_function_plugin( cf )
+	: transfer_function_plugin( cf, cp )
 	{
 		m_filename_Tk = pcf_->get_value<std::string>("cosmology","transfer_file");
 		
@@ -127,20 +127,20 @@ public:
 		acc_dtot = gsl_interp_accel_alloc();
 		acc_dcdm = gsl_interp_accel_alloc();
 		acc_dbaryon = gsl_interp_accel_alloc();
-		acc_vcdm = gsl_interp_accel_alloc();
+		acc_theta_cdm = gsl_interp_accel_alloc();
-		acc_vbaryon = gsl_interp_accel_alloc();
+		acc_theta_baryon = gsl_interp_accel_alloc();
 		
 		spline_dtot = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
 		spline_dcdm = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
 		spline_dbaryon = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
-		spline_vcdm = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
+		spline_theta_cdm = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
-		spline_vbaryon = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
+		spline_theta_baryon = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
 		
 		gsl_spline_init (spline_dtot, &m_tab_k[0], &m_tab_Tk_tot[0], m_tab_k.size() );
 		gsl_spline_init (spline_dcdm, &m_tab_k[0], &m_tab_Tk_cdm[0], m_tab_k.size() );
 		gsl_spline_init (spline_dbaryon, &m_tab_k[0], &m_tab_Tk_baryon[0], m_tab_k.size() );
-		gsl_spline_init (spline_vcdm, &m_tab_k[0], &m_tab_Tvk_cdm[0], m_tab_k.size() );
+		gsl_spline_init (spline_theta_cdm, &m_tab_k[0], &m_tab_Tvk_cdm[0], m_tab_k.size() );
-		gsl_spline_init (spline_vbaryon, &m_tab_k[0], &m_tab_Tvk_baryon[0], m_tab_k.size() );
+		gsl_spline_init (spline_theta_baryon, &m_tab_k[0], &m_tab_Tvk_baryon[0], m_tab_k.size() );
 		
 		tf_distinct_ = true;
 		tf_withvel_  = true;
@@ -151,14 +151,14 @@ public:
 		gsl_spline_free (spline_dtot);
 		gsl_spline_free (spline_dcdm);
 		gsl_spline_free (spline_dbaryon);
-		gsl_spline_free (spline_vcdm);
+		gsl_spline_free (spline_theta_cdm);
-		gsl_spline_free (spline_vbaryon);
+		gsl_spline_free (spline_theta_baryon);
 		
 		gsl_interp_accel_free (acc_dtot);
 		gsl_interp_accel_free (acc_dcdm);
 		gsl_interp_accel_free (acc_dbaryon);
-		gsl_interp_accel_free (acc_vcdm);
+		gsl_interp_accel_free (acc_theta_cdm);
-		gsl_interp_accel_free (acc_vbaryon);
+		gsl_interp_accel_free (acc_theta_baryon);
 	}
 	
 	inline double extrap_left( double k, const tf_type& type ) 
@@ -169,23 +169,23 @@ public:
 		double v1(1.0), v2(1.0);
 		switch( type )
 		{
-			case cdm:
+			case delta_cdm:
 				v1 = m_tab_Tk_cdm[0];
 				v2 = m_tab_Tk_cdm[1];
 				break;
-			case baryon:
+			case delta_baryon:
 				v1 = m_tab_Tk_baryon[0];
 				v2 = m_tab_Tk_baryon[1];
 				break;
-			case vcdm:
+			case theta_cdm:
 				v1 = m_tab_Tvk_cdm[0];
 				v2 = m_tab_Tvk_cdm[1];
 				break;
-			case vbaryon:
+			case theta_baryon:
 				v1 = m_tab_Tvk_baryon[0];
 				v2 = m_tab_Tvk_baryon[1];
 				break;
-			case total: 
+			case delta_matter: 
 				v1 = m_tab_Tk_tot[0];
 				v2 = m_tab_Tk_tot[1];
 				break;
@@ -209,23 +209,23 @@ public:
 		int n=m_tab_k.size()-1, n1=n-1;
 		switch( type )
 		{
-			case cdm:
+			case delta_cdm:
 				v1 = m_tab_Tk_cdm[n1];
 				v2 = m_tab_Tk_cdm[n];
 				break;
-			case baryon:
+			case delta_baryon:
 				v1 = m_tab_Tk_baryon[n1];
 				v2 = m_tab_Tk_baryon[n];
 				break;
-			case vcdm:
+			case theta_cdm:
 				v1 = m_tab_Tvk_cdm[n1];
 				v2 = m_tab_Tvk_cdm[n];
 				break;
-			case vbaryon:
+			case theta_baryon:
 				v1 = m_tab_Tvk_baryon[n1];
 				v2 = m_tab_Tvk_baryon[n];
 				break;
-			case total: 
+			case delta_matter: 
 				v1 = m_tab_Tk_tot[n1];
 				v2 = m_tab_Tk_tot[n];
 				break;
@@ -261,15 +261,15 @@ public:
 		
 		switch( type )
 		{
-			case cdm:
+			case delta_cdm:
 				return pow(10.0, gsl_spline_eval (spline_dcdm, lk, acc_dcdm) );
-			case baryon:
+			case delta_baryon:
 				return pow(10.0, gsl_spline_eval (spline_dbaryon, lk, acc_dbaryon) );
-			case vcdm:
+			case theta_cdm:
-				return pow(10.0, gsl_spline_eval (spline_vcdm, lk, acc_vcdm) );
+				return pow(10.0, gsl_spline_eval (spline_theta_cdm, lk, acc_theta_cdm) );
-			case vbaryon:
+			case theta_baryon:
-				return pow(10.0, gsl_spline_eval (spline_vbaryon, lk, acc_vbaryon) );
+				return pow(10.0, gsl_spline_eval (spline_theta_baryon, lk, acc_theta_baryon) );
-			case total: 
+			case delta_matter: 
 				return pow(10.0, gsl_spline_eval (spline_dtot, lk, acc_dtot) );
 				
 			default:

src/random.hh

@@ -95,7 +95,7 @@ protected:
 
 public:
 	//! constructor
-	random_number_generator(config_file &cf, transfer_function *ptf = NULL)
+	explicit random_number_generator(config_file &cf)
 			: pcf_(&cf) //, prefh_( &refh )
 	{
 		levelmin_ = pcf_->get_value<int>("setup", "levelmin");

src/transfer_function.cc

@@ -30,17 +30,14 @@ void print_transfer_function_plugins()
 		if( (*it).second )
 			std::cout << "\t\'" << (*it).first << "\'\n";
 		++it;
-	}
+	}	
-	
-	
 }
 
-transfer_function_plugin *select_transfer_function_plugin( config_file& cf )
+std::unique_ptr<transfer_function_plugin> select_transfer_function_plugin( config_file& cf, const cosmology::parameters& cp )
 {
 	std::string tfname = cf.get_value<std::string>( "cosmology", "transfer" );
 	
-	transfer_function_plugin_creator *the_transfer_function_plugin_creator 
+	transfer_function_plugin_creator *the_transfer_function_plugin_creator = get_transfer_function_plugin_map()[ tfname ];
-	= get_transfer_function_plugin_map()[ tfname ];
 	
 	if( !the_transfer_function_plugin_creator )
 	{	
@@ -55,9 +52,6 @@ transfer_function_plugin *select_transfer_function_plugin( config_file& cf )
 		music::ulog.Print("Selecting transfer function plug-in  : %s",tfname.c_str() );
 	}
 	
-	transfer_function_plugin *the_transfer_function_plugin 
+	return the_transfer_function_plugin_creator->create( cf, cp );
-	= the_transfer_function_plugin_creator->create( cf );
-	
-	return the_transfer_function_plugin;
 }
 

src/transfer_function.hh

@@ -1,15 +1,13 @@
 /*
- 
- transfer_function.hh - This file is part of MUSIC -
- a code to generate multi-scale initial conditions 
- for cosmological simulations 
- 
- Copyright (C) 2010  Oliver Hahn
- 
-*/
 
-#ifndef __TRANSFERFUNCTION_HH
+ transfer_function.hh - This file is part of MUSIC -
-#define __TRANSFERFUNCTION_HH
+ a code to generate multi-scale initial conditions
+ for cosmological simulations
+
+ Copyright (C) 2010  Oliver Hahn
+
+*/
+#pragma once
 
 #include <vector>
 #include <sstream>
@@ -23,12 +21,26 @@
 #include <gsl/gsl_spline.h>
 #include <gsl/gsl_sf_gamma.h>
 
-#include "Numerics.hh"
+#include <Numerics.hh>
-#include "config_file.hh"
+#include <config_file.hh>
+#include <cosmology_parameters.hh>
 
-
+enum tf_type
-enum tf_type{
+{
-	total, cdm, baryon, vtotal, vcdm, vbaryon, total0
+	delta_matter,
+	delta_cdm,
+	delta_baryon,
+	theta_matter,
+	theta_cdm,
+	theta_baryon,
+	delta_bc,
+	theta_bc,
+	delta_matter0,
+	delta_cdm0,
+	delta_baryon0,
+	theta_matter0,
+	theta_cdm0,
+	theta_baryon0,
 };
 
 #define GSL_INTEGRATION_ERR 1e-5
@@ -37,97 +49,98 @@ enum tf_type{
 /*!
  This class implements a purely virtual interface that can be
  used to derive instances implementing various transfer functions.
- */ 
+ */
-class transfer_function_plugin{
+class transfer_function_plugin
+{
 public:
-	Cosmology cosmo_;		//!< cosmological parameter, read from config_file
+	config_file *pcf_;													//!< pointer to config_file from which to read parameters
-	config_file *pcf_;		//!< pointer to config_file from which to read parameters
+	const cosmology::parameters &cosmo_params_; //!< cosmological parameters are stored here
-	bool tf_distinct_;		//!< bool if density transfer function is distinct for baryons and DM
+	bool tf_distinct_;													//!< bool if density transfer function is distinct for baryons and DM
-	bool tf_withvel_;		//!< bool if also have velocity transfer functions
+	bool tf_withvel_;														//!< bool if also have velocity transfer functions
-	bool tf_withtotal0_;	//!< have the z=0 spectrum for normalisation purposes
+	bool tf_withtotal0_;												//!< have the z=0 spectrum for normalisation purposes
-	bool tf_velunits_;		//!< velocities are in velocity units (km/s)
+	bool tf_velunits_;													//!< velocities are in velocity units (km/s)
-public:
+	bool tf_isnormalised_;											//!< assume that transfer functions come already correctly normalised and need be re-normalised to a specified value
-	
+
 	//! constructor
-	transfer_function_plugin( config_file& cf ) 
+	transfer_function_plugin(config_file &cf,const cosmology::parameters &cosmo_params)
-	: pcf_( &cf ), tf_distinct_(false), tf_withvel_(false), tf_withtotal0_(false), tf_velunits_(false)
+			: pcf_(&cf), cosmo_params_(cosmo_params), tf_distinct_(false), tf_withvel_(false), tf_withtotal0_(false), tf_velunits_(false), tf_isnormalised_(false)
 	{
-		real_t zstart;
-		zstart				= pcf_->get_value<real_t>( "setup", "zstart" );
-		cosmo_.astart		= 1.0/(1.0+zstart);
-		cosmo_.Omega_b		= pcf_->get_value<real_t>( "cosmology", "Omega_b" );
-		cosmo_.Omega_m		= pcf_->get_value<real_t>( "cosmology", "Omega_m" );
-		cosmo_.Omega_DE		= pcf_->get_value<real_t>( "cosmology", "Omega_L" );
-		cosmo_.H0			= pcf_->get_value<real_t>( "cosmology", "H0" );
-		cosmo_.sigma8		= pcf_->get_value<real_t>( "cosmology", "sigma_8" );
-		cosmo_.nspect		= pcf_->get_value<real_t>( "cosmology", "nspec" );
 	}
-	
+
 	//! destructor
-	virtual ~transfer_function_plugin(){ };
+	virtual ~transfer_function_plugin(){};
-	
+
+	//! initialise, i.e. prepare data for later usage
+	virtual void intialise(void) {}
+
 	//! compute value of transfer function at waven umber
-	virtual double compute( double k, tf_type type) = 0;
+	virtual double compute(double k, tf_type type) = 0;
 
 	//! return maximum wave number allowed
-	virtual double get_kmax( void ) = 0;
+	virtual double get_kmax(void) = 0;
-	
-	//! return minimum wave number allowed
-	virtual double get_kmin( void ) = 0;
-	
-	//! return if density transfer function is distinct for baryons and DM
-	bool tf_is_distinct( void )
-	{	return tf_distinct_;	}
-	
-	//! return if we also have velocity transfer functions
-	bool tf_has_velocities( void )
-	{	return tf_withvel_; }
-    
-    //! return if we also have a z=0 transfer function for normalisation
-    bool tf_has_total0( void )
-    {   return tf_withtotal0_; }
-	
-	//! return if velocity returned is in velocity or in displacement units
-	bool tf_velocity_units( void )
-	{	return tf_velunits_; }
-};
 
+	//! return minimum wave number allowed
+	virtual double get_kmin(void) = 0;
+
+	//! return if density transfer function is distinct for baryons and DM
+	bool tf_is_distinct(void) const
+	{
+		return tf_distinct_;
+	}
+
+	//! return if we also have velocity transfer functions
+	bool tf_has_velocities(void) const
+	{
+		return tf_withvel_;
+	}
+
+	//! return if we also have a z=0 transfer function for normalisation
+	bool tf_has_total0(void) const
+	{
+		return tf_withtotal0_;
+	}
+
+	//! return if velocity returned is in velocity or in displacement units
+	bool tf_velocity_units(void) const
+	{
+		return tf_velunits_;
+	}
+};
 
 //! Implements abstract factory design pattern for transfer function plug-ins
 struct transfer_function_plugin_creator
 {
 	//! create an instance of a transfer function plug-in
-	virtual transfer_function_plugin * create( config_file& cf ) const = 0;
+	virtual std::unique_ptr<transfer_function_plugin> create(config_file &cf, const cosmology::parameters& cp) const = 0;
-	
+
-	//! destroy an instance of a plug-in 
+	//! destroy an instance of a plug-in
-	virtual ~transfer_function_plugin_creator() { }
+	virtual ~transfer_function_plugin_creator() {}
 };
 
 //! Write names of registered transfer function plug-ins to stdout
-std::map< std::string, transfer_function_plugin_creator *>& get_transfer_function_plugin_map();
+std::map<std::string, transfer_function_plugin_creator *> &get_transfer_function_plugin_map();
-void print_transfer_function_plugins( void );
+void print_transfer_function_plugins(void);
 
 //! Concrete factory pattern for transfer function plug-ins
-template< class Derived >
+template <class Derived>
 struct transfer_function_plugin_creator_concrete : public transfer_function_plugin_creator
 {
-	//! register the plug-in by its name 
+	//! register the plug-in by its name
-	transfer_function_plugin_creator_concrete( const std::string& plugin_name )
+	transfer_function_plugin_creator_concrete(const std::string &plugin_name)
 	{
-		get_transfer_function_plugin_map()[ plugin_name ] = this;
+		get_transfer_function_plugin_map()[plugin_name] = this;
 	}
-	
+
-	//! create an instance of the plug-in 
+	//! create an instance of the plug-in
-	transfer_function_plugin * create( config_file& cf ) const
+	std::unique_ptr<transfer_function_plugin> create(config_file &cf, const cosmology::parameters& cp) const
 	{
-		return new Derived( cf );
+		return std::make_unique<Derived>(cf,cp);
 	}
 };
 
 typedef transfer_function_plugin transfer_function;
 
-transfer_function_plugin *select_transfer_function_plugin( config_file& cf );
+std::unique_ptr<transfer_function_plugin> select_transfer_function_plugin(config_file &cf, const cosmology::parameters &cp);
-
 
 /**********************************************************************/
 /**********************************************************************/
@@ -141,561 +154,78 @@ public:
 	static real_t nspec_;
 	double pnorm_, sqrtpnorm_;
 	static tf_type type_;
-	
+
-	TransferFunction_k( tf_type type, transfer_function *tf, real_t nspec, real_t pnorm )
+	TransferFunction_k(tf_type type, transfer_function *tf, real_t nspec, real_t pnorm)
-	: pnorm_(pnorm)
+			: pnorm_(pnorm)
 	{
 		ptf_ = tf;
 		nspec_ = nspec;
-		sqrtpnorm_ = sqrt( pnorm_ );
+		sqrtpnorm_ = sqrt(pnorm_);
 		type_ = type;
 
 		std::string fname("input_powerspec.txt");
-		if( type == cdm || type == total )
+		if (type == delta_cdm || type == delta_matter)
-		  {
-		    std::ofstream ofs(fname.c_str());
-		    double kmin=log10(tf->get_kmin()), kmax= log10(tf->get_kmax());
-		    double dk=(kmax-kmin)/300.;
-		    
-		    ofs << "# The power spectrum definition is smaller than CAMB by a factor 8 pi^3."
-		        << std::endl;
-
-		    if( tf->tf_is_distinct() )
-		      {
-			ofs << "#"
-			    << std::setw(15) << "k [h/Mpc]"
-			    << std::setw(16) << "P_cdm"
-			    << std::setw(16) << "P_vcdm"
-			    << std::setw(16) << "P_bar"
-			    << std::setw(16) << "P_vbar"
-			    << std::setw(16) << "P_total"
-                << std::setw(16) << "P_vtotal"
-                << std::endl;
-			
-			for( int i=0; i<300; ++i )
-			{ 
-				double k = pow(10.0,kmin+i*dk);
-				ofs << std::setw(16) << k 
-				    << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,cdm),2)
-				    << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,vcdm),2)
-				    << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,baryon),2)
-				    << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,vbaryon),2)
-				    << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,total),2)
-                    << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,vtotal),2)
-				    << std::endl;
-			}
-		      }
-		    else
-		      {
-                  ofs << "#"
-                  << std::setw(16) << "k [h/Mpc]"
-                  << std::setw(16) << "P_cdm"
-                  << std::setw(16) << "P_vcdm"
-                  << std::setw(16) << "P_total"
-                  << std::endl;
-                  
-			for( int i=0; i<300; ++i )
-			{
-                double k = pow(10.0,kmin+i*dk);
-				ofs << std::setw(16) << k 
-				    << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,cdm),2)
-				    << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,vcdm),2)
-				    << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,total),2)
-				    << std::endl;
-			}
-		      }
-		  }
-
-
-	}
-	
-	inline real_t compute( real_t k ) const
-	{
-		return sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,type_);
-	}
-};
-
-
-/**********************************************************************/
-/**********************************************************************/
-/**********************************************************************/
-
-#define NZERO_Q
-typedef std::complex<double> complex;
-class TransferFunction_real
-{
-	
-public:
-	
-	double Tr0_;
-	real_t Tmin_, Tmax_, Tscale_;
-	real_t rneg_, rneg2_, kny_;
-	static transfer_function *ptf_;
-	static real_t nspec_;
-	
-protected:
-	
-	real_t krgood( real_t mu, real_t q, real_t dlnr, real_t kr )
-	{
-		double krnew = kr;
-		complex zm, zp;
-		double arg, iarg, xneg, xpos, y;
-		gsl_sf_result g_a, g_p;
-		
-		xpos = 0.5*(mu+1.0+q);
-		xneg = 0.5*(mu+1.0-q);
-		y = M_PI/(2.0*dlnr);
-		zp=complex(xpos,y);
-		zm=complex(xneg,y);
-		
-		gsl_sf_lngamma_complex_e (zp.real(), zp.imag(), &g_a, &g_p);
-		zp=std::polar(exp(g_a.val),g_p.val);
-		real_t zpa = g_p.val;
-
-		gsl_sf_lngamma_complex_e (zm.real(), zm.imag(), &g_a, &g_p);
-		zm=std::polar(exp(g_a.val),g_p.val);
-		real_t zma = g_p.val;
-		
-		arg=log(2.0/kr)/dlnr+(zpa+zma)/M_PI;
-		iarg=(real_t)((int)(arg + 0.5));
-		
-		if( arg!=iarg )
-			krnew=kr*exp((arg-iarg)*dlnr);
-		
-		return krnew;
-	}
-	
-	void transform( real_t pnorm, unsigned N, real_t q, std::vector<double>& rr, std::vector<double>& TT )
-	{
-		const double mu = 0.5;
-		double qmin = 1.0e-7, qmax = 1.0e+7;
-		
-		q = 0.0;
-		
-		//N = 16384;
-		N = 1<<12;
-		
-#ifdef NZERO_Q
-		q=0.4;
-		//q=-0.1;
-#endif
-		
-		double kmin = qmin, kmax=qmax;
-		double rmin = qmin, rmax = qmax;
-		double k0 = exp(0.5*(log(kmax)+log(kmin)));
-		double r0 = exp(0.5*(log(rmax)+log(rmin)));
-		double L = log(rmax)-log(rmin);
-		double k0r0 = k0*r0;
-		double dlnk = L/N, dlnr = L/N;
-		
-		double sqrtpnorm = sqrt(pnorm);
-		
-		double dir = 1.0;
-		
-		double fftnorm = 1.0/N;
-		
-		complex_t *in, *out;
-		
-		in = new complex_t[N];
-		out = new complex_t[N];
-		
-		//... perform anti-ringing correction from Hamilton (2000)
-		k0r0 = krgood( mu, q, dlnr, k0r0 );
-		
-		std::string ofname;
-		switch( type_ )
 		{
-			case cdm:		
-				ofname = "input_powerspec_cdm.txt"; break;
-			case baryon:	
-				ofname = "input_powerspec_baryon.txt"; break;
-			case total:		
-				ofname = "input_powerspec_total.txt"; break;
-			case vcdm:		
-				ofname = "input_powerspec_vcdm.txt"; break;
-			case vbaryon:	
-				ofname = "input_powerspec_vbaryon.txt"; break;
-			default:
-				throw std::runtime_error("Unknown transfer function type in TransferFunction_real::transform");
-		}
-			
-		
-		std::ofstream ofsk(ofname.c_str());
-		
-		ofsk << "# The power spectrum definition is smaller than CAMB by a factor 8 pi^3."
-		    << std::endl;
-
-		for( unsigned i=0; i<N; ++i )
-		{
-			double k = k0*exp(((int)i - (int)N/2+1) * dlnk);
-			double T = ptf_->compute( k, type_ );
-			double del = sqrtpnorm*T*pow(k,0.5*nspec_);
-			
-			RE(in[i]) = del*pow(k,1.5-q);
-			IM(in[i]) = 0.0;
-			
-			ofsk << std::setw(16) << k <<std::setw(16) << del*del << std::setw(16) << T << std::endl;
-		}
-		ofsk.close();
-		
-		fftw_plan_t p,ip;
-		p = FFTW_API(plan_dft_1d)(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
-		ip = FFTW_API(plan_dft_1d)(N, out, in, FFTW_BACKWARD, FFTW_ESTIMATE);
-		FFTW_API(execute)(p);
-		
-		//... compute the Hankel transform by convolution with the Bessel function
-		for( unsigned i=0; i<N; ++i )
-		{
-			int ii=i;
-			if( ii > (int)N/2 )
-				ii -= N;
-			
-#ifndef NZERO_Q
-			double y=ii*M_PI/L;
-			complex zp((mu+1.0)*0.5,y);
-			gsl_sf_result g_a, g_p;
-			gsl_sf_lngamma_complex_e(zp.real(), zp.imag(), &g_a, &g_p);
-			
-			double arg = 2.0*(log(2.0/k0r0)*y+g_p.val);
-			//complex cu = complex(out[i].re,out[i].im)*std::polar(1.0,arg);
-			//out[i].re = cu.real()*fftnorm;
-			//out[i].im = cu.imag()*fftnorm;
-			
-			complex cu = complex( RE(out[i]), IM(out[i]) ) * std::polar(1.0,arg);
-			RE(out[i]) = cu.real()*fftnorm;
-			IM(out[i]) = cu.imag()*fftnorm;
-			
-			
-#else		
-			complex x(dir*q, (double)ii*2.0*M_PI/L);
-			gsl_sf_result g_a, g_p;
-			
-			complex g1, g2, garg, U, phase;						
-			complex twotox = pow(complex(2.0,0.0),x);
-			
-			/////////////////////////////////////////////////////////
-			//.. evaluate complex Gamma functions
-			
-			garg = 0.5*(mu+1.0+x);
-			gsl_sf_lngamma_complex_e (garg.real(), garg.imag(), &g_a, &g_p);
-			g1 = std::polar(exp(g_a.val),g_p.val);
-
-			
-			garg = 0.5*(mu+1.0-x);
-			gsl_sf_lngamma_complex_e (garg.real(), garg.imag(), &g_a, &g_p);
-			g2 = std::polar(exp(g_a.val),g_p.val);
-
-			/////////////////////////////////////////////////////////
-			//.. compute U
-			
-			if( (fabs(g2.real()) < 1e-19 && fabs(g2.imag()) < 1e-19) )
-			{
-				//std::cerr << "Warning : encountered possible singularity in TransferFunction_real::transform!\n";
-				g1 = 1.0; g2 = 1.0;
-			}
-			
-			
-			U = twotox * g1 / g2;
-			phase = pow(complex(k0r0,0.0),complex(0.0,2.0*M_PI*(double)ii/L));
-			
-			complex cu = complex(RE(out[i]),IM(out[i]))*U*phase*fftnorm;
-			
-			RE(out[i]) = cu.real();
-			IM(out[i]) = cu.imag();
-
-			/*if( (RE(out[i]) != RE(out[i]))||(IM(out[i]) != IM(out[i])) )
-			{	std::cerr << "NaN @ i=" << i << ", U= " << U << ", phase = " << phase << ", g1 = " << g1 << ", g2 = " << g2 << std::endl;
-				std::cerr << "mu+1+q = " << mu+1.0+q << std::endl;
-				//break;
-			}*/
-			
-#endif
-
-		}
-
-		FFTW_API(execute)(ip);
-		
-		rr.assign(N,0.0);
-		TT.assign(N,0.0);
-		
-		r0 = k0r0/k0;
-		
-		for( unsigned i=0; i<N; ++i )
-		{
-			int ii = i;
-			ii -= N/2-1;
-			double r = r0*exp(-ii*dlnr);
-			
-			rr[N-i-1] = r;
-			TT[N-i-1] = 4.0*M_PI* sqrt(M_PI/2.0) *  RE(in[i]) * pow(r,-(1.5+q));
-		}
-		
-		
-		
-		{
-			std::string fname;
-			if(type_==total) fname = "transfer_real_total.txt";
-			if(type_==cdm) fname = "transfer_real_cdm.txt";
-			if(type_==baryon) fname = "transfer_real_baryon.txt";
-			if(type_==vcdm) fname = "transfer_real_vcdm.txt";
-			if(type_==vbaryon) fname = "transfer_real_vbaryon.txt";
-			
 			std::ofstream ofs(fname.c_str());
-							  
+			double kmin = log10(tf->get_kmin()), kmax = log10(tf->get_kmax());
-			for( unsigned i=0; i<N; ++i )
+			double dk = (kmax - kmin) / 300.;
-			{
-				int ii = i;
-				ii -= N/2-1;
-				
-				double r = r0*exp(-ii*dlnr);//r0*exp(ii*dlnr);
 
-				double T = 4.0*M_PI* sqrt(M_PI/2.0) *  RE(in[i]) * pow(r,-(1.5+q));
+			ofs << "# The power spectrum definition is smaller than CAMB by a factor 8 pi^3."
-				ofs << r << "\t\t" << T << "\t\t" << IM(in[i]) << std::endl;				
+					<< std::endl;
+
+			if (tf->tf_is_distinct())
+			{
+				ofs << "#"
+						<< std::setw(15) << "k [h/Mpc]"
+						<< std::setw(16) << "P_cdm"
+						<< std::setw(16) << "P_theta_cdm"
+						<< std::setw(16) << "P_bar"
+						<< std::setw(16) << "P_vbar"
+						<< std::setw(16) << "P_total"
+						<< std::setw(16) << "P_vtotal"
+						<< std::endl;
+
+				for (int i = 0; i < 300; ++i)
+				{
+					double k = pow(10.0, kmin + i * dk);
+					ofs << std::setw(16) << k
+							<< std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, delta_cdm), 2)
+							<< std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, theta_cdm), 2)
+							<< std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, delta_baryon), 2)
+							<< std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, theta_baryon), 2)
+							<< std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, delta_matter), 2)
+							<< std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, theta_matter), 2)
+							<< std::endl;
+				}
+			}
+			else
+			{
+				ofs << "#"
+						<< std::setw(16) << "k [h/Mpc]"
+						<< std::setw(16) << "P_cdm"
+						<< std::setw(16) << "P_theta_cdm"
+						<< std::setw(16) << "P_total"
+						<< std::endl;
+
+				for (int i = 0; i < 300; ++i)
+				{
+					double k = pow(10.0, kmin + i * dk);
+					ofs << std::setw(16) << k
+							<< std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, delta_cdm), 2)
+							<< std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, theta_cdm), 2)
+							<< std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, delta_matter), 2)
+							<< std::endl;
+				}
 			}
 		}
-		
-		delete[] in;
-		delete[] out;
-
-		FFTW_API(destroy_plan)(p);
-		FFTW_API(destroy_plan)(ip);
 	}
-	std::vector<real_t> m_xtable,m_ytable,m_dytable;
+
-	double m_xmin, m_xmax, m_dx, m_rdx;
+	inline real_t compute(real_t k) const
-	static tf_type type_;
-	
-	
-public:
-	
-	TransferFunction_real( double boxlength, int nfull, tf_type type, transfer_function *tf, 
-						   real_t nspec, real_t pnorm, real_t rmin, real_t rmax, real_t knymax, unsigned nr )
 	{
-		real_t q = 0.8;
+		return sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, type_);
-		
-		ptf_	= tf;
-		nspec_	= nspec;
-		type_	= type;
-		kny_	= knymax;
-		
-		
-		/*****************************************************************/
-		//... compute the FFTlog transform of the k^n T(k) kernel
-
-		std::vector<double> r,T;
-		transform( pnorm, nr, q, r, T );
-		
-		gsl_set_error_handler_off ();
-		
-		
-		/*****************************************************************/
-		//... compute T(r=0) by 3D k-space integration
-		{
-			const double REL_PRECISION=1.e-5;
-			
-			gsl_integration_workspace * ws = gsl_integration_workspace_alloc (1000);
-			gsl_function ff;
-			
-			double kmin = 2.0*M_PI/boxlength;
-			double kmax = nfull*M_PI/boxlength;
-
-			//... integrate 0..kmax
-			double a[6];
-			a[3] = 0.1*kmin;
-			a[4] = kmax;
-			  
-			ff.function = &call_x;
-			ff.params = reinterpret_cast<void*> (a);
-			double res, err, res2, err2;
-			
-			gsl_integration_qags( &ff, a[3], a[4], 0.0, GSL_INTEGRATION_ERR, 1000, ws, &res, &err );
-			
-			if( err/res > REL_PRECISION )
-				std::cerr << " - Warning: no convergence in \'TransferFunction_real\', rel. error=" << err/res << std::endl;
-			
-			//... integrate 0..kmin
-			a[3] = 0.1*kmin;
-			a[4] = kmin;
-			gsl_integration_qags( &ff, a[3], a[4], 0.0, GSL_INTEGRATION_ERR, 1000, ws, &res2, &err2 );
-			
-			if( err2/res2 > 10*REL_PRECISION )
-				std::cerr << " - Warning: no convergence in \'TransferFunction_real\', rel. error=" << err2/res2 << std::endl;
-
-			gsl_integration_workspace_free ( ws );
-
-			//.. get kmin..kmax
-			res -= res2;
-			//.. *8 because we only integrated one octant
-			res *= 8.0*sqrt(pnorm);
-			Tr0_ = res;
-		}
-		
-		/*****************************************************************/
-		//... store as table for spline interpolation
-		
-		gsl_interp_accel *accp;
-		gsl_spline *splinep;
-		
-		std::vector<double> xsp, ysp;
-		
-		for( unsigned i=0; i<r.size(); ++i )
-		{
-			if( r[i] > rmin/512. && r[i] < rmax )
-			{
-				xsp.push_back( log10(r[i]) );
-				ysp.push_back( T[i]*r[i]*r[i] );
-			}
-			
-		}
-		
-		accp = gsl_interp_accel_alloc ();
-
-		//... spline interpolation is only marginally slower here
-		splinep = gsl_spline_alloc (gsl_interp_akima, xsp.size() );
-
-		//... set up everything for spline interpolation
-		gsl_spline_init (splinep, &xsp[0], &ysp[0], xsp.size() );
-		
-		//.. build lookup table using spline interpolation
-		m_xmin = log10(rmin);
-		m_xmax = log10(rmax);
-		m_dx   = (m_xmax-m_xmin)/nr;
-		m_rdx  = 1.0/m_dx;
-		
-		for(unsigned i=0; i<nr; ++i )
-		{
-			m_xtable.push_back( m_xmin+i*m_dx );
-			m_ytable.push_back( gsl_spline_eval(splinep, (m_xtable.back()), accp) );
-		}
-		
-		for(unsigned i=0; i<nr-1; ++i )
-		{
-			real_t dy,dr;
-			dy = m_ytable[i+1]/pow(m_xtable[i+1],2)-m_ytable[i]/pow(m_xtable[i],2);
-			dr = pow(10.0,m_xtable[i+1])-pow(10.0,m_xtable[i]);
-			m_dytable.push_back(dy/dr);
-		}
-		
-		gsl_spline_free (splinep);
-		gsl_interp_accel_free (accp);
-	}
-	
-	
-	static double call_wrapper( double k, void *arg )
-	{
-		double *a = (double*)arg;
-
-		double T = ptf_->compute( k, type_ );
-		
-		return 4.0*M_PI*a[0]*T*pow(k,0.5*nspec_)*k*k;
-	}
-	
-	
-	
-	static double call_x( double kx, void *arg )
-	{
-		gsl_integration_workspace * wx = gsl_integration_workspace_alloc (1000);
-		
-		double *a = (double*)arg;
-		double kmin = a[3], kmax = a[4];
-		
-		a[0] = kx;
-		
-		gsl_function FX;
-		FX.function = &call_y;
-		FX.params = reinterpret_cast<void*> (a);
-		
-		double resx, errx;
-		gsl_integration_qags( &FX, kmin, kmax, 0.0, GSL_INTEGRATION_ERR, 1000, wx, &resx, &errx );
-							 
-		gsl_integration_workspace_free (wx);
-		
-		return resx;
-	}
-	
-	static double call_y( double ky, void *arg )
-	{
-		gsl_integration_workspace * wy = gsl_integration_workspace_alloc (1000);
-		
-		double *a = (double*)arg;
-		double kmin = a[3], kmax = a[4];
-		
-		a[1] = ky;
-		
-		gsl_function FY;
-		FY.function = &call_z;
-		FY.params = reinterpret_cast<void*> (a);
-		
-		double resy, erry;
-		gsl_integration_qags( &FY, kmin, kmax, 0.0, GSL_INTEGRATION_ERR, 1000, wy, &resy, &erry );
-		
-		gsl_integration_workspace_free (wy);
-		
-		return resy;
-	}
-	
-	static double call_z( double kz, void *arg )
-	{
-		double *a = (double*)arg;
-		double kx = a[0], ky = a[1];
-		
-		double kk = sqrt(kx*kx+ky*ky+kz*kz);
-		double T = ptf_->compute( kk, type_ );
-		
-		return pow(kk,0.5*nspec_)*T;
-		
-	}
-	
-	~TransferFunction_real()
-	{ }
-
-	inline real_t get_grad( real_t r2 ) const
-	{
-		double r = 0.5*fast_log10(r2);
-		double ii = (r-m_xmin)*m_rdx;
-		int i = (int)ii;
-		
-		i=std::max(0,i);
-		i=std::min(i, (int)m_xtable.size()-2);
-		
-		return (real_t)m_dytable[i];
-	}
-	
-	//... fast version
-	inline real_t compute_real( real_t r2 ) const
-	{
-		const double EPS = 1e-8;
-		const double Reps2 = EPS*EPS;
-		
-		if( r2 <Reps2 )
-			return Tr0_;
-		
-		//double r = 0.5*log10(r2);
-		double r = 0.5*fast_log10(r2);
-		
-		double ii = (r-m_xmin)*m_rdx;
-		int i = (int)ii;
-		
-		i=std::max(0,i);
-		i=std::min(i, (int)m_xtable.size()-2);
-		
-		double y1,y2;
-		y1 = m_ytable[i];
-		y2 = m_ytable[i+1];
-		
-		//divide by r**2 because r^2 T is tabulated
-		//return (real_t)((y1 + (y2-y1)*(ii-(double)i))/r2);
-        
-        real_t retval = (real_t)((y1 + (y2-y1)*(ii-(double)i))/r2);
-        
-        if( retval != retval ){
-            std::cerr << "FAILURE FAILURE FAILURE" << std::endl;
-            fprintf(stderr,"r2 = %f, r = %f, i = %d, y1 = %f, y2 = %f, xtable[i]=%f",r2,r,i,y1,y2, m_xtable[i]);
-            abort();
-        }
-        
-        return retval;
 	}
 };
 
-
+/**********************************************************************/
-#endif
+/**********************************************************************/
+/**********************************************************************/