cooling_rates.c File Reference

#include "cooling_rates.h"
#include "cooling_qso_lightup.h"
#include "cosmology.h"
#include <omp.h>
#include <math.h>
#include <mpi.h>
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_interp.h>
#include "physconst.h"
#include "utils/endrun.h"
#include "utils/paramset.h"
#include "utils/mymalloc.h"

Include dependency graph for cooling_rates.c:

Go to the source code of this file.

Classes
struct	itp_type
struct	he_ions

Macros
#define	NGRAY   6
#define	NRECOMBTAB   1000
#define	RECOMBTMAX   log(1e9)
#define	RECOMBTMIN   0
#define	MAXITER   1000
#define	ITERCONV   1e-6

Functions
static void	init_itp_type (double *xarr, struct itp_type *Gamma, int Nelem)
static double	load_tree_value (char **saveptr)
static void	load_treecool (const char *TreeCoolFile)
static double	get_photo_rate (double redshift, struct itp_type *Gamma_tab)
static double	self_shield_dens (double redshift, const struct UVBG *uvbg)
struct UVBG	get_global_UVBG (double redshift)
static double	self_shield_corr (double nh, double logt, double ssdens)
static double	_Verner96Fit (double temp, double aa, double bb, double temp0, double temp1)
double	recomb_alphaHp (double temp)
static double	_Verner96alphaHep (double temp)
static double	recomb_alphaHep (double temp)
static double	recomb_alphad (double temp)
static double	recomb_alphaHepd (double temp)
static double	recomb_alphaHepp (double temp)
static double	_Voronov96Fit (double temp, double dE, double PP, double AA, double XX, double KK)
static double	recomb_GammaeH0 (double temp)
static double	recomb_GammaeHe0 (double temp)
static double	recomb_GammaeHep (double temp)
static double	get_interpolated_recomb (double logt, double *rec_tab, double rec_func(double))
static double	nH0_internal (double logt, double ne, const struct UVBG *uvbg, double photofac)
static double	nHp_internal (double nH0)
static struct he_ions	nHe_internal (double nh, double logt, double ne, const struct UVBG *uvbg, double photofac)
static double	get_temp_internal (double nebynh, double ienergy, double helium)
static double	ne_internal (double nh, double ienergy, double ne, double helium, double *logt, const struct UVBG *uvbg)
static double	scipy_optimize_fixed_point (double ne_init, double nh, double ienergy, double helium, double *logt, const struct UVBG *uvbg)
double	get_equilib_ne (double density, double ienergy, double helium, double *logt, const struct UVBG *uvbg, double ne_init)
double	get_ne_by_nh (double density, double ienergy, double helium, const struct UVBG *uvbg, double ne_init)
static double	_t5 (double temp)
static double	cool_CollisionalExciteH0 (double temp)
static double	cool_CollisionalExciteHeP (double temp)
static double	cool_CollisionalExciteHe0 (double temp)
static double	cool_CollisionalIonizeH0 (double temp)
static double	cool_CollisionalIonizeHe0 (double temp)
static double	cool_CollisionalIonizeHeP (double temp)
static double	cool_CollisionalH0 (double temp)
static double	cool_CollisionalHe0 (double temp)
static double	cool_CollisionalHeP (double temp)
static double	cool_RecombHp (double temp)
static double	cool_RecombDielect (double temp)
static double	cool_RecombHeP (double temp)
static double	cool_RecombHePP (double temp)
static double	cool_FreeFree (double temp, int zz)
static double	cool_FreeFree1 (double temp)
static double	cool_InverseCompton (double temp, double redshift)
static double	cool_he_reion_factor (double nHcgs, double helium, double redshift)
void	set_coolpar (struct cooling_params cp)
void	set_cooling_params (ParameterSet *ps)
void	init_cooling_rates (const char *TreeCoolFile, const char *MetalCoolFile, Cosmology *CP)
double	get_compton_cooling (double density, double ienergy, double helium, double redshift, double nebynh)
double	get_individual_cooling (enum CoolProcess process, double density, double ienergy, double helium, const struct UVBG *uvbg, double *ne_equilib)
double	get_heatingcooling_rate (double density, double ienergy, double helium, double redshift, double metallicity, const struct UVBG *uvbg, double *ne_equilib)
double	get_temp (double density, double ienergy, double helium, const struct UVBG *uvbg, double *ne_init)
double	get_neutral_fraction_phys_cgs (double density, double ienergy, double helium, const struct UVBG *uvbg, double *ne_init)
double	get_helium_ion_phys_cgs (int ion, double density, double ienergy, double helium, const struct UVBG *uvbg, double ne_init)

Variables
static struct cooling_params	CoolingParams
static gsl_interp *	GrayOpac
static const double	GrayOpac_ydata [NGRAY] = { 2.59e-18, 2.37e-18, 2.27e-18, 2.15e-18, 2.02e-18, 1.94e-18}
static const double	GrayOpac_zz [NGRAY] = {0, 1, 2, 3, 4, 5}
static int	NTreeCool
static double *	Gamma_log1z
static struct itp_type Gamma_HI Gamma_HeI	Gamma_HeII
static struct itp_type Eps_HI Eps_HeI	Eps_HeII
static double *	temp_tab
static double *	rec_alphaHp
static double *	rec_alphaHep
static double *	rec_alphaHepp
static double *	rec_GammaH0
static double *	rec_GammaHe0
static double *	rec_GammaHep
static double *	cool_collisH0
static double *	cool_collisHe0
static double *	cool_collisHeP
static double *	cool_recombHp
static double *	cool_recombHeP
static double *	cool_recombHePP
static double *	cool_freefree1

Macro Definition Documentation

ITERCONV

#define ITERCONV   1e-6

Definition at line 636 of file cooling_rates.c.

MAXITER

#define MAXITER   1000

Definition at line 633 of file cooling_rates.c.

NGRAY

#define NGRAY   6

Definition at line 73 of file cooling_rates.c.

NRECOMBTAB

#define NRECOMBTAB   1000

Definition at line 95 of file cooling_rates.c.

RECOMBTMAX

#define RECOMBTMAX   log(1e9)

Definition at line 96 of file cooling_rates.c.

RECOMBTMIN

#define RECOMBTMIN   0

Definition at line 97 of file cooling_rates.c.

Function Documentation

_t5()

static double _t5 ( double  temp )

static

Definition at line 746 of file cooling_rates.c.

{
    double t0;
    if(CoolingParams.cooling == KWH92)
        t0 = 1e5;
    /* The original Cen 92 correction is implemented in order to achieve the
     * correct asymptotic behaviour at high temperatures. However, the rates he is correcting
     * should be good up to 10^7 K, so he is being too aggressive.*/
    else
        t0 = 5e7;
    return 1+sqrt(temp/t0);
}

References cooling_params::cooling, CoolingParams, and KWH92.

Referenced by cool_CollisionalExciteH0(), cool_CollisionalExciteHe0(), and cool_CollisionalExciteHeP().

Here is the caller graph for this function:

_Verner96alphaHep()

static double _Verner96alphaHep ( double  temp )

static

Definition at line 367 of file cooling_rates.c.

{
    /*Accurate to ~2% for T < 10^6 and 5% for T< 10^10.
     * VF96 give two rates. The first is more accurate for T < 10^6, the second is valid up to T = 10^10.
     * We use the most accurate allowed. See lines 2 and 3 of Table 1 of VF96.*/
    double lowTfit = _Verner96Fit(temp, 3.294e-11, 0.6910, 1.554e+01, 3.676e+07);
    double highTfit = _Verner96Fit(temp, 9.356e-10, 0.7892, 4.266e-02, 4.677e+06);
    /*Note that at 10^6K the two fits differ by ~10%. This may lead one to disbelieve the quoted accuracies!
    * We thus switch over at a slightly lower temperature. The two fits cross at T ~ 3e5K.*/
    double swtmp = 7e5;
    double deltat = 1e5;
    double upper = swtmp + deltat;
    double lower = swtmp - deltat;
    //In order to avoid a sharp feature at 10^6 K, we linearly interpolate between the two fits around 10^6 K.
    double interpfit = (lowTfit * (upper - temp) + highTfit * (temp - lower))/(2*deltat);
    return (temp < lower)*lowTfit + (temp > upper)*highTfit + (upper > temp)*(temp > lower)*interpfit;
}

References _Verner96Fit().

Referenced by recomb_alphaHep().

Here is the call graph for this function:

Here is the caller graph for this function:

_Verner96Fit()

static double _Verner96Fit ( double  temp, double  aa, double  bb, double  temp0, double  temp1  )

static

Definition at line 338 of file cooling_rates.c.

{
        double sqrttt0 = sqrt(temp/temp0);
        double sqrttt1 = sqrt(temp/temp1);
        return aa / ( sqrttt0 * pow(1 + sqrttt0, 1-bb)*pow(1+sqrttt1, 1+bb) );
}

Referenced by _Verner96alphaHep(), recomb_alphaHep(), recomb_alphaHepp(), and recomb_alphaHp().

Here is the caller graph for this function:

_Voronov96Fit()

static double _Voronov96Fit ( double  temp, double  dE, double  PP, double  AA, double  XX, double  KK  )

static

Definition at line 451 of file cooling_rates.c.

{
    double UU = dE / (BOLEVK * temp);
    return AA * (1 + PP * sqrt(UU))/(XX+UU) * pow(UU, KK) * exp(-UU);
}

References BOLEVK.

Referenced by recomb_GammaeH0(), recomb_GammaeHe0(), and recomb_GammaeHep().

Here is the caller graph for this function:

cool_CollisionalExciteH0()

static double cool_CollisionalExciteH0 ( double  temp )

static

Definition at line 761 of file cooling_rates.c.

{
    return 7.5e-19 * exp(-118348.0/temp) / _t5(temp);
}

References _t5().

Referenced by cool_CollisionalH0().

Here is the call graph for this function:

Here is the caller graph for this function:

cool_CollisionalExciteHe0()

static double cool_CollisionalExciteHe0 ( double  temp )

static

Definition at line 775 of file cooling_rates.c.

{
    return 9.1e-27 * pow(temp, -0.1687) * exp(-473638/temp) / _t5(temp);
}

References _t5().

Referenced by cool_CollisionalHe0().

Here is the call graph for this function:

Here is the caller graph for this function:

cool_CollisionalExciteHeP()

static double cool_CollisionalExciteHeP ( double  temp )

static

Definition at line 768 of file cooling_rates.c.

{
    return 5.54e-17 * pow(temp, -0.397)*exp(-473638./temp)/_t5(temp);
}

References _t5().

Referenced by cool_CollisionalHeP().

Here is the call graph for this function:

Here is the caller graph for this function:

cool_CollisionalH0()

static double cool_CollisionalH0 ( double  temp )

static

Definition at line 804 of file cooling_rates.c.

{
    if(CoolingParams.cooling == Enzo2Nyx) {
        /*Formula from Eq. 23, Table 4 of Scholz & Walters, claimed good to 0.45 %.
        Note though that they have two datasets which differ by a factor of two.
        Differs from Cen 92 by a factor of two.*/
        //Technically only good for T > 2000.
        double y = log(temp);
        //Constant is 0.75/k_B in Rydberg
        double Ryd = 2.1798741e-11;
        double tot = -0.75/BOLTZMANN *Ryd/temp;
        double coeffslowT[6] = {213.7913, 113.9492, 25.06062, 2.762755, 0.1515352, 3.290382e-3};
        double coeffshighT[6] = {271.25446, 98.019455, 14.00728, 0.9780842, 3.356289e-2, 4.553323e-4};
        int j;
        for(j=0; j < 6; j++)
            tot += ((temp < 1e5)*coeffslowT[j]+(temp >=1e5)*coeffshighT[j])*pow(-y, j);
        return 1e-20 * exp(tot);
    }
    else /*Everyone else splits collisional from excitation, because collisional is mostly exact*/
        return cool_CollisionalExciteH0(temp) + cool_CollisionalIonizeH0(temp);
}

References BOLTZMANN, cool_CollisionalExciteH0(), cool_CollisionalIonizeH0(), cooling_params::cooling, CoolingParams, and Enzo2Nyx.

Referenced by get_heatingcooling_rate(), get_individual_cooling(), and init_cooling_rates().

Here is the call graph for this function:

Here is the caller graph for this function:

cool_CollisionalHe0()

static double cool_CollisionalHe0 ( double  temp )

static

Definition at line 828 of file cooling_rates.c.

{
    return cool_CollisionalExciteHe0(temp) + cool_CollisionalIonizeHe0(temp);
}

References cool_CollisionalExciteHe0(), and cool_CollisionalIonizeHe0().

Referenced by get_heatingcooling_rate(), get_individual_cooling(), and init_cooling_rates().

Here is the call graph for this function:

Here is the caller graph for this function:

cool_CollisionalHeP()

static double cool_CollisionalHeP ( double  temp )

static

Definition at line 835 of file cooling_rates.c.

{
    return cool_CollisionalExciteHeP(temp) + cool_CollisionalIonizeHeP(temp);
}

References cool_CollisionalExciteHeP(), and cool_CollisionalIonizeHeP().

Referenced by get_heatingcooling_rate(), get_individual_cooling(), and init_cooling_rates().

Here is the call graph for this function:

Here is the caller graph for this function:

cool_CollisionalIonizeH0()

static double cool_CollisionalIonizeH0 ( double  temp )

static

Definition at line 782 of file cooling_rates.c.

{
    /*H ionization potential*/
    return 13.5984 * eVinergs * recomb_GammaeH0(temp);
}

References eVinergs, and recomb_GammaeH0().

Referenced by cool_CollisionalH0().

Here is the call graph for this function:

Here is the caller graph for this function:

cool_CollisionalIonizeHe0()

static double cool_CollisionalIonizeHe0 ( double  temp )

static

Definition at line 790 of file cooling_rates.c.

{
    return 24.5874 * eVinergs * recomb_GammaeHe0(temp);
}

References eVinergs, and recomb_GammaeHe0().

Referenced by cool_CollisionalHe0().

Here is the call graph for this function:

Here is the caller graph for this function:

cool_CollisionalIonizeHeP()

static double cool_CollisionalIonizeHeP ( double  temp )

static

Definition at line 797 of file cooling_rates.c.

{
    return 54.417760 * eVinergs * recomb_GammaeHep(temp);
}

References eVinergs, and recomb_GammaeHep().

Referenced by cool_CollisionalHeP().

Here is the call graph for this function:

Here is the caller graph for this function:

cool_FreeFree()

static double cool_FreeFree ( double  temp, int  zz  )

static

Definition at line 881 of file cooling_rates.c.

{
    double gff;
    /*Formula for the Gaunt factor from Shapiro & Kang 1987. ZZ is 1 for H+ and He+ and 2 for He++.
      This is almost identical to the KWH rate but not continuous.*/
    if(CoolingParams.cooling == Enzo2Nyx) {
        double lt = 2 * log10(temp/zz);
        if(lt <= log10(3.2e5))
            gff = (0.79464 + 0.1243*lt);
        else
            gff = (2.13164 - 0.1240*lt);
    }
    else {
        /*Formula for the Gaunt factor. KWH takes this from Spitzer 1978.*/
        gff = 1.1+0.34*exp(-pow(5.5 - log10(temp),2) /3.);
    }
    return 1.426e-27*sqrt(temp)* pow(zz,2) * gff;
}

References cooling_params::cooling, CoolingParams, and Enzo2Nyx.

Referenced by cool_FreeFree1(), get_heatingcooling_rate(), get_individual_cooling(), and init_cooling_rates().

Here is the caller graph for this function:

cool_FreeFree1()

static double cool_FreeFree1 ( double  temp )

static

Definition at line 901 of file cooling_rates.c.

{
    return cool_FreeFree(temp, 1);
}

References cool_FreeFree().

Referenced by get_heatingcooling_rate(), and get_individual_cooling().

Here is the call graph for this function:

Here is the caller graph for this function:

cool_he_reion_factor()

static double cool_he_reion_factor ( double  nHcgs, double  helium, double  redshift  )

static

Definition at line 924 of file cooling_rates.c.

{
  if(!CoolingParams.HeliumHeatOn)
      return 1.;
  const double rho = PROTONMASS * nHcgs / (1 - helium);
  double overden = rho/(CoolingParams.rho_crit_baryon * pow(1+redshift,3.0));
  if (overden >= CoolingParams.HeliumHeatThresh)
      overden = CoolingParams.HeliumHeatThresh;
  return CoolingParams.HeliumHeatAmp*pow(overden, CoolingParams.HeliumHeatExp);
}

References CoolingParams, cooling_params::HeliumHeatAmp, cooling_params::HeliumHeatExp, cooling_params::HeliumHeatOn, cooling_params::HeliumHeatThresh, PROTONMASS, and cooling_params::rho_crit_baryon.

Referenced by get_heatingcooling_rate().

Here is the caller graph for this function:

cool_InverseCompton()

static double cool_InverseCompton ( double  temp, double  redshift  )

static

Definition at line 910 of file cooling_rates.c.

{
        double tcmb_red = CoolingParams.CMBTemperature * (1+redshift);
        return 4 * THOMPSON * RAD_CONST / (ELECTRONMASS * LIGHTCGS ) * pow(tcmb_red, 4) * BOLTZMANN * (temp - tcmb_red);
}

References BOLTZMANN, cooling_params::CMBTemperature, CoolingParams, ELECTRONMASS, LIGHTCGS, RAD_CONST, and THOMPSON.

Referenced by get_compton_cooling(), and get_heatingcooling_rate().

Here is the caller graph for this function:

cool_RecombDielect()

static double cool_RecombDielect ( double  temp )

static

Definition at line 853 of file cooling_rates.c.

{
    /*What is this magic number?*/
    return 6.526e-11* recomb_alphad(temp);
}

References recomb_alphad().

Referenced by cool_RecombHeP().

Here is the call graph for this function:

Here is the caller graph for this function:

cool_RecombHeP()

static double cool_RecombHeP ( double  temp )

static

Definition at line 861 of file cooling_rates.c.

{
    return 0.75 * BOLTZMANN * temp * recomb_alphaHep(temp) + cool_RecombDielect(temp);
}

References BOLTZMANN, cool_RecombDielect(), and recomb_alphaHep().

Referenced by get_heatingcooling_rate(), get_individual_cooling(), and init_cooling_rates().

Here is the call graph for this function:

Here is the caller graph for this function:

cool_RecombHePP()

static double cool_RecombHePP ( double  temp )

static

Definition at line 868 of file cooling_rates.c.

{
    if(CoolingParams.cooling == Enzo2Nyx) {
        /*Recombination cooling rate from Black 1981: these increase rapidly for T > 5e5K. */
        return 1.140e-26 * sqrt(temp) * (6.607 - 0.5 * log(temp) + 7.459e-3 * pow(temp, 1./3));
    }
    return 0.75 * BOLTZMANN * temp * recomb_alphaHepp(temp);
}

References BOLTZMANN, cooling_params::cooling, CoolingParams, Enzo2Nyx, and recomb_alphaHepp().

Referenced by get_heatingcooling_rate(), get_individual_cooling(), and init_cooling_rates().

Here is the call graph for this function:

Here is the caller graph for this function:

cool_RecombHp()

static double cool_RecombHp ( double  temp )

static

Definition at line 842 of file cooling_rates.c.

{
    if(CoolingParams.cooling == Enzo2Nyx) {
        /*Recombination cooling rate from Black 1981: these increase rapidly for T > 5e5K. */
        return 2.851e-27 * sqrt(temp) * (5.914 - 0.5 * log(temp) + 0.01184 * pow(temp, 1./3));
    }
    return 0.75 * BOLTZMANN * temp * recomb_alphaHp(temp);
}

References BOLTZMANN, cooling_params::cooling, CoolingParams, Enzo2Nyx, and recomb_alphaHp().

Referenced by get_heatingcooling_rate(), get_individual_cooling(), and init_cooling_rates().

Here is the call graph for this function:

Here is the caller graph for this function:

get_compton_cooling()

double get_compton_cooling	(	double	density,
		double	ienergy,
		double	helium,
		double	redshift,
		double	nebynh
	)

Definition at line 1034 of file cooling_rates.c.

{
    double nh = density * (1 - helium);
    double temp = get_temp_internal(nebynh, ienergy, helium);
    /*Compton cooling in erg/s cm^3*/
    double LambdaCmptn = -1*nebynh * cool_InverseCompton(temp, redshift) / nh;
    return LambdaCmptn * pow(1 - helium, 2) * density / PROTONMASS;
}

References cool_InverseCompton(), density(), get_temp_internal(), and PROTONMASS.

Here is the call graph for this function:

get_equilib_ne()

double get_equilib_ne	(	double	density,
		double	ienergy,
		double	helium,
		double *	logt,
		const struct UVBG *	uvbg,
		double	ne_init
	)

Definition at line 683 of file cooling_rates.c.

{
    /*Get hydrogen number density*/
    double nh = density * (1-helium);
    /* Avoid getting stuck in an alternate solution
     * where there is never any heating in the presence of roundoff.*/
    if(ne_init <= 0)
        ne_init = 1.0;
    return scipy_optimize_fixed_point(ne_init, nh, ienergy, helium, logt, uvbg);
}

References density(), and scipy_optimize_fixed_point().

Referenced by get_heatingcooling_rate(), get_helium_ion_phys_cgs(), get_individual_cooling(), get_ne_by_nh(), get_neutral_fraction_phys_cgs(), get_temp(), and test_rate_network().

Here is the call graph for this function:

Here is the caller graph for this function:

get_global_UVBG()

struct UVBG get_global_UVBG

(

double

redshift

)

Definition at line 244 of file cooling_rates.c.

{
    struct UVBG GlobalUVBG = {0};
 
    if(!CoolingParams.PhotoIonizationOn)
        return GlobalUVBG;
 
    /* Set the homogeneous reionization redshift as the point when the UVB switches on.*/
    GlobalUVBG.zreion = pow(10,Gamma_log1z[NTreeCool - 1])-1;
 
    if(CoolingParams.UVRedshiftThreshold >= 0.)
        GlobalUVBG.zreion = CoolingParams.UVRedshiftThreshold;
 
    /* if a threshold is set, disable UV bg above that redshift */
    if(CoolingParams.UVRedshiftThreshold >= 0. && redshift > CoolingParams.UVRedshiftThreshold)
        return GlobalUVBG;
 
 
    GlobalUVBG.gJH0 = get_photo_rate(redshift, &Gamma_HI);
    GlobalUVBG.gJHe0 = get_photo_rate(redshift, &Gamma_HeI);
    GlobalUVBG.gJHep = get_photo_rate(redshift, &Gamma_HeII);
 
    GlobalUVBG.epsH0 = get_photo_rate(redshift, &Eps_HI);
    GlobalUVBG.epsHe0 = get_photo_rate(redshift, &Eps_HeI);
    /* During helium reionization we have a model for the inhomogeneous non-equilibrium heating.
     * To avoid double counting, remove the heating in the existing UVB*/
    if(during_helium_reionization(redshift))
        GlobalUVBG.epsHep = 0;
    else
        GlobalUVBG.epsHep = get_photo_rate(redshift, &Eps_HeII);
    GlobalUVBG.self_shield_dens = self_shield_dens(redshift, &GlobalUVBG);
    return GlobalUVBG;
}

References CoolingParams, cooling_params::fBar, UVBG::gJH0, GrayOpac, GrayOpac_ydata, GrayOpac_zz, and NGRAY.

Referenced by cooling_and_starformation(), get_helium_neutral_fraction_sfreff(), get_neutral_fraction_sfreff(), sfreff_on_eeqos(), test_DoCooling(), test_heatingcooling_rate(), test_rate_network(), and test_uvbg_loader().

Here is the caller graph for this function:

get_heatingcooling_rate()

double get_heatingcooling_rate	(	double	density,
		double	ienergy,
		double	helium,
		double	redshift,
		double	metallicity,
		const struct UVBG *	uvbg,
		double *	ne_equilib
	)

Definition at line 1104 of file cooling_rates.c.

{
    double logt;
    double ne = get_equilib_ne(density, ienergy, helium, &logt, uvbg, *ne_equilib);
    double nh = density * (1 - helium);
    double nebynh = ne/nh;
    /*Faster than running the exp.*/
    double temp = get_temp_internal(nebynh, ienergy, helium);
    double photofac = self_shield_corr(nh, logt, uvbg->self_shield_dens);
 
    /*The helium number fraction*/
    double yy = helium / 4 / (1 - helium);
 
    double nH0 = nH0_internal(logt, ne, uvbg, photofac);
    double nHp = nHp_internal(nH0);
    struct he_ions He = nHe_internal(nh, logt, ne, uvbg, photofac);
    /*Put the abundances in units of nH to avoid underflows*/
    He.nHep*= yy/nh;
    He.nHe0*= yy/nh;
    He.nHepp*= yy/nh;
    /*Collisional ionization and excitation rate*/
    double LambdaCollis = nebynh * (get_interpolated_recomb(logt, cool_collisH0, cool_CollisionalH0) * nH0 +
            get_interpolated_recomb(logt, cool_collisHe0, cool_CollisionalHe0) * He.nHe0 +
            get_interpolated_recomb(logt, cool_collisHeP, cool_CollisionalHeP) * He.nHep);
    double LambdaRecomb = nebynh * (get_interpolated_recomb(logt, cool_recombHp, cool_RecombHp) * nHp +
            get_interpolated_recomb(logt, cool_recombHeP, cool_RecombHeP) * He.nHep +
            get_interpolated_recomb(logt, cool_recombHePP, cool_RecombHePP) * He.nHepp);
    /*Free-free cooling rate*/
    double LambdaFF = 0;
 
    double cff = get_interpolated_recomb(logt, cool_freefree1, cool_FreeFree1);
 
    if(CoolingParams.cooling == Enzo2Nyx) {
        LambdaFF = nebynh * (cff * (nHp + He.nHep) + cool_FreeFree(temp, 2) * He.nHepp);
    } else {
        /*The factor of (zz=2)^2 has been pulled out, so if we use the Spitzer gaunt factor we don't need
         * to call the FreeFree function again.*/
        LambdaFF = nebynh * (cff * (nHp + He.nHep) + 4 * cff * He.nHepp);
    }
    /*Compton cooling in erg/s cm^3*/
    double LambdaCmptn = nebynh * cool_InverseCompton(temp, redshift) / nh;
    /*Total cooling rate per proton in erg/s cm^3*/
    double Lambda = LambdaCollis + LambdaRecomb + LambdaFF + LambdaCmptn;
 
    /*Total heating rate per proton in erg/s cm^3*/
    double Heat = (nH0 * uvbg->epsH0 + He.nHe0 * uvbg->epsHe0 + He.nHep * uvbg->epsHep)/nh;
 
    Heat *= cool_he_reion_factor(density, helium, redshift);
    /*Set external equilibrium electron density*/
    *ne_equilib = nebynh;
 
    /*Apply metal cooling. Does nothing if metal cooling is disabled*/
    double MetalCooling = metallicity * TableMetalCoolingRate(redshift, temp, nh);
 
    double LambdaNet = Heat - Lambda - MetalCooling;
 
    //message(1, "Heat = %g Lambda = %g MetalCool = %g LC = %g LR = %g LFF = %g LCmptn = %g, ne = %g, nH0 = %g, nHp = %g, nHe0 = %g, nHep = %g, nHepp = %g, nh=%g, temp=%g, ienergy=%g\n", Heat, Lambda, MetalCooling, LambdaCollis, LambdaRecomb, LambdaFF, LambdaCmptn, nebynh, nH0, nHp, nHe0, nHep, nHepp, nh, temp, ienergy);
 
    /* LambdaNet in erg cm^3 /s, Density in protons/cm^3, PROTONMASS in protons/g.
     * Convert to erg/s/g*/
    return LambdaNet * pow(1 - helium, 2) * density / PROTONMASS;
}

References cool_collisH0, cool_collisHe0, cool_collisHeP, cool_CollisionalH0(), cool_CollisionalHe0(), cool_CollisionalHeP(), cool_FreeFree(), cool_freefree1, cool_FreeFree1(), cool_he_reion_factor(), cool_InverseCompton(), cool_recombHeP, cool_RecombHeP(), cool_recombHePP, cool_RecombHePP(), cool_recombHp, cool_RecombHp(), cooling_params::cooling, CoolingParams, density(), Enzo2Nyx, UVBG::epsH0, UVBG::epsHe0, UVBG::epsHep, get_equilib_ne(), get_interpolated_recomb(), get_temp_internal(), nH0_internal(), he_ions::nHe0, nHe_internal(), he_ions::nHep, he_ions::nHepp, nHp_internal(), PROTONMASS, self_shield_corr(), UVBG::self_shield_dens, and TableMetalCoolingRate().

Referenced by get_lambdanet(), GetCoolingTime(), and test_heatingcooling_rate().

Here is the call graph for this function:

Here is the caller graph for this function:

get_helium_ion_phys_cgs()

double get_helium_ion_phys_cgs	(	int	ion,
		double	density,
		double	ienergy,
		double	helium,
		const struct UVBG *	uvbg,
		double	ne_init
	)

Definition at line 1202 of file cooling_rates.c.

{
    double logt;
    double ne = get_equilib_ne(density, ienergy, helium, &logt, uvbg, ne_init);
    double yy = helium / 4 / (1 - helium);
    double nh = density * (1-helium);
    double photofac = self_shield_corr(nh, logt, uvbg->self_shield_dens);
    struct he_ions He = nHe_internal(nh, logt, ne, uvbg, photofac);
    if(ion == 0)
        return yy * He.nHe0 / nh;
    else if (ion == 1)
        return yy * He.nHep / nh;
    else
        return yy * He.nHepp / nh;
}

References density(), get_equilib_ne(), he_ions::nHe0, nHe_internal(), he_ions::nHep, he_ions::nHepp, self_shield_corr(), and UVBG::self_shield_dens.

Referenced by GetHeliumIonFraction().

Here is the call graph for this function:

Here is the caller graph for this function:

get_individual_cooling()

double get_individual_cooling	(	enum CoolProcess	process,
		double	density,
		double	ienergy,
		double	helium,
		const struct UVBG *	uvbg,
		double *	ne_equilib
	)

Definition at line 1046 of file cooling_rates.c.

{
    double logt;
    double ne = get_equilib_ne(density, ienergy, helium, &logt, uvbg, *ne_equilib);
    double nh = density * (1 - helium);
    double nebynh = ne/nh;
    /*Faster than running the exp.*/
    double temp = get_temp_internal(nebynh, ienergy, helium);
    double photofac = self_shield_corr(nh, logt, uvbg->self_shield_dens);
    double nH0 = nH0_internal(logt, ne, uvbg, photofac);
    double nHp = nHp_internal(nH0);
    /*The helium number fraction*/
    double yy = helium / 4 / (1 - helium);
    struct he_ions He = nHe_internal(nh, logt, ne, uvbg, photofac);
    /*Put the abundances in units of nH to avoid underflows*/
    He.nHep*= yy/nh;
    He.nHe0*= yy/nh;
    He.nHepp*= yy/nh;
 
    /*Compton cooling in erg/s cm^3*/
    double Lambda = 0;
 
    if(process == FREEFREE) {
        double cff = get_interpolated_recomb(logt, cool_freefree1, cool_FreeFree1);
        if(CoolingParams.cooling == Enzo2Nyx) {
            Lambda = -1*nebynh * (cff * (nHp + He.nHep) + cool_FreeFree(temp, 2) * He.nHepp);
        } else {
            /*The factor of (zz=2)^2 has been pulled out, so if we use the Spitzer gaunt factor we don't need
            * to call the FreeFree function again.*/
            Lambda = -1*nebynh * (cff * (nHp + He.nHep) + 4 * cff * He.nHepp);
        }
    } else if(process == HEAT) {
            /*Total heating rate per proton in erg/s cm^3*/
            Lambda = (nH0 * uvbg->epsH0 + He.nHe0 * uvbg->epsHe0 + He.nHep * uvbg->epsHep)/nh;
    }
    else if(process == RECOMB) {
        Lambda = -1*nebynh * (get_interpolated_recomb(logt, cool_recombHp, cool_RecombHp) * nHp +
            get_interpolated_recomb(logt, cool_recombHeP, cool_RecombHeP) * He.nHep +
            get_interpolated_recomb(logt, cool_recombHePP, cool_RecombHePP) * He.nHepp);
    }
    else if(process == COLLIS) {
        Lambda = -1*nebynh * (get_interpolated_recomb(logt, cool_collisH0, cool_CollisionalH0) * nH0 +
            get_interpolated_recomb(logt, cool_collisHe0, cool_CollisionalHe0) * He.nHe0 +
            get_interpolated_recomb(logt, cool_collisHeP, cool_CollisionalHeP) * He.nHep);
    }
 
    return Lambda * pow(1 - helium, 2) * density / PROTONMASS;
}

References COLLIS, cool_collisH0, cool_collisHe0, cool_collisHeP, cool_CollisionalH0(), cool_CollisionalHe0(), cool_CollisionalHeP(), cool_FreeFree(), cool_freefree1, cool_FreeFree1(), cool_recombHeP, cool_RecombHeP(), cool_recombHePP, cool_RecombHePP(), cool_recombHp, cool_RecombHp(), cooling_params::cooling, CoolingParams, density(), Enzo2Nyx, UVBG::epsH0, UVBG::epsHe0, UVBG::epsHep, FREEFREE, get_equilib_ne(), get_interpolated_recomb(), get_temp_internal(), HEAT, nH0_internal(), he_ions::nHe0, nHe_internal(), he_ions::nHep, he_ions::nHepp, nHp_internal(), PROTONMASS, RECOMB, self_shield_corr(), and UVBG::self_shield_dens.

Here is the call graph for this function:

get_interpolated_recomb()

static double get_interpolated_recomb ( double  logt, double *  rec_tab, double   rec_funcdouble  )

static

Definition at line 510 of file cooling_rates.c.

{
    /*Find the index to use in our temperature table.*/
    double dind = (logt - RECOMBTMIN) / (RECOMBTMAX - RECOMBTMIN) * NRECOMBTAB;
    int index = (int) dind;
    /*Just call the function directly if we are out of interpolation range*/
    if(index < 0 || index >= NRECOMBTAB-1)
        return rec_func(exp(logt));
    //if (temp_tab[index] > logt || temp_tab[index+1] < logt || index < 0 || index >= NRECOMBTAB)
    //    endrun(2, "Incorrect indexing of recombination array\n");
    return rec_tab[index + 1] * (dind - index) + rec_tab[index] * (1 - (dind - index));
}

References NRECOMBTAB, RECOMBTMAX, and RECOMBTMIN.

Referenced by get_heatingcooling_rate(), get_individual_cooling(), and nH0_internal().

Here is the caller graph for this function:

get_ne_by_nh()

double get_ne_by_nh	(	double	density,
		double	ienergy,
		double	helium,
		const struct UVBG *	uvbg,
		double	ne_init
	)

Definition at line 696 of file cooling_rates.c.

{
    double logt;
    return get_equilib_ne(density, ienergy, helium, &logt, uvbg, ne_init)/(density*(1-helium));
}

References density(), and get_equilib_ne().

Here is the call graph for this function:

get_neutral_fraction_phys_cgs()

double get_neutral_fraction_phys_cgs	(	double	density,
		double	ienergy,
		double	helium,
		const struct UVBG *	uvbg,
		double *	ne_init
	)

Definition at line 1186 of file cooling_rates.c.

{
    double logt;
    double ne = get_equilib_ne(density, ienergy, helium, &logt, uvbg, *ne_init);
    double nh = density * (1-helium);
    double photofac = self_shield_corr(nh, logt, uvbg->self_shield_dens);
    *ne_init = ne/nh;
    return nH0_internal(logt, ne, uvbg, photofac);
}

References density(), get_equilib_ne(), nH0_internal(), self_shield_corr(), and UVBG::self_shield_dens.

Referenced by GetNeutralFraction(), and test_rate_network().

Here is the call graph for this function:

Here is the caller graph for this function:

get_photo_rate()

static double get_photo_rate ( double  redshift, struct itp_type *  Gamma_tab  )

static

Definition at line 215 of file cooling_rates.c.

{
    if(!CoolingParams.PhotoIonizationOn)
        return 0;
    double log1z = log10(1+redshift);
    double photo_rate;
    if (NTreeCool < 1 || log1z >= Gamma_log1z[NTreeCool - 1])
        return 0;
    else if (log1z < Gamma_log1z[0])
        photo_rate = Gamma_tab->ydata[0];
    else {
        photo_rate = gsl_interp_eval(Gamma_tab->intp, Gamma_log1z, Gamma_tab->ydata, log1z, NULL);
    }
    return pow(10, photo_rate) * CoolingParams.PhotoIonizeFactor;
}

References CoolingParams, Gamma_log1z, itp_type::intp, NTreeCool, cooling_params::PhotoIonizationOn, cooling_params::PhotoIonizeFactor, and itp_type::ydata.

get_temp()

double get_temp	(	double	density,
		double	ienergy,
		double	helium,
		const struct UVBG *	uvbg,
		double *	ne_init
	)

Definition at line 1172 of file cooling_rates.c.

{
    double logt;
    double ne = get_equilib_ne(density, ienergy, helium, &logt, uvbg, *ne_init);
    double nh = density * (1 - helium);
    *ne_init = ne/nh;
    return get_temp_internal(ne/nh, ienergy, helium);
}

References density(), get_equilib_ne(), and get_temp_internal().

Referenced by test_rate_network().

Here is the call graph for this function:

Here is the caller graph for this function:

get_temp_internal()

static double get_temp_internal ( double  nebynh, double  ienergy, double  helium  )

static

Definition at line 601 of file cooling_rates.c.

{
    /*convert U (erg/g) to T (K) : U = N k T / (γ - 1)
    T = U (γ-1) μ m_P / k_B
    where k_B is the Boltzmann constant
    γ is 5/3, the perfect gas constant
    m_P is the proton mass
    μ is 1 / (mean no. molecules per unit atomic weight) calculated in loop.
    Internal energy units are 10^-10 erg/g*/
    double hy_mass = 1 - helium;
    double muienergy = 4 / (hy_mass * (3 + 4*nebynh) + 1)*ienergy;
    /*So for T in K, Boltzmann in erg/K, internal energy has units of erg/g*/
    double temp = GAMMA_MINUS1 * PROTONMASS / BOLTZMANN * muienergy;
    if(temp < CoolingParams.MinGasTemp)
        return CoolingParams.MinGasTemp;
    return temp;
}

References BOLTZMANN, CoolingParams, GAMMA_MINUS1, cooling_params::MinGasTemp, and PROTONMASS.

Referenced by get_compton_cooling(), get_heatingcooling_rate(), get_individual_cooling(), get_temp(), ne_internal(), and scipy_optimize_fixed_point().

Here is the caller graph for this function:

init_cooling_rates()

void init_cooling_rates	(	const char *	TreeCoolFile,
		const char *	MetalCoolFile,
		Cosmology *	CP
	)

Definition at line 970 of file cooling_rates.c.

{
    CoolingParams.fBar = CP->OmegaBaryon / CP->OmegaCDM;
    CoolingParams.rho_crit_baryon = CP->OmegaBaryon * 3.0 * pow(CP->HubbleParam*HUBBLE,2.0) /(8.0*M_PI*GRAVITY);
 
    /* Initialize the interpolation for the self-shielding module as a function of redshift.
     * A crash has been observed in GSL with a cspline interpolator. */
    GrayOpac = gsl_interp_alloc(gsl_interp_linear,NGRAY);
    gsl_interp_init(GrayOpac,GrayOpac_zz,GrayOpac_ydata, NGRAY);
 
    if(!TreeCoolFile || strnlen(TreeCoolFile,100) == 0) {
        CoolingParams.PhotoIonizationOn = 0;
        message(0, "No TreeCool file is provided. Cooling is broken. OK for DM only runs. \n");
    }
    else {
        message(0, "Using uniform UVB from file %s\n", TreeCoolFile);
        /* Load the TREECOOL into Gamma_HI->ydata, and initialise the interpolators*/
        load_treecool(TreeCoolFile);
    }
 
    /*Initialize the recombination tables*/
    temp_tab = (double *) mymalloc("Recombination_tables", NRECOMBTAB * sizeof(double) * 14);
 
    rec_GammaH0 = temp_tab + NRECOMBTAB;
    rec_GammaHe0 = temp_tab + 2 * NRECOMBTAB;
    rec_GammaHep = temp_tab + 3 * NRECOMBTAB;
    rec_alphaHp = temp_tab + 4 * NRECOMBTAB;
    rec_alphaHep = temp_tab + 5 * NRECOMBTAB;
    rec_alphaHepp = temp_tab + 6 * NRECOMBTAB;
    cool_collisH0 = temp_tab + 7 * NRECOMBTAB;
    cool_collisHe0 = temp_tab + 8 * NRECOMBTAB;
    cool_collisHeP = temp_tab + 9 * NRECOMBTAB;
    cool_recombHp = temp_tab + 10 * NRECOMBTAB;
    cool_recombHeP = temp_tab + 11 * NRECOMBTAB;
    cool_recombHePP = temp_tab + 12 * NRECOMBTAB;
    cool_freefree1 = temp_tab + 13 * NRECOMBTAB;
 
    int i;
    for(i = 0 ; i < NRECOMBTAB; i++)
    {
        temp_tab[i] = RECOMBTMIN + (RECOMBTMAX - RECOMBTMIN) * i / NRECOMBTAB;
        double tt = exp(temp_tab[i]);
        rec_GammaH0[i] = recomb_GammaeH0(tt);
        rec_GammaHe0[i] = recomb_GammaeHe0(tt);
        rec_GammaHep[i] = recomb_GammaeHep(tt);
        rec_alphaHp[i] = recomb_alphaHp(tt);
        /* Note this includes dielectronic recombination*/
        rec_alphaHep[i] = recomb_alphaHepd(tt);
        rec_alphaHepp[i] = recomb_alphaHepp(tt);
        cool_collisH0[i] = cool_CollisionalH0(tt);
        cool_collisHe0[i] = cool_CollisionalHe0(tt);
        cool_collisHeP[i] = cool_CollisionalHeP(tt);
        cool_recombHp[i] = cool_RecombHp(tt);
        cool_recombHeP[i] = cool_RecombHeP(tt);
        cool_recombHePP[i] = cool_RecombHePP(tt);
        cool_freefree1[i] = cool_FreeFree(tt, 1);
    }
 
    /*Initialize the metal cooling table*/
    InitMetalCooling(MetalCoolFile);
}

References cool_collisH0, cool_collisHe0, cool_collisHeP, cool_CollisionalH0(), cool_CollisionalHe0(), cool_CollisionalHeP(), cool_FreeFree(), cool_freefree1, cool_recombHeP, cool_RecombHeP(), cool_recombHePP, cool_RecombHePP(), cool_recombHp, cool_RecombHp(), CoolingParams, CP, cooling_params::fBar, GRAVITY, GrayOpac, GrayOpac_ydata, GrayOpac_zz, HUBBLE, Cosmology::HubbleParam, InitMetalCooling(), load_treecool(), message(), mymalloc, NGRAY, NRECOMBTAB, Cosmology::OmegaBaryon, Cosmology::OmegaCDM, cooling_params::PhotoIonizationOn, rec_alphaHep, rec_alphaHepp, rec_alphaHp, rec_GammaH0, rec_GammaHe0, rec_GammaHep, recomb_alphaHepd(), recomb_alphaHepp(), recomb_alphaHp(), recomb_GammaeH0(), recomb_GammaeHe0(), recomb_GammaeHep(), RECOMBTMAX, RECOMBTMIN, cooling_params::rho_crit_baryon, temp_tab, and tt.

Referenced by init_cooling(), test_heatingcooling_rate(), test_rate_network(), test_recomb_rates(), and test_uvbg_loader().

Here is the call graph for this function:

Here is the caller graph for this function:

init_itp_type()

static void init_itp_type ( double *  xarr, struct itp_type *  Gamma, int  Nelem  )

static

Definition at line 107 of file cooling_rates.c.

{
    Gamma->intp = gsl_interp_alloc(gsl_interp_linear,Nelem);
    gsl_interp_init(Gamma->intp, xarr, Gamma->ydata, Nelem);
}

References itp_type::intp, and itp_type::ydata.

Referenced by load_treecool().

Here is the caller graph for this function:

load_tree_value()

static double load_tree_value ( char **  saveptr )

static

Definition at line 115 of file cooling_rates.c.

{
    double data = atof(strtok_r(NULL, " \t", saveptr));
    if(data > 0)
        return log10(data);
    return -9000;
}

Referenced by load_treecool().

Here is the caller graph for this function:

load_treecool()

static void load_treecool ( const char *  TreeCoolFile )

static

Definition at line 130 of file cooling_rates.c.

{
    if(!CoolingParams.PhotoIonizationOn)
        return;
    FILE * fd = fopen(TreeCoolFile, "r");
    if(!fd)
        endrun(456, "Could not open photon background (TREECOOL) file at: '%s'\n", TreeCoolFile);
 
    /*Find size of file*/
    NTreeCool = 0;
    do
    {
        char buffer[1024];
        char * retval = fgets(buffer, 1024, fd);
        /*Happens on end of file*/
        if(!retval)
            break;
        retval = strtok(buffer, " \t");
        /*Discard comments*/
        if(!retval || retval[0] == '#')
            continue;
        NTreeCool++;
    }
    while(1);
    rewind(fd);
 
    MPI_Bcast(&(NTreeCool), 1, MPI_INT, 0, MPI_COMM_WORLD);
 
    if(NTreeCool<= 2)
        endrun(1, "Photon background contains: %d entries, not enough.\n", NTreeCool);
 
    /*Allocate memory for the photon background table.*/
    Gamma_log1z = (double *) mymalloc("TreeCoolTable", 7 * NTreeCool * sizeof(double));
    Gamma_HI.ydata = Gamma_log1z + NTreeCool;
    Gamma_HeI.ydata = Gamma_log1z + 2 * NTreeCool;
    Gamma_HeII.ydata = Gamma_log1z + 3 * NTreeCool;
    Eps_HI.ydata = Gamma_log1z + 4 * NTreeCool;
    Eps_HeI.ydata = Gamma_log1z + 5 * NTreeCool;
    Eps_HeII.ydata = Gamma_log1z + 6 * NTreeCool;
 
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    if(ThisTask == 0)
    {
        int i = 0;
        while(i < NTreeCool)
        {
            char buffer[1024];
            char * saveptr;
            char * line = fgets(buffer, 1024, fd);
            /*Happens on end of file*/
            if(!line)
                break;
            char * retval = strtok_r(line, " \t", &saveptr);
            if(!retval || retval[0] == '#')
                continue;
            Gamma_log1z[i] = atof(retval);
            /*Get the rest*/
            Gamma_HI.ydata[i]   = load_tree_value(&saveptr);
            Gamma_HeI.ydata[i]  = load_tree_value(&saveptr);
            Gamma_HeII.ydata[i] = load_tree_value(&saveptr);
            Eps_HI.ydata[i]     = load_tree_value(&saveptr)+ CoolingParams.HydrogenHeatAmp;
            Eps_HeI.ydata[i]    = load_tree_value(&saveptr);
            Eps_HeII.ydata[i]   = load_tree_value(&saveptr);
            i++;
        }
 
        fclose(fd);
    }
 
    /*Broadcast data to other processors*/
    MPI_Bcast(Gamma_log1z, 7 * NTreeCool, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    /*Initialize the UVB redshift interpolation: reticulate the splines*/
    init_itp_type(Gamma_log1z, &Gamma_HI, NTreeCool);
    init_itp_type(Gamma_log1z, &Gamma_HeI, NTreeCool);
    init_itp_type(Gamma_log1z, &Gamma_HeII, NTreeCool);
    init_itp_type(Gamma_log1z, &Eps_HI, NTreeCool);
    init_itp_type(Gamma_log1z, &Eps_HeI, NTreeCool);
    init_itp_type(Gamma_log1z, &Eps_HeII, NTreeCool);
 
    message(0, "Read %d lines z = %g - %g from file %s\n", NTreeCool, pow(10, Gamma_log1z[0])-1, pow(10, Gamma_log1z[NTreeCool-1])-1, TreeCoolFile);
}

References CoolingParams, endrun(), Eps_HeII, Gamma_HeII, Gamma_log1z, cooling_params::HydrogenHeatAmp, init_itp_type(), load_tree_value(), message(), mymalloc, NTreeCool, cooling_params::PhotoIonizationOn, ThisTask, and itp_type::ydata.

Referenced by init_cooling_rates().

Here is the call graph for this function:

Here is the caller graph for this function:

ne_internal()

static double ne_internal ( double  nh, double  ienergy, double  ne, double  helium, double *  logt, const struct UVBG *  uvbg  )

static

Definition at line 621 of file cooling_rates.c.

{
    double yy = helium / 4 / (1 - helium);
    *logt = log(get_temp_internal(ne/nh, ienergy, helium));
    double photofac = self_shield_corr(nh, *logt, uvbg->self_shield_dens);
    double nH0 = nH0_internal(*logt, ne, uvbg, photofac);
    double nHp = nHp_internal(nH0);
    struct he_ions He = nHe_internal(nh, *logt, ne, uvbg, photofac);
    return nh * nHp + yy * He.nHep + 2 * yy * He.nHepp;
}

References get_temp_internal(), nH0_internal(), nHe_internal(), he_ions::nHep, he_ions::nHepp, nHp_internal(), self_shield_corr(), and UVBG::self_shield_dens.

Referenced by scipy_optimize_fixed_point().

Here is the call graph for this function:

Here is the caller graph for this function:

nH0_internal()

static double nH0_internal ( double  logt, double  ne, const struct UVBG *  uvbg, double  photofac  )

static

Definition at line 526 of file cooling_rates.c.

{
    double alphaHp = get_interpolated_recomb(logt, rec_alphaHp, &recomb_alphaHp);
    double GammaeH0 = get_interpolated_recomb(logt, rec_GammaH0, &recomb_GammaeH0);
    /*Be careful when there is no ionization.*/
    double photorate = 0;
    if(uvbg->gJH0 > 0. && ne > 1e-50)
        photorate = uvbg->gJH0/ne * photofac;
    return alphaHp/ (alphaHp + GammaeH0 + photorate);
}

References get_interpolated_recomb(), UVBG::gJH0, rec_alphaHp, rec_GammaH0, recomb_alphaHp(), and recomb_GammaeH0().

Referenced by get_heatingcooling_rate(), get_individual_cooling(), get_neutral_fraction_phys_cgs(), and ne_internal().

Here is the call graph for this function:

Here is the caller graph for this function:

nHe_internal()

static struct he_ions nHe_internal ( double  nh, double  logt, double  ne, const struct UVBG *  uvbg, double  photofac  )

static

Definition at line 539 of file cooling_rates.c.

{
    double alphaHep = get_interpolated_recomb(logt, rec_alphaHep, &recomb_alphaHepd);
    double alphaHepp = get_interpolated_recomb(logt, rec_alphaHepp, &recomb_alphaHepp);
    double GammaHe0 = get_interpolated_recomb(logt, rec_GammaHe0, &recomb_GammaeHe0);
    double GammaHep = get_interpolated_recomb(logt, rec_GammaHep, &recomb_GammaeHep);
    struct he_ions He;
    /*Be careful when there is no ionization.*/
    if(uvbg->gJHe0 > 0. && ne > 1e-50) {
        GammaHe0 += uvbg->gJHe0/ne * photofac;
        GammaHep += uvbg->gJHep/ne * photofac;
    }
    /*Deal with the case where there is no ionization separately to avoid NaN.*/
    if(GammaHe0 > 1e-50) {
        He.nHep = nh / (1 + alphaHep / GammaHe0 + GammaHep/alphaHepp);
        He.nHe0 = He.nHep * alphaHep / GammaHe0;
        He.nHepp = He.nHep * GammaHep / alphaHepp;
    }
    else {
        He.nHep = 0;
        He.nHe0 = 1;
        He.nHepp = 0;
    }
    return He;
}

Referenced by get_heatingcooling_rate(), get_helium_ion_phys_cgs(), get_individual_cooling(), and ne_internal().

Here is the caller graph for this function:

nHp_internal()

static double nHp_internal ( double  nH0 )

static

Definition at line 539 of file cooling_rates.c.

{
    double nHp = 1. - nH0;
    if (nHp < 0)
        return 0;
    return nHp;
}

Referenced by get_heatingcooling_rate(), get_individual_cooling(), and ne_internal().

Here is the caller graph for this function:

recomb_alphad()

static double recomb_alphad ( double  temp )

static

Definition at line 406 of file cooling_rates.c.

{
    switch(CoolingParams.recomb)
    {
        case Cen92:
            return 1.9e-3 / pow(temp,1.5) * exp(-4.7e5/temp)*(1+0.3*exp(-9.4e4/temp));
        case Verner96:
        case Badnell06:
             /* Private communication from Avery Meiksin:
                The rate in Black (1981) is wrong by a factor of 0.65. He took the rate fit from Aldrovandi & Pequignot (1973),
                who based it on Burgess (1965), but was unaware of the correction factor in Burgess & Tworkowski (1976, ApJ 205:L105-L107, fig 1).
                The correct dielectronic cooling rate should have the coefficient 0.813e-13 instead of 1.24e-13 as in Black's table 3.
             */
            return 1.23e-3 / pow(temp,1.5) * exp(-4.72e5/temp)*(1+0.3*exp(-9.4e4/temp));
        default:
            endrun(3, "Recombination rate not supported\n");
    }
}

References Badnell06, Cen92, CoolingParams, endrun(), cooling_params::recomb, and Verner96.

Referenced by cool_RecombDielect(), and recomb_alphaHepd().

Here is the call graph for this function:

Here is the caller graph for this function:

recomb_alphaHep()

static double recomb_alphaHep ( double  temp )

static

Definition at line 387 of file cooling_rates.c.

{
    switch(CoolingParams.recomb)
    {
        case Cen92:
            return 1.5e-10 / pow(temp,0.6353);
        case Verner96:
            return _Verner96alphaHep(temp);
        case Badnell06:
            return _Verner96Fit(temp, 1.818E-10, 0.7492, 10.17, 2.786e6);
        default:
            endrun(3, "Recombination rate not supported\n");
    }
}

References _Verner96alphaHep(), _Verner96Fit(), Badnell06, Cen92, CoolingParams, endrun(), cooling_params::recomb, and Verner96.

Referenced by cool_RecombHeP(), and recomb_alphaHepd().

Here is the call graph for this function:

Here is the caller graph for this function:

recomb_alphaHepd()

static double recomb_alphaHepd ( double  temp )

static

Definition at line 427 of file cooling_rates.c.

{
    return recomb_alphad(temp) + recomb_alphaHep(temp);
}

References recomb_alphad(), and recomb_alphaHep().

Referenced by init_cooling_rates().

Here is the call graph for this function:

Here is the caller graph for this function:

recomb_alphaHepp()

static double recomb_alphaHepp ( double  temp )

static

Definition at line 434 of file cooling_rates.c.

{
    switch(CoolingParams.recomb)
    {
        case Cen92:
            return 4 * recomb_alphaHp(temp);
        case Verner96:
            return _Verner96Fit(temp, 1.891e-10, 0.7524, 9.370, 2.774e6);
        case Badnell06:
            return _Verner96Fit(temp, 5.235E-11, 0.6988 + 0.0829*exp(-1.682e5/temp), 7.301, 4.475e6);
        default:
            endrun(3, "Recombination rate not supported\n");
    }
}

References _Verner96Fit(), Badnell06, Cen92, CoolingParams, endrun(), cooling_params::recomb, recomb_alphaHp(), and Verner96.

Referenced by cool_RecombHePP(), and init_cooling_rates().

Here is the call graph for this function:

Here is the caller graph for this function:

recomb_alphaHp()

double recomb_alphaHp

(

double

temp

)

Definition at line 347 of file cooling_rates.c.

{
    switch(CoolingParams.recomb)
    {
        case Cen92:
            return 8.4e-11 / sqrt(temp) / pow(temp/1000, 0.2) / (1+ pow(temp/1e6, 0.7));
        case Verner96:
            /*The V&F 96 fitting formula is accurate to < 1% in the worst case.
            See line 1 of V&F96 table 1.*/
            return _Verner96Fit(temp, 7.982e-11, 0.748, 3.148, 7.036e+05);
        case Badnell06:
            /*Slightly different fit coefficients*/
            return _Verner96Fit(temp, 8.318e-11, 0.7472, 2.965, 7.001e5);
        default:
            endrun(3, "Recombination rate not supported\n");
    }
}

References _Verner96Fit(), Badnell06, Cen92, CoolingParams, endrun(), cooling_params::recomb, and Verner96.

Referenced by cool_RecombHp(), init_cooling_rates(), nH0_internal(), recomb_alphaHepp(), and test_recomb_rates().

Here is the call graph for this function:

Here is the caller graph for this function:

recomb_GammaeH0()

static double recomb_GammaeH0 ( double  temp )

static

Definition at line 459 of file cooling_rates.c.

{
    switch(CoolingParams.recomb)
    {
        case Cen92:
            return 5.85e-11 * sqrt(temp) * exp(-157809.1/temp) / (1 + sqrt(temp/1e5));
        case Verner96:
        case Badnell06:
            //Voronov 97 Table 1
            return _Voronov96Fit(temp, 13.6, 0, 0.291e-07, 0.232, 0.39);
        default:
            endrun(3, "Collisional rate not supported\n");
    }
}

References _Voronov96Fit(), Badnell06, Cen92, CoolingParams, endrun(), cooling_params::recomb, and Verner96.

Referenced by cool_CollisionalIonizeH0(), init_cooling_rates(), and nH0_internal().

Here is the call graph for this function:

Here is the caller graph for this function:

recomb_GammaeHe0()

static double recomb_GammaeHe0 ( double  temp )

static

Definition at line 476 of file cooling_rates.c.

{
    switch(CoolingParams.recomb)
    {
        case Cen92:
            return 2.38e-11 * sqrt(temp) * exp(-285335.4/temp) / (1+ sqrt(temp/1e5));
        case Verner96:
        case Badnell06:
            //Voronov 97 Table 1
            return _Voronov96Fit(temp, 24.6, 0, 0.175e-07, 0.180, 0.35);
        default:
            endrun(3, "Collisional rate not supported\n");
    }
}

References _Voronov96Fit(), Badnell06, Cen92, CoolingParams, endrun(), cooling_params::recomb, and Verner96.

Referenced by cool_CollisionalIonizeHe0(), and init_cooling_rates().

Here is the call graph for this function:

Here is the caller graph for this function:

recomb_GammaeHep()

static double recomb_GammaeHep ( double  temp )

static

Definition at line 493 of file cooling_rates.c.

{
    switch(CoolingParams.recomb)
    {
        case Cen92:
            return 5.68e-12 * sqrt(temp) * exp(-631515.0/temp) / (1+ sqrt(temp/1e5));
        case Verner96:
        case Badnell06:
            //Voronov 97 Table 1
            return _Voronov96Fit(temp, 54.4, 1, 0.205e-08, 0.265, 0.25);
        default:
            endrun(3, "Collisional rate not supported\n");
    }
}

References _Voronov96Fit(), Badnell06, Cen92, CoolingParams, endrun(), cooling_params::recomb, and Verner96.

Referenced by cool_CollisionalIonizeHeP(), and init_cooling_rates().

Here is the call graph for this function:

Here is the caller graph for this function:

scipy_optimize_fixed_point()

static double scipy_optimize_fixed_point ( double  ne_init, double  nh, double  ienergy, double  helium, double *  logt, const struct UVBG *  uvbg  )

static

Definition at line 645 of file cooling_rates.c.

{
    int i;
    double ne0 = ne_init;
    for(i = 0; i < MAXITER; i++)
    {
        double logt1;
        double ne1 = ne_internal(nh, ienergy, ne0*nh, helium, &logt1, uvbg) / nh;
 
        if(fabs(ne1 - ne0) < ITERCONV) {
            *logt = logt1;
            ne0 = ne1;
            break;
        }
 
        double ne2 = ne_internal(nh, ienergy, ne1*nh, helium, &logt1, uvbg) / nh;
        double d = ne0 + ne2 - 2.0 * ne1;
        double pp = ne2;
        /*This is del^2*/
        if (d > 1e-15 || d < -1e-15)
            pp = ne0 - (ne1 - ne0)*(ne1 - ne0) / d;
        ne0 = pp;
        /*Enforce positivity*/
        if(ne0 < 0)
            ne0 = 0;
    }
    if (!isfinite(ne0) || i == MAXITER)
        endrun(1, "Ionization rate network failed to converge for nh = %g temp = %g helium=%g ienergy=%g: last ne = %g (init=%g)\n", nh, get_temp_internal(ne0, ienergy, helium), helium, ienergy, ne0, ne_init);
    return ne0 * nh;
}

References endrun(), get_temp_internal(), ITERCONV, MAXITER, and ne_internal().

Referenced by get_equilib_ne().

Here is the call graph for this function:

Here is the caller graph for this function:

self_shield_corr()

static double self_shield_corr ( double  nh, double  logt, double  ssdens  )

static

Definition at line 304 of file cooling_rates.c.

{
    /* Turn off self-shielding for low-density gas.
     * If such gas becomes very cold, this is not strictly what they find,
     * but I think it is more physical*/
    if(!CoolingParams.SelfShieldingOn || nh < ssdens * 0.01)
        return 1;
    /* T/1e4**0.17*/
    double T4 = exp(0.17 * (logt - log(1e4)));
    double nSSh = 1.003*ssdens*T4;
    return 0.98*pow(1+pow(nh/nSSh,1.64),-2.28)+0.02*pow(1+nh/nSSh, -0.84);
}

References CoolingParams, and cooling_params::SelfShieldingOn.

Referenced by get_heatingcooling_rate(), get_helium_ion_phys_cgs(), get_individual_cooling(), get_neutral_fraction_phys_cgs(), and ne_internal().

Here is the caller graph for this function:

self_shield_dens()

static double self_shield_dens ( double  redshift, const struct UVBG *  uvbg  )

static

Definition at line 244 of file cooling_rates.c.

{
    /*Before the UVBG switches on, no need for self-shielding*/
    if(uvbg->gJH0 == 0)
        return 1e10;
    double G12 = uvbg->gJH0/1e-12;
    double greyopac;
    if (redshift <= GrayOpac_zz[0])
        greyopac = GrayOpac_ydata[0];
    else if (redshift >= GrayOpac_zz[NGRAY-1])
        greyopac = GrayOpac_ydata[NGRAY-1];
    else {
        greyopac = gsl_interp_eval(GrayOpac, GrayOpac_zz, GrayOpac_ydata,redshift, NULL);
    }
    return 6.73e-3 * pow(greyopac / 2.49e-18, -2./3)*pow(G12, 2./3)*pow(CoolingParams.fBar/0.17,-1./3);
}

set_cooling_params()

void set_cooling_params

(

ParameterSet *

ps

)

Definition at line 942 of file cooling_rates.c.

{
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    if(ThisTask == 0) {
        /*Cooling rate network parameters*/
        CoolingParams.CMBTemperature = param_get_double(ps, "CMBTemperature");
        CoolingParams.cooling = (enum CoolingType) param_get_enum(ps, "CoolingRates"); // Sherwood;
        CoolingParams.recomb = (enum RecombType) param_get_enum(ps, "RecombRates"); // Verner96;
        CoolingParams.SelfShieldingOn = param_get_int(ps, "SelfShieldingOn");
        CoolingParams.PhotoIonizeFactor = param_get_double(ps, "PhotoIonizeFactor");
        CoolingParams.PhotoIonizationOn = param_get_int(ps, "PhotoIonizationOn");
        CoolingParams.MinGasTemp = param_get_double(ps, "MinGasTemp");
        CoolingParams.UVRedshiftThreshold = param_get_double(ps, "UVRedshiftThreshold");
        CoolingParams.HydrogenHeatAmp = log10(param_get_double(ps, "HydrogenHeatAmp"));
 
        /*Helium model parameters*/
        CoolingParams.HeliumHeatOn = param_get_int(ps, "HeliumHeatOn");
        CoolingParams.HeliumHeatThresh = param_get_double(ps, "HeliumHeatThresh");
        CoolingParams.HeliumHeatAmp = param_get_double(ps, "HeliumHeatAmp");
        CoolingParams.HeliumHeatExp = param_get_double(ps, "HeliumHeatExp");
    }
    MPI_Bcast(&CoolingParams, sizeof(struct cooling_params), MPI_BYTE, 0, MPI_COMM_WORLD);
}

References cooling_params::CMBTemperature, cooling_params::cooling, CoolingParams, cooling_params::HeliumHeatAmp, cooling_params::HeliumHeatExp, cooling_params::HeliumHeatOn, cooling_params::HeliumHeatThresh, cooling_params::HydrogenHeatAmp, cooling_params::MinGasTemp, param_get_double(), param_get_enum(), param_get_int(), cooling_params::PhotoIonizationOn, cooling_params::PhotoIonizeFactor, cooling_params::recomb, cooling_params::SelfShieldingOn, ThisTask, and cooling_params::UVRedshiftThreshold.

Referenced by read_parameter_file().

Here is the call graph for this function:

Here is the caller graph for this function:

set_coolpar()

void set_coolpar

(

struct cooling_params

cp

)

Definition at line 936 of file cooling_rates.c.

{
    CoolingParams = cp;
}

References CoolingParams.

Referenced by test_DoCooling(), test_heatingcooling_rate(), test_rate_network(), test_recomb_rates(), and test_uvbg_loader().

Here is the caller graph for this function:

Variable Documentation

cool_collisH0

double* cool_collisH0

static

Definition at line 101 of file cooling_rates.c.

Referenced by get_heatingcooling_rate(), get_individual_cooling(), and init_cooling_rates().

cool_collisHe0

double * cool_collisHe0

static

Definition at line 101 of file cooling_rates.c.

Referenced by get_heatingcooling_rate(), get_individual_cooling(), and init_cooling_rates().

cool_collisHeP

double * cool_collisHeP

static

Definition at line 101 of file cooling_rates.c.

Referenced by get_heatingcooling_rate(), get_individual_cooling(), and init_cooling_rates().

cool_freefree1

double* cool_freefree1

static

Definition at line 104 of file cooling_rates.c.

Referenced by get_heatingcooling_rate(), get_individual_cooling(), and init_cooling_rates().

cool_recombHeP

double * cool_recombHeP

static

Definition at line 102 of file cooling_rates.c.

Referenced by get_heatingcooling_rate(), get_individual_cooling(), and init_cooling_rates().

cool_recombHePP

double * cool_recombHePP

static

Definition at line 102 of file cooling_rates.c.

Referenced by get_heatingcooling_rate(), get_individual_cooling(), and init_cooling_rates().

cool_recombHp

double* cool_recombHp

static

Definition at line 102 of file cooling_rates.c.

Referenced by get_heatingcooling_rate(), get_individual_cooling(), and init_cooling_rates().

CoolingParams

struct cooling_params CoolingParams

static

Definition at line 1 of file cooling_rates.c.

Referenced by _t5(), cool_CollisionalH0(), cool_FreeFree(), cool_he_reion_factor(), cool_InverseCompton(), cool_RecombHePP(), cool_RecombHp(), get_global_UVBG(), get_heatingcooling_rate(), get_individual_cooling(), get_photo_rate(), get_temp_internal(), init_cooling_rates(), load_treecool(), recomb_alphad(), recomb_alphaHep(), recomb_alphaHepp(), recomb_alphaHp(), recomb_GammaeH0(), recomb_GammaeHe0(), recomb_GammaeHep(), self_shield_corr(), set_cooling_params(), and set_coolpar().

Eps_HeII

struct itp_type Eps_HI Eps_HeI Eps_HeII

static

Definition at line 88 of file cooling_rates.c.

Referenced by load_treecool().

Gamma_HeII

struct itp_type Gamma_HI Gamma_HeI Gamma_HeII

static

Definition at line 88 of file cooling_rates.c.

Referenced by load_treecool().

Gamma_log1z

double* Gamma_log1z

static

Definition at line 88 of file cooling_rates.c.

Referenced by get_photo_rate(), and load_treecool().

GrayOpac

gsl_interp* GrayOpac

static

Definition at line 70 of file cooling_rates.c.

Referenced by get_global_UVBG(), and init_cooling_rates().

GrayOpac_ydata

const double GrayOpac_ydata[NGRAY] = { 2.59e-18, 2.37e-18, 2.27e-18, 2.15e-18, 2.02e-18, 1.94e-18}

static

Definition at line 75 of file cooling_rates.c.

Referenced by get_global_UVBG(), and init_cooling_rates().

GrayOpac_zz

const double GrayOpac_zz[NGRAY] = {0, 1, 2, 3, 4, 5}

static

Definition at line 76 of file cooling_rates.c.

Referenced by get_global_UVBG(), and init_cooling_rates().

NTreeCool

int NTreeCool

static

Definition at line 86 of file cooling_rates.c.

Referenced by get_photo_rate(), and load_treecool().

rec_alphaHep

double * rec_alphaHep

static

Definition at line 99 of file cooling_rates.c.

Referenced by init_cooling_rates().

rec_alphaHepp

double * rec_alphaHepp

static

Definition at line 99 of file cooling_rates.c.

Referenced by init_cooling_rates().

rec_alphaHp

double* rec_alphaHp

static

Definition at line 99 of file cooling_rates.c.

Referenced by init_cooling_rates(), and nH0_internal().

rec_GammaH0

double* rec_GammaH0

static

Definition at line 100 of file cooling_rates.c.

Referenced by init_cooling_rates(), and nH0_internal().

rec_GammaHe0

double * rec_GammaHe0

static

Definition at line 100 of file cooling_rates.c.

Referenced by init_cooling_rates().

rec_GammaHep

double * rec_GammaHep

static

Definition at line 100 of file cooling_rates.c.

Referenced by init_cooling_rates().

temp_tab

double* temp_tab

static

Definition at line 98 of file cooling_rates.c.

Referenced by init_cooling_rates().