power.c File Reference

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stddef.h>
#include <mpi.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_interp.h>
#include <bigfile-mpi.h>
#include <libgadget/cosmology.h>
#include <libgadget/utils.h>
#include <libgadget/physconst.h>
#include "power.h"
#include "proto.h"

Include dependency graph for power.c:

Go to the source code of this file.

Classes
struct	table

Macros
#define	NUGGET   1e-30
#define	MAXCOLS   9

Typedefs
typedef void(*	_parse_fn) (int i, double k, char *line, struct table *, int *InputInLog10, const double InitTime, int NumCol)

Functions
static double	Delta_EH (double k)
static double	Delta_Tabulated (double k, enum TransferType Type)
static double	sigma2_int (double k, void *params)
static double	TopHatSigma2 (double R)
static double	tk_eh (double k)
double	DeltaSpec (double k, enum TransferType Type)
static double	get_Tabulated (double k, enum TransferType Type, double oobval)
double	dlogGrowth (double kmag, enum TransferType Type)
static void	save_transfer (BigFile *bf, int ncol, struct table *ttable, const char *bname, int ThisTask, const char *colnames[])
void	save_all_transfer_tables (BigFile *bf, int ThisTask)
void	parse_power (int i, double k, char *line, struct table *out_tab, int *InputInLog10, const double InitTime, int NumCol)
void	parse_transfer (int i, double k, char *line, struct table *out_tab, int *InputInLog10, const double InitTime, int NumCol)
void	read_power_table (int ThisTask, const char *inputfile, const int ncols, struct table *out_tab, const double InitTime, _parse_fn parse_line)
int	init_transfer_table (int ThisTask, double InitTime, const struct power_params *const ppar)
int	init_powerspectrum (int ThisTask, double InitTime, double UnitLength_in_cm_in, Cosmology *CPin, struct power_params *ppar)

Variables
static double	Norm
static int	WhichSpectrum
static double	PrimordialIndex
static double	UnitLength_in_cm
static Cosmology *	CP
static struct table	power_table
static struct table	transfer_table
static const char *	tnames [MAXCOLS] = {"DELTA_BAR", "DELTA_CDM", "DELTA_NU", "DELTA_CB", "VEL_BAR", "VEL_CDM", "VEL_NU", "VEL_CB", "VEL_TOT"}

Macro Definition Documentation

MAXCOLS

#define MAXCOLS   9

Definition at line 32 of file power.c.

NUGGET

#define NUGGET   1e-30

Definition at line 31 of file power.c.

Typedef Documentation

_parse_fn

typedef void(* _parse_fn) (int i, double k, char *line, struct table *, int *InputInLog10, const double InitTime, int NumCol)

Definition at line 43 of file power.c.

Function Documentation

Delta_EH()

double Delta_EH ( double  k )

static

Definition at line 432 of file power.c.

{
  return sqrt(k * pow(tk_eh(k), 2)* pow(k, PrimordialIndex - 1.0));
}

References PrimordialIndex, and tk_eh().

Referenced by DeltaSpec().

Here is the call graph for this function:

Here is the caller graph for this function:

Delta_Tabulated()

double Delta_Tabulated ( double  k, enum TransferType  Type  )

static

Definition at line 93 of file power.c.

{
    if(Type >= VEL_BAR && Type <= VEL_TOT)
        endrun(1, "Velocity Type %d passed to Delta_Tabulated\n", Type);
 
    return get_Tabulated(k, Type, 0);
}

References endrun(), get_Tabulated(), VEL_BAR, and VEL_TOT.

Referenced by DeltaSpec().

Here is the call graph for this function:

Here is the caller graph for this function:

DeltaSpec()

double DeltaSpec	(	double	k,
		enum TransferType	Type
	)

Definition at line 52 of file power.c.

{
  double power;
 
  if(WhichSpectrum == 2)
      power = Delta_Tabulated(k, Type);
  else
      power = Delta_EH(k);
 
  /*Normalise the power spectrum*/
  power *= Norm;
 
  return power;
}

References Delta_EH(), Delta_Tabulated(), Norm, and WhichSpectrum.

Referenced by density_transfer(), disp_transfer(), print_spec(), sigma2_int(), test_growth_numerical(), test_read_no_rescale(), and test_read_rescale_sigma8().

Here is the call graph for this function:

Here is the caller graph for this function:

dlogGrowth()

double dlogGrowth	(	double	kmag,
		enum TransferType	Type
	)

Definition at line 101 of file power.c.

{
    /*Default to total growth: type 3 is cdm + baryons.*/
    if(Type < DELTA_BAR || Type > DELTA_CB)
        Type = VEL_TOT;
    else
        /*Type should be an offset from the first velocity*/
        Type = (enum TransferType) ((int) VEL_BAR + ((int) Type - (int) DELTA_BAR));
    return get_Tabulated(kmag, Type, 1);
}

References DELTA_BAR, DELTA_CB, get_Tabulated(), VEL_BAR, and VEL_TOT.

Referenced by disp_transfer(), and test_growth_numerical().

Here is the call graph for this function:

Here is the caller graph for this function:

get_Tabulated()

static double get_Tabulated ( double  k, enum TransferType  Type, double  oobval  )

static

Definition at line 69 of file power.c.

{
    /*Convert k to Mpc/h*/
    const double scale = (CM_PER_MPC / UnitLength_in_cm);
    const double logk = log10(k*scale);
 
    if(logk < power_table.logk[0] || logk > power_table.logk[power_table.Nentry - 1])
      return oobval;
 
    double logD = gsl_interp_eval(power_table.mat_intp[0], power_table.logk, power_table.logD[0], logk, NULL);
    double trans = 1;
    /*Transfer table stores (T_type(k) / T_tot(k))*/
    if(transfer_table.Nentry > 0)
       if(Type >= DELTA_BAR && Type < DELTA_TOT)
          trans = gsl_interp_eval(transfer_table.mat_intp[Type], transfer_table.logk, transfer_table.logD[Type], logk, NULL);
 
    /*Convert delta from (Mpc/h)^3/2 to kpc/h^3/2*/
    logD += 1.5 * log10(scale);
    double delta = (pow(10.0, logD)-NUGGET) * trans;
    if(!isfinite(delta))
        endrun(1,"infinite delta or growth: %g for k = %g, Type = %d (tk = %g, logD = %g)\n",delta, k, Type, trans, logD);
    return delta;
}

References CM_PER_MPC, DELTA_BAR, DELTA_TOT, endrun(), table::logD, table::logk, table::mat_intp, table::Nentry, NUGGET, power_table, transfer_table, and UnitLength_in_cm.

Referenced by Delta_Tabulated(), and dlogGrowth().

Here is the call graph for this function:

Here is the caller graph for this function:

init_powerspectrum()

int init_powerspectrum	(	int	ThisTask,
		double	InitTime,
		double	UnitLength_in_cm_in,
		Cosmology *	CPin,
		struct power_params *	ppar
	)

Definition at line 394 of file power.c.

{
    WhichSpectrum = ppar->WhichSpectrum;
    /*Used only for tk_eh*/
    PrimordialIndex = ppar->PrimordialIndex;
    UnitLength_in_cm = UnitLength_in_cm_in;
    CP = CPin;
 
    if(ppar->WhichSpectrum == 2) {
        read_power_table(ThisTask, ppar->FileWithInputSpectrum, 1, &power_table, InitTime, parse_power);
        /*Initialise the interpolation*/
        gsl_interp_init(power_table.mat_intp[0],power_table.logk, power_table.logD[0],power_table.Nentry);
        transfer_table.Nentry = 0;
        if(ppar->DifferentTransferFunctions || ppar->ScaleDepVelocity) {
            init_transfer_table(ThisTask, InitTime, ppar);
        }
    }
 
    Norm = 1.0;
    if(ppar->InputPowerRedshift >= 0 || ppar->Sigma8 > 0) {
        const double R8 = 8 * (CM_PER_MPC / UnitLength_in_cm);  /* 8 Mpc/h */
        if(ppar->Sigma8 > 0) {
            double res = TopHatSigma2(R8);
            if(isfinite(res) && res > 0)
                Norm = ppar->Sigma8 / sqrt(res);
            else
                endrun(1, "Could not normalize P(k) to Sigma8=%g! Measured Sigma8^2 is %g\n", ppar->Sigma8, res);
        }
        double Dplus = GrowthFactor(CP, InitTime, 1/(1+ppar->InputPowerRedshift));
        if(ppar->InputPowerRedshift >= 0) {
            Norm *= Dplus;
            message(0,"Growth factor from z=%g (InputPowerRedshift) to z=%g (Init): %g \n", ppar->InputPowerRedshift, 1. / InitTime - 1, Dplus);
        }
        message(0, "Normalization adjusted to  Sigma8=%g (at z=0)  (Normfac=%g). \n", sqrt(TopHatSigma2(R8))/Dplus, Norm);
    }
    return power_table.Nentry;
}

References CM_PER_MPC, CP, power_params::DifferentTransferFunctions, endrun(), power_params::FileWithInputSpectrum, GrowthFactor(), init_transfer_table(), power_params::InputPowerRedshift, table::logD, table::logk, table::mat_intp, message(), table::Nentry, Norm, parse_power(), power_table, PrimordialIndex, power_params::PrimordialIndex, read_power_table(), power_params::ScaleDepVelocity, power_params::Sigma8, ThisTask, TopHatSigma2(), transfer_table, UnitLength_in_cm, WhichSpectrum, and power_params::WhichSpectrum.

Referenced by main(), test_growth_numerical(), test_read_no_rescale(), and test_read_rescale_sigma8().

Here is the call graph for this function:

Here is the caller graph for this function:

init_transfer_table()

int init_transfer_table	(	int	ThisTask,
		double	InitTime,
		const struct power_params *const	ppar
	)

Definition at line 322 of file power.c.

{
    int i, t;
    const int nnu = (CP->MNu[0] > 0) + (CP->MNu[1] > 0) + (CP->MNu[2] > 0);
    if(strlen(ppar->FileWithTransferFunction) > 0) {
        read_power_table(ThisTask, ppar->FileWithTransferFunction, MAXCOLS, &transfer_table, InitTime, parse_transfer);
    }
    if(transfer_table.Nentry == 0) {
        endrun(1, "Could not read transfer table at: '%s'\n",ppar->FileWithTransferFunction);
    }
    /*Normalise the transfer functions.*/
 
    /*Now normalise the velocity transfer functions: divide by a * hubble, where hubble is the hubble function in Mpc^-1, so H0/c*/
    const double fac = InitTime * hubble_function(CP, InitTime)/CP->Hubble * 100 * CP->HubbleParam/(LIGHTCGS / 1e5);
    const double onu = get_omega_nu(&CP->ONu, InitTime)*pow(InitTime,3);
    double meangrowth[VEL_TOT-VEL_BAR+1] = {0};
    /* At this point the transfer table contains: (3,4,5) t_b, 0.5 * h_prime, t_ncdm.
     * After, if t_cdm = 0.5 h_prime / (a H(a) / H0 /c) we need:
     * (t_b + t_cdm) / d_b, t_cdm/d_cdm, (t_ncdm + t_cdm) / d_ncdm*/
    for(i=0; i< transfer_table.Nentry; i++) {
        /* Now row 4 is t_cdm*/
        transfer_table.logD[VEL_CDM][i] /= fac;
        transfer_table.logD[VEL_BAR][i] += transfer_table.logD[VEL_CDM][i];
        transfer_table.logD[VEL_NU][i] += transfer_table.logD[VEL_CDM][i];
 
        /*Set up the CDM + baryon rows*/
        transfer_table.logD[DELTA_CB][i] = CP->OmegaBaryon * transfer_table.logD[DELTA_BAR][i] + CP->OmegaCDM * transfer_table.logD[DELTA_CDM][i];
        transfer_table.logD[VEL_CB][i] = CP->OmegaBaryon * transfer_table.logD[VEL_BAR][i] + CP->OmegaCDM * transfer_table.logD[VEL_CDM][i];
        /*total growth and delta: start as CDM + baryon and then add nu if needed.*/
        transfer_table.logD[VEL_TOT][i] = transfer_table.logD[VEL_CB][i];
        double T_tot = transfer_table.logD[DELTA_CB][i];
        /*Normalise the cb rows*/
        double Omega0a3 = CP->OmegaBaryon + CP->OmegaCDM;
        transfer_table.logD[DELTA_CB][i] /= Omega0a3;
        transfer_table.logD[VEL_CB][i] /= Omega0a3;
        /*Total matter density in T_tot. Neutrinos may be slightly relativistic, so
         * Omega0a3 >= CP->Omega0 if neutrinos are massive.*/
        if(nnu > 0) {
            /*Add neutrino growth to total growth*/
            transfer_table.logD[VEL_TOT][i] += onu *  transfer_table.logD[VEL_NU][i];
            T_tot += onu * transfer_table.logD[DELTA_NU][i];
            Omega0a3 += onu;
        }
        /*Normalise the totals now we have neutrinos*/
        transfer_table.logD[VEL_TOT][i] /= Omega0a3;
        T_tot /= Omega0a3;
        /*Normalize growth_i and delta_i by delta_tot */
        for(t = DELTA_BAR; t <= VEL_TOT; t++) {
            transfer_table.logD[t][i] /= T_tot;
        }
    }
 
    /*Now compute mean growths*/
    for(t = VEL_BAR; t <= VEL_TOT; t++) {
        int nmean=0;
        for(i=0; i< transfer_table.Nentry; i++)
            if(transfer_table.logk[i] > power_table.logk[0]) {
                meangrowth[t-VEL_BAR] += transfer_table.logD[t][i];
                nmean++;
            }
        if(nmean > 0)
            meangrowth[t-VEL_BAR]/= nmean;
    }
    /*Initialise the interpolation*/
    for(t = 0; t < MAXCOLS; t++)
        gsl_interp_init(transfer_table.mat_intp[t],transfer_table.logk, transfer_table.logD[t],transfer_table.Nentry);
 
    message(0,"Scale-dependent growth calculated. Mean = %g %g %g %g %g\n",meangrowth[0], meangrowth[1], meangrowth[2], meangrowth[3], meangrowth[4]);
    message(0, "Power spectrum rows: %d, Transfer: %d (%g -> %g)\n", power_table.Nentry, transfer_table.Nentry, transfer_table.logD[DELTA_BAR][0],transfer_table.logD[DELTA_BAR][transfer_table.Nentry-1]);
    return transfer_table.Nentry;
}

References CP, DELTA_BAR, DELTA_CB, DELTA_CDM, DELTA_NU, endrun(), power_params::FileWithTransferFunction, get_omega_nu(), Cosmology::Hubble, hubble_function(), Cosmology::HubbleParam, LIGHTCGS, table::logD, table::logk, table::mat_intp, MAXCOLS, message(), Cosmology::MNu, table::Nentry, Cosmology::OmegaBaryon, Cosmology::OmegaCDM, Cosmology::ONu, parse_transfer(), power_table, read_power_table(), ThisTask, transfer_table, VEL_BAR, VEL_CB, VEL_CDM, VEL_NU, and VEL_TOT.

Referenced by init_powerspectrum().

Here is the call graph for this function:

Here is the caller graph for this function:

parse_power()

void parse_power	(	int	i,
		double	k,
		char *	line,
		struct table *	out_tab,
		int *	InputInLog10,
		const double	InitTime,
		int	NumCol
	)

Definition at line 158 of file power.c.

{
    char * retval;
    if((*InputInLog10) == 0) {
        if(k < 0) {
            message(1, "some input k is negative, guessing the file is in log10 units\n");
            *InputInLog10 = 1;
        }
        else
            k = log10(k);
    }
    out_tab->logk[i] = k;
    retval = strtok(NULL, " \t");
    if(!retval)
        endrun(1,"Incomplete line in power spectrum: %s\n",line);
    double p = atof(retval);
    if ((*InputInLog10) == 0)
        p = log10(p+NUGGET);
    /*Store delta, square root of power*/
    out_tab->logD[0][i] = p/2;
}

References endrun(), table::logD, table::logk, message(), and NUGGET.

Referenced by init_powerspectrum().

Here is the call graph for this function:

Here is the caller graph for this function:

parse_transfer()

void parse_transfer	(	int	i,
		double	k,
		char *	line,
		struct table *	out_tab,
		int *	InputInLog10,
		const double	InitTime,
		int	NumCol
	)

Definition at line 180 of file power.c.

{
    int j;
    int ncols = NumCol - 1; /* The first column k is already read in read_power_table. */
    int nnu = round((ncols - 15)/2);
    double * transfers = (double *) mymalloc("transfers", sizeof(double) * ncols);
    k = log10(k);
    out_tab->logk[i] = k;
    /* Note: the ncdm entries change depending on the number of neutrino species. The first row, k,
     * is read in read_power_table and passed as a parameter.
     * but only the first is used for particles.
     * 1:k (h/Mpc)              2:d_g                    3:d_b                    4:d_cdm                  5:d_ur
     * 6:d_ncdm[0]              7:d_ncdm[1]              8:d_ncdm[2]              9:d_tot                 10:phi
     * 11:psi                   12:h                     13:h_prime               14:eta                   15:eta_prime
     * 16:t_g                   17:t_b                   18:t_ur
     * 19:t_ncdm[0]             20:t_ncdm[1]             21:t_ncdm[2]             22:t_tot*/
    for(j = 0; j< ncols; j++) {
        char * retval = strtok(NULL, " \t");
        /*This happens if we do not have as many neutrino species as we expect, or we don't find h_prime.*/
        if(!retval)
            endrun(1,"Incomplete line in power spectrum: only %d columns, expecting %d. Did you remember to set extra metric transfer functions=y?\n",j, ncols);
        transfers[j] = atof(retval);
    }
    /*Order of the transfer table matches the particle types:
     * 0 is baryons, 1 is CDM, 2 is massive neutrinos (we use the most massive neutrino species).
     * 3 is growth function for baryons, 4 is growth function for CDM, 5 is growth function for massive neutrinos.
     * We use the formulae for velocities from fastpm in synchronous gauge:
     * https://github.com/rainwoodman/fastpm-python/blob/02ce2ff87897f713c7b9204630f4e0257d703784/fastpm/multi.py#L185
     * CDM = - h_prime / 2 / d_cdm
     * bar = - (h_prime / 2  + t_b) / d_b
     * nu = - (h_prime / 2 + t_ncdm) / d_ncdm
     * and there is a normalisation factor of (1+z)/ hubble applied later on.
     *
     * See http://adsabs.harvard.edu/abs/1995ApJ...455....7M
     * Eq. 42 (cdm), not using 49(ur), eq. 66 (baryons, and ncdm, truncation at the same order as baryon).
     * These are in CLASS units: they are converted to Gadget units in init_transfer_table().
     * */
    out_tab->logD[DELTA_BAR][i] = -1*transfers[1];
    out_tab->logD[DELTA_CDM][i] = -1*transfers[2];
    /*This should be the weighted average sum of the three neutrino species*/
    const _omega_nu * Onu = &CP->ONu;
    out_tab->logD[DELTA_NU][i] = 0;
    for(j=0; j < nnu; j++)
        out_tab->logD[DELTA_NU][i] = -1*transfers[4+j] * omega_nu_single(Onu, InitTime, j);
    const double onu = get_omega_nu(&CP->ONu, InitTime);
    /*Should be weighted by omega_nu*/
    out_tab->logD[DELTA_NU][i] /= onu;
    /*h_prime is entry 8 + nnu. t_b is 12 + nnu, t_ncdm[2] is 13 + nnu * 2.*/
    out_tab->logD[VEL_BAR][i] = transfers[12+nnu];
    out_tab->logD[VEL_CDM][i] = transfers[8+nnu] * 0.5;
    out_tab->logD[VEL_NU][i] = 0;
    for(j=0; j < nnu; j++)
        out_tab->logD[VEL_NU][i] = transfers[13 + nnu + j] * omega_nu_single(Onu, InitTime, j);
    /*Should be weighted by omega_nu*/
    out_tab->logD[VEL_NU][i] /= onu;
    myfree(transfers);
}

References CP, DELTA_BAR, DELTA_CDM, DELTA_NU, endrun(), get_omega_nu(), table::logD, table::logk, myfree, mymalloc, omega_nu_single(), Cosmology::ONu, VEL_BAR, VEL_CDM, and VEL_NU.

Referenced by init_transfer_table().

Here is the call graph for this function:

Here is the caller graph for this function:

read_power_table()

void read_power_table	(	int	ThisTask,
		const char *	inputfile,
		const int	ncols,
		struct table *	out_tab,
		const double	InitTime,
		_parse_fn	parse_line
	)

Definition at line 238 of file power.c.

{
    FILE *fd = NULL;
    int j;
    int InputInLog10 = 0;
 
    if(ThisTask == 0) {
        if(!(fd = fopen(inputfile, "r")))
            endrun(1, "can't read input spectrum in file '%s' on task %d\n", inputfile, ThisTask);
 
        out_tab->Nentry = 0;
        do
        {
            char buffer[1024];
            char * retval = fgets(buffer, 1024, fd);
            /*Happens on end of file*/
            if(!retval)
                break;
            retval = strtok(buffer, " \t");
            if(!retval || retval[0] == '#')
                continue;
            out_tab->Nentry++;
        }
        while(1);
        rewind(fd);
    }
    MPI_Bcast(&(out_tab->Nentry), 1, MPI_INT, 0, MPI_COMM_WORLD);
 
    if(out_tab->Nentry < 2)
        endrun(1, "Input spectrum too short\n");
    out_tab->logk = (double *) mymalloc("Powertable", (ncols+1)*out_tab->Nentry * sizeof(double));
    for(j=0; j<ncols; j++)
        out_tab->logD[j] = out_tab->logk + (j+1)*out_tab->Nentry;
 
    if(ThisTask == 0)
    {
        /* detect the columns of the input file */
        char line1[1024];
 
        while(fgets(line1,1024,fd))
        {
            char * content = strtok(line1, " \t");
            if(content[0] != '#') /*Find the first line*/
                break;
        }
        int Ncolumns = 0;
        char *c;
        do
        {
            Ncolumns++;
            c = strtok(NULL," \t");
        }
        while(c != NULL);
 
        rewind(fd);
        message(0, "Detected %d columns in file '%s'. \n", Ncolumns, inputfile);
 
        int i = 0;
        do
        {
            char buffer[1024];
            char * line = fgets(buffer, 1024, fd);
            /*Happens on end of file*/
            if(!line)
                break;
            char * retval = strtok(line, " \t");
            if(!retval || retval[0] == '#')
                continue;
            double k = atof(retval);
            parse_line(i, k, line, out_tab, &InputInLog10, InitTime, Ncolumns);
            i++;
        }
        while(1);
 
        fclose(fd);
    }
 
    MPI_Bcast(out_tab->logk, (ncols+1)*out_tab->Nentry, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    for(j=0; j<ncols; j++) {
        out_tab->mat_intp[j] = gsl_interp_alloc(gsl_interp_cspline,out_tab->Nentry);
    }
}

References endrun(), table::logD, table::logk, table::mat_intp, message(), mymalloc, table::Nentry, and ThisTask.

Referenced by init_powerspectrum(), and init_transfer_table().

Here is the call graph for this function:

Here is the caller graph for this function:

save_all_transfer_tables()

void save_all_transfer_tables	(	BigFile *	bf,
		int	ThisTask
	)

Definition at line 150 of file power.c.

{
    const char * pname = "DELTA_MAT";
    save_transfer(bf, 1, &power_table, "ICPower", ThisTask, &pname);
    save_transfer(bf, MAXCOLS, &transfer_table, "ICTransfers", ThisTask, tnames);
}

References MAXCOLS, power_table, save_transfer(), ThisTask, tnames, and transfer_table.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

save_transfer()

static void save_transfer ( BigFile *  bf, int  ncol, struct table *  ttable, const char *  bname, int  ThisTask, const char *  colnames[]  )

static

Definition at line 113 of file power.c.

{
    BigBlock btransfer;
    int i;
    if(0 != big_file_mpi_create_block(bf, &btransfer, bname, NULL, 0, 0, 0, MPI_COMM_WORLD)) {
        endrun(0, "failed to create block %s:%s", bname,
                big_file_get_error_message());
    }
 
    if ( (0 != big_block_set_attr(&btransfer, "Nentry", &(ttable->Nentry), "u8", 1)) ) {
        endrun(0, "Failed to write table size %s\n",
                    big_file_get_error_message());
    }
    if(0 != big_block_mpi_close(&btransfer, MPI_COMM_WORLD)) {
        endrun(0, "Failed to close block %s\n",
                    big_file_get_error_message());
    }
    size_t dims[2] = {0 , 1};
    /*The transfer state is shared between all processors,
     *so only write on master task*/
    if(ThisTask == 0) {
        dims[0] = ttable->Nentry;
    }
 
    char buf[100];
    snprintf(buf, 100, "%s/logk", bname);
 
    _bigfile_utils_create_block_from_c_array(bf, ttable->logk, buf, "f8", dims, sizeof(double), 1, 1, MPI_COMM_WORLD);
 
    for(i = 0; i < ncol; i++)
    {
        snprintf(buf, 100, "%s/%s", bname, colnames[i]);
        _bigfile_utils_create_block_from_c_array(bf, ttable->logD[i], buf, "f8", dims, sizeof(double), 1, 1, MPI_COMM_WORLD);
    }
}

References _bigfile_utils_create_block_from_c_array(), endrun(), table::logD, table::logk, table::Nentry, and ThisTask.

Referenced by save_all_transfer_tables().

Here is the call graph for this function:

Here is the caller graph for this function:

sigma2_int()

double sigma2_int ( double  k, void *  params  )

static

Definition at line 486 of file power.c.

{
  double w, x;
 
  double r_tophat = *(double *) params;
  const double kr = r_tophat * k;
  const double kr2 = kr * kr;
 
  /*Series expansion; actually good until kr~1*/
  if(kr < 1e-3)
      w = 1./3. - kr2/30. +kr2*kr2/840.;
  else
      w = 3 * (sin(kr) / kr - cos(kr)) / kr2;
  x = 4 * M_PI / (2 * M_PI * 2 * M_PI * 2 * M_PI) * k * k * w * w * pow(DeltaSpec(k, DELTA_TOT),2);
 
  return x;
 
}

References DELTA_TOT, and DeltaSpec().

Referenced by TopHatSigma2().

Here is the call graph for this function:

Here is the caller graph for this function:

tk_eh()

double tk_eh ( double  k )

static

Definition at line 438 of file power.c.

{
  double q, theta, ommh2, a, s, gamma, L0, C0;
  double tmp;
  double omegam, ombh2, hubble;
 
  /* other input parameters */
  hubble = CP->HubbleParam;
 
  omegam = CP->Omega0;
  ombh2 = CP->OmegaBaryon * CP->HubbleParam * CP->HubbleParam;
 
  if(CP->OmegaBaryon == 0)
    ombh2 = 0.044 * CP->HubbleParam * CP->HubbleParam;
 
  k *= (CM_PER_MPC / UnitLength_in_cm); /* convert to h/Mpc */
 
  theta = 2.728 / 2.7;
  ommh2 = omegam * hubble * hubble;
  s = 44.5 * log(9.83 / ommh2) / sqrt(1. + 10. * exp(0.75 * log(ombh2))) * hubble;
  a = 1. - 0.328 * log(431. * ommh2) * ombh2 / ommh2
    + 0.380 * log(22.3 * ommh2) * (ombh2 / ommh2) * (ombh2 / ommh2);
  gamma = a + (1. - a) / (1. + exp(4 * log(0.43 * k * s)));
  gamma *= omegam * hubble;
  q = k * theta * theta / gamma;
  L0 = log(2. * exp(1.) + 1.8 * q);
  C0 = 14.2 + 731. / (1. + 62.5 * q);
  tmp = L0 / (L0 + C0 * q * q);
  return (tmp);
}

References CM_PER_MPC, CP, Cosmology::HubbleParam, Cosmology::Omega0, Cosmology::OmegaBaryon, and UnitLength_in_cm.

Referenced by Delta_EH().

Here is the caller graph for this function:

TopHatSigma2()

double TopHatSigma2 ( double  R )

static

Definition at line 470 of file power.c.

{
  gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);
  double result,abserr;
  gsl_function F;
  F.function = &sigma2_int;
  F.params = &R;
 
  /* note: 500/R is here chosen as integration boundary (infinity) */
  gsl_integration_qags (&F, 0, 500. / R, 0, 1e-4,1000,w,&result, &abserr);
/*   printf("gsl_integration_qng in TopHatSigma2. Result %g, error: %g, intervals: %lu\n",result, abserr,w->size); */
  gsl_integration_workspace_free (w);
  return result;
}

References sigma2_int().

Referenced by init_powerspectrum().

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

CP

Cosmology* CP

static

Definition at line 27 of file power.c.

Referenced by apply_half_kick(), apply_PM_half_kick(), atime_integ(), atime_to_myr(), blackhole(), blackhole_accretion_postprocess(), check_density_entropy(), check_omega(), check_units(), compute_mass(), cooling_and_starformation(), delta_nu_from_power(), density(), displacement_fields(), do_density_test(), do_force_test(), do_heiii_reionization(), drift_all_particles(), drift_integ(), dump_snapshot(), exact_drift_factor(), F_Omega(), find_timesteps(), fof_save_groups(), fof_save_particles(), fof_write_header(), force_tree_build(), fslength(), fslength_int(), get_delta_nu(), get_delta_nu_combined(), get_exact_drift_factor(), get_exact_factor(), get_exact_gravkick_factor(), get_exact_hydrokick_factor(), get_long_range_timestep_dloga(), get_PM_timestep_ti(), gravkick_integ(), gravpm_force(), growth(), growth_ode(), GrowthFactor(), hubble_function(), hybrid_nu_tracer(), hydro_force(), hydrokick_integ(), init(), init_cooling(), init_cooling_and_star_formation(), init_cooling_rates(), init_cosmology(), init_powerspectrum(), init_transfer_table(), kernel(), lightcone_compute(), lightcone_init(), lightcone_init_entry(), main(), metal_return(), metal_return_init(), OmegaFLD(), parse_transfer(), petaio_read_header_internal(), petaio_read_snapshot(), petaio_save_snapshot(), petaio_write_header(), print_spec(), read_parameterfile(), run_gravity_test(), saveheader(), setup_cosmology(), setup_density_indep_entropy(), setup_smoothinglengths(), test_cosmology(), test_DoCooling(), test_drift_factor(), test_fslength(), test_growth_numerical(), test_heatingcooling_rate(), test_rate_network(), test_read_no_rescale(), test_read_rescale_sigma8(), test_recomb_rates(), test_uvbg_loader(), tk_eh(), turn_on_quasars(), and write_checkpoint().

Norm

double Norm

static

Definition at line 22 of file power.c.

Referenced by DeltaSpec(), and init_powerspectrum().

power_table

struct table power_table

static

Definition at line 43 of file power.c.

Referenced by get_Tabulated(), init_powerspectrum(), init_transfer_table(), and save_all_transfer_tables().

PrimordialIndex

double PrimordialIndex

static

Definition at line 25 of file power.c.

Referenced by Delta_EH(), and init_powerspectrum().

tnames

const char* tnames[MAXCOLS] = {"DELTA_BAR", "DELTA_CDM", "DELTA_NU", "DELTA_CB", "VEL_BAR", "VEL_CDM", "VEL_NU", "VEL_CB", "VEL_TOT"}

static

Definition at line 50 of file power.c.

Referenced by save_all_transfer_tables().

transfer_table

struct table transfer_table

static

Definition at line 43 of file power.c.

Referenced by get_Tabulated(), init_powerspectrum(), init_transfer_table(), and save_all_transfer_tables().

UnitLength_in_cm

double UnitLength_in_cm

static

Definition at line 26 of file power.c.

Referenced by get_local_UVBG(), get_Tabulated(), glass_evolve(), gravpm_force(), init_neutrinos_lra(), init_powerspectrum(), lightcone_init(), lightcone_init_entry(), setup_glass(), test_DoCooling(), test_heatingcooling_rate(), and tk_eh().

WhichSpectrum

int WhichSpectrum

static

Definition at line 23 of file power.c.

Referenced by DeltaSpec(), and init_powerspectrum().