Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F87024257
cooling_with_metals.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
Thu, Oct 10, 01:49
Size
41 KB
Mime Type
text/x-c
Expires
Sat, Oct 12, 01:49 (2 d)
Engine
blob
Format
Raw Data
Handle
21524444
Attached To
rPNBODY pNbody
cooling_with_metals.c
View Options
#include <Python.h>
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <numpy/arrayobject.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
#define TO_DOUBLE(a) ( (PyArrayObject*) PyArray_CastToType(a, PyArray_DescrFromType(NPY_DOUBLE) ,0) )
/*
****************************************************
these variables are already defined in Gadget
****************************************************
*/
#ifdef DOUBLEPRECISION
/*!< If defined, the variable type FLOAT is set to "double", otherwise to FLOAT */
#define FLOAT double
#else
#define FLOAT float
#endif
#define MAXLEN_FILENAME 100
#define HYDROGEN_MASSFRAC 0.76
#define GAMMA (5.0/3)
#define GAMMA_MINUS1 (GAMMA-1)
/* Some physical constants in cgs units */
#define GRAVITY 6.672e-8
/*!< Gravitational constant (in cgs units) */
#define SOLAR_MASS 1.989e33
#define SOLAR_LUM 3.826e33
#define RAD_CONST 7.565e-15
#define AVOGADRO 6.0222e23
#define BOLTZMANN 1.3806e-16
#define GAS_CONST 8.31425e7
#define C 2.9979e10
#define PLANCK 6.6262e-27
#define CM_PER_MPC 3.085678e24
#define PROTONMASS 1.6726e-24
#define ELECTRONMASS 9.10953e-28
#define THOMPSON 6.65245e-25
#define ELECTRONCHARGE 4.8032e-10
#define HUBBLE 3.2407789e-18
/* in h/sec */
#define YEAR_IN_SECOND 31536000.0
/* year in sec */
#define METALLICITY_SOLAR 0.02
#define FEH_SOLAR 0.001771
/* 0.00181 */
#define MGH_SOLAR 0.00091245
#define PI 3.1415926535897931
#define TWOPI 6.2831853071795862
#define FE 0
struct
global_data_all_processes
{
double
InitGasMetallicity
;
char
CoolingFile
[
MAXLEN_FILENAME
];
double
*
logT
;
double
*
logL
;
gsl_interp_accel
*
acc_cooling_spline
;
gsl_spline
*
cooling_spline
;
/*
new metal dependent cooling
*/
double
CoolingParameters_zmin
;
double
CoolingParameters_zmax
;
double
CoolingParameters_slz
;
double
CoolingParameters_tmin
;
double
CoolingParameters_tmax
;
double
CoolingParameters_slt
;
double
CoolingParameters_FeHSolar
;
double
CoolingParameters_cooling_data_max
;
double
CoolingParameters_cooling_data
[
9
][
171
];
int
CoolingParameters_p
;
int
CoolingParameters_q
;
double
Boltzmann
;
double
ProtonMass
;
double
mumh
;
double
UnitLength_in_cm
;
double
UnitMass_in_g
;
double
UnitVelocity_in_cm_per_s
;
double
UnitTime_in_s
;
double
UnitEnergy_in_cgs
;
double
Timebase_interval
;
double
MinEgySpec
;
}
All
;
int
ThisTask
=
0
;
/****************************************************************************************/
/*
/*
/*
/* GADGET COOLING PART
/*
/*
/*
/****************************************************************************************/
/*! 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, 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
;
}
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
;
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
);
#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
;
}
/****************************************************************************************/
/*
/*
/*
/* OTHER C FUNCTIONS
/*
/*
/*
/****************************************************************************************/
double
set_default_parameters
()
{
double
meanweight
;
strcpy
(
All
.
CoolingFile
,
"/home/epfl/revaz/.pNbody/cooling_with_metals.dat"
);
All
.
UnitLength_in_cm
=
3.085e+21
;
All
.
UnitMass_in_g
=
1.989e+43
;
All
.
UnitVelocity_in_cm_per_s
=
20725573.785998672
;
All
.
Timebase_interval
=
0.001
;
All
.
MinEgySpec
=
0
;
}
double
set_parameters
(
PyObject
*
params_dict
)
{
PyObject
*
key
;
PyObject
*
value
;
Py_ssize_t
pos
=
0
;
int
ivalue
;
float
fvalue
;
double
dvalue
;
while
(
PyDict_Next
(
params_dict
,
&
pos
,
&
key
,
&
value
))
{
if
(
PyString_Check
(
key
))
{
if
(
strcmp
(
PyString_AsString
(
key
),
"UnitLength_in_cm"
)
==
0
)
if
(
PyInt_Check
(
value
)
||
PyLong_Check
(
value
)
||
PyFloat_Check
(
value
))
{
dvalue
=
PyFloat_AsDouble
(
value
);
All
.
UnitLength_in_cm
=
dvalue
;
}
if
(
strcmp
(
PyString_AsString
(
key
),
"UnitMass_in_g"
)
==
0
)
if
(
PyInt_Check
(
value
)
||
PyLong_Check
(
value
)
||
PyFloat_Check
(
value
))
{
dvalue
=
PyFloat_AsDouble
(
value
);
All
.
UnitMass_in_g
=
dvalue
;
}
if
(
strcmp
(
PyString_AsString
(
key
),
"UnitVelocity_in_cm_per_s"
)
==
0
)
if
(
PyInt_Check
(
value
)
||
PyLong_Check
(
value
)
||
PyFloat_Check
(
value
))
{
dvalue
=
PyFloat_AsDouble
(
value
);
All
.
UnitVelocity_in_cm_per_s
=
dvalue
;
}
if
(
strcmp
(
PyString_AsString
(
key
),
"CoolingFile"
)
==
0
)
if
(
PyString_Check
(
value
))
{
strcpy
(
All
.
CoolingFile
,
PyString_AsString
(
value
));
}
}
}
}
double
set_units
()
{
double
meanweight
;
All
.
UnitTime_in_s
=
All
.
UnitLength_in_cm
/
All
.
UnitVelocity_in_cm_per_s
;
All
.
UnitEnergy_in_cgs
=
All
.
UnitMass_in_g
*
pow
(
All
.
UnitLength_in_cm
,
2
)
/
pow
(
All
.
UnitTime_in_s
,
2
);
meanweight
=
4.0
/
(
1
+
3
*
HYDROGEN_MASSFRAC
);
/* note: we assume neutral gas here */
/*meanweight = 4 / (8 - 5 * (1 - HYDROGEN_MASSFRAC));*/
/* note: we assume FULL ionized gas here */
All
.
Boltzmann
=
BOLTZMANN
/
All
.
UnitEnergy_in_cgs
;
All
.
ProtonMass
=
PROTONMASS
/
All
.
UnitMass_in_g
;
All
.
mumh
=
All
.
ProtonMass
*
meanweight
;
}
double
lambda_from_Density_Entropy_FeH
(
FLOAT
Density
,
FLOAT
Entropy
,
FLOAT
FeH
)
{
/*
Density and Entropy in code units
return lambda in code units
The function corresponds to
double lambda(FLOAT Density,FLOAT Entropy,int phase,int i)
in Gadget
*/
double
nH
,
nH2
,
T
,
l
;
nH
=
HYDROGEN_MASSFRAC
*
Density
/
All
.
ProtonMass
;
nH2
=
nH
*
nH
;
T
=
All
.
mumh
/
All
.
Boltzmann
*
Entropy
*
pow
(
Density
,
GAMMA_MINUS1
);
l
=
cooling_function_with_metals
(
T
,
FeH
);
/* convert lambda' in user units */
l
=
l
/
All
.
UnitEnergy_in_cgs
/
pow
(
All
.
UnitLength_in_cm
,
3
)
*
All
.
UnitTime_in_s
;
l
=
l
*
nH2
;
return
l
;
}
double
lambda_from_Density_EnergyInt_FeH
(
FLOAT
Density
,
FLOAT
EnergyInt
,
FLOAT
FeH
)
{
/*
Density and EnergyInt in code units
return lambda in code units
*/
double
nH
,
nH2
,
T
,
l
;
nH
=
HYDROGEN_MASSFRAC
*
Density
/
All
.
ProtonMass
;
nH2
=
nH
*
nH
;
T
=
GAMMA_MINUS1
*
All
.
mumh
/
All
.
Boltzmann
*
EnergyInt
;
l
=
cooling_function_with_metals
(
T
,
FeH
);
/* convert lambda' in user units */
l
=
l
/
All
.
UnitEnergy_in_cgs
/
pow
(
All
.
UnitLength_in_cm
,
3
)
*
All
.
UnitTime_in_s
;
l
=
l
*
nH2
;
return
l
;
}
double
lambda_from_Density_Temperature_FeH
(
FLOAT
Density
,
FLOAT
Temperature
,
FLOAT
FeH
)
{
/*
Density in code units
return lambda in code units
*/
double
nH
,
nH2
,
T
,
l
;
nH
=
HYDROGEN_MASSFRAC
*
Density
/
All
.
ProtonMass
;
nH2
=
nH
*
nH
;
T
=
Temperature
;
l
=
cooling_function_with_metals
(
T
,
FeH
);
/* convert lambda' in user units */
l
=
l
/
All
.
UnitEnergy_in_cgs
/
pow
(
All
.
UnitLength_in_cm
,
3
)
*
All
.
UnitTime_in_s
;
l
=
l
*
nH2
;
return
l
;
}
double
lambda_normalized_from_Temperature_FeH
(
FLOAT
Temperature
,
FLOAT
FeH
)
{
/*
return lambda normalized (in mks)
*/
double
T
,
l
;
T
=
Temperature
;
l
=
cooling_function_with_metals
(
T
,
FeH
);
return
l
;
}
double
cooling_time_from_Density_Temperature_FeH
(
FLOAT
Density
,
FLOAT
Temperature
,
FLOAT
FeH
)
{
/*
Cooling time in code units
*/
double
u
,
l
,
dudt
,
tc
;
/* energy int */
u
=
Temperature
/
GAMMA_MINUS1
/
All
.
mumh
*
All
.
Boltzmann
;
/* lambda in user units */
l
=
lambda_from_Density_Temperature_FeH
(
Density
,
Temperature
,
FeH
);
dudt
=
l
/
Density
;
tc
=
u
/
dudt
;
return
tc
;
}
double
cooling_time_from_Density_EnergyInt_FeH
(
FLOAT
Density
,
FLOAT
EnergyInt
,
FLOAT
FeH
)
{
/*
Cooling time in code units
*/
double
u
,
l
,
dudt
,
tc
;
/* energy int */
u
=
EnergyInt
;
/* lambda in user units */
l
=
lambda_from_Density_EnergyInt_FeH
(
Density
,
EnergyInt
,
FeH
);
dudt
=
l
/
Density
;
tc
=
u
/
dudt
;
return
tc
;
}
/****************************************************************************************/
/*
/*
/*
/* NEW C FUNCTIONS FOR COOLING INTEGRATION
/*
/*
/*
/****************************************************************************************/
double
a3inv
=
1.
;
double
hubble_a
=
1
;
double
lambda
(
FLOAT
Density
,
FLOAT
Entropy
,
FLOAT
Metal
)
{
double
l
;
l
=
lambda_from_Density_Entropy_FeH
(
Density
,
Entropy
,
Metal
);
return
l
;
}
/*! returns the maximum of two double
*/
double
dmax
(
double
x
,
double
y
)
{
if
(
x
>
y
)
return
x
;
else
return
y
;
}
/*! returns the minimum of two double
*/
double
dmin
(
double
x
,
double
y
)
{
if
(
x
<
y
)
return
x
;
else
return
y
;
}
/*! returns the maximum of two integers
*/
int
imax
(
int
x
,
int
y
)
{
if
(
x
>
y
)
return
x
;
else
return
y
;
}
/*! returns the minimum of two integers
*/
int
imin
(
int
x
,
int
y
)
{
if
(
x
<
y
)
return
x
;
else
return
y
;
}
float
cooling1
(
int
tstart
,
int
tend
,
FLOAT
Entropy
,
FLOAT
Density
,
FLOAT
Metal
,
FLOAT
DtEntropy
,
FLOAT
*
newEntropy
,
FLOAT
*
newDtEntropy
)
{
double
dt
,
dadt
,
tcool
,
dt_entr
;
double
minentropy
;
double
MinSizeTimestep
,
ErrTolIntAccuracy
;
int
ti_current
,
istep
;
int
ti_step
;
FLOAT
oldEntropy
,
oldDtEntropy
,
DEntropy
;
double
lmbd
=
0
;
if
(
All
.
MinEgySpec
)
minentropy
=
All
.
MinEgySpec
*
GAMMA_MINUS1
/
pow
(
Density
*
a3inv
,
GAMMA_MINUS1
);
oldEntropy
=
Entropy
;
oldDtEntropy
=
DtEntropy
;
dt_entr
=
(
tend
-
tstart
)
*
All
.
Timebase_interval
;
ErrTolIntAccuracy
=
0.02
;
MinSizeTimestep
=
0.01
*
dt_entr
;
/***************************************/
/* integrate with adaptative timesteps */
/***************************************/
ti_current
=
tstart
;
istep
=
0
;
#ifdef CHIMIE_THERMAL_FEEDBACK
if
(
SphP
[
i
].
ThermalTime
>
(
All
.
Time
-
All
.
ChimieThermalTime
))
{
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
printf
(
"%g %g %g
\n
"
,
ti_current
*
All
.
Timebase_interval
,
Entropy
,
lmbd
);
while
(
ti_current
<
tend
)
{
/* compute lambda */
lmbd
=
lambda
(
Density
*
a3inv
,
Entropy
,
Metal
);
/* compute da/dt */
dadt
=
fabs
(
-
GAMMA_MINUS1
*
pow
(
Density
*
a3inv
,
-
GAMMA
)
*
lmbd
/
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
;
/* normal integration of Entropy */
Entropy
+=
DtEntropy
*
dt
;
/* viscosity */
Entropy
+=
-
GAMMA_MINUS1
*
pow
(
Density
*
a3inv
,
-
GAMMA
)
*
lmbd
/
hubble_a
*
dt
;
/* cooling */
/* avoid Entropy to be less than minentropy */
if
(
All
.
MinEgySpec
)
if
(
Entropy
<
minentropy
)
{
Entropy
=
minentropy
;
break
;
}
printf
(
"%g %g %g
\n
"
,
ti_current
*
All
.
Timebase_interval
,
Entropy
,
lmbd
);
ti_current
+=
ti_step
;
istep
=
istep
+
1
;
}
printf
(
"%g %g %g
\n
"
,
ti_current
*
All
.
Timebase_interval
,
Entropy
,
lmbd
);
/* entropy only due to cooling */
DEntropy
=
Entropy
-
oldEntropy
-
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);
printf
(
"Entropy Final %g
\n
"
,
Entropy
);
/* update particle */
*
newEntropy
=
Entropy
;
*
newDtEntropy
=
DtEntropy
+
oldDtEntropy
;
/* need to return */
/* Entropy, DtEntropy */
return
0
;
}
struct
cooling_solver_params
{
double
Entropy
;
double
Density
;
double
Metal
;
double
DtEntropy
;
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
;
double
Metal
=
p
->
Metal
;
double
DtEntropyVisc
=
p
->
DtEntropy
;
double
dt
=
p
->
dt
;
double
hubble_a
=
p
->
hubble_a
;
double
DtEntropyRadSph
=
0
;
DtEntropyRadSph
=
-
GAMMA_MINUS1
*
pow
(
Density
,
-
GAMMA
)
*
lambda
(
Density
,
EntropyVar
,
Metal
)
/
hubble_a
;
return
Entropy
+
(
DtEntropyVisc
+
DtEntropyRadSph
)
*
dt
-
EntropyVar
;
};
float
cooling2
(
int
tstart
,
int
tend
,
FLOAT
Entropy
,
FLOAT
Density
,
FLOAT
Metal
,
FLOAT
DtEntropy
,
FLOAT
*
newEntropy
,
FLOAT
*
newDtEntropy
)
{
double
dt
,
dadt
,
tcool
,
dt_entr
;
double
minentropy
;
double
MinSizeTimestep
,
ErrTolIntAccuracy
;
int
ti_current
,
istep
;
int
ti_step
;
FLOAT
oldEntropy
,
oldDtEntropy
,
DEntropy
;
double
lmbd
=
0
;
if
(
All
.
MinEgySpec
)
minentropy
=
All
.
MinEgySpec
*
GAMMA_MINUS1
/
pow
(
Density
*
a3inv
,
GAMMA_MINUS1
);
oldEntropy
=
Entropy
;
oldDtEntropy
=
DtEntropy
;
/***********************************/
double
EntropyNew
;
double
Entropy_lo
=
0
,
Entropy_hi
=
0
;
double
lo
,
hi
;
int
status
;
int
iter
=
0
;
int
max_iter
=
100
;
dt
=
(
tend
-
tstart
)
*
All
.
Timebase_interval
;
const
gsl_root_fsolver_type
*
T
;
gsl_root_fsolver
*
s
;
gsl_function
F
;
struct
cooling_solver_params
params
=
{(
double
)
Entropy
,(
double
)
Density
,(
double
)
Metal
,(
double
)
DtEntropy
,(
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
);
printf
(
"----->%g
\n
"
,
Entropy
);
do
{
iter
++
;
status
=
gsl_root_fsolver_iterate
(
s
);
EntropyNew
=
gsl_root_fsolver_root
(
s
);
printf
(
"----->%g
\n
"
,
EntropyNew
);
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);
}
printf
(
"--------->%g
\n
"
,
EntropyNew
);
/***********************************/
/* update particle */
*
newEntropy
=
Entropy
;
*
newDtEntropy
=
DtEntropy
+
oldDtEntropy
;
/* need to return */
/* Entropy, DtEntropy */
return
0
;
}
/****************************************************************************************/
/*
/*
/*
/* PYTHON INTERFACE
/*
/*
/*
/****************************************************************************************/
static
PyObject
*
cooling_get_lambda_from_Density_EnergyInt_FeH
(
self
,
args
)
PyObject
*
self
;
PyObject
*
args
;
{
double
density
,
energyint
,
feh
;
PyObject
*
densityO
,
*
energyintO
,
*
fehO
;
PyArrayObject
*
densityA
,
*
energyintA
,
*
fehA
;
double
l
;
PyArrayObject
*
ls
;
int
n
,
i
;
if
(
!
PyArg_ParseTuple
(
args
,
"OOO"
,
&
densityO
,
&
energyintO
,
&
fehO
))
return
PyString_FromString
(
"error"
);
/* a scalar */
if
(
PyArray_IsAnyScalar
(
densityO
)
&&
PyArray_IsAnyScalar
(
energyintO
)
&&
PyArray_IsAnyScalar
(
fehO
))
{
l
=
lambda_from_Density_EnergyInt_FeH
(
PyFloat_AsDouble
(
densityO
),
PyFloat_AsDouble
(
energyintO
),
PyFloat_AsDouble
(
fehO
));
return
Py_BuildValue
(
"d"
,
l
);
}
/* an array scalar */
if
(
PyArray_Check
(
densityO
)
&&
PyArray_Check
(
energyintO
)
&&
PyArray_Check
(
fehO
))
{
/* convert into array */
densityA
=
(
PyArrayObject
*
)
densityO
;
energyintA
=
(
PyArrayObject
*
)
energyintO
;
fehA
=
(
PyArrayObject
*
)
fehO
;
/* convert arrays to double */
densityA
=
TO_DOUBLE
(
densityA
);
energyintA
=
TO_DOUBLE
(
energyintA
);
fehA
=
TO_DOUBLE
(
fehA
);
/* check dimension and size */
if
(
(
densityA
->
nd
!=
1
)
||
(
energyintA
->
nd
!=
1
)
||
(
fehA
->
nd
!=
1
)
)
{
PyErr_SetString
(
PyExc_ValueError
,
"array objects must be of dimention 1."
);
return
NULL
;
}
n
=
densityA
->
dimensions
[
0
];
if
(
(
energyintA
->
dimensions
[
0
]
!=
n
)
||
(
fehA
->
dimensions
[
0
]
!=
n
)
)
{
PyErr_SetString
(
PyExc_ValueError
,
"arrays must have the same size."
);
return
NULL
;
}
/* create output array */
/*ls = (PyArrayObject *) PyArray_SimpleNewFromDescr(densityA->nd,densityA->dimensions,ld,PyArray_DescrFromType(NPY_DOUBLE));*/
ls
=
(
PyArrayObject
*
)
PyArray_SimpleNew
(
densityA
->
nd
,
densityA
->
dimensions
,
NPY_DOUBLE
);
for
(
i
=
0
;
i
<
n
;
i
++
)
{
*
(
double
*
)(
ls
->
data
+
(
i
)
*
(
ls
->
strides
[
0
]))
=
lambda_from_Density_EnergyInt_FeH
(
*
(
double
*
)(
densityA
->
data
+
(
i
)
*
(
densityA
->
strides
[
0
])),
*
(
double
*
)(
energyintA
->
data
+
(
i
)
*
(
energyintA
->
strides
[
0
])),
*
(
double
*
)(
fehA
->
data
+
(
i
)
*
(
fehA
->
strides
[
0
]))
);
}
return
Py_BuildValue
(
"O"
,
ls
);
}
/* something else */
PyErr_SetString
(
PyExc_ValueError
,
"parameters must be either scalars or array objects."
);
return
NULL
;
}
static
PyObject
*
cooling_get_lambda_from_Density_Entropy_FeH
(
self
,
args
)
PyObject
*
self
;
PyObject
*
args
;
{
double
density
,
entropy
,
fhe
;
PyObject
*
densityO
,
*
entropyO
,
*
fehO
;
PyArrayObject
*
densityA
,
*
entropyA
,
*
fehA
;
double
l
;
PyArrayObject
*
ls
;
int
n
,
i
;
if
(
!
PyArg_ParseTuple
(
args
,
"OOO"
,
&
densityO
,
&
entropyO
,
&
fehO
))
return
PyString_FromString
(
"error"
);
/* a scalar */
if
(
PyArray_IsAnyScalar
(
densityO
)
&&
PyArray_IsAnyScalar
(
entropyO
)
&&
PyArray_IsAnyScalar
(
fehO
))
{
l
=
lambda_from_Density_Entropy_FeH
(
PyFloat_AsDouble
(
densityO
),
PyFloat_AsDouble
(
entropyO
),
PyFloat_AsDouble
(
fehO
));
return
Py_BuildValue
(
"d"
,
l
);
}
/* an array scalar */
if
(
PyArray_Check
(
densityO
)
&&
PyArray_Check
(
entropyO
)
&&
PyArray_Check
(
fehO
))
{
/* convert into array */
densityA
=
(
PyArrayObject
*
)
densityO
;
entropyA
=
(
PyArrayObject
*
)
entropyO
;
fehA
=
(
PyArrayObject
*
)
fehO
;
/* convert arrays to double */
densityA
=
TO_DOUBLE
(
densityA
);
entropyA
=
TO_DOUBLE
(
entropyA
);
fehA
=
TO_DOUBLE
(
fehA
);
/* check dimension and size */
if
(
(
densityA
->
nd
!=
1
)
||
(
entropyA
->
nd
!=
1
)
||
(
fehA
->
nd
!=
1
)
)
{
PyErr_SetString
(
PyExc_ValueError
,
"array objects must be of dimention 1."
);
return
NULL
;
}
n
=
densityA
->
dimensions
[
0
];
if
(
(
entropyA
->
dimensions
[
0
]
!=
n
)
||
(
fehA
->
dimensions
[
0
]
!=
n
)
)
{
PyErr_SetString
(
PyExc_ValueError
,
"arrays must have the same size."
);
return
NULL
;
}
/* create output array */
/*ls = (PyArrayObject *) PyArray_SimpleNewFromDescr(densityA->nd,densityA->dimensions,ld,PyArray_DescrFromType(NPY_DOUBLE));*/
ls
=
(
PyArrayObject
*
)
PyArray_SimpleNew
(
densityA
->
nd
,
densityA
->
dimensions
,
NPY_DOUBLE
);
for
(
i
=
0
;
i
<
n
;
i
++
)
{
*
(
double
*
)(
ls
->
data
+
(
i
)
*
(
ls
->
strides
[
0
]))
=
lambda_from_Density_Entropy_FeH
(
*
(
double
*
)(
densityA
->
data
+
(
i
)
*
(
densityA
->
strides
[
0
])),
*
(
double
*
)(
entropyA
->
data
+
(
i
)
*
(
entropyA
->
strides
[
0
])),
*
(
double
*
)(
fehA
->
data
+
(
i
)
*
(
fehA
->
strides
[
0
]))
);
}
return
Py_BuildValue
(
"O"
,
ls
);
}
/* something else */
PyErr_SetString
(
PyExc_ValueError
,
"parameters must be either scalars or array objects."
);
return
NULL
;
}
static
PyObject
*
cooling_get_lambda_from_Density_Temperature_FeH
(
self
,
args
)
PyObject
*
self
;
PyObject
*
args
;
{
double
density
,
temperature
,
fhe
;
PyObject
*
densityO
,
*
temperatureO
,
*
fehO
;
PyArrayObject
*
densityA
,
*
temperatureA
,
*
fehA
;
double
l
;
PyArrayObject
*
ls
;
int
n
,
i
;
if
(
!
PyArg_ParseTuple
(
args
,
"OOO"
,
&
densityO
,
&
temperatureO
,
&
fehO
))
return
PyString_FromString
(
"error"
);
/* a scalar */
if
(
PyArray_IsAnyScalar
(
densityO
)
&&
PyArray_IsAnyScalar
(
temperatureO
)
&&
PyArray_IsAnyScalar
(
fehO
))
{
l
=
lambda_from_Density_Temperature_FeH
(
PyFloat_AsDouble
(
densityO
),
PyFloat_AsDouble
(
temperatureO
),
PyFloat_AsDouble
(
fehO
));
return
Py_BuildValue
(
"d"
,
l
);
}
/* an array scalar */
if
(
PyArray_Check
(
densityO
)
&&
PyArray_Check
(
temperatureO
)
&&
PyArray_Check
(
fehO
))
{
/* convert into array */
densityA
=
(
PyArrayObject
*
)
densityO
;
temperatureA
=
(
PyArrayObject
*
)
temperatureO
;
fehA
=
(
PyArrayObject
*
)
fehO
;
/* convert arrays to double */
densityA
=
TO_DOUBLE
(
densityA
);
temperatureA
=
TO_DOUBLE
(
temperatureA
);
fehA
=
TO_DOUBLE
(
fehA
);
/* check dimension and size */
if
(
(
densityA
->
nd
!=
1
)
||
(
temperatureA
->
nd
!=
1
)
||
(
fehA
->
nd
!=
1
)
)
{
PyErr_SetString
(
PyExc_ValueError
,
"array objects must be of dimention 1."
);
return
NULL
;
}
n
=
densityA
->
dimensions
[
0
];
if
(
(
temperatureA
->
dimensions
[
0
]
!=
n
)
||
(
fehA
->
dimensions
[
0
]
!=
n
)
)
{
PyErr_SetString
(
PyExc_ValueError
,
"arrays must have the same size."
);
return
NULL
;
}
/* create output array */
/*ls = (PyArrayObject *) PyArray_SimpleNewFromDescr(densityA->nd,densityA->dimensions,ld,PyArray_DescrFromType(NPY_DOUBLE));*/
ls
=
(
PyArrayObject
*
)
PyArray_SimpleNew
(
densityA
->
nd
,
densityA
->
dimensions
,
NPY_DOUBLE
);
for
(
i
=
0
;
i
<
n
;
i
++
)
{
*
(
double
*
)(
ls
->
data
+
(
i
)
*
(
ls
->
strides
[
0
]))
=
lambda_from_Density_Temperature_FeH
(
*
(
double
*
)(
densityA
->
data
+
(
i
)
*
(
densityA
->
strides
[
0
])),
*
(
double
*
)(
temperatureA
->
data
+
(
i
)
*
(
temperatureA
->
strides
[
0
])),
*
(
double
*
)(
fehA
->
data
+
(
i
)
*
(
fehA
->
strides
[
0
]))
);
}
return
Py_BuildValue
(
"O"
,
ls
);
}
/* something else */
PyErr_SetString
(
PyExc_ValueError
,
"parameters must be either scalars or array objects."
);
return
NULL
;
}
static
PyObject
*
cooling_get_lambda_normalized_from_Temperature_FeH
(
self
,
args
)
PyObject
*
self
;
PyObject
*
args
;
{
double
density
,
temperature
,
fhe
;
PyObject
*
temperatureO
,
*
fehO
;
PyArrayObject
*
temperatureA
,
*
fehA
;
double
l
;
PyArrayObject
*
ls
;
int
n
,
i
;
if
(
!
PyArg_ParseTuple
(
args
,
"OO"
,
&
temperatureO
,
&
fehO
))
return
PyString_FromString
(
"error"
);
/* a scalar */
if
(
PyArray_IsAnyScalar
(
temperatureO
)
&&
PyArray_IsAnyScalar
(
fehO
))
{
l
=
lambda_normalized_from_Temperature_FeH
(
PyFloat_AsDouble
(
temperatureO
),
PyFloat_AsDouble
(
fehO
));
return
Py_BuildValue
(
"d"
,
l
);
}
/* an array scalar */
if
(
PyArray_Check
(
temperatureO
)
&&
PyArray_Check
(
fehO
))
{
/* convert into array */
temperatureA
=
(
PyArrayObject
*
)
temperatureO
;
fehA
=
(
PyArrayObject
*
)
fehO
;
/* convert arrays to double */
temperatureA
=
TO_DOUBLE
(
temperatureA
);
fehA
=
TO_DOUBLE
(
fehA
);
/* check dimension and size */
if
((
temperatureA
->
nd
!=
1
)
||
(
fehA
->
nd
!=
1
)
)
{
PyErr_SetString
(
PyExc_ValueError
,
"array objects must be of dimention 1."
);
return
NULL
;
}
n
=
temperatureA
->
dimensions
[
0
];
if
(
fehA
->
dimensions
[
0
]
!=
n
)
{
PyErr_SetString
(
PyExc_ValueError
,
"arrays must have the same size."
);
return
NULL
;
}
/* create output array */
/*ls = (PyArrayObject *) PyArray_SimpleNewFromDescr(densityA->nd,densityA->dimensions,ld,PyArray_DescrFromType(NPY_DOUBLE));*/
ls
=
(
PyArrayObject
*
)
PyArray_SimpleNew
(
temperatureA
->
nd
,
temperatureA
->
dimensions
,
NPY_DOUBLE
);
for
(
i
=
0
;
i
<
n
;
i
++
)
{
*
(
double
*
)(
ls
->
data
+
(
i
)
*
(
ls
->
strides
[
0
]))
=
lambda_normalized_from_Temperature_FeH
(
*
(
double
*
)(
temperatureA
->
data
+
(
i
)
*
(
temperatureA
->
strides
[
0
])),
*
(
double
*
)(
fehA
->
data
+
(
i
)
*
(
fehA
->
strides
[
0
]))
);
}
return
Py_BuildValue
(
"O"
,
ls
);
}
/* something else */
PyErr_SetString
(
PyExc_ValueError
,
"parameters must be either scalars or array objects."
);
return
NULL
;
}
static
PyObject
*
cooling_get_cooling_time_from_Density_Temperature_FeH
(
self
,
args
)
PyObject
*
self
;
PyObject
*
args
;
{
double
density
,
temperature
,
fhe
;
PyObject
*
densityO
,
*
temperatureO
,
*
fehO
;
PyArrayObject
*
densityA
,
*
temperatureA
,
*
fehA
;
double
tc
;
PyArrayObject
*
tcs
;
int
n
,
i
;
if
(
!
PyArg_ParseTuple
(
args
,
"OOO"
,
&
densityO
,
&
temperatureO
,
&
fehO
))
return
PyString_FromString
(
"error"
);
/* a scalar */
if
(
PyArray_IsAnyScalar
(
densityO
)
&&
PyArray_IsAnyScalar
(
temperatureO
)
&&
PyArray_IsAnyScalar
(
fehO
))
{
tc
=
cooling_time_from_Density_Temperature_FeH
(
PyFloat_AsDouble
(
densityO
),
PyFloat_AsDouble
(
temperatureO
),
PyFloat_AsDouble
(
fehO
));
return
Py_BuildValue
(
"d"
,
tc
);
}
/* an array scalar */
if
(
PyArray_Check
(
densityO
)
&&
PyArray_Check
(
temperatureO
)
&&
PyArray_Check
(
fehO
))
{
/* convert into array */
densityA
=
(
PyArrayObject
*
)
densityO
;
temperatureA
=
(
PyArrayObject
*
)
temperatureO
;
fehA
=
(
PyArrayObject
*
)
fehO
;
/* convert arrays to double */
densityA
=
TO_DOUBLE
(
densityA
);
temperatureA
=
TO_DOUBLE
(
temperatureA
);
fehA
=
TO_DOUBLE
(
fehA
);
/* check dimension and size */
if
(
(
densityA
->
nd
!=
1
)
||
(
temperatureA
->
nd
!=
1
)
||
(
fehA
->
nd
!=
1
)
)
{
PyErr_SetString
(
PyExc_ValueError
,
"array objects must be of dimention 1."
);
return
NULL
;
}
n
=
densityA
->
dimensions
[
0
];
if
(
(
temperatureA
->
dimensions
[
0
]
!=
n
)
||
(
fehA
->
dimensions
[
0
]
!=
n
)
)
{
PyErr_SetString
(
PyExc_ValueError
,
"arrays must have the same size."
);
return
NULL
;
}
/* create output array */
/*tcs = (PyArrayObject *) PyArray_SimpleNewFromDescr(densityA->nd,densityA->dimensions,ld,PyArray_DescrFromType(NPY_DOUBLE));*/
tcs
=
(
PyArrayObject
*
)
PyArray_SimpleNew
(
densityA
->
nd
,
densityA
->
dimensions
,
NPY_DOUBLE
);
for
(
i
=
0
;
i
<
n
;
i
++
)
{
*
(
double
*
)(
tcs
->
data
+
(
i
)
*
(
tcs
->
strides
[
0
]))
=
cooling_time_from_Density_Temperature_FeH
(
*
(
double
*
)(
densityA
->
data
+
(
i
)
*
(
densityA
->
strides
[
0
])),
*
(
double
*
)(
temperatureA
->
data
+
(
i
)
*
(
temperatureA
->
strides
[
0
])),
*
(
double
*
)(
fehA
->
data
+
(
i
)
*
(
fehA
->
strides
[
0
]))
);
}
return
Py_BuildValue
(
"O"
,
tcs
);
}
/* something else */
PyErr_SetString
(
PyExc_ValueError
,
"parameters must be either scalars or array objects."
);
return
NULL
;
}
static
PyObject
*
cooling_get_cooling_time_from_Density_EnergyInt_FeH
(
self
,
args
)
PyObject
*
self
;
PyObject
*
args
;
{
double
density
,
energyint
,
fhe
;
PyObject
*
densityO
,
*
energyintO
,
*
fehO
;
PyArrayObject
*
densityA
,
*
energyintA
,
*
fehA
;
double
tc
;
PyArrayObject
*
tcs
;
int
n
,
i
;
if
(
!
PyArg_ParseTuple
(
args
,
"OOO"
,
&
densityO
,
&
energyintO
,
&
fehO
))
return
PyString_FromString
(
"error"
);
/* a scalar */
if
(
PyArray_IsAnyScalar
(
densityO
)
&&
PyArray_IsAnyScalar
(
energyintO
)
&&
PyArray_IsAnyScalar
(
fehO
))
{
tc
=
cooling_time_from_Density_EnergyInt_FeH
(
PyFloat_AsDouble
(
densityO
),
PyFloat_AsDouble
(
energyintO
),
PyFloat_AsDouble
(
fehO
));
return
Py_BuildValue
(
"d"
,
tc
);
}
/* an array scalar */
if
(
PyArray_Check
(
densityO
)
&&
PyArray_Check
(
energyintO
)
&&
PyArray_Check
(
fehO
))
{
/* convert into array */
densityA
=
(
PyArrayObject
*
)
densityO
;
energyintA
=
(
PyArrayObject
*
)
energyintO
;
fehA
=
(
PyArrayObject
*
)
fehO
;
/* convert arrays to double */
densityA
=
TO_DOUBLE
(
densityA
);
energyintA
=
TO_DOUBLE
(
energyintA
);
fehA
=
TO_DOUBLE
(
fehA
);
/* check dimension and size */
if
(
(
densityA
->
nd
!=
1
)
||
(
energyintA
->
nd
!=
1
)
||
(
fehA
->
nd
!=
1
)
)
{
PyErr_SetString
(
PyExc_ValueError
,
"array objects must be of dimention 1."
);
return
NULL
;
}
n
=
densityA
->
dimensions
[
0
];
if
(
(
energyintA
->
dimensions
[
0
]
!=
n
)
||
(
fehA
->
dimensions
[
0
]
!=
n
)
)
{
PyErr_SetString
(
PyExc_ValueError
,
"arrays must have the same size."
);
return
NULL
;
}
/* create output array */
/*tcs = (PyArrayObject *) PyArray_SimpleNewFromDescr(densityA->nd,densityA->dimensions,ld,PyArray_DescrFromType(NPY_DOUBLE));*/
tcs
=
(
PyArrayObject
*
)
PyArray_SimpleNew
(
densityA
->
nd
,
densityA
->
dimensions
,
NPY_DOUBLE
);
for
(
i
=
0
;
i
<
n
;
i
++
)
{
*
(
double
*
)(
tcs
->
data
+
(
i
)
*
(
tcs
->
strides
[
0
]))
=
cooling_time_from_Density_EnergyInt_FeH
(
*
(
double
*
)(
densityA
->
data
+
(
i
)
*
(
densityA
->
strides
[
0
])),
*
(
double
*
)(
energyintA
->
data
+
(
i
)
*
(
energyintA
->
strides
[
0
])),
*
(
double
*
)(
fehA
->
data
+
(
i
)
*
(
fehA
->
strides
[
0
]))
);
}
return
Py_BuildValue
(
"O"
,
tcs
);
}
/* something else */
PyErr_SetString
(
PyExc_ValueError
,
"parameters must be either scalars or array objects."
);
return
NULL
;
}
/*********************************/
/* */
/*********************************/
static
PyObject
*
cooling_init_cooling
(
self
,
args
)
PyObject
*
self
;
PyObject
*
args
;
{
PyObject
*
params_dict
=
PyDict_New
();
set_default_parameters
();
if
(
!
PyArg_ParseTuple
(
args
,
"|O"
,
&
params_dict
))
{
PyErr_SetString
(
PyExc_AttributeError
,
"error."
);
return
NULL
;
}
if
(
!
PyDict_Check
(
params_dict
))
{
PyErr_SetString
(
PyExc_AttributeError
,
"argument is not a dictionary."
);
return
NULL
;
}
if
(
PyDict_Size
(
params_dict
))
set_parameters
(
params_dict
);
set_units
();
init_cooling_with_metals
();
return
Py_BuildValue
(
"d"
,
0
);
}
/*********************************/
/* */
/*********************************/
static
PyObject
*
cooling_PrintParameters
(
self
,
args
)
PyObject
*
self
;
PyObject
*
args
;
{
printf
(
"UnitLength_in_cm = %g
\n
"
,
All
.
UnitLength_in_cm
);
printf
(
"UnitMass_in_g = %g
\n
"
,
All
.
UnitMass_in_g
);
printf
(
"UnitVelocity_in_cm_per_s = %g
\n
"
,
All
.
UnitVelocity_in_cm_per_s
);
printf
(
"CoolingFile = %s
\n
"
,
All
.
CoolingFile
);
return
Py_BuildValue
(
"i"
,
0
);
}
/*********************************/
/* */
/*********************************/
static
PyObject
*
cooling_integrate1
(
self
,
args
)
PyObject
*
self
;
PyObject
*
args
;
{
double
tstart
,
tend
,
Density
,
Entropy
,
Metal
;
double
DtEntropy
=
0
;
int
itstart
,
itend
;
FLOAT
newEntropy
;
FLOAT
newDtEntropy
;
if
(
!
PyArg_ParseTuple
(
args
,
"|ddddd"
,
&
tstart
,
&
tend
,
&
Density
,
&
Entropy
,
&
Metal
))
{
PyErr_SetString
(
PyExc_AttributeError
,
"error."
);
return
NULL
;
}
itstart
=
tstart
/
All
.
Timebase_interval
;
itend
=
tend
/
All
.
Timebase_interval
;
cooling1
(
itstart
,
itend
,
(
FLOAT
)
Entropy
,
(
FLOAT
)
Density
,
(
FLOAT
)
Metal
,
DtEntropy
,
&
newEntropy
,
&
newDtEntropy
)
;
return
Py_BuildValue
(
"i"
,
0
);
}
/*********************************/
/* */
/*********************************/
static
PyObject
*
cooling_integrate2
(
self
,
args
)
PyObject
*
self
;
PyObject
*
args
;
{
double
tstart
,
tend
,
Density
,
Entropy
,
Metal
;
double
DtEntropy
=
0
;
int
itstart
,
itend
;
FLOAT
newEntropy
;
FLOAT
newDtEntropy
;
if
(
!
PyArg_ParseTuple
(
args
,
"|ddddd"
,
&
tstart
,
&
tend
,
&
Density
,
&
Entropy
,
&
Metal
))
{
PyErr_SetString
(
PyExc_AttributeError
,
"error."
);
return
NULL
;
}
itstart
=
tstart
/
All
.
Timebase_interval
;
itend
=
tend
/
All
.
Timebase_interval
;
cooling2
(
itstart
,
itend
,
(
FLOAT
)
Entropy
,
(
FLOAT
)
Density
,
(
FLOAT
)
Metal
,
DtEntropy
,
&
newEntropy
,
&
newDtEntropy
)
;
return
Py_BuildValue
(
"i"
,
0
);
}
/* definition of the method table */
static
PyMethodDef
coolingMethods
[]
=
{
{
"init_cooling"
,
cooling_init_cooling
,
METH_VARARGS
,
"Init cooling."
},
{
"get_lambda_from_Density_EnergyInt_FeH"
,
cooling_get_lambda_from_Density_EnergyInt_FeH
,
METH_VARARGS
,
"Get the lambda function in user units."
},
{
"get_lambda_from_Density_Entropy_FeH"
,
cooling_get_lambda_from_Density_Entropy_FeH
,
METH_VARARGS
,
"Get the lambda function in user units."
},
{
"get_lambda_from_Density_Temperature_FeH"
,
cooling_get_lambda_from_Density_Temperature_FeH
,
METH_VARARGS
,
"Get the lambda function in user units."
},
{
"get_lambda_normalized_from_Temperature_FeH"
,
cooling_get_lambda_normalized_from_Temperature_FeH
,
METH_VARARGS
,
"Get the normalized lambda function in mks."
},
{
"get_cooling_time_from_Density_Temperature_FeH"
,
cooling_get_cooling_time_from_Density_Temperature_FeH
,
METH_VARARGS
,
"Get the cooling time in user units."
},
{
"get_cooling_time_from_Density_EnergyInt_FeH"
,
cooling_get_cooling_time_from_Density_EnergyInt_FeH
,
METH_VARARGS
,
"Get the cooling time in user units."
},
{
"PrintParameters"
,
cooling_PrintParameters
,
METH_VARARGS
,
"Print parameters."
},
{
"integrate1"
,
cooling_integrate1
,
METH_VARARGS
,
"Integrate cooling during a timestep using first scheme of integration."
},
{
"integrate2"
,
cooling_integrate2
,
METH_VARARGS
,
"Integrate cooling during a timestep using second scheme of integration."
},
{
NULL
,
NULL
,
0
,
NULL
}
/* Sentinel */
};
void
initcooling_with_metals
(
void
)
{
(
void
)
Py_InitModule
(
"cooling_with_metals"
,
coolingMethods
);
import_array
();
}
Event Timeline
Log In to Comment