thermal.c File Reference

#include <gsl/gsl_integration.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_rng.h>
#include <assert.h>
#include "thermal.h"
#include <libgadget/physconst.h>
#include <libgadget/utils.h>

Include dependency graph for thermal.c:

Go to the source code of this file.

Macros
#define	BOLEVK   8.61734e-5

Functions
double	NU_V0 (const double Time, const double kBTNubyMNu, const double UnitVelocity_in_cm_per_s)
double	WDM_V0 (const double Time, const double WDM_therm_mass, const double Omega_CDM, const double HubbleParam, const double UnitVelocity_in_cm_per_s)
static double	fermi_dirac_kernel (double x, void *params)
double	init_thermalvel (struct thermalvel *thermals, const double v_amp, double max_fd, const double min_fd)
unsigned int *	init_rng (int Seed, int Nmesh)
void	add_thermal_speeds (struct thermalvel *thermals, gsl_rng *g_rng, float Vel[])

Macro Definition Documentation

BOLEVK

#define BOLEVK   8.61734e-5

Definition at line 11 of file thermal.c.

Function Documentation

add_thermal_speeds()

void add_thermal_speeds	(	struct thermalvel *	thermals,
		gsl_rng *	g_rng,
		float	Vel[]
	)

Definition at line 109 of file thermal.c.

{
    const double p = gsl_rng_uniform (g_rng);
    /*m_vamp multiples by the dimensional factor to get a velocity again.*/
    const double v = thermals->m_vamp * gsl_interp_eval(thermals->fd_intp,thermals->fermi_dirac_cumprob, thermals->fermi_dirac_vel, p, thermals->fd_intp_acc);
 
    /*Random phase*/
    const double phi = 2 * M_PI * gsl_rng_uniform (g_rng);
    const double theta = acos(2 * gsl_rng_uniform (g_rng) - 1);
 
    Vel[0] += v * sin(theta) * cos(phi);
    Vel[1] += v * sin(theta) * sin(phi);
    Vel[2] += v * cos(theta);
}

References thermalvel::fd_intp, thermalvel::fd_intp_acc, thermalvel::fermi_dirac_cumprob, thermalvel::fermi_dirac_vel, and thermalvel::m_vamp.

Referenced by main(), and test_thermal_vel().

Here is the caller graph for this function:

fermi_dirac_kernel()

static double fermi_dirac_kernel ( double  x, void *  params  )

static

Definition at line 40 of file thermal.c.

{
  return x * x / (exp(x) + 1);
}

Referenced by init_thermalvel().

Here is the caller graph for this function:

init_rng()

unsigned int* init_rng	(	int	Seed,
		int	Nmesh
	)

Definition at line 90 of file thermal.c.

{
    unsigned int * seedtable = (unsigned int *) mymalloc("randseeds", Nmesh*Nmesh*sizeof(unsigned int));
    gsl_rng * rng = gsl_rng_alloc(gsl_rng_ranlxd1);
    gsl_rng_set(rng, Seed);
 
    int i, j;
    for(i = 0; i < Nmesh; i++)
        for(j=0; j < Nmesh; j++)
        {
            seedtable[i+Nmesh*j] = gsl_rng_get(rng);
        }
    gsl_rng_free(rng);
    return seedtable;
}

References mymalloc.

Referenced by main().

Here is the caller graph for this function:

init_thermalvel()

double init_thermalvel	(	struct thermalvel *	thermals,
		const double	v_amp,
		double	max_fd,
		const double	min_fd
	)

Definition at line 47 of file thermal.c.

{
    int i;
    if(max_fd <= min_fd)
        endrun(1,"Thermal vel called with negative interval: %g <= %g\n", max_fd, min_fd);
 
    if(max_fd > MAX_FERMI_DIRAC)
        max_fd = MAX_FERMI_DIRAC;
    thermals->m_vamp = v_amp;
 
    /*These functions are so smooth that we don't need much space*/
    gsl_integration_workspace * w = gsl_integration_workspace_alloc (100);
    double abserr;
    gsl_function F;
    F.function = &fermi_dirac_kernel;
    F.params = NULL;
    for(i = 0; i < LENGTH_FERMI_DIRAC_TABLE; i++) {
        thermals->fermi_dirac_vel[i] = min_fd+(max_fd-min_fd)* i / (LENGTH_FERMI_DIRAC_TABLE - 1.0);
        gsl_integration_qag (&F, min_fd, thermals->fermi_dirac_vel[i], 0, 1e-6,100,GSL_INTEG_GAUSS61, w,&(thermals->fermi_dirac_cumprob[i]), &abserr);
    //       printf("gsl_integration_qng in fermi_dirac_init_nu. Result %g, error: %g, intervals: %lu\n",fermi_dirac_cumprob[i], abserr,w->size);
    }
    /*Save the largest cum. probability, pre-normalisation,
     * divided by the total F-D probability (which is 3 Zeta(3)/2 ~ 1.8 if MAX_FERMI_DIRAC is large enough*/
    double total_fd;
    gsl_integration_qag (&F, 0, MAX_FERMI_DIRAC, 0, 1e-6,100,GSL_INTEG_GAUSS61, w,&(total_fd), &abserr);
    assert(total_fd > 1.8);
 
    gsl_integration_workspace_free (w);
 
    double total_frac = thermals->fermi_dirac_cumprob[LENGTH_FERMI_DIRAC_TABLE-1]/total_fd;
    //Normalise total integral to unity
    for(i = 0; i < LENGTH_FERMI_DIRAC_TABLE; i++)
        thermals->fermi_dirac_cumprob[i] /= thermals->fermi_dirac_cumprob[LENGTH_FERMI_DIRAC_TABLE - 1];
 
    /*Initialise the GSL table*/
    thermals->fd_intp = gsl_interp_alloc(gsl_interp_cspline,LENGTH_FERMI_DIRAC_TABLE);
    thermals->fd_intp_acc = gsl_interp_accel_alloc();
    gsl_interp_init(thermals->fd_intp,thermals->fermi_dirac_cumprob, thermals->fermi_dirac_vel,LENGTH_FERMI_DIRAC_TABLE);
    return total_frac;
}

References endrun(), thermalvel::fd_intp, thermalvel::fd_intp_acc, thermalvel::fermi_dirac_cumprob, fermi_dirac_kernel(), thermalvel::fermi_dirac_vel, LENGTH_FERMI_DIRAC_TABLE, thermalvel::m_vamp, and MAX_FERMI_DIRAC.

Referenced by main(), and test_thermal_vel().

Here is the call graph for this function:

Here is the caller graph for this function:

NU_V0()

double NU_V0	(	const double	Time,
		const double	kBTNubyMNu,
		const double	UnitVelocity_in_cm_per_s
	)

Definition at line 24 of file thermal.c.

{
    return kBTNubyMNu / Time * (LIGHTCGS / UnitVelocity_in_cm_per_s);
}

References LIGHTCGS.

Referenced by main(), and test_mean_velocity().

Here is the caller graph for this function:

WDM_V0()

double WDM_V0	(	const double	Time,
		const double	WDM_therm_mass,
		const double	Omega_CDM,
		const double	HubbleParam,
		const double	UnitVelocity_in_cm_per_s
	)

Definition at line 30 of file thermal.c.

{
        //Not actually sure where this equation comes from: the fiducial values are from Bode, Ostriker & Turok 2001.
        double WDM_V0 = 0.012 / Time * pow(Omega_CDM / 0.3, 1.0 / 3) * pow(HubbleParam / 0.65, 2.0 / 3) * pow(1.0 /WDM_therm_mass,4.0 / 3);
        WDM_V0 *= 1.0e5 / UnitVelocity_in_cm_per_s;
        return WDM_V0;
}

Referenced by main().

Here is the caller graph for this function: