Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F83525171
cooling.c
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Sep 17, 15:04
Size
46 KB
Mime Type
text/x-c
Expires
Thu, Sep 19, 15:04 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
20852220
Attached To
R195 SCITAS application benchmarks
cooling.c
View Options
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
#include "allvars.h"
#include "proto.h"
#ifdef COOLING
/*! initialize cooling function (the metallicity is fixed)
*
* T = temperature
* L0 = m000 primordial metallicity
* L1 = m-30
* L2 = m-20
* L3 = m-10
* L4 = m-05
* L5 = m-00 solar metallicity
* L6 = m+05
*/
int
init_cooling
(
FLOAT
metallicity
)
{
FILE
*
fd
;
int
n
,
i
;
char
line
[
72
];
float
T
,
L0
,
L1
,
L2
,
L3
,
L4
,
L5
,
L6
;
int
MetallicityIndex
=
4
;
/* find the right index */
if
(
All
.
InitGasMetallicity
<-
3
)
MetallicityIndex
=
0
;
else
{
if
(
All
.
InitGasMetallicity
<-
2
)
MetallicityIndex
=
1
;
else
{
if
(
All
.
InitGasMetallicity
<-
1
)
MetallicityIndex
=
2
;
else
{
if
(
All
.
InitGasMetallicity
<-
0.5
)
MetallicityIndex
=
3
;
else
{
if
(
All
.
InitGasMetallicity
<
0
)
MetallicityIndex
=
4
;
else
{
MetallicityIndex
=
5
;
}
}
}
}
}
fd
=
fopen
(
All
.
CoolingFile
,
"r"
);
fscanf
(
fd
,
"# %6d
\n
"
,
&
n
);
/* allocate memory */
All
.
logT
=
malloc
(
n
*
sizeof
(
double
));
All
.
logL
=
malloc
(
n
*
sizeof
(
double
));
/* read empty line */
fgets
(
line
,
sizeof
(
line
),
fd
);
/* read file */
for
(
i
=
0
;
i
<
n
;
i
++
){
fscanf
(
fd
,
"%f %f %f %f %f %f %f %f
\n
"
,
&
T
,
&
L0
,
&
L1
,
&
L2
,
&
L3
,
&
L4
,
&
L5
,
&
L6
);
//printf("%8.3f %8.3f\n",T,L0);
/* keep only solar values */
All
.
logT
[
i
]
=
(
double
)
T
;
switch
(
MetallicityIndex
)
{
case
0
:
All
.
logL
[
i
]
=
(
double
)
L0
;
break
;
case
1
:
All
.
logL
[
i
]
=
(
double
)
L1
;
break
;
case
2
:
All
.
logL
[
i
]
=
(
double
)
L2
;
break
;
case
3
:
All
.
logL
[
i
]
=
(
double
)
L3
;
break
;
case
4
:
All
.
logL
[
i
]
=
(
double
)
L4
;
break
;
case
5
:
All
.
logL
[
i
]
=
(
double
)
L5
;
break
;
case
6
:
All
.
logL
[
i
]
=
(
double
)
L6
;
break
;
}
}
fclose
(
fd
);
/* init interpolation */
All
.
acc_cooling_spline
=
gsl_interp_accel_alloc
();
All
.
cooling_spline
=
gsl_spline_alloc
(
gsl_interp_cspline
,
n
);
gsl_spline_init
(
All
.
cooling_spline
,
All
.
logT
,
All
.
logL
,
n
);
#ifdef OUTPUT_COOLING_FUNCTION
/* test cooling */
double
logT
;
double
l
;
logT
=
1.
;
while
(
logT
<
8
)
{
T
=
pow
(
10
,
logT
);
l
=
log10
(
cooling_function
(
T
));
if
(
ThisTask
==
0
)
printf
(
"%8.3f %8.3f
\n
"
,
logT
,
l
);
logT
=
logT
+
0.05
;
}
#endif
return
0
;
}
/*! This function return the normalized cooling function (no metallicity dependency)
*/
double
cooling_function
(
double
temperature
)
{
double
logT
;
if
(
temperature
>=
All
.
CutofCoolingTemperature
)
{
logT
=
log10
(
temperature
);
if
(
logT
>
8.5
)
logT
=
8.5
;
return
pow
(
10
,
gsl_spline_eval
(
All
.
cooling_spline
,
logT
,
All
.
acc_cooling_spline
));
}
else
return
1e-100
;
}
/***************************************************************************
METALLICITY DEPENDENT COOLING
**************************************************************************/
int
init_cooling_with_metals
()
{
/*
zmin zmax slz
tmin tmax slt
FeHSolar
p k
*/
FILE
*
fd
;
int
p
,
k
,
i
,
j
;
float
zmin
,
zmax
,
slz
,
tmin
,
tmax
,
slt
,
FeHSolar
;
float
lbd
;
if
(
ThisTask
==
0
)
{
fd
=
fopen
(
All
.
CoolingFile
,
"r"
);
fscanf
(
fd
,
"%f %f %f
\n
"
,
&
zmin
,
&
zmax
,
&
slz
);
fscanf
(
fd
,
"%f %f %f
\n
"
,
&
tmin
,
&
tmax
,
&
slt
);
fscanf
(
fd
,
"%f
\n
"
,
&
FeHSolar
);
fscanf
(
fd
,
"%d %d
\n
"
,
&
p
,
&
k
);
All
.
CoolingParameters_zmin
=
zmin
;
All
.
CoolingParameters_zmax
=
zmax
;
All
.
CoolingParameters_slz
=
slz
;
All
.
CoolingParameters_tmin
=
tmin
;
All
.
CoolingParameters_tmax
=
tmax
;
All
.
CoolingParameters_slt
=
slt
;
//All.CoolingParameters_FeHSolar = FEH_SOLAR; /* instead of FeHSolar*/
All
.
CoolingParameters_cooling_data_max
=
k
-
2
;
for
(
i
=
0
;
i
<
p
;
i
++
)
for
(
j
=
0
;
j
<
k
;
j
++
)
{
fscanf
(
fd
,
"%f
\n
"
,
&
lbd
);
All
.
CoolingParameters_cooling_data
[
i
][
j
]
=
lbd
;
}
fclose
(
fd
);
}
/* now broadcast */
/*
This is quite bad to do it like this, however, there is no other solution, as
All has already be sent.
The other solution will be to create a structure devoted to the cooling (like Cps in chimie.c)
avoiding to link the parameters to All.
*/
MPI_Bcast
(
&
All
.
CoolingParameters_zmin
,
1
,
MPI_DOUBLE
,
0
,
MPI_COMM_WORLD
);
MPI_Bcast
(
&
All
.
CoolingParameters_zmax
,
1
,
MPI_DOUBLE
,
0
,
MPI_COMM_WORLD
);
MPI_Bcast
(
&
All
.
CoolingParameters_slz
,
1
,
MPI_DOUBLE
,
0
,
MPI_COMM_WORLD
);
MPI_Bcast
(
&
All
.
CoolingParameters_tmin
,
1
,
MPI_DOUBLE
,
0
,
MPI_COMM_WORLD
);
MPI_Bcast
(
&
All
.
CoolingParameters_tmax
,
1
,
MPI_DOUBLE
,
0
,
MPI_COMM_WORLD
);
MPI_Bcast
(
&
All
.
CoolingParameters_slt
,
1
,
MPI_DOUBLE
,
0
,
MPI_COMM_WORLD
);
MPI_Bcast
(
&
All
.
CoolingParameters_cooling_data_max
,
1
,
MPI_DOUBLE
,
0
,
MPI_COMM_WORLD
);
MPI_Bcast
(
&
All
.
CoolingParameters_cooling_data
,
(
COOLING_NMETALICITIES
*
COOLING_NTEMPERATURES
),
MPI_DOUBLE
,
0
,
MPI_COMM_WORLD
);
#ifdef OUTPUT_COOLING_FUNCTION
/* test cooling */
double
logT
,
T
;
double
l
;
double
metal
;
logT
=
1.
;
metal
=
(
pow
(
10
,
All
.
InitGasMetallicity
)
-
1e-10
)
*
All
.
CoolingParameters_FeHSolar
;
while
(
logT
<
8
)
{
T
=
pow
(
10
,
logT
);
l
=
log10
(
cooling_function_with_metals
(
T
,
metal
));
if
(
ThisTask
==
0
)
printf
(
"%8.3f %8.3f
\n
"
,
logT
,
l
);
logT
=
logT
+
0.05
;
}
#endif
return
0
;
}
/*! This function return the normalized cooling function, that depends on metallicity
*/
double
cooling_function_with_metals
(
double
temperature
,
double
metal
)
{
double
cooling
;
double
T
,
Z
;
double
rt
,
rz
,
ft
,
fz
,
v1
,
v2
,
v
;
int
it
,
iz
,
itp
,
izp
;
double
zmin
,
zmax
,
slz
,
tmin
,
tmax
,
slt
,
FeHSolar
,
cooling_data_max
;
zmin
=
All
.
CoolingParameters_zmin
;
zmax
=
All
.
CoolingParameters_zmax
;
slz
=
All
.
CoolingParameters_slz
;
tmin
=
All
.
CoolingParameters_tmin
;
tmax
=
All
.
CoolingParameters_tmax
;
slt
=
All
.
CoolingParameters_slt
;
FeHSolar
=
All
.
CoolingParameters_FeHSolar
;
cooling_data_max
=
All
.
CoolingParameters_cooling_data_max
;
cooling
=
0.0
;
T
=
log10
(
temperature
);
Z
=
log10
(
metal
/
FeHSolar
+
1.e-10
);
if
(
Z
>
zmax
)
{
/*print *,'Warning: Z>Zmax for',i*/
Z
=
zmax
;
}
if
(
Z
<
zmin
)
{
rt
=
(
T
-
tmin
)
/
slt
;
it
=
(
int
)
rt
;
if
(
it
<
cooling_data_max
)
it
=
(
int
)
rt
;
else
it
=
cooling_data_max
;
itp
=
it
+
1
;
ft
=
rt
-
it
;
fz
=
(
10.
+
Z
)
/
(
10.
+
zmin
);
v1
=
ft
*
(
All
.
CoolingParameters_cooling_data
[
1
][
itp
]
-
All
.
CoolingParameters_cooling_data
[
1
][
it
]
)
+
All
.
CoolingParameters_cooling_data
[
1
][
it
];
v2
=
ft
*
(
All
.
CoolingParameters_cooling_data
[
0
][
itp
]
-
All
.
CoolingParameters_cooling_data
[
0
][
it
]
)
+
All
.
CoolingParameters_cooling_data
[
0
][
it
];
v
=
v2
+
fz
*
(
v1
-
v2
);
}
else
{
rt
=
(
T
-
tmin
)
/
slt
;
rz
=
(
Z
-
zmin
)
/
slz
+
1.0
;
it
=
(
int
)
rt
;
if
(
it
<
cooling_data_max
)
it
=
(
int
)
rt
;
else
it
=
cooling_data_max
;
iz
=
(
int
)
rz
;
itp
=
it
+
1
;
izp
=
iz
+
1
;
ft
=
rt
-
it
;
fz
=
rz
-
iz
;
v1
=
ft
*
(
All
.
CoolingParameters_cooling_data
[
izp
][
itp
]
-
All
.
CoolingParameters_cooling_data
[
izp
][
it
])
+
All
.
CoolingParameters_cooling_data
[
izp
][
it
];
v2
=
ft
*
(
All
.
CoolingParameters_cooling_data
[
iz
][
itp
]
-
All
.
CoolingParameters_cooling_data
[
iz
][
it
])
+
All
.
CoolingParameters_cooling_data
[
iz
][
it
];
v
=
v2
+
fz
*
(
v1
-
v2
);
}
cooling
=
pow
(
10
,
v
);
return
cooling
;
}
/***************************************************************************
END OF METALLICITY DEPENDENT COOLING
**************************************************************************/
/*! \file cooling.c
* \brief Compute gas cooling
*
*/
static
double
hubble_a
,
a3inv
;
static
double
eV
=
1.6022000e-12
;
static
double
normfacJ0
=
0.74627
;
static
double
J0min
=
1.e-29
;
static
double
alpha
=
1.0
;
static
int
Norderweinberg
=
7
;
/* polynom order+1 */
static
double
coefweinberg
[
7
][
6
];
static
double
z
;
static
double
J0
;
static
double
Cte_G_gHI
;
static
double
Cte_G_gHeI
;
static
double
Cte_G_gHeII
;
static
double
Cte_heating_radiative_HI
;
static
double
Cte_heating_radiative_HeI
;
static
double
Cte_heating_radiative_HeII
;
/*
* init some variables that depends only on redshift
*/
void
init_from_new_redshift
(
double
Redshift
)
{
/* init weinberg coeff */
coefweinberg
[
0
][
0
]
=
-
0.31086729929951613e+002
;
coefweinberg
[
1
][
0
]
=
0.34803667059463761e+001
;
coefweinberg
[
2
][
0
]
=
-
0.15145716066316397e+001
;
coefweinberg
[
3
][
0
]
=
0.54649951450632972e+000
;
coefweinberg
[
4
][
0
]
=
-
0.16395924120387340e+000
;
coefweinberg
[
5
][
0
]
=
0.25197466148524143e-001
;
coefweinberg
[
6
][
0
]
=
-
0.15352763785487806e-002
;
coefweinberg
[
0
][
1
]
=
-
0.31887274113252204e+002
;
coefweinberg
[
1
][
1
]
=
0.44178493140927095e+001
;
coefweinberg
[
2
][
1
]
=
-
0.20158132553082293e+001
;
coefweinberg
[
3
][
1
]
=
0.64080497292269134e+000
;
coefweinberg
[
4
][
1
]
=
-
0.15981267091909040e+000
;
coefweinberg
[
5
][
1
]
=
0.22056900050237707e-001
;
coefweinberg
[
6
][
1
]
=
-
0.12837570029562849e-002
;
coefweinberg
[
0
][
2
]
=
-
0.35693331167978656e+002
;
coefweinberg
[
1
][
2
]
=
0.20207245722165794e+001
;
coefweinberg
[
2
][
2
]
=
-
0.76856976101363744e-001
;
coefweinberg
[
3
][
2
]
=
-
0.75691470654320359e-001
;
coefweinberg
[
4
][
2
]
=
-
0.54502220282734729e-001
;
coefweinberg
[
5
][
2
]
=
0.20633345104660583e-001
;
coefweinberg
[
6
][
2
]
=
-
0.18410307456285177e-002
;
coefweinberg
[
0
][
3
]
=
-
0.56967559787460921e+002
;
coefweinberg
[
1
][
3
]
=
0.38601174525546353e+001
;
coefweinberg
[
2
][
3
]
=
-
0.18318926655684415e+001
;
coefweinberg
[
3
][
3
]
=
0.67360594266440688e+000
;
coefweinberg
[
4
][
3
]
=
-
0.18983466813215341e+000
;
coefweinberg
[
5
][
3
]
=
0.27768907786915147e-001
;
coefweinberg
[
6
][
3
]
=
-
0.16330066969315893e-002
;
coefweinberg
[
0
][
4
]
=
-
0.56977907250821026e+002
;
coefweinberg
[
1
][
4
]
=
0.38686249565302266e+001
;
coefweinberg
[
2
][
4
]
=
-
0.13330942368518774e+001
;
coefweinberg
[
3
][
4
]
=
0.33988839029092172e+000
;
coefweinberg
[
4
][
4
]
=
-
0.98997915675929332e-001
;
coefweinberg
[
5
][
4
]
=
0.16781612113050747e-001
;
coefweinberg
[
6
][
4
]
=
-
0.11514328893746039e-002
;
coefweinberg
[
0
][
5
]
=
-
0.59825233828609278e+002
;
coefweinberg
[
1
][
5
]
=
0.21898162706563347e+001
;
coefweinberg
[
2
][
5
]
=
-
0.42982055888598525e+000
;
coefweinberg
[
3
][
5
]
=
0.50312144291614215e-001
;
coefweinberg
[
4
][
5
]
=
-
0.61550639239553132e-001
;
coefweinberg
[
5
][
5
]
=
0.18017109270959387e-001
;
coefweinberg
[
6
][
5
]
=
-
0.15438891584271634e-002
;
z
=
Redshift
;
J0
=
J_0
();
/* here, we initialize the ctes that uses J_nu(z) */
/* Tessier */
/*
Cte_G_gHI = G_gHI();
Cte_G_gHeI = G_gHeI();
Cte_G_gHeII = G_gHeII();
Cte_heating_radiative_HI = heating_radiative_HI();
Cte_heating_radiative_HeI = heating_radiative_HeI();
Cte_heating_radiative_HeII = heating_radiative_HeII();
*/
/* Theuns */
/*
Cte_G_gHI = G_gHI_t(J0);
Cte_G_gHeI = G_gHeI_t(J0);
Cte_G_gHeII = G_gHeII_t(J0);
Cte_heating_radiative_HI = heating_radiative_HI_t(J0);
Cte_heating_radiative_HeI = heating_radiative_HeI_t(J0);
Cte_heating_radiative_HeII = heating_radiative_HeII_t(J0);
*/
/* Weinberg */
Cte_G_gHI
=
G_gHI_w
();
Cte_G_gHeI
=
G_gHeI_w
();
Cte_G_gHeII
=
G_gHeII_w
();
Cte_heating_radiative_HI
=
heating_radiative_HI_w
();
Cte_heating_radiative_HeI
=
heating_radiative_HeI_w
();
Cte_heating_radiative_HeII
=
heating_radiative_HeII_w
();
}
/*
* J0
*/
double
J_0
()
{
double
Fz
;
if
(
z
>
6
)
Fz
=
0
;
else
{
if
(
z
>
3
)
Fz
=
4
/
(
z
+
1
);
else
{
if
(
z
>
2
)
Fz
=
1
;
else
Fz
=
pow
(((
1
+
z
)
/
3.
),
3
);
}
}
return
1.0e-22
*
Fz
;
}
/*
* UV background intensity
*/
double
J_nu
(
double
e
)
{
double
e_L
;
e_L
=
13.598
*
eV
;
return
(
e_L
/
e
)
*
J_0
();
}
/*
* sigma_rad
*/
double
sigma_rad_HI
(
double
e
)
{
double
xxx
,
alph
,
e_i
;
e_i
=
13.598
*
eV
;
xxx
=
e
/
e_i
;
alph
=
sqrt
(
xxx
-
1.0
);
return
6.30e-18
/
pow
(
xxx
,
4
)
*
exp
(
4.0
-
4.0
*
atan
(
alph
)
/
alph
)
/
(
1.0
-
exp
(
-
TWOPI
/
alph
));
}
double
sigma_rad_HeI
(
double
e
)
{
double
xxx
,
alph
,
e_i
;
e_i
=
24.587
*
eV
;
xxx
=
e
/
e_i
;
alph
=
sqrt
(
xxx
-
1.0
);
return
7.42e-18
*
(
1.660
/
pow
(
xxx
,
2.050
)
-
0.660
/
pow
(
xxx
,
3.050
));
}
double
sigma_rad_HeII
(
double
e
)
{
double
xxx
,
alph
,
e_i
;
e_i
=
54.416
*
eV
;
xxx
=
e
/
e_i
;
alph
=
sqrt
(
xxx
-
1.0
);
return
1.58e-18
/
pow
(
xxx
,
4
)
*
exp
(
4.0
-
4.0
*
atan
(
alph
)
/
alph
)
/
(
1.0
-
exp
(
-
TWOPI
/
alph
));
}
/*
* cooling rates
*/
/* Bremstrahlung */
double
cooling_bremstrahlung_HI
(
double
T
)
{
return
1.42e-27
*
sqrt
(
T
)
*
(
1.10
+
0.340
*
exp
(
-
pow
((
5.50
-
log10
(
T
)),
2
)
/
3.0
));
}
double
cooling_bremstrahlung_HeI
(
double
T
)
{
return
1.42e-27
*
sqrt
(
T
)
*
(
1.10
+
0.340
*
exp
(
-
pow
((
5.50
-
log10
(
T
)),
2
)
/
3.0
));
}
double
cooling_bremstrahlung_HeII
(
double
T
)
{
return
5.68e-27
*
sqrt
(
T
)
*
(
1.10
+
0.340
*
exp
(
-
pow
((
5.50
-
log10
(
T
)),
2
)
/
3.0
));
}
/* Ionization */
double
cooling_ionization_HI
(
double
T
)
{
double
T5
;
T5
=
T
/
1e5
;
return
2.54e-21
*
sqrt
(
T
)
*
exp
(
-
157809.1
/
T
)
/
(
1
+
sqrt
(
T5
));
}
double
cooling_ionization_HeI
(
double
T
)
{
double
T5
;
T5
=
T
/
1e5
;
return
1.88e-21
*
sqrt
(
T
)
*
exp
(
-
285335.4
/
T
)
/
(
1
+
sqrt
(
T5
));
}
double
cooling_ionization_HeII
(
double
T
)
{
double
T5
;
T5
=
T
/
1e5
;
return
9.90e-22
*
sqrt
(
T
)
*
exp
(
-
631515.0
/
T
)
/
(
1
+
sqrt
(
T5
));
}
/* Recombination */
double
cooling_recombination_HI
(
double
T
)
{
double
T3
,
T6
;
T3
=
T
/
1e3
;
T6
=
T
/
1e6
;
return
8.70e-27
*
sqrt
(
T
)
/
pow
(
T3
,
0.2
)
/
(
1.0
+
pow
(
T6
,
0.7
));
}
double
cooling_recombination_HeI
(
double
T
)
{
return
1.55e-26
*
pow
(
T
,
0.3647
);
}
double
cooling_recombination_HeII
(
double
T
)
{
double
T3
,
T6
;
T3
=
T
/
1e3
;
T6
=
T
/
1e6
;
return
3.48e-26
*
sqrt
(
T
)
/
pow
(
T3
,
0.2
)
/
(
1.0
+
pow
(
T6
,
0.7
));
}
/* Dielectric Recombination */
double
cooling_dielectric_recombination
(
double
T
)
{
return
1.24e-13
*
pow
(
T
,
-
1.5
)
*
exp
(
-
470000.0
/
T
)
*
(
1.0
+
0.3
*
exp
(
-
94000.0
/
T
));
}
/* Ecitation cooling (line cooling) */
double
cooling_excitation_HI
(
double
T
)
{
double
T5
;
T5
=
T
/
1e5
;
return
7.50e-19
*
exp
(
-
118348.0
/
T
)
/
(
1
+
sqrt
(
T5
));
}
double
cooling_excitation_HII
(
double
T
)
{
double
T5
;
T5
=
T
/
1e5
;
return
5.54e-17
/
pow
(
T
,
0.397
)
*
exp
(
-
473638.0
/
T
)
/
(
1
+
sqrt
(
T5
));
}
/* Compton cooling */
double
cooling_compton
(
double
T
)
{
return
5.406e-36
*
(
T
-
2.7
*
(
1
+
z
))
*
pow
((
1
+
z
),
4
);
}
/*
* recombination rates (taux_rec)
*/
double
A_HII
(
double
T
)
{
double
T3
,
T6
;
T3
=
T
/
1e3
;
T6
=
T
/
1e6
;
return
6.30e-11
/
sqrt
(
T
)
/
pow
(
T3
,
0.2
)
/
(
1
+
pow
(
T6
,
0.7
));
}
double
A_HeIId
(
double
T
)
{
return
1.9e-3
/
pow
(
T
,
1.50
)
*
exp
(
-
470000.0
/
T
)
*
(
1.0
+
0.30
*
exp
(
-
94000.0
/
T
));
}
double
A_HeII
(
double
T
)
{
return
1.5e-10
/
pow
(
T
,
0.6353
)
+
A_HeIId
(
T
);
}
double
A_HeIII
(
double
T
)
{
double
T3
,
T6
;
T3
=
T
/
1e3
;
T6
=
T
/
1e6
;
return
3.36e-10
/
sqrt
(
T
)
/
pow
(
T3
,
0.2
)
/
(
1.0
+
pow
(
T6
,
0.7
));
}
/*
* collisional rates (taux_ion)
*/
double
G_HI
(
double
T
)
{
double
T5
;
T5
=
T
/
1e5
;
return
1.17e-10
*
sqrt
(
T
)
*
exp
(
-
157809.1
/
T
)
/
(
1.0
+
sqrt
(
T5
));
}
double
G_HeI
(
double
T
)
{
double
T5
;
T5
=
T
/
1e5
;
return
2.38e-11
*
sqrt
(
T
)
*
exp
(
-
285335.4
/
T
)
/
(
1.0
+
sqrt
(
T5
));
}
double
G_HeII
(
double
T
)
{
double
T5
;
T5
=
T
/
1e5
;
return
5.68e-12
*
sqrt
(
T
)
*
exp
(
-
631515.0
/
T
)
/
(
1.0
+
sqrt
(
T5
));
}
/*
* photoionisation rates (depend only on z)
*/
double
G_gHI
()
{
double
e_i
,
integ
,
e
,
de
,
error
;
e_i
=
13.598
*
eV
;
integ
=
0.0
;
e
=
e_i
;
de
=
e
/
100.0
;
error
=
1.0
;
while
(
error
>
1.e-6
)
{
e
=
e
+
de
;
de
=
e
/
100.0
;
error
=
2
*
TWOPI
*
J_nu
(
e
)
*
sigma_rad_HI
(
e
)
*
de
/
e
;
integ
=
integ
+
error
;
error
=
error
/
fabs
(
integ
);
}
return
integ
/
PLANCK
;
}
double
G_gHeI
()
{
double
e_i
,
integ
,
e
,
de
,
error
;
e_i
=
24.587
*
eV
;
integ
=
0.0
;
e
=
e_i
;
de
=
e
/
100.0
;
error
=
1.0
;
while
(
error
>
1.e-6
)
{
e
=
e
+
de
;
de
=
e
/
100.0
;
error
=
2
*
TWOPI
*
J_nu
(
e
)
*
sigma_rad_HeI
(
e
)
*
de
/
e
;
integ
=
integ
+
error
;
error
=
error
/
fabs
(
integ
);
}
return
integ
/
PLANCK
;
}
double
G_gHeII
()
{
double
e_i
,
integ
,
e
,
de
,
error
;
e_i
=
54.416
*
eV
;
integ
=
0.0
;
e
=
e_i
;
de
=
e
/
100.0
;
error
=
1.0
;
while
(
error
>
1.e-6
)
{
e
=
e
+
de
;
de
=
e
/
100.0
;
error
=
2
*
TWOPI
*
J_nu
(
e
)
*
sigma_rad_HeII
(
e
)
*
de
/
e
;
integ
=
integ
+
error
;
error
=
error
/
fabs
(
integ
);
}
return
integ
/
PLANCK
;
}
double
G_gHI_t
(
double
J0
)
{
return
1.26e10
*
J0
/
(
3.0
+
alpha
);
}
double
G_gHeI_t
(
double
J0
)
{
return
1.48e10
*
J0
*
pow
(
0.5530
,
alpha
)
*
(
1.660
/
(
alpha
+
2.050
)
-
0.660
/
(
alpha
+
3.050
));
}
double
G_gHeII_t
(
double
J0
)
{
return
3.34e9
*
J0
*
pow
(
0.2490
,
alpha
)
/
(
3.0
+
alpha
);
}
double
G_gHI_w
()
{
double
taux_rad_weinbergint
;
double
hh
,
tt
,
zz
;
int
i
;
if
(
z
<
8.50
)
{
hh
=
0.0
;
zz
=
dmax
(
z
,
1.0e-15
);
for
(
i
=
0
;
i
<
Norderweinberg
;
i
++
)
hh
=
hh
+
coefweinberg
[
i
][
0
]
*
pow
(
zz
,
i
);
taux_rad_weinbergint
=
normfacJ0
*
exp
(
hh
);
}
else
taux_rad_weinbergint
=
0.0
;
tt
=
G_gHI_t
(
J0min
);
if
(
taux_rad_weinbergint
<
tt
)
taux_rad_weinbergint
=
tt
;
return
taux_rad_weinbergint
;
}
double
G_gHeI_w
()
{
double
taux_rad_weinbergint
;
double
hh
,
tt
,
zz
;
int
i
;
if
(
z
<
8.50
)
{
hh
=
0.0
;
zz
=
dmax
(
z
,
1.0e-15
);
for
(
i
=
0
;
i
<
Norderweinberg
;
i
++
)
hh
=
hh
+
coefweinberg
[
i
][
1
]
*
pow
(
zz
,
i
);
taux_rad_weinbergint
=
normfacJ0
*
exp
(
hh
);
}
else
taux_rad_weinbergint
=
0.0
;
tt
=
G_gHeI_t
(
J0min
);
if
(
taux_rad_weinbergint
<
tt
)
taux_rad_weinbergint
=
tt
;
return
taux_rad_weinbergint
;
}
double
G_gHeII_w
()
{
double
taux_rad_weinbergint
;
double
hh
,
tt
,
zz
;
int
i
;
if
(
z
<
8.50
)
{
hh
=
0.0
;
zz
=
dmax
(
z
,
1.0e-15
);
for
(
i
=
0
;
i
<
Norderweinberg
;
i
++
)
hh
=
hh
+
coefweinberg
[
i
][
2
]
*
pow
(
zz
,
i
);
taux_rad_weinbergint
=
normfacJ0
*
exp
(
hh
);
}
else
taux_rad_weinbergint
=
0.0
;
tt
=
G_gHeII_t
(
J0min
);
if
(
taux_rad_weinbergint
<
tt
)
taux_rad_weinbergint
=
tt
;
return
taux_rad_weinbergint
;
}
/*
* heating rates (depend only on z)
*/
double
heating_radiative_HI
()
/* use J_nu */
{
double
e_i
,
integ
,
e
,
de
,
error
;
e_i
=
13.598
*
eV
;
integ
=
0.0
;
e
=
e_i
;
de
=
e
/
100.0
;
error
=
1.0
;
while
(
error
>
1.e-6
)
{
e
=
e
+
de
;
de
=
e
/
100.0
;
error
=
2.0
*
TWOPI
*
J_nu
(
e
)
*
sigma_rad_HI
(
e
)
*
(
e
/
e_i
-
1.0
)
*
de
/
e
;
integ
=
integ
+
error
;
error
=
error
/
fabs
(
integ
);
}
return
integ
/
PLANCK
*
e_i
;
}
double
heating_radiative_HeI
()
/* use J_nu */
{
double
e_i
,
integ
,
e
,
de
,
error
;
e_i
=
24.587
*
eV
;
integ
=
0.0
;
e
=
e_i
;
de
=
e
/
100.0
;
error
=
1.0
;
while
(
error
>
1.e-6
)
{
e
=
e
+
de
;
de
=
e
/
100.0
;
error
=
2.0
*
TWOPI
*
J_nu
(
e
)
*
sigma_rad_HeI
(
e
)
*
(
e
/
e_i
-
1.0
)
*
de
/
e
;
integ
=
integ
+
error
;
error
=
error
/
fabs
(
integ
);
}
return
integ
/
PLANCK
*
e_i
;
}
double
heating_radiative_HeII
()
/* use J_nu */
{
double
e_i
,
integ
,
e
,
de
,
error
;
e_i
=
54.416
*
eV
;
integ
=
0.0
;
e
=
e_i
;
de
=
e
/
100.0
;
error
=
1.0
;
while
(
error
>
1.e-6
)
{
e
=
e
+
de
;
de
=
e
/
100.0
;
error
=
2.0
*
TWOPI
*
J_nu
(
e
)
*
sigma_rad_HeII
(
e
)
*
(
e
/
e_i
-
1.0
)
*
de
/
e
;
integ
=
integ
+
error
;
error
=
error
/
fabs
(
integ
);
}
return
integ
/
PLANCK
*
e_i
;
}
double
heating_radiative_HI_t
(
double
J0
)
/* use Theuns */
{
return
(
2.91e-1
*
J0
/
(
2.0
+
alpha
))
/
(
3.0
+
alpha
);
}
double
heating_radiative_HeI_t
(
double
J0
)
/* use Theuns */
{
return
5.84e-1
*
J0
*
pow
(
0.5530
,
alpha
)
*
(
1.660
/
(
alpha
+
1.050
)
-
2.320
/
(
alpha
+
2.050
)
+
0.660
/
(
alpha
+
3.050
));
}
double
heating_radiative_HeII_t
(
double
J0
)
/* use Theuns */
{
return
(
2.92e-1
*
J0
*
pow
(
0.2490
,
alpha
)
/
(
2.0
+
alpha
))
/
(
3.0
+
alpha
);
}
double
heating_radiative_HI_w
()
/* use weinberg coeff */
{
double
heat_rad_weinbergint
;
double
hh
,
tt
,
zz
;
int
i
;
if
(
z
<
8.50
)
{
hh
=
0.0
;
zz
=
dmax
(
z
,
1.0e-15
);
for
(
i
=
0
;
i
<
Norderweinberg
;
i
++
)
hh
=
hh
+
coefweinberg
[
i
][
3
]
*
pow
(
zz
,
i
);
heat_rad_weinbergint
=
normfacJ0
*
exp
(
hh
);
}
else
heat_rad_weinbergint
=
0.0
;
tt
=
heating_radiative_HI_t
(
J0min
);
if
(
heat_rad_weinbergint
<
tt
)
heat_rad_weinbergint
=
tt
;
return
heat_rad_weinbergint
;
}
double
heating_radiative_HeI_w
()
/* use weinberg coeff */
{
double
heat_rad_weinbergint
;
double
hh
,
tt
,
zz
;
int
i
;
if
(
z
<
8.50
)
{
hh
=
0.0
;
zz
=
dmax
(
z
,
1.0e-15
);
for
(
i
=
0
;
i
<
Norderweinberg
;
i
++
)
hh
=
hh
+
coefweinberg
[
i
][
4
]
*
pow
(
zz
,
i
);
heat_rad_weinbergint
=
normfacJ0
*
exp
(
hh
);
}
else
heat_rad_weinbergint
=
0.0
;
tt
=
heating_radiative_HeI_t
(
J0min
);
if
(
heat_rad_weinbergint
<
tt
)
heat_rad_weinbergint
=
tt
;
return
heat_rad_weinbergint
;
}
double
heating_radiative_HeII_w
()
/* use weinberg coeff */
{
double
heat_rad_weinbergint
;
double
hh
,
tt
,
zz
;
int
i
;
if
(
z
<
8.50
)
{
hh
=
0.0
;
zz
=
dmax
(
z
,
1.0e-15
);
for
(
i
=
0
;
i
<
Norderweinberg
;
i
++
)
hh
=
hh
+
coefweinberg
[
i
][
5
]
*
pow
(
zz
,
i
);
heat_rad_weinbergint
=
normfacJ0
*
exp
(
hh
);
}
else
heat_rad_weinbergint
=
0.0
;
tt
=
heating_radiative_HeII_t
(
J0min
);
if
(
heat_rad_weinbergint
<
tt
)
heat_rad_weinbergint
=
tt
;
return
heat_rad_weinbergint
;
}
double
heating_compton
()
{
/* Abel, Tom; Haehnelt, Martin G.Apj 520 */
//return 5.406e-36*2.726*pow((1+z),5); /* from Ramses */
//if (z>6)
// return 0;
//else
// return 1.25e-31*pow((1+z),13/3.);
return
0
;
}
void
compute_densities
(
double
T
,
double
X
,
double
*
pn_H
,
double
*
pn_HI
,
double
*
pn_HII
,
double
*
pn_HEI
,
double
*
pn_HEII
,
double
*
pn_HEIII
,
double
*
pn_E
,
double
*
pmu
)
{
double
Y
,
yy
,
x1
;
double
t_rad_HI
,
t_rec_HI
,
t_ion_HI
;
double
t_rad_HEI
,
t_rec_HEI
,
t_ion_HEI
;
double
t_rad_HEII
,
t_rec_HEII
,
t_ion_HEII
;
double
t_ion2_HI
,
t_ion2_HEI
,
t_ion2_HEII
;
double
err_nE
;
double
n_T
;
double
n_H
,
n_HI
,
n_HII
,
n_HEI
,
n_HEII
,
n_HEIII
,
n_E
,
mu
;
Y
=
1
-
X
;
yy
=
Y
/
(
4
-
4
*
Y
);
t_rad_HI
=
Cte_G_gHI
;
t_rec_HI
=
A_HII
(
T
);
t_ion_HI
=
G_HI
(
T
);
t_rad_HEI
=
Cte_G_gHeI
;
t_rec_HEI
=
A_HeII
(
T
);
t_ion_HEI
=
G_HeI
(
T
);
t_rad_HEII
=
Cte_G_gHeII
;
t_rec_HEII
=
A_HeIII
(
T
);
t_ion_HEII
=
G_HeII
(
T
);
n_H
=
*
pn_H
;
n_E
=
n_H
;
err_nE
=
1.
;
while
(
err_nE
>
1.e-8
)
{
/* compute densities (Ramses implementation) */
t_ion2_HI
=
t_ion_HI
+
t_rad_HI
/
dmax
(
n_E
,
1e-15
*
n_H
);
t_ion2_HEI
=
t_ion_HEI
+
t_rad_HEI
/
dmax
(
n_E
,
1e-15
*
n_H
);
t_ion2_HEII
=
t_ion_HEII
+
t_rad_HEII
/
dmax
(
n_E
,
1e-15
*
n_H
);
n_HI
=
t_rec_HI
/
(
t_ion2_HI
+
t_rec_HI
)
*
n_H
;
n_HII
=
t_ion2_HI
/
(
t_ion2_HI
+
t_rec_HI
)
*
n_H
;
x1
=
(
t_rec_HEII
*
t_rec_HEI
+
t_ion2_HEI
*
t_rec_HEII
+
t_ion2_HEII
*
t_ion2_HEI
);
n_HEIII
=
yy
*
t_ion2_HEII
*
t_ion2_HEI
/
x1
*
n_H
;
n_HEII
=
yy
*
t_ion2_HEI
*
t_rec_HEII
/
x1
*
n_H
;
n_HEI
=
yy
*
t_rec_HEII
*
t_rec_HEI
/
x1
*
n_H
;
err_nE
=
fabs
((
n_E
-
(
n_HII
+
n_HEII
+
2.
*
n_HEIII
))
/
n_H
);
n_E
=
0.5
*
n_E
+
0.5
*
(
n_HII
+
n_HEII
+
2.
*
n_HEIII
);
}
n_T
=
n_HI
+
n_HII
+
n_HEI
+
n_HEII
+
n_HEIII
+
n_E
;
mu
=
n_H
/
X
/
n_T
;
*
pn_H
=
n_H
;
*
pn_HI
=
n_HI
;
*
pn_HII
=
n_HII
;
*
pn_HEI
=
n_HEI
;
*
pn_HEII
=
n_HEII
;
*
pn_HEIII
=
n_HEIII
;
*
pn_E
=
n_E
;
*
pmu
=
mu
;
}
void
print_cooling
(
double
T
,
double
c1
,
double
c2
,
double
c3
,
double
c4
,
double
c5
,
double
c6
,
double
c7
,
double
c8
,
double
c9
,
double
c10
,
double
c11
,
double
c12
,
double
c13
,
double
h1
,
double
h2
,
double
h3
,
double
h4
)
{
double
ctot
,
htot
,
chtot
;
ctot
=
c1
+
c2
+
c3
+
c4
+
c5
+
c6
+
c7
+
c8
+
c9
+
c10
+
c11
+
c12
+
c13
;
htot
=
h1
+
h2
+
h3
+
h4
;
chtot
=
ctot
-
htot
;
printf
(
"%g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g
\n
"
,
T
,
c1
,
c2
,
c3
,
c4
,
c5
,
c6
,
c7
,
c8
,
c9
,
c10
,
c11
,
c12
,
c13
,
h1
,
h2
,
h3
,
h4
,
ctot
,
htot
,
chtot
);
}
void
compute_cooling_from_T_and_Nh
(
double
T
,
double
X
,
double
n_H
,
double
*
c1
,
double
*
c2
,
double
*
c3
,
double
*
c4
,
double
*
c5
,
double
*
c6
,
double
*
c7
,
double
*
c8
,
double
*
c9
,
double
*
c10
,
double
*
c11
,
double
*
c12
,
double
*
c13
,
double
*
h1
,
double
*
h2
,
double
*
h3
,
double
*
h4
)
{
double
n_HI
,
n_HII
,
n_HEI
,
n_HEII
,
n_HEIII
,
n_E
,
mu
;
double
nH2
;
//double c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13;
//double h1,h2,h3,h4;
compute_densities
(
T
,
X
,
&
n_H
,
&
n_HI
,
&
n_HII
,
&
n_HEI
,
&
n_HEII
,
&
n_HEIII
,
&
n_E
,
&
mu
);
nH2
=
n_H
*
n_H
;
/*
* compute cooling
*/
/* Bremstrahlung (cool_bre) */
*
c1
=
cooling_bremstrahlung_HI
(
T
)
*
n_E
*
n_HII
/
nH2
;
*
c2
=
cooling_bremstrahlung_HeI
(
T
)
*
n_E
*
n_HEII
/
nH2
;
*
c3
=
cooling_bremstrahlung_HeII
(
T
)
*
n_E
*
n_HEIII
/
nH2
;
/* Ionization cooling (cool_ion) */
*
c4
=
cooling_ionization_HI
(
T
)
*
n_E
*
n_HI
/
nH2
;
*
c5
=
cooling_ionization_HeI
(
T
)
*
n_E
*
n_HEI
/
nH2
;
*
c6
=
cooling_ionization_HeII
(
T
)
*
n_E
*
n_HEII
/
nH2
;
/* Recombination cooling (cool_rec) */
*
c7
=
cooling_recombination_HI
(
T
)
*
n_E
*
n_HII
/
nH2
;
*
c8
=
cooling_recombination_HeI
(
T
)
*
n_E
*
n_HEII
/
nH2
;
*
c9
=
cooling_recombination_HeII
(
T
)
*
n_E
*
n_HEIII
/
nH2
;
/* Dielectric recombination cooling (cool_die) */
*
c10
=
cooling_dielectric_recombination
(
T
)
*
n_E
*
n_HEII
/
nH2
;
/* Line cooling (cool_exc) */
*
c11
=
cooling_excitation_HI
(
T
)
*
n_E
*
n_HI
/
nH2
;
*
c12
=
cooling_excitation_HII
(
T
)
*
n_E
*
n_HEII
/
nH2
;
/* Compton cooling (cool_com) */
*
c13
=
cooling_compton
(
T
)
*
n_E
/
nH2
;
/* !! dep on z */
/*
* compute heating
*/
/* Radiative heating (h_rad_spec) */
*
h1
=
Cte_heating_radiative_HI
*
n_HI
/
nH2
;
*
h2
=
Cte_heating_radiative_HeI
*
n_HEI
/
nH2
;
*
h3
=
Cte_heating_radiative_HeII
*
n_HEII
/
nH2
;
/* Compton heating (heat_com) */
*
h4
=
heating_compton
()
*
n_E
/
nH2
;
/* !! dep on z */
}
double
compute_cooling_from_Egyspec_and_Density
(
double
Egyspec
,
double
Density
,
double
*
MeanWeight
)
{
double
T
,
mu
,
n_H
;
double
n_HI
,
n_HII
,
n_HEI
,
n_HEII
,
n_HEIII
,
n_E
;
double
err_mu
,
mu_left
,
mu_right
,
mu_old
;
int
niter
;
double
c1
,
c2
,
c3
,
c4
,
c5
,
c6
,
c7
,
c8
,
c9
,
c10
,
c11
,
c12
,
c13
;
double
h1
,
h2
,
h3
,
h4
;
double
nH2
;
/* Hydrogen density (cgs) */
n_H
=
HYDROGEN_MASSFRAC
*
Density
/
PROTONMASS
;
/* itterate to find the right mu and T */
err_mu
=
1.
;
mu_left
=
0.5
;
mu_right
=
1.3
;
niter
=
0
;
while
(
(
err_mu
>
1.e-4
)
&&
(
niter
<=
50
)
)
{
mu_old
=
0.5
*
(
mu_left
+
mu_right
);
/* compute temperature */
T
=
GAMMA_MINUS1
*
mu_old
*
PROTONMASS
/
BOLTZMANN
*
Egyspec
;
/* compute all */
compute_densities
(
T
,
HYDROGEN_MASSFRAC
,
&
n_H
,
&
n_HI
,
&
n_HII
,
&
n_HEI
,
&
n_HEII
,
&
n_HEIII
,
&
n_E
,
&
mu
);
err_mu
=
(
mu
-
mu_old
)
/
mu_old
;
if
(
err_mu
>
0.
)
{
mu_left
=
0.5
*
(
mu_left
+
mu_right
);
mu_right
=
mu_right
;
}
else
{
mu_left
=
mu_left
;
mu_right
=
0.5
*
(
mu_left
+
mu_right
);
}
err_mu
=
fabs
(
err_mu
);
niter
=
niter
+
1
;
}
if
(
niter
>
50
)
printf
(
"ERROR : too many iterations."
);
*
MeanWeight
=
0.5
*
(
mu_left
+
mu_right
);
/* now, compute cooling */
nH2
=
n_H
*
n_H
;
/*
* compute cooling
*/
/* Bremstrahlung (cool_bre) */
c1
=
cooling_bremstrahlung_HI
(
T
)
*
n_E
*
n_HII
/
nH2
;
c2
=
cooling_bremstrahlung_HeI
(
T
)
*
n_E
*
n_HEII
/
nH2
;
c3
=
cooling_bremstrahlung_HeII
(
T
)
*
n_E
*
n_HEIII
/
nH2
;
/* Ionization cooling (cool_ion) */
c4
=
cooling_ionization_HI
(
T
)
*
n_E
*
n_HI
/
nH2
;
c5
=
cooling_ionization_HeI
(
T
)
*
n_E
*
n_HEI
/
nH2
;
c6
=
cooling_ionization_HeII
(
T
)
*
n_E
*
n_HEII
/
nH2
;
/* Recombination cooling (cool_rec) */
c7
=
cooling_recombination_HI
(
T
)
*
n_E
*
n_HII
/
nH2
;
c8
=
cooling_recombination_HeI
(
T
)
*
n_E
*
n_HEII
/
nH2
;
c9
=
cooling_recombination_HeII
(
T
)
*
n_E
*
n_HEIII
/
nH2
;
/* Dielectric recombination cooling (cool_die) */
c10
=
cooling_dielectric_recombination
(
T
)
*
n_E
*
n_HEII
/
nH2
;
/* Line cooling (cool_exc) */
c11
=
cooling_excitation_HI
(
T
)
*
n_E
*
n_HI
/
nH2
;
c12
=
cooling_excitation_HII
(
T
)
*
n_E
*
n_HEII
/
nH2
;
/* Compton cooling (cool_com) */
c13
=
cooling_compton
(
T
)
*
n_E
/
nH2
;
/* !! dep on z */
/*
* compute heating
*/
/* Radiative heating (h_rad_spec) */
h1
=
Cte_heating_radiative_HI
*
n_HI
/
nH2
;
h2
=
Cte_heating_radiative_HeI
*
n_HEI
/
nH2
;
h3
=
Cte_heating_radiative_HeII
*
n_HEII
/
nH2
;
/* Compton heating (heat_com) */
h4
=
heating_compton
()
*
n_E
/
nH2
;
/* !! dep on z */
/* output info */
//print_cooling(T,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,h1,h2,h3,h4);
c1
=
dmax
(
c1
,
0
);
c2
=
dmax
(
c2
,
0
);
c3
=
dmax
(
c3
,
0
);
c4
=
dmax
(
c4
,
0
);
c5
=
dmax
(
c5
,
0
);
c6
=
dmax
(
c6
,
0
);
c7
=
dmax
(
c7
,
0
);
c8
=
dmax
(
c8
,
0
);
c9
=
dmax
(
c9
,
0
);
c10
=
dmax
(
c10
,
0
);
c11
=
dmax
(
c11
,
0
);
c12
=
dmax
(
c12
,
0
);
c13
=
dmax
(
c13
,
0
);
h1
=
dmax
(
h1
,
0
);
h2
=
dmax
(
h2
,
0
);
h3
=
dmax
(
h3
,
0
);
h4
=
dmax
(
h4
,
0
);
return
(
c1
+
c2
+
c3
+
c4
+
c5
+
c6
+
c7
+
c8
+
c9
+
c10
+
c11
+
c12
+
c13
)
-
(
h1
+
h2
+
h3
+
h4
);
}
struct
cooling_solver_params
{
double
Entropy
;
double
Density
;
int
Phase
;
int
i
;
double
DtEntropyVisc
;
double
dt
;
double
hubble_a
;
};
double
cooling_solver_function
(
double
EntropyVar
,
void
*
params
)
{
struct
cooling_solver_params
*
p
=
(
struct
cooling_solver_params
*
)
params
;
double
Entropy
=
p
->
Entropy
;
double
Density
=
p
->
Density
;
int
Phase
=
p
->
Phase
;
int
i
=
p
->
i
;
double
DtEntropyVisc
=
p
->
DtEntropyVisc
;
double
dt
=
p
->
dt
;
double
hubble_a
=
p
->
hubble_a
;
double
DtEntropyRadSph
=
0
;
#ifdef MULTIPHASE
switch
(
Phase
)
{
case
GAS_SPH
:
DtEntropyRadSph
=
-
GAMMA_MINUS1
*
pow
(
Density
,
-
GAMMA
)
*
lambda
(
Density
,
EntropyVar
,
Phase
,
i
)
/
hubble_a
;
break
;
case
GAS_STICKY
:
case
GAS_DARK
:
DtEntropyRadSph
=
-
1
/
(
Density
*
a3inv
)
*
lambda
(
Density
,
EntropyVar
,
Phase
,
i
)
/
hubble_a
;
break
;
}
#else
DtEntropyRadSph
=
-
GAMMA_MINUS1
*
pow
(
Density
,
-
GAMMA
)
*
lambda
(
Density
,
EntropyVar
,
Phase
,
i
)
/
hubble_a
;
#endif
return
Entropy
+
(
DtEntropyVisc
+
DtEntropyRadSph
)
*
dt
-
EntropyVar
;
};
/*! This function compute the new Entropy due to isochoric cooling
* using an implicit iteration scheme
*
* !!! here Density is already expressed in comobile coord
*
*/
double
DoCooling
(
FLOAT
Density
,
FLOAT
Entropy
,
int
Phase
,
int
i
,
FLOAT
DtEntropyVisc
,
double
dt
,
double
hubble_a
)
{
double
EntropyNew
;
double
Entropy_lo
=
0
,
Entropy_hi
=
0
;
double
lo
,
hi
;
int
status
;
int
iter
=
0
;
int
max_iter
=
100
;
const
gsl_root_fsolver_type
*
T
;
gsl_root_fsolver
*
s
;
gsl_function
F
;
struct
cooling_solver_params
params
=
{(
double
)
Entropy
,(
double
)
Density
,(
int
)
Phase
,(
int
)
i
,(
double
)
DtEntropyVisc
,(
double
)
dt
,(
double
)
hubble_a
};
F
.
function
=
&
cooling_solver_function
;
F
.
params
=
&
params
;
T
=
gsl_root_fsolver_brent
;
s
=
gsl_root_fsolver_alloc
(
T
);
Entropy_lo
=
0.5
*
Entropy
;
Entropy_hi
=
1.1
*
Entropy
;
lo
=
cooling_solver_function
(
Entropy_lo
,
&
params
);
hi
=
cooling_solver_function
(
Entropy_hi
,
&
params
);
if
(
lo
*
hi
>
0
)
{
do
{
Entropy_hi
=
2
*
Entropy_hi
;
Entropy_lo
=
0.5
*
Entropy_lo
;
lo
=
cooling_solver_function
(
Entropy_lo
,
&
params
);
hi
=
cooling_solver_function
(
Entropy_hi
,
&
params
);
//printf("here, we need to iterate...\n");
}
while
(
lo
*
hi
>
0
);
}
gsl_root_fsolver_set
(
s
,
&
F
,
Entropy_lo
,
Entropy_hi
);
do
{
iter
++
;
status
=
gsl_root_fsolver_iterate
(
s
);
EntropyNew
=
gsl_root_fsolver_root
(
s
);
Entropy_lo
=
gsl_root_fsolver_x_lower
(
s
);
Entropy_hi
=
gsl_root_fsolver_x_upper
(
s
);
status
=
gsl_root_test_interval
(
Entropy_lo
,
Entropy_hi
,
0
,
0.001
);
}
while
(
status
==
GSL_CONTINUE
&&
iter
<
max_iter
);
gsl_root_fsolver_free
(
s
);
if
(
status
!=
GSL_SUCCESS
)
{
printf
(
"WARNING, HERE WE DO NOT CONVERGE...%g %g
\n
"
,
Entropy_lo
,
Entropy_hi
);
endrun
(
3737
);
}
return
EntropyNew
;
}
/*! This function computes the entropy variation due to the cooling.
* Cooling is computed only for sph active particles.
*/
void
cooling
()
{
int
i
;
double
dt
=
0
;
double
EntropyNew
;
/* set the right Redshift and compute value indep of Temperature */
if
(
All
.
CoolingType
==
1
)
init_from_new_redshift
(
1.0
/
(
All
.
Time
)
-
1
);
if
(
All
.
ComovingIntegrationOn
)
{
hubble_a
=
All
.
Omega0
/
(
All
.
Time
*
All
.
Time
*
All
.
Time
)
+
(
1
-
All
.
Omega0
-
All
.
OmegaLambda
)
/
(
All
.
Time
*
All
.
Time
)
+
All
.
OmegaLambda
;
hubble_a
=
All
.
Hubble
*
sqrt
(
hubble_a
);
a3inv
=
1
/
(
All
.
Time
*
All
.
Time
*
All
.
Time
);
}
else
a3inv
=
hubble_a
=
1
;
for
(
i
=
0
;
i
<
NumPart
;
i
++
)
{
if
(
P
[
i
].
Ti_endstep
==
All
.
Ti_Current
)
/* active particles */
{
if
(
P
[
i
].
Type
==
0
)
/* SPH stuff */
{
//DtEntropyRadSph=0.;
//SphP[i].DtEntropyRadSph = 0;
#ifdef MULTIPHASE
if
(
SphP
[
i
].
Phase
==
GAS_SPH
)
{
#endif
/* note : SphP[i].DtEntropyRadSph should not be necessary */
dt
=
(
All
.
Ti_Current
-
P
[
i
].
Ti_begstep
)
*
All
.
Timebase_interval
;
//SphP[i].DtEntropyRadSph = -GAMMA_MINUS1*pow(SphP[i].Density * a3inv,-GAMMA)*lambda(SphP[i].Density *a3inv,SphP[i].Entropy,0,i)/hubble_a;
//if (fabs((SphP[i].DtEntropyRadSph+SphP[i].DtEntropy)*dt) > 0.1*fabs(SphP[i].Entropy))
{
/* do implicit isochoric cooling */
EntropyNew
=
DoCooling
(
SphP
[
i
].
Density
*
a3inv
,
SphP
[
i
].
Entropy
,
0
,
i
,
SphP
[
i
].
DtEntropy
,
dt
,
hubble_a
);
if
(
dt
>
0
)
SphP
[
i
].
DtEntropy
=
(
EntropyNew
-
SphP
[
i
].
Entropy
)
/
dt
;
}
//else
{
// SphP[i].DtEntropy += SphP[i].DtEntropyRadSph;
}
//SphP[i].DtEgySpecRadSph = - 1/GAMMA_MINUS1 * pow(SphP[i].Density * a3inv,GAMMA_MINUS1) * (SphP[i].DtEntropyRadSph);
#ifdef MULTIPHASE
}
else
/* STICKY OR DARK */
{
//SphP[i].DtEntropyRadSph = -1/(SphP[i].Density * a3inv)*lambda(SphP[i].Density *a3inv,SphP[i].Entropy,SphP[i].Phase,i)/hubble_a;
//SphP[i].DtEntropy += SphP[i].DtEntropyRadSph;
/* do implicit isochoric cooling */
dt
=
(
All
.
Ti_Current
-
P
[
i
].
Ti_begstep
)
*
All
.
Timebase_interval
;
EntropyNew
=
DoCooling
(
SphP
[
i
].
Density
*
a3inv
,
SphP
[
i
].
Entropy
,
SphP
[
i
].
Phase
,
i
,
SphP
[
i
].
DtEntropy
,
dt
,
hubble_a
);
if
(
dt
>
0
)
SphP
[
i
].
DtEntropy
=
(
EntropyNew
-
SphP
[
i
].
Entropy
)
/
dt
;
/* !!! here, we do not take into account the energy variation !!! */
/* SphP[i].DtEgySpecRadSph = SphP[i].DtEntropy, no ?*/
}
#endif
/* finally sum to the entropy variation */
/* WARNING : we do not compute DtEntropy here, it is updated in timestep.c */
/* no, because, it is updated juste above, no ? */
//SphP[i].DtEntropy += SphP[i].DtEntropyRadSph;
//SphP[i].DtEntropy += SphP[i].DtEgySpecRadSph / (-1/GAMMA_MINUS1 * pow(SphP[i].Density * a3inv,GAMMA_MINUS1));
}
}
}
}
/*! This function computes the new entropy due to the cooling,
* between step t0 and t1.
*/
void
CoolingForOne
(
int
i
,
int
tstart
,
int
tend
,
double
a3inv
,
double
hubble_a
)
{
double
dt
,
dadt
,
tcool
,
dt_entr
;
double
MinSizeTimestep
,
ErrTolIntAccuracy
;
int
ti_current
,
istep
;
int
ti_step
;
double
minentropy
;
FLOAT
Entropy
,
DEntropy
,
DtEntropy
,
DtEgySpec
;
if
(
All
.
MinEgySpec
)
minentropy
=
All
.
MinEgySpec
*
GAMMA_MINUS1
/
pow
(
SphP
[
i
].
Density
*
a3inv
,
GAMMA_MINUS1
);
/* compute dt */
/* here we use the convention of Gadget */
/* this assume that DtEntropy = dA/dt/hubble_a */
/* and not only dA/dt */
dt_entr
=
(
tend
-
tstart
)
*
All
.
Timebase_interval
;
ErrTolIntAccuracy
=
0.02
;
MinSizeTimestep
=
0.01
*
dt_entr
;
/* compute da/dt */
//dadt = fabs( -GAMMA_MINUS1*pow(SphP[i].Density * a3inv,-GAMMA)*lambda(SphP[i].Density *a3inv,SphP[i].Entropy,0,i)/hubble_a );
/* compute cooling time */
//tcool = SphP[i].Entropy / dadt;
//if (ErrTolIntAccuracy*tcool/dt_entr < 1)
// printf("** %g %g\n",ErrTolIntAccuracy*tcool,dt_entr); /* --> verifier le cooling time */
/***************************************/
/* integrate with adaptative timesteps */
/***************************************/
Entropy
=
SphP
[
i
].
Entropy
;
ti_current
=
tstart
;
istep
=
0
;
#ifdef CHIMIE_THERMAL_FEEDBACK
int
no_cooling_SNII
,
no_cooling_SNIa
,
no_cooling
;
int
Tis
,
Tic
;
double
td
;
no_cooling
=
0
;
no_cooling_SNIa
=
0
;
no_cooling_SNII
=
0
;
/* check if we are in an adiabatic phase or not */
if
(
All
.
ComovingIntegrationOn
)
{
Tic
=
All
.
Ti_Current
;
if
(
SphP
[
i
].
SNIaThermalTime
>
0
)
/* only if the time has been set at least once, it is negative instead (see init.c) */
{
Tis
=
log
(
SphP
[
i
].
SNIaThermalTime
/
All
.
TimeBegin
)
/
All
.
Timebase_interval
;
td
=
get_cosmictime_difference
(
Tis
,
Tic
);
if
(
td
<
All
.
ChimieSNIaThermalTime
)
no_cooling_SNIa
=
1
;
}
if
(
SphP
[
i
].
SNIIThermalTime
>
0
)
/* only if the time has been set at least once, it is negative instead (see init.c) */
{
Tis
=
log
(
SphP
[
i
].
SNIIThermalTime
/
All
.
TimeBegin
)
/
All
.
Timebase_interval
;
td
=
get_cosmictime_difference
(
Tis
,
Tic
);
if
(
td
<
All
.
ChimieSNIIThermalTime
)
no_cooling_SNII
=
1
;
}
}
else
{
if
(
SphP
[
i
].
SNIaThermalTime
>
0
)
/* only if the time has been set at least once, it is negative instead (see init.c) */
if
((
All
.
Time
-
SphP
[
i
].
SNIaThermalTime
)
<
All
.
ChimieSNIaThermalTime
)
no_cooling_SNIa
=
1
;
if
(
SphP
[
i
].
SNIIThermalTime
>
0
)
/* only if the time has been set at least once, it is negative instead (see init.c) */
if
((
All
.
Time
-
SphP
[
i
].
SNIIThermalTime
)
<
All
.
ChimieSNIIThermalTime
)
no_cooling_SNII
=
1
;
}
no_cooling
=
no_cooling_SNIa
+
no_cooling_SNII
;
if
(
no_cooling
)
{
Entropy
=
Entropy
+
SphP
[
i
].
DtEntropy
*
dt_entr
;
/* avoid Entropy to be less than minentropy */
if
(
All
.
MinEgySpec
)
if
(
Entropy
<
minentropy
)
Entropy
=
minentropy
;
/* update particle */
SphP
[
i
].
DtEntropy
=
(
Entropy
-
SphP
[
i
].
Entropy
)
/
dt_entr
;
SphP
[
i
].
Entropy
=
Entropy
;
return
;
}
#endif
while
(
ti_current
<
tend
)
{
/* compute da/dt */
dadt
=
fabs
(
-
GAMMA_MINUS1
*
pow
(
SphP
[
i
].
Density
*
a3inv
,
-
GAMMA
)
*
lambda
(
SphP
[
i
].
Density
*
a3inv
,
Entropy
,
0
,
i
)
/
hubble_a
);
/* compute cooling time */
/* this is similar in comobile integraction */
tcool
=
Entropy
/
dadt
;
/* find dt */
dt
=
dmax
(
MinSizeTimestep
,
tcool
*
ErrTolIntAccuracy
);
dt
=
dmin
(
dt
,
dt_entr
);
ti_step
=
dt
/
All
.
Timebase_interval
;
ti_step
=
imax
(
1
,
ti_step
);
ti_step
=
imin
(
ti_step
,
tend
-
ti_current
);
dt
=
ti_step
*
All
.
Timebase_interval
;
#ifndef IMPLICIT_COOLING_INTEGRATION
/* normal integration of Entropy */
Entropy
+=
SphP
[
i
].
DtEntropy
*
dt
;
/* viscosity */
Entropy
+=
-
GAMMA_MINUS1
*
pow
(
SphP
[
i
].
Density
*
a3inv
,
-
GAMMA
)
*
lambda
(
SphP
[
i
].
Density
*
a3inv
,
Entropy
,
0
,
i
)
/
hubble_a
*
dt
;
/* cooling */
#else
/* or use implicit integration of Entropy */
/* need this if there is also heating like UV */
if
(
All
.
ComovingIntegrationOn
)
{
printf
(
"CoolingForOne : this must be checked !
\n
"
);
endrun
(
123321
);
}
Entropy
=
DoCooling
(
SphP
[
i
].
Density
*
a3inv
,
Entropy
,
0
,
i
,
SphP
[
i
].
DtEntropy
,
dt
,
hubble_a
);
#endif
/* avoid Entropy to be less than minentropy */
if
(
All
.
MinEgySpec
)
if
(
Entropy
<
minentropy
)
{
Entropy
=
minentropy
;
break
;
}
ti_current
+=
ti_step
;
istep
=
istep
+
1
;
}
/* entropy only due to cooling */
DEntropy
=
Entropy
-
SphP
[
i
].
Entropy
-
SphP
[
i
].
DtEntropy
*
dt_entr
;
DtEntropy
=
DEntropy
/
dt_entr
;
//if (istep>1)
// printf("* %d %d %d %g %g\n",istep,ti_current,tend,SphP[i].DtEntropy,DtEntropy);
/* update particle */
SphP
[
i
].
Entropy
=
Entropy
;
SphP
[
i
].
DtEntropy
=
SphP
[
i
].
DtEntropy
+
DtEntropy
;
DtEgySpec
=
-
1
/
GAMMA_MINUS1
*
pow
(
SphP
[
i
].
Density
*
a3inv
,
GAMMA_MINUS1
)
*
(
DtEntropy
);
LocalSysState
.
RadiatedEnergy
+=
DtEgySpec
*
dt_entr
*
P
[
i
].
Mass
;
#ifdef CHECK_ENTROPY_SIGN
if
(
SphP
[
i
].
Entropy
<
0
)
{
printf
(
"
\n
task=%d: entropy less than zero in CoolingForOne !
\n
"
,
ThisTask
);
printf
(
"ID=%d Entropy=%g EntropyPred=%g DtEntropy=%g
\n
"
,
P
[
i
].
ID
,
SphP
[
i
].
Entropy
,
SphP
[
i
].
EntropyPred
,
SphP
[
i
].
DtEntropy
);
fflush
(
stdout
);
endrun
(
444001
);
}
#endif
}
/*! cooling function
*
*/
double
lambda
(
FLOAT
Density
,
FLOAT
Entropy
,
int
phase
,
int
i
)
{
/*
* These function returns the Lambda (not the Lambda_n)
* Here, we assume that Lambda may also contain the heating term.
*
* Here, the Entropy and Density are physical, but in h units
*
* The function is used only in cooling.c
*
*/
double
EgySpec
;
double
MeanWeight
;
double
T
=
0
,
nH
=
0
,
nH2
=
0
,
l
=
0
;
double
nHcgs
=
0
,
nH2cgs
=
0
;
#ifdef HEATING
double
Gpe
=
0
;
double
X
,
XTne
,
eps
,
ne
,
flux_in_cgs
;
#endif
/* number of Hydrogen atoms per unit volume (user units, not corrected from h : [nH] = h2/cm^3 ) */
#ifndef DO_NO_USE_HYDROGEN_MASSFRAC_IN_COOLING
nH
=
HYDROGEN_MASSFRAC
*
Density
/
All
.
ProtonMass
;
#else
nH
=
1
*
Density
/
All
.
ProtonMass
;
#endif
nH2
=
nH
*
nH
;
/* in cgs, corrected from h */
nHcgs
=
nH
/
pow
(
All
.
UnitLength_in_cm
,
3
)
*
(
All
.
HubbleParam
*
All
.
HubbleParam
);
nH2cgs
=
nHcgs
*
nHcgs
;
/* compute temperature */
#ifdef MULTIPHASE
switch
(
phase
)
{
case
GAS_SPH
:
T
=
All
.
mumh
/
All
.
Boltzmann
*
Entropy
*
pow
(
Density
,
GAMMA_MINUS1
);
break
;
case
GAS_STICKY
:
case
GAS_DARK
:
T
=
All
.
mumh
/
All
.
Boltzmann
*
GAMMA_MINUS1
*
Entropy
;
break
;
}
#else
T
=
All
.
mumh
/
All
.
Boltzmann
*
Entropy
*
pow
(
Density
,
GAMMA_MINUS1
);
#endif
/*******************
* * * COOLING * * *
*******************/
if
(
All
.
CoolingType
==
0
||
All
.
CoolingType
==
2
)
{
/**************/
/* Sutherland */
/**************/
#ifdef MULTIPHASE
switch
(
phase
)
{
case
GAS_SPH
:
if
(
T
>
All
.
CutofCoolingTemperature
)
if
(
All
.
CoolingType
==
0
)
l
=
cooling_function
(
T
);
else
#ifdef CHIMIE
l
=
cooling_function_with_metals
(
T
,
SphP
[
i
].
Metal
[
FE
]);
#else
l
=
cooling_function_with_metals
(
T
,(
pow
(
10
,
All
.
InitGasMetallicity
)
-
1e-10
)
*
All
.
CoolingParameters_FeHSolar
);
#endif
else
l
=
0
;
break
;
case
GAS_STICKY
:
case
GAS_DARK
:
if
(
T
>
All
.
CutofCoolingTemperature
)
if
(
All
.
CoolingType
==
0
)
l
=
cooling_function
(
T
);
else
#ifdef CHIMIE
l
=
cooling_function_with_metals
(
T
,
SphP
[
i
].
Metal
[
FE
]);
#else
l
=
cooling_function_with_metals
(
T
,(
pow
(
10
,
All
.
InitGasMetallicity
)
-
1e-10
)
*
All
.
CoolingParameters_FeHSolar
);
#endif
else
l
=
0
;
break
;
}
#else
/* here, lambda' is in erg*cm^3/s = kg*m^5/s^3 */
if
(
T
>
All
.
CutofCoolingTemperature
)
if
(
All
.
CoolingType
==
0
)
l
=
cooling_function
(
T
);
else
#ifdef CHIMIE
l
=
cooling_function_with_metals
(
T
,
SphP
[
i
].
Metal
[
FE
]);
#else
l
=
cooling_function_with_metals
(
T
,(
pow
(
10
,
All
.
InitGasMetallicity
)
-
1e-10
)
*
All
.
CoolingParameters_FeHSolar
);
#endif
else
l
=
0
;
#endif
}
else
{
/******************************/
/* cooling with UV background */
/******************************/
/* get the right density and egyspec in cgs */
/* entropy and density are already physical */
#ifdef MULTIPHASE
/* WARNING, HERE, WE MUST DIFERENCIATE ACORDING TO THE PHASE... */
printf
(
"WARNING, HERE, WE MUST DIFERENCIATE ACORDING TO THE PHASE...
\n
"
);
exit
(
0
);
// if (phase == GAS_SPH)
// EgySpec = Entropy / GAMMA_MINUS1 * pow(Density, GAMMA_MINUS1);
// else
// EgySpec = Entropy;
#else
EgySpec
=
Entropy
/
GAMMA_MINUS1
*
pow
(
Density
,
GAMMA_MINUS1
);
#endif
/* into cgs, corrected from h */
EgySpec
*=
All
.
UnitEnergy_in_cgs
/
All
.
UnitMass_in_g
;
Density
*=
All
.
UnitDensity_in_cgs
;
//if(All.ComovingIntegrationOn)
// Density *= (All.HubbleParam*All.HubbleParam);
/* compute cooling from EnergySpec and Density */
l
=
compute_cooling_from_Egyspec_and_Density
(
EgySpec
,
Density
,
&
MeanWeight
);
/* compute temperature */
/*
Temperature = GAMMA_MINUS1 *MeanWeight*PROTONMASS/BOLTZMANN *EgySpec;
//printf("%g %g %g\n",Temperature,MeanWeight,Lambda);
logT = log10(Temperature);
*/
}
/*******************
* * * HEATING * * *
*******************/
#ifdef HEATING
#ifdef HEATING_PE
/**************************/
/* Photo-electric heating */
/* all must be in cgs */
/**************************/
X
=
0
;
#ifdef STELLAR_FLUX
flux_in_cgs
=
SphP
[
i
].
EnergyFlux
*
All
.
UnitEnergy_in_cgs
/
All
.
UnitTime_in_s
/
pow
(
All
.
UnitLength_in_cm
,
2
);
X
=
X
+
flux_in_cgs
/
C
/
All
.
HeatingPeSolarEnergyDensity
;
#endif
#ifdef EXTERNAL_FLUX
X
=
X
+
All
.
HeatingExternalFLuxEnergyDensity
/
All
.
HeatingPeSolarEnergyDensity
;
#endif
ne
=
nHcgs
*
All
.
HeatingPeElectronFraction
;
XTne
=
X
*
sqrt
(
T
)
/
ne
;
eps
=
4.87e-2
/
(
1
+
4e-3
*
pow
(
XTne
,
0.73
))
+
3.65e-2
*
(
T
/
1e4
)
/
(
1
+
2e-4
*
XTne
);
Gpe
=
(
1e-24
*
eps
*
X
*
nHcgs
)
/
nH2cgs
;
l
=
l
-
Gpe
;
#endif
/*HEATING_PE*/
#endif
/**********************************
* * * final unit conversions * * *
***********************************/
/* convert lambda' in user units */
l
=
l
/
All
.
UnitEnergy_in_cgs
/
pow
(
All
.
UnitLength_in_cm
,
3
)
*
All
.
UnitTime_in_s
;
/* in unit with h */
l
=
l
*
All
.
HubbleParam
;
/* correct from h */
/*
* [ Lambda / H / rho_p ] = [u] = cm^2/s^2
*
* [H] = h/s
* [rho_p] = g/cm^3 * h^2
* [Lambda_n] = g * m^5 / s^3
* [n] = h^2/m^5
*
* => Lambda_n must be multiplied by h (in order to remove one h !! not a unit !!)
*
*/
//if(All.ComovingIntegrationOn)
// l = l * All.HubbleParam;
/* get the final lambda by multiplying lambda' by nH2 (all in user units) */
l
=
l
*
nH2
;
return
l
;
}
#endif
Event Timeline
Log In to Comment