Homec4science
Diffusion pNbody (master) Browse

pNbody/src/cooling_with_metals/sto.cooling_with_metals.c
4fbc3430dc96master

Last Change
Raw File

sto.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>
#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;
} 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_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 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;
}
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;
}
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;
}
/****************************************************************************************/
/*
/*
/*
/* 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_init_cooling(self, args)
PyObject *self;
PyObject *args;
{
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.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;
init_cooling_with_metals();
//printf("----> %g\n",lambda_from_Density_EnergyInt_FeH(3.0364363e-06,0.0015258887,0.));
double nH,nH2,T,l;
double Density,EnergyInt,FeH;
Density = 3.0364363e-06;
EnergyInt = 0.0015258887;
FeH = 0.;
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;
printf("----> %g\n",l);
printf("----> %g\n",lambda_from_Density_EnergyInt_FeH(3.0364363e-06,0.0015258887,0.));
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);
}
/* 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."},
{"PrintParameters", cooling_PrintParameters, METH_VARARGS,
"Print parameters."},
{NULL, NULL, 0, NULL} /* Sentinel */
};
void initcooling_with_metals(void)
{
(void) Py_InitModule("cooling_with_metals", coolingMethods);
import_array();
}
c4science · Help