Data Structures
class	HaloModel
	Class for computing the standard halomodel predictions for \(P(k,z)\) plus \(n(M,z)\) or just doing spherical collapse. More...

Typedefs
using	Timer = FML::UTILS::Timings

using	DVector = FML::SOLVERS::ODESOLVER::DVector

using	DVector2D = FML::SOLVERS::ODESOLVER::DVector2D

using	ODEFunction = FML::SOLVERS::ODESOLVER::ODEFunction

using	ODEFunctionJacobian = FML::SOLVERS::ODESOLVER::ODEFunctionJacobian

using	ODESolver = FML::SOLVERS::ODESOLVER::ODESolver

using	Spline = FML::INTERPOLATION::SPLINE::Spline

using	Spline2D = FML::INTERPOLATION::SPLINE::Spline2D

using	Func1 = std::function< double(double)>

using	Func2 = std::function< double(double, double)>

using	SphericalCollapse = FML::COSMOLOGY::SPHERICALCOLLAPSE::SphericalCollapse

using	SphericalCollapseModel = FML::COSMOLOGY::SPHERICALCOLLAPSE::SphericalCollapseModel

using	RhoFunction = std::function< double(double, double, HaloModel *hm)>

using	cofMFunction = std::function< double(double, HaloModel *hm)>

using	HaloMultiplicityFunction = std::function< double(double, HaloModel *hm)>

using	HaloBiasFunction = std::function< double(double, HaloModel *hm)>

using	WindowFunction = std::function< double(double)>

Functions
void	powerlaw_extend_logY_of_logx_spline (Spline &logY_of_logx_spline, double xmin_new, double xmax_new, int npts_new, std::string label="logY(log(x)) extrapolated")

void	read_pofk_file (const std::string filename_pofk, Spline &logDelta_of_logk_spline, double normalization, int col_k, int col_pofk)

std::pair< double, double >	compute_neff (const Spline &logsigma_of_logR_spline, const Spline &lognu_of_logR_spline)

double	compute_sigmav (const Func1 logDelta_of_logk, const WindowFunction window_of_kR, double R, const DVector *fiducial_logk_array)

void	compute_sigma (const Func1 logDelta_of_logk, const WindowFunction window_of_kR, Spline &logsigma_of_logR_spline, const DVector fiducial_logR_array, const DVector fiducial_logk_array)

double	RvirofM (double M, double DeltaVir, double OmegaM, double x)

double	RofM (double M, double OmegaM, double x)

double	MofRvir (double Rvir, double DeltaVir, double OmegaM, double x)

double	MofR (double R, double OmegaM, double x)

DVector	compute_rho_fourier_single (const Func2 &rho_function, const DVector &k_array, const double M, const double Rvir)

void	compute_rho_fourier (const Func2 &rho_function, Spline2D &y_of_logR_and_logk_spline, const double DeltaVir, const double OmegaM, const double x, const DVector fiducial_logR_array, const DVector fiducial_logk_array)

void	compute_lognu_of_logR (const double deltac, const Spline &logsigma_of_logR_spline, Spline &lognu_of_logR_spline, Spline &logR_of_lognu_spline, const DVector *fiducial_logR_array)

void	compute_massfunction (const Spline &logR_of_lognu_spline, const Spline &lognu_of_logR_spline, Spline &dndlogM_of_logM_spline, Spline &n_of_logM_spline, const Func1 &pdf_of_nu, const double OmegaM, const double x)

void	compute_one_and_two_halo_terms (const Spline &logR_of_lognu_spline, const Spline2D &y_of_logR_and_logk_spline, Spline &logDelta_onehalo_of_logk_spline, Spline &logDelta_twohalo_of_logk_spline, Spline &logDelta_full_of_logk_spline, const Func1 &logDeltaLin_of_logk, const Func1 &pdf_of_nu, const Func1 &bias_of_nu, const double eta_hmcode, const double kmin, const double kmax, const int npts_k)

void	compute_formation_redshift_factor (const Spline &growthfactor_of_x_spline, const Spline &logsigma_of_logR_spline, Spline &formationredshift_of_logM_spline, double deltac, double OmegaM, double xcollapse, double massfraction)

HaloMultiplicityFunction	get_sheth_tormen_halo_multiplicity_function ()

HaloBiasFunction	get_sheth_tormen_halo_bias_function ()

WindowFunction	get_tophat_window_fourier ()

cofMFunction	get_nfw_halo_mass_concentration_relation ()

RhoFunction	get_nfw_halo_density_profile ()

Variables
int	_ACCURACY_BOOST = 1

Typedef Documentation

◆ cofMFunction

using FML::COSMOLOGY::HALOMODEL::cofMFunction = typedef std::function<double(double,HaloModel *hm)>

Definition at line 38 of file Halomodel.h.

◆ DVector

using FML::COSMOLOGY::HALOMODEL::DVector = typedef FML::SOLVERS::ODESOLVER::DVector

Definition at line 24 of file Halomodel.h.

◆ DVector2D

using FML::COSMOLOGY::HALOMODEL::DVector2D = typedef FML::SOLVERS::ODESOLVER::DVector2D

Definition at line 25 of file Halomodel.h.

◆ Func1

using FML::COSMOLOGY::HALOMODEL::Func1 = typedef std::function<double(double)>

Definition at line 31 of file Halomodel.h.

◆ Func2

using FML::COSMOLOGY::HALOMODEL::Func2 = typedef std::function<double(double,double)>

Definition at line 32 of file Halomodel.h.

◆ HaloBiasFunction

using FML::COSMOLOGY::HALOMODEL::HaloBiasFunction = typedef std::function<double(double,HaloModel *hm)>

Definition at line 40 of file Halomodel.h.

◆ HaloMultiplicityFunction

using FML::COSMOLOGY::HALOMODEL::HaloMultiplicityFunction = typedef std::function<double(double,HaloModel *hm)>

Definition at line 39 of file Halomodel.h.

◆ ODEFunction

using FML::COSMOLOGY::HALOMODEL::ODEFunction = typedef FML::SOLVERS::ODESOLVER::ODEFunction

Definition at line 26 of file Halomodel.h.

◆ ODEFunctionJacobian

using FML::COSMOLOGY::HALOMODEL::ODEFunctionJacobian = typedef FML::SOLVERS::ODESOLVER::ODEFunctionJacobian

Definition at line 27 of file Halomodel.h.

◆ ODESolver

using FML::COSMOLOGY::HALOMODEL::ODESolver = typedef FML::SOLVERS::ODESOLVER::ODESolver

Definition at line 28 of file Halomodel.h.

◆ RhoFunction

using FML::COSMOLOGY::HALOMODEL::RhoFunction = typedef std::function<double(double,double,HaloModel *hm)>

Definition at line 37 of file Halomodel.h.

◆ SphericalCollapse

using FML::COSMOLOGY::HALOMODEL::SphericalCollapse = typedef FML::COSMOLOGY::SPHERICALCOLLAPSE::SphericalCollapse

Definition at line 33 of file Halomodel.h.

◆ SphericalCollapseModel

using FML::COSMOLOGY::HALOMODEL::SphericalCollapseModel = typedef FML::COSMOLOGY::SPHERICALCOLLAPSE::SphericalCollapseModel

Definition at line 34 of file Halomodel.h.

◆ Spline

using FML::COSMOLOGY::HALOMODEL::Spline = typedef FML::INTERPOLATION::SPLINE::Spline

Definition at line 29 of file Halomodel.h.

◆ Spline2D

using FML::COSMOLOGY::HALOMODEL::Spline2D = typedef FML::INTERPOLATION::SPLINE::Spline2D

Definition at line 30 of file Halomodel.h.

◆ Timer

using FML::COSMOLOGY::HALOMODEL::Timer = typedef FML::UTILS::Timings

Definition at line 23 of file Halomodel.h.

◆ WindowFunction

using FML::COSMOLOGY::HALOMODEL::WindowFunction = typedef std::function<double(double)>

Definition at line 41 of file Halomodel.h.

Function Documentation

◆ compute_formation_redshift_factor()

void FML::COSMOLOGY::HALOMODEL::compute_formation_redshift_factor	(	const Spline &	growthfactor_of_x_spline,
		const Spline &	logsigma_of_logR_spline,
		Spline &	formationredshift_of_logM_spline,
		double	deltac,
		double	OmegaM,
		double	xcollapse,
		double	massfraction
	)

Definition at line 580 of file Halomodel.cpp.

                              {
 
        // Solve D(zf)/D(z) sigma(massfraction * M, z) == deltac(z) for zf
        // where z is collapse redshift
        auto compute_formation_redshift = [&](double M){
          const double Mform = massfraction * M;
          const double R = RofM(Mform, OmegaM, xcollapse);
          auto condition = [&](double x){
            const double xminspline = growthfactor_of_x_spline.get_xrange().first;
            double growthratio = growthfactor_of_x_spline(x) / growthfactor_of_x_spline(xcollapse);
            if(x < xminspline){
              growthratio = growthfactor_of_x_spline(xminspline) / growthfactor_of_x_spline(xcollapse);
              growthratio *= std::exp(x-xminspline);
            }
            const double sigma = std::exp(logsigma_of_logR_spline(std::log(R)));
            return growthratio * sigma - deltac;
          };
          double xformation = xcollapse;
          try{
            xformation = FML::MATH::find_root_bisection(condition, {std::log(1e-10), xcollapse}, 1e-8);
          } catch(...){
          }
          return xformation;
        };
 
        //========================================
        // Make R array
        //========================================
        const double Rmin = 1e-8, Rmax = 1e8;
        const int npts_R = 10000;
        const auto logR_array = FML::MATH::linspace(std::log(Rmin), std::log(Rmax), npts_R);
 
        // Compute formation-factor (1+zf(m))/(1+zcollapse) for all radii
        DVector formationtime_array;
        DVector logM_array;
        for(const auto logR : logR_array){
          const double M = MofR(std::exp(logR), OmegaM, xcollapse);
          const double xformation = compute_formation_redshift(M);
          formationtime_array.push_back(std::exp(-xformation)-1.0);
          logM_array.push_back(std::log(M));
        }
        formationredshift_of_logM_spline = Spline(logM_array, formationtime_array, "zf(logM)");
      }

◆ compute_lognu_of_logR()

void FML::COSMOLOGY::HALOMODEL::compute_lognu_of_logR	(	const double	deltac,
		const Spline &	logsigma_of_logR_spline,
		Spline &	lognu_of_logR_spline,
		Spline &	logR_of_lognu_spline,
		const DVector *	fiducial_logR_array
	)

Definition at line 376 of file Halomodel.cpp.

                                             {
 
        //================================================
        // Radius array
        //================================================
        const double logRmin = logsigma_of_logR_spline.get_xrange().first;
        const double logRmax = logsigma_of_logR_spline.get_xrange().second;
        const int npts_R = 500;
        auto logR_array = FML::MATH::linspace(logRmin, logRmax, npts_R);
        if(fiducial_logR_array != nullptr)
          logR_array = *fiducial_logR_array;
 
        DVector lognu_array;
        for(const auto logR : logR_array){
          const double sigma = std::exp(logsigma_of_logR_spline(logR));
          const double nu = deltac/sigma;
          lognu_array.push_back(std::log(nu));
        }
 
        lognu_of_logR_spline = Spline(logR_array, lognu_array, "lognu(logR)");
        logR_of_lognu_spline = Spline(lognu_array, logR_array, "logR(lognu)");
      }

◆ compute_massfunction()

void FML::COSMOLOGY::HALOMODEL::compute_massfunction	(	const Spline &	logR_of_lognu_spline,
		const Spline &	lognu_of_logR_spline,
		Spline &	dndlogM_of_logM_spline,
		Spline &	n_of_logM_spline,
		const Func1 &	pdf_of_nu,
		const double	OmegaM,
		const double	x
	)

Definition at line 408 of file Halomodel.cpp.

                         {
 
        // Units of (h/Mpc)^3 (same as 1/R^3)
        auto dndlogM_of_lognu = [&](double lognu){
          const double nu = std::exp(lognu);
          const double logR = logR_of_lognu_spline(lognu);
          const double R = std::exp(logR);
          const double f = pdf_of_nu(nu);
          const double dnudlogR = nu * lognu_of_logR_spline.deriv_x(logR);
          return f * std::abs(dnudlogR) / (4.0 * M_PI * R * R * R);
        };
 
        //================================================
        // M array
        //================================================
        const double Mmin = 1e5, Mmax = 1e18;
        const int npts_M = 500;
        const auto logM_array = FML::MATH::linspace(std::log(Mmin),std::log(Mmax),npts_M);
        const double dlogM = logM_array[1]-logM_array[0];
 
        DVector dndlogM_array(npts_M), n_array(npts_M, 0.0), lognu_array(npts_M,0.0);
        for(int i = npts_M-1; i >= 0; i--){
          const double logM = logM_array[i];
          const double M = std::exp(logM);
          const double R = RofM(M, OmegaM, x);
          const double logR = std::log(R);
          const double lognu = lognu_of_logR_spline(logR);
          const double dndlogM = dndlogM_of_lognu(lognu);
          dndlogM_array[i] = dndlogM;
          n_array[i] = dndlogM * dlogM;
          if(i < npts_M-1)
            n_array[i] += n_array[i+1];
        }
        dndlogM_of_logM_spline = Spline(logM_array, dndlogM_array, "dndlogM(logM)");
        n_of_logM_spline = Spline(logM_array, n_array, "n(logM)");
      }

◆ compute_neff()

std::pair< double, double > FML::COSMOLOGY::HALOMODEL::compute_neff	(	const Spline &	logsigma_of_logR_spline,
		const Spline &	lognu_of_logR_spline
	)

Definition at line 101 of file Halomodel.cpp.

                                              {
        const double logRmin = std::log(1e-4);
        const double logRmax = std::log(1e4);
        const double epsilon = 1e-10;
 
        const double lognu_goal = std::log(1.0);
        auto condition = [&](double logR){
          return lognu_of_logR_spline(logR) - lognu_goal;
        };
        const double logR = FML::MATH::find_root_bisection(condition, {logRmin, logRmax}, epsilon);
        const double neff = -3.0 - 2.0 * logsigma_of_logR_spline.deriv_x(logR);
        return std::pair(std::exp(logR), neff);
      }

◆ compute_one_and_two_halo_terms()

void FML::COSMOLOGY::HALOMODEL::compute_one_and_two_halo_terms	(	const Spline &	logR_of_lognu_spline,
		const Spline2D &	y_of_logR_and_logk_spline,
		Spline &	logDelta_onehalo_of_logk_spline,
		Spline &	logDelta_twohalo_of_logk_spline,
		Spline &	logDelta_full_of_logk_spline,
		const Func1 &	logDeltaLin_of_logk,
		const Func1 &	pdf_of_nu,
		const Func1 &	bias_of_nu,
		const double	eta_hmcode,
		const double	kmin,
		const double	kmax,
		const int	npts_k
	)

Definition at line 458 of file Halomodel.cpp.

                           {
 
        int accuracy_boost = _ACCURACY_BOOST;
 
        //================================================
        // nu array
        //================================================
        const double numin = 1e-10, numax = 10;
        const int npts_nu = 200 * accuracy_boost;
        const auto lognu_array = FML::MATH::linspace(std::log(numin),std::log(numax),npts_nu);
 
        //================================================
        // k array
        //================================================
        const auto logk_array = FML::MATH::linspace(std::log(kmin), std::log(kmax), npts_k);
        auto k_array = logk_array;
        for(auto & k : k_array) k = std::exp(k);
 
        auto integrand_onehalo = [&logR_of_lognu_spline,&y_of_logR_and_logk_spline,&eta_hmcode,&pdf_of_nu](double lognu, double logk){
          const double nu = std::exp(lognu);
          const double logR = logR_of_lognu_spline(lognu);
          const double y = y_of_logR_and_logk_spline(logR, logk + lognu*eta_hmcode);
          const double f = pdf_of_nu(nu);
          const double Moverrhomean = 4.0 * M_PI / 3.0 * std::exp(3.0*logR);
          return Moverrhomean * f * (y*y) * nu;
        };
 
        auto integrand_twohalo = [&eta_hmcode,&logR_of_lognu_spline,&y_of_logR_and_logk_spline,&pdf_of_nu,&bias_of_nu](double lognu, double logk){
          const double nu = std::exp(lognu);
          const double logR = logR_of_lognu_spline(lognu);
          const double y = y_of_logR_and_logk_spline(logR, logk + lognu*eta_hmcode);
          const double f = pdf_of_nu(nu);
          const double bias = bias_of_nu(nu);
          return bias * f * y * nu;
        };
 
        // This should integrate to 1
        auto integrand_norm_onehalo = [&pdf_of_nu](double lognu){
          const double nu = std::exp(lognu);
          const double f = pdf_of_nu(nu);
          return f * nu;
        };
 
        // This should integrate to 1
        auto integrand_norm_twohalo = [&pdf_of_nu,&bias_of_nu](double lognu){
          const double nu = std::exp(lognu);
          const double f = pdf_of_nu(nu);
          const double bias = bias_of_nu(nu);
          return bias * f * nu;
        };
 
        // Compute Int fdnu and Int bf dnu
        // Both which should be unity
        double normone = 0.0;
        double normtwo = 0.0;
        for(const auto lognu : lognu_array){
          normone += integrand_norm_onehalo(lognu);
          normtwo += integrand_norm_twohalo(lognu);
        }
        normone *= (lognu_array[1]-lognu_array[0]);
        normtwo *= (lognu_array[1]-lognu_array[0]);
        // Sanity check Int f dnu = 1 and Int f b dnu = 1
        if(std::abs(normone-1.0) > 0.001 or std::abs(normtwo-1.0) > 0.001){
          // Can happen normally if massfunction is not normalized but in any case Int bf / Int f == 1
          std::cout << "Warning: Unity check failed for Int f(nu)dnu = " << normone << " Int b(nu)f(nu)dnu = " << normtwo << "\n";
          std::cout << "Could be that massfunction is not normalized. In that case this should be unity: " << normtwo/normone << "\n";
        }
 
        DVector logDelta_onehalo_array;
        DVector logDelta_twohalo_array;
        DVector logDelta_full_array;
        for(const auto logk : logk_array){
          const double k = std::exp(logk);
 
          // Integrate up one and two halo integrals
          double onehalo = 0.0;
          double twohalo = 0.0;
          double full = 0.0;
          for(const auto lognu : lognu_array){
            onehalo += integrand_onehalo(lognu, logk);
            twohalo += integrand_twohalo(lognu, logk);
          }
          onehalo *= (lognu_array[1]-lognu_array[0]);
          twohalo *= (lognu_array[1]-lognu_array[0]);
 
          // Normalize
          onehalo *= (k*k*k)/(2.0*M_PI*M_PI) / normone;
          twohalo *= 1.0 / normtwo;
          twohalo = twohalo * twohalo * std::exp(logDeltaLin_of_logk(logk));
          full = onehalo + twohalo;
 
          logDelta_onehalo_array.push_back(std::log(onehalo));
          logDelta_twohalo_array.push_back(std::log(twohalo));
          logDelta_full_array.push_back(std::log(full));
 
        }
 
        logDelta_onehalo_of_logk_spline = Spline(logk_array, logDelta_onehalo_array, "logDelta1h(logk)");
        logDelta_twohalo_of_logk_spline = Spline(logk_array, logDelta_twohalo_array, "logDelta2h(logk)");
        logDelta_full_of_logk_spline = Spline(logk_array, logDelta_full_array, "logDelta(logk)");
      }

◆ compute_rho_fourier()

void FML::COSMOLOGY::HALOMODEL::compute_rho_fourier	(	const Func2 &	rho_function,
		Spline2D &	y_of_logR_and_logk_spline,
		const double	DeltaVir,
		const double	OmegaM,
		const double	x,
		const DVector *	fiducial_logR_array,
		const DVector *	fiducial_logk_array
	)

Definition at line 325 of file Halomodel.cpp.

                                             {
 
        int accuracy_boost = _ACCURACY_BOOST;
 
        //================================================
        // k array
        //================================================
        const double kmin = 1e-8, kmax = 1e8;
        const int npts_k = 75 * accuracy_boost;
        auto logk_array = FML::MATH::linspace(std::log(kmin), std::log(kmax), npts_k);
        if(fiducial_logk_array != nullptr)
          logk_array = *fiducial_logk_array;
        auto k_array = logk_array;
        for(auto & k : k_array) k = std::exp(k);
 
        //================================================
        // R array
        //================================================
        const double Rmin = 1e-8, Rmax = 1e8;
        const int npts_R = 75 * accuracy_boost;
        auto logR_array = FML::MATH::linspace(std::log(Rmin), std::log(Rmax), npts_R);
        if(fiducial_logR_array != nullptr)
          logR_array = *fiducial_logR_array;
        auto R_array = logR_array;
        for(auto & R : R_array) R = std::exp(R);
 
        DVector2D yofkandM_array(R_array.size(), DVector(k_array.size()));
#ifdef USE_OMP
#pragma omp parallel for schedule(dynamic,1)
#endif
        for(size_t i = 0; i < R_array.size(); i++){
          const double R = R_array[i];
          const double M = MofR(R, OmegaM, x);
          const double Rvir = R / std::pow(DeltaVir,1.0/3.0);
          const auto result = compute_rho_fourier_single(rho_function, k_array, M, Rvir);
          yofkandM_array[i] = result;
        }
        y_of_logR_and_logk_spline = Spline2D(logR_array, logk_array, yofkandM_array, "y(logR,logk)");
      }

◆ compute_rho_fourier_single()

DVector FML::COSMOLOGY::HALOMODEL::compute_rho_fourier_single	(	const Func2 &	rho_function,
		const DVector &	k_array,
		const double	M,
		const double	Rvir
	)

Definition at line 257 of file Halomodel.cpp.

                            { 
 
        int accuracy_boost = _ACCURACY_BOOST;
 
        //================================================
        // Determine lower integration limits
        // by requiring rho(rmin) < 1e-6 rho(Rvir/10)
        //================================================
        double rho_start =  rho_function(Rvir/10.0,M); 
        int power_min = 1;
        while(power_min < 5) {
          ++power_min;
          double rho = rho_function(Rvir/std::pow(10,power_min),M);
          if(rho < 1e-4 * rho_start) break;
        };
 
        //================================================
        // Determine upper integration limit
        //================================================
        int power_max = 0;
        while(power_max < 7) {
          ++power_max;
          double rho = rho_function(Rvir*std::pow(2,power_max),M);
          if(rho == 0.0) break;
        }
 
        //================================================
        // Radius integration array
        //================================================
        const double rmin = Rvir/std::pow(10,power_min), rmax = Rvir*std::pow(2,power_max);
        const int npts_r = (power_max+power_min) * 50 * accuracy_boost;
        const auto logr_array = FML::MATH::linspace(std::log(rmin), std::log(rmax), npts_r);
        auto r_array = logr_array;
        for(auto & r : r_array) r = std::exp(r);
 
        DVector y_array;
        for(const auto k : k_array){
 
          auto integrand = [&rho_function](double k, double r, double M){
            return rho_function(r, M) * std::sin(k*r)/(k*r);
          };
 
          // Integrate up to get the dimensionless y = rho(k,M)/M
          double integral = 0.0;
          for(const auto r : r_array){
            integral += integrand(k, r, M);
          }
          integral *= (logr_array[1] - logr_array[0]);
          y_array.push_back(integral);
        }
 
        // Ensure normalization in case rho is not properly normalized
        // y(k,M) -> 1 as k -> 0
        const double norm = y_array[0];
        for(auto & y : y_array) y /= norm;
 
        return y_array;
      }

◆ compute_sigma()

void FML::COSMOLOGY::HALOMODEL::compute_sigma	(	const Func1	logDelta_of_logk,
		const WindowFunction	window_of_kR,
		Spline &	logsigma_of_logR_spline,
		const DVector *	fiducial_logR_array,
		const DVector *	fiducial_logk_array
	)

Definition at line 149 of file Halomodel.cpp.

                                             {
 
        //================================================
        // Radius array
        // If not provided then we generate it
        //================================================
        const double Rmin = 1e-8, Rmax = 1e8;
        const int npts_R = 500;
        auto logR_array = FML::MATH::linspace(std::log(Rmin), std::log(Rmax), npts_R);
        if(fiducial_logR_array != nullptr)
          logR_array = *fiducial_logR_array;
        auto R_array = logR_array;
        for(auto & R : R_array) R = std::exp(R);
 
        //================================================
        // k integration array. Linear in logk
        // If not provided then we generate it
        //================================================
        const double kmin = 1e-8, kmax = 1e8;
        const int npts_k = 500;
        auto logk_array = FML::MATH::linspace(std::log(kmin), std::log(kmax), npts_k);
        if(fiducial_logk_array != nullptr)
          logk_array = *fiducial_logk_array;
        auto k_array = logk_array;
        for(auto & k : k_array) k = std::exp(k);
 
        //================================================
        // Compute and spline sigma(R)
        //================================================
        DVector logsigma_array;
        for(const auto R : R_array){
 
          auto integrand = [&logDelta_of_logk,&window_of_kR](double logk, double R){
            const double k = std::exp(logk);
            const double W = window_of_kR(k*R);
            return std::exp(logDelta_of_logk(logk)) * W * W;
          };
 
          double integral = 0.0;
          for(const auto logk : logk_array){
            integral += integrand(logk, R);
          }
          integral *= (logk_array[1]-logk_array[0]);
          logsigma_array.push_back(0.5*std::log(integral));
        }
        logsigma_of_logR_spline = Spline(logR_array, logsigma_array, "sigma(R)");
      }

◆ compute_sigmav()

double FML::COSMOLOGY::HALOMODEL::compute_sigmav	(	const Func1	logDelta_of_logk,
		const WindowFunction	window_of_kR,
		double	R,
		const DVector *	fiducial_logk_array
	)

Definition at line 121 of file Halomodel.cpp.

                                             {
        auto integrand = [&](double logk){
          const double W = window_of_kR(std::exp(logk)*R);
          return W*W*std::exp(logDelta_of_logk(logk))*std::exp(-2*logk)/3.0;
        };
 
        const double kmin = 1e-8, kmax = 1e4;
        const int npts_k = 1000;
        auto logk_array = FML::MATH::linspace(std::log(kmin), std::log(kmax), npts_k);
        if(fiducial_logk_array != nullptr)
          logk_array = *fiducial_logk_array;
        auto k_array = logk_array;
        for(auto & k : k_array) k = std::exp(k);
        double integral = 0.0;
        for(const auto logk : logk_array){
          integral += integrand(logk);
        }
        integral *= (logk_array[1]-logk_array[0]);
        return std::sqrt(integral);
      }

◆ get_nfw_halo_density_profile()

RhoFunction FML::COSMOLOGY::HALOMODEL::get_nfw_halo_density_profile ( )

Definition at line 1017 of file Halomodel.cpp.

                                                {
        return [&](double r, double M, HaloModel *hm){
          const double DeltaVir = hm->DeltaVir;
          const double OmegaM = hm->OmegaM;
          const double x = hm->xcollapse;
          const double c = hm->halo_mass_concentration_relation(M,hm);
          const double R = RofM(M, OmegaM, x);
          double Rvir = R / std::pow(DeltaVir,1.0/3.0);
 
          // The modification Rvir *= nu^eta is for the halo bloating effect
          // Fiducial value is 0.0 which has no effect
          // This can be included *either* here or when doing the 1h 2h integrals 
          // Rvir *= std::pow(std::exp(hm->lognu_of_logR_spline(std::log(R))), hm->hmcode_eta);
          // but not both. Currently we do this in the y-integral
 
          const double rs = Rvir / c;
          if(r > Rvir) return 0.0;
          return std::pow(r/rs, 3) / ( (r/rs) * (1 + r/rs) * (1 + r/rs) ) / (std::log(1 + c) - c/(1 + c));
        };
      }

◆ get_nfw_halo_mass_concentration_relation()

cofMFunction FML::COSMOLOGY::HALOMODEL::get_nfw_halo_mass_concentration_relation ( )

Definition at line 1001 of file Halomodel.cpp.

                                                             {
        return [&](double M, HaloModel *hm){
          if(hm->cofM_model == "Bullock2001"){
            const double zf = hm->formationredshift_of_logM_spline(std::log(M));
            const double z  = std::exp(-hm->xcollapse)-1;
            return hm->hmcode_cmin * std::max(1.0, (1.0+zf)/(1+z) );
          }
          // Simple c = c0/(1+z) * (M/M0)^(-beta) power-law model
          return hm->cofM_0 * std::pow(M / hm->cofM_M0, -hm->cofM_beta) * std::exp(hm->xcollapse);
        };
      }

◆ get_sheth_tormen_halo_bias_function()

HaloBiasFunction FML::COSMOLOGY::HALOMODEL::get_sheth_tormen_halo_bias_function ( )

Definition at line 665 of file Halomodel.cpp.

                                                            {
        return [&](double nu, HaloModel *hm){
          const double deltac = hm->deltac;
          const double p_ST = hm->p_ST;
          const double a_ST = hm->a_ST;
          return 1.0 + (a_ST*nu*nu-1.0)/deltac + 2.0*p_ST/deltac/(1.0+std::pow(a_ST*nu*nu,p_ST));
        };
      }

◆ get_sheth_tormen_halo_multiplicity_function()

HaloMultiplicityFunction FML::COSMOLOGY::HALOMODEL::get_sheth_tormen_halo_multiplicity_function ( )

Definition at line 657 of file Halomodel.cpp.

                                                                            {
        return [&](double nu, HaloModel *hm){
          const double A_ST = hm->A_ST;
          const double p_ST = hm->p_ST;
          const double a_ST = hm->a_ST;
          return A_ST * (1.0 + std::pow(a_ST*nu*nu,-p_ST)) * std::exp(-a_ST*nu*nu/2.0);
        };
      }

◆ get_tophat_window_fourier()

WindowFunction FML::COSMOLOGY::HALOMODEL::get_tophat_window_fourier ( )

Definition at line 677 of file Halomodel.cpp.

                                                {
        return [](double kR){
          if(kR < 1e-6) return 1.0 - kR*kR/10.;
          return 3.0/(kR*kR*kR)*(std::sin(kR)-kR*std::cos(kR));
        };
      }

◆ MofR()

double FML::COSMOLOGY::HALOMODEL::MofR	(	double	R,
		double	OmegaM,
		double	x
	)

Definition at line 241 of file Halomodel.cpp.

                   {
        const double Mpl2overMsunH0overh = 2.493e21;
        const double H0overh_in_hmpc = 1.0 / 2997.92458;
        const double rho_mean_in_msunh_over_mpch3 = 3.0 * OmegaM * std::exp(-3.0*x) 
          * std::pow(H0overh_in_hmpc,3) * Mpl2overMsunH0overh;
        return 4.0 * M_PI / 3.0 * rho_mean_in_msunh_over_mpch3 * std::pow(R, 3);
      }

◆ MofRvir()

double FML::COSMOLOGY::HALOMODEL::MofRvir	(	double	Rvir,
		double	DeltaVir,
		double	OmegaM,
		double	x
	)

Definition at line 230 of file Halomodel.cpp.

                   {
        return MofR(Rvir, OmegaM, x) * DeltaVir;
     }

◆ powerlaw_extend_logY_of_logx_spline()

void FML::COSMOLOGY::HALOMODEL::powerlaw_extend_logY_of_logx_spline	(	Spline &	logY_of_logx_spline,
		double	xmin_new,
		double	xmax_new,
		int	npts_new,
		std::string	label = `"logY(log(x)) extrapolated"`
	)

Definition at line 22 of file Halomodel.cpp.

                                                        {
 
        // Fetch logx-array from spline
        const auto logx_array = logY_of_logx_spline.get_x_data();
        const double logxmin = logx_array.front();
        const double logxmax = logx_array.back();
 
        // Compute slope near end-points
        double slope_min  = logY_of_logx_spline.deriv_x(logxmin*0.99);
        double slope_max = logY_of_logx_spline.deriv_x(logxmax*1.01);
 
        // Extrapolation function
        auto logY_extrapolated = [&](double logx){
          if(logx < logxmin){
            return logY_of_logx_spline(logxmin) + (logx-logxmin) * slope_min;
          } else if(logx > logxmax){
            return logY_of_logx_spline(logxmax) + (logx-logxmax) * slope_max;
          }
          return logY_of_logx_spline(logx);
        };
 
        // Extend spline to one covering a larger region
        const auto logx_array_new = FML::MATH::linspace(std::log(xmin_new), std::log(xmax_new), npts_new);
        DVector logY_array_new;
        for(const auto logx : logx_array_new){
          logY_array_new.push_back( logY_extrapolated(logx) );
        }
        logY_of_logx_spline = Spline(logx_array_new, logY_array_new, label);
      }

◆ read_pofk_file()

void FML::COSMOLOGY::HALOMODEL::read_pofk_file	(	const std::string	filename_pofk,
		Spline &	logDelta_of_logk_spline,
		double	normalization,
		int	col_k,
		int	col_pofk
	)

Definition at line 61 of file Halomodel.cpp.

                       {
 
        // Read file
        auto data = FML::FILEUTILS::loadtxt(filename_pofk);
 
        // Extract data, transform to Delta, normalize and spline
        DVector logk_array, logDelta_array;
        for(size_t i = 0; i < data.size(); i++){
          double k = data[i][col_k];
          double pofk = data[i][col_pofk];
          double Delta_over_pofk = (k*k*k) / (2.0*M_PI*M_PI);
          logk_array.push_back( std::log(k) );
          logDelta_array.push_back( std::log( pofk * Delta_over_pofk * normalization ) );
        }
 
        // Make spline
        logDelta_of_logk_spline = Spline(logk_array, logDelta_array, "logDelta(logk)");
 
        // Extrapolate assuming ~ k^n at each end
        // and remake spline
        const double kmin = 1e-12;
        const double kmax = 1e12;
        const int npts_spline = 1000;
        powerlaw_extend_logY_of_logx_spline(
            logDelta_of_logk_spline,
            kmin,
            kmax,
            npts_spline,
            "logDelta(logk) extrapolated");
      }

◆ RofM()

double FML::COSMOLOGY::HALOMODEL::RofM	(	double	M,
		double	OmegaM,
		double	x
	)

Definition at line 216 of file Halomodel.cpp.

                   {
        const double Mpl2overMsunH0overh = 2.493e21;
        const double H0overh_in_hmpc = 1.0 / 2997.92458;
        const double rho_mean_in_msunh_over_mpch3 = 3.0 * OmegaM * std::exp(-3.0*x) 
          * std::pow(H0overh_in_hmpc,3) * Mpl2overMsunH0overh;
        return std::pow(M / (4.0 * M_PI / 3.0 * rho_mean_in_msunh_over_mpch3), 1.0/3.0);
      }

◆ RvirofM()

double FML::COSMOLOGY::HALOMODEL::RvirofM	(	double	M,
		double	DeltaVir,
		double	OmegaM,
		double	x
	)

Definition at line 205 of file Halomodel.cpp.

                   {
        return RofM(M, OmegaM, x) * std::pow(DeltaVir, -1.0/3.0);
      }

Variable Documentation

◆ _ACCURACY_BOOST

int FML::COSMOLOGY::HALOMODEL::_ACCURACY_BOOST = 1

Definition at line 16 of file Halomodel.cpp.

Data Structures

Typedefs

Functions

Variables