-
Notifications
You must be signed in to change notification settings - Fork 0
/
halo.cpp
155 lines (142 loc) · 4.77 KB
/
halo.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
/*
* halo.cpp
*
* Created on: Feb 2, 2012
* Author: cgiocoli
*/
#include <halo.h>
/** \ingroup cosmolib
* \brief Constructor initializing a NFW-halo
*/
HALOCalculator::HALOCalculator (
COSMOLOGY *cos /// pointer to a COSMOLOGY
,double mass /// halo mass in M_sun/h TODO Carlo: Is this really M_sun/h
,double redshift /// redshift
)
: co(cos),m(mass),z(redshift){
Set_Parameters();
}
void HALOCalculator::Set_Parameters(){
Omz=co->Omegam(z);
Omo=co->getOmega_matter();
Oml=co->getOmega_lambda();
sigma2M=co->TopHatVariance(m);
deltac0=1.68647;
if(Omo<1 && Oml==0) deltac0*=pow(Omz,0.0185);
if(Omo+Oml==1) deltac0*=pow(Omz,0.0055);
}
HALOCalculator::~HALOCalculator () {}
/** \ingroup cosmolib
* \brief Reset halo mass and redshift
*/
void HALOCalculator::reset(double mr,double zr){
m = mr;
z = zr;
Set_Parameters();
}
/** \ingroup cosmolib
* \brief Virial radius of the halo in Mpc
*/
double HALOCalculator:: getRvir(
int caseunit /// by default uses the Brayan and Norman fit, if equal to 1 uses the fit by Felix and Stoer
){
double d=co->DeltaVir(z,caseunit)*co->rho_crit(z);
return pow( 3*m/(4*M_PI*d), 0.3333 )/(1+z);
}
/** \ingroup cosmolib
* \brief Radius in Mpc at which the enclosed density reach 200 times the critical value
*/
double HALOCalculator:: getR200(){
double d=200.0*co->rho_crit(z);
return pow( 3*m/(4*M_PI*d), 0.3333 )/(1+z);
}
/** \ingroup cosmolib
* \brief The median redshift at which the main halo progenitor assembles
* 50% of the halo mass, if 0<f<1 is given, is returned the redshift at which
* this fraction is assembled
*/
double HALOCalculator:: getFormationRedshift(double f){
double sigma2fM=co->TopHatVariance(f*m);
double alphaf = 0.815*exp(-2*f*f*f)/pow(f,0.707);
double wmed = sqrt(2*log(alphaf+1));
double deltacz = deltac0/co->Dgrowth(z)+wmed*sqrt(sigma2fM-sigma2M);
return co->getZfromDeltaC(deltacz);
}
/** \ingroup cosmolib
* \brief The median time at which the main halo progenitor assembles
* 50% of the halo mass, if 0<f<1 is given, is returned the redshift at which
* this fraction is assembled
*/
double HALOCalculator:: getFormationTime(double f){
double sigma2fM=co->TopHatVariance(f*m);
double alphaf = 0.815*exp(-2*f*f*f)/pow(f,0.707);
double wmed = sqrt(2*log(alphaf+1));
double deltacz = deltac0/co->Dgrowth(z)+wmed*sqrt(sigma2fM-sigma2M);
return co->getTimefromDeltaC(deltacz);
}
/** \ingroup cosmolib
* \brief The halo concentration according to Zhao et al. 2009 is returned,
* if the integer 1 is given the concentration according to Munoz-Cuartas et al. 2011 is returned
* if 2 the concentration according to Giocoli et al. 2012 is returned
* is returned
*/
double HALOCalculator:: getConcentration(
int caseunit /// set relation used - (0) Zhao et al. 2009, (1) Munoz-Cuartas et al. 2011 (2) Giocoli et al. 2012 (3) power-law c-m relation
,double alpha0 /// power-law slope if caseunit == 3
){
double w=0.029;
double mu=0.097;
double alpha=-110.001;
double beta=2469.720;
double gamma=16.885;
double a,b,logc;
double t004,t05,t0;
double h0,hz;
switch (caseunit){
case (1): // Munoz-Cuartas et al. 2011
a = w*z-mu;
b=alpha/(z+gamma)+beta/(z+gamma)/(z+gamma);
logc=a*log10(m*co->gethubble())+b;
return pow(10.,logc);
break;
case (-1): // Munoz-Cuartas et al. 2011, as the previous but M is assumet in M/h
a = w*z-mu;
b=alpha/(z+gamma)+beta/(z+gamma)/(z+gamma);
logc=a*log10(m)+b;
return pow(10.,logc);
break;
case (2): // Giocoli et al. 2012
t004=getFormationTime(0.04);
t05=getFormationTime(0.5);
t0=co->time(z);
return 0.45*(4.23+pow(t0/t004,1.15)+pow(t0/t05,2.3));
break;
case (3): // power-law c-m relation
alpha0 = alpha0/(1+z);
h0 = co->gethubble();
hz = h0/co->drdz(1+z);
return 10.*pow(m*h0/1.e+12,alpha0)*pow(h0/hz,2./3.);
break;
case (-3): // power-law c-m relation, as the previous but M is assumed in M/h
alpha0 = alpha0/(1+z);
h0 = co->gethubble();
hz = h0/co->drdz(1+z);
return 10.*pow(m/1.e+12,alpha0)*pow(h0/hz,2./3.);
break;
default: // Zhao et al. 2009
t004=getFormationTime(0.04);
t0=co->time(z);
return 4*pow(1+pow(t0/3.75/t004,8.4),1./8.);
}
return 0;
}
/** \brief returns the fraction of halo+galaxy in stars according to from Moster et al. 2010ApJ...710..903M
*/
double HALOCalculator::MosterStellarMassFraction(
double Mtotal /// total mass of halo+galaxy in solar masses
){
static const double mo=7.3113e10,M1=2.8575e10,gam1=7.17,gam2=0.201,be=0.557;
// from Moster et al. 2010ApJ...710..903M
return mo*pow(Mtotal/M1,gam1)
/pow(1+pow(Mtotal/M1,be),(gam1-gam2)/be)/Mtotal;
}