Homec4science
Diffusion Gear (master) Browse

Gear/src/cooling_wiersma.c
c1eaec5e7047master

Last Change
Raw File

cooling_wiersma.c
View Options

#ifdef COOLING
#ifdef COOLING_WIERSMA
#ifdef PY_INTERFACE
#include <Python.h>
#else
#include <stdio.h>
#include <string.h>
#endif
#include "lt3d.h"
#include "hdf5io.h"
#include <stdlib.h>
#include <math.h>
#include <mpi.h>
#include <dirent.h>
#include <hdf5.h>
#include <math.h>
#include <mpi.h>
#define HEADER_GRP "/Header"
#define TOTAL_METAL_GRP "/Total_Metals"
#define METAL_FREE_GRP "/Metal_free"
#define SOLAR_GRP "/Solar"
#define TABLE_DIR "/home/epfl/revaz/code/gear/PyCool/tables_wiersma/coolingtables/"
static int locThisTask=0;
/* some global variables */
static int RedshiftInterpolation=0;
struct RedshiftFileStruct
{
double Redshift;
char * filename;
}
*RedshiftFiles;
static int nFiles;
double currentRedshift;
struct CoolingTableStruct
{
char * filename;
int index;
double Redshift;
int Number_of_temperature_bins;
int Number_of_density_bins;
int Number_of_helium_fractions;
int Number_of_species;
//double *Temperature_bins;
//double *Hydrogen_density_bins;
//double *Helium_mass_fraction_bins;
struct lt3d_RegularTable3D lambda_metal_free;
struct lt3d_RegularTable2D lambda_metal;
struct lt3d_RegularTable3D ne_n_H;
struct lt3d_RegularTable2D ne_n_H_solar;
}
cTable0,cTable1;
/*! This is a comparison kernel used in a sort routine.
*/
int redshift_compare_key(const void *a, const void *b)
{
if(((struct RedshiftFileStruct *) a)->Redshift < (((struct RedshiftFileStruct *) b)->Redshift))
return -1;
if(((struct RedshiftFileStruct *) a)->Redshift > (((struct RedshiftFileStruct *) b)->Redshift))
return +1;
return 0;
}
double i2Redshift(int i)
{
double xmin=0;
double dx=1/48.;
return pow(10,xmin+i*dx)-1;
}
int Redshift2i(double Redshift)
{
double xmin=0;
double dx=1/48.;
int idx;
idx = (int)round( (log10(Redshift+1)-xmin)/(dx) );
if (idx<0)
return 0;
if (idx>48)
return 48;
return idx;
}
int Redshift2ij(double Redshift, int *idx0, int *idx1)
{
double xmin=0;
double dx=1/48.;
*idx1 = (int)round( (log10(Redshift+1)-xmin)/(dx) );
if (*idx1<0)
*idx1 =0;
if (*idx1>48)
*idx1= 48;
*idx0 = *idx1-1;
if (*idx0<0)
*idx0 = 0;
if (*idx1>48)
*idx0 = 48;
return 0;
}
int findRedshiftij(double Redshift, int *idx0, int *idx1)
{
int i;
/* more than the max */
if (Redshift>=RedshiftFiles[nFiles-1].Redshift)
{
*idx0 = nFiles-1;
*idx1 = nFiles-1;
return 0;
}
i=0;
for (i=0;i<nFiles-2;i++)
{
if (Redshift<=RedshiftFiles[i].Redshift)
break;
}
if (i<=0)
{
*idx0 = 0;
*idx1 = 0;
}
else
{
*idx0 = i-1;
*idx1 = i;
}
return 0;
}
double i2Helium_mass_fraction(int i)
{
double xmin=0.238;
double dx=0.01;
return xmin + i*dx;
}
int Helium_mass_fraction2i(double Helium_mass_fraction)
{
double xmin=0.238;
double dx=0.01;
int idx;
idx = (int)round( (Helium_mass_fraction - xmin)/dx );
if (idx<0)
return 0;
if (idx>6)
return 6;
return idx;
}
double i2Hydrogen_density(int i)
{
double xmin=-8.0;
double dx=0.1;
return pow( 10.0, xmin + i*dx);
}
int Hydrogen_density2i(double Hydrogen_density)
{
double xmin=-8.0;
double dx=0.1;
int idx;
idx = (int)round( (log10(Hydrogen_density) - xmin)/dx );
if (idx<0)
return 0;
if (idx>80)
return 80;
return idx;
}
double i2Temperature(int i)
{
double xmin=2.0;
double xmax=8.98185062;
double n = 352;
double dx;
dx = (xmax-xmin)/(n-1);
return pow(10,(xmin + i*dx));
}
int Temperature2i(double Temperature)
{
double xmin=2.0;
double xmax=8.98185062;
double n = 352;
double dx;
int idx;
dx = (xmax-xmin)/(n-1);
idx = (int)round( (log10(Temperature) - xmin)/dx );
if (idx<0)
return 0;
if (idx>n-1)
return n-1;
return idx;
}
/*! determine if the string str ends by the string suffix
* function used here to determine the extension of file given their name
*
*/
int endsWith(const char *str, const char *suffix) {
if (!str || !suffix)
return 0;
size_t lenstr = strlen(str);
size_t lensuffix = strlen(suffix);
if (lensuffix > lenstr)
return 0;
return strncmp(str + lenstr - lensuffix, suffix, lensuffix) == 0;
}
int SetRedshiftInterpolationOn()
{
RedshiftInterpolation=1;
return 0;
}
int SetRedshiftInterpolationOff()
{
RedshiftInterpolation=0;
return 0;
}
int InitWiersmaCooling(char * TablesDirectory)
{
DIR *dir;
struct dirent* content;
int hdf5 = 0;
int i;
MPI_Comm_rank(MPI_COMM_WORLD, &locThisTask);
SetRedshiftInterpolationOn();
if (locThisTask==0)
{
/* some general inits */
cTable0.index=-1;
cTable1.index=-1;
/* open directory */
dir = opendir(TablesDirectory);
if(dir == NULL){
printf("an error occured while opening hdf5 directory %s\n", TablesDirectory);
return 1;
}
/* scan directory */
nFiles=0;
while((content = readdir(dir)) != NULL)
{
hdf5 = endsWith(content->d_name, ".hdf5");
if(hdf5 == 1)
nFiles++;
}
closedir(dir);
printf("%d hdf5 files\n",nFiles);
RedshiftFiles = malloc(nFiles*sizeof(struct RedshiftFileStruct));
i=0;
dir = opendir(TablesDirectory);
while((content = readdir(dir)) != NULL)
{
hdf5 = endsWith(content->d_name, ".hdf5");
if(hdf5 == 1)
{
/* read the file and record redshift */
char file_path[(int)strlen(TablesDirectory)+(int)strlen(content->d_name)];
strcpy(file_path, TablesDirectory);
strcat(file_path, content->d_name);
//printf("reading: %s\n",file_path);
hid_t table;
hid_t group;
hid_t dset;
herr_t status;
float fRedshift;
table = H5Fopen(file_path, H5F_ACC_RDONLY, H5P_DEFAULT);
group = H5Gopen( table , HEADER_GRP, H5P_DEFAULT);
dset = H5Dopen( group , "Redshift", H5P_DEFAULT);
status = H5Dread(dset,H5T_IEEE_F32LE,H5S_ALL, H5S_ALL, H5P_DEFAULT,&fRedshift);
status = H5Dclose (dset);
status = H5Gclose(group);
status = H5Fclose(table);
//printf("%g\n",fRedshift);
/* copy to the file list */
RedshiftFiles[i].Redshift = (double)fRedshift;
RedshiftFiles[i].filename = (char *) malloc((strlen(file_path)+1)*sizeof(char));
strcpy(RedshiftFiles[i].filename,file_path);
//printf("%8.3f %s\n",RedshiftFiles[i].Redshift,RedshiftFiles[i].filename);
i++;
}
}
/* sort according to the redshift */
qsort(RedshiftFiles, nFiles, sizeof(struct RedshiftFileStruct), redshift_compare_key);
//for(i=0;i<nFiles;i++)
// printf("%d %8.5f %8.5f %d %s\n",i,RedshiftFiles[i].Redshift, i2Redshift(i) , Redshift2i(RedshiftFiles[i].Redshift) ,RedshiftFiles[i].filename);
/* check */
//for (i=0;i<100;i++)
// printf("i=%03d z=%8.5f idx=%d z=%8.5f\n",i, i*0.1 , Redshift2i(i*0.1), RedshiftFiles[Redshift2i(i*0.1)].Redshift );
}
return 0;
}
int CoolingTableFree(struct CoolingTableStruct ctab)
{
/* tables */
//lt3d_FreeRegularTable3D(ctab.lambda_metal_free);
//lt3d_FreeRegularTable2D(ctab.lambda_metal);
//lt3d_FreeRegularTable3D(ctab.ne_n_H);
//lt3d_FreeRegularTable2D(ctab.ne_n_H_solar);
//free(ctab.filename);
/* bins */
//free(ctab.Temperature_bins);
//free(ctab.Hydrogen_density_bins);
//free(ctab.Helium_mass_fraction_bins);
return 0;
}
int Free()
{
CoolingTableFree(cTable0);
CoolingTableFree(cTable1);
return 0;
}
int readCoolingTable(char *table_name,struct CoolingTableStruct *ctab)
{
hid_t table;
hid_t group;
herr_t status;
printf("reading %s...\n",table_name);
table = H5Fopen(table_name, H5F_ACC_RDONLY, H5P_DEFAULT);
/* store the name */
ctab->filename = malloc(strlen(table_name)*sizeof(char)+1);
strcpy(ctab->filename, table_name);
/* read header */
group = H5Gopen( table , HEADER_GRP, H5P_DEFAULT);
ctab->Redshift = readDataSetAsDouble(group,"Redshift");
ctab->Number_of_temperature_bins = readDataSetAsInt(group,"Number_of_temperature_bins");
ctab->Number_of_density_bins = readDataSetAsInt(group,"Number_of_density_bins");
//ctab->Number_of_helium_fractions = readDataSetAsInt(group,"Number_of_helium_fractions");
ctab->Number_of_helium_fractions = 7; /* seems to be a bug in the table */
ctab->Number_of_species = readDataSetAsInt(group,"Number_of_species");
status = H5Gclose (group);
/*****************/
/* read the bins */
/*****************/
group = H5Gopen( table , METAL_FREE_GRP, H5P_DEFAULT);
/* not needed (only for tests) */
//ctab->Temperature_bins = malloc( ctab->Number_of_temperature_bins*sizeof(double));
//ctab->Hydrogen_density_bins = malloc( ctab->Number_of_density_bins*sizeof(double));
//ctab->Helium_mass_fraction_bins = malloc( ctab->Number_of_helium_fractions*sizeof(double));
//ctab->Temperature_bins = readDatasetAsArrayDouble(group,"Temperature_bins");
//ctab->Hydrogen_density_bins = readDatasetAsArrayDouble(group,"Hydrogen_density_bins");
//ctab->Helium_mass_fraction_bins = readDatasetAsArrayDouble(group,"Helium_mass_fraction_bins");
/* use lt3d objects */
/* Temperature */
struct lt3d_RegularBinning logTemperature_bins;
logTemperature_bins.name = "logTemperature";
logTemperature_bins.n = ctab->Number_of_temperature_bins;
logTemperature_bins.x0 = 2.0;
logTemperature_bins.dx = 0.019891312307692307;
/* Hydrogen */
struct lt3d_RegularBinning logHydrogen_density_bins;
logHydrogen_density_bins.name = "logHydrogen_density";
logHydrogen_density_bins.n = ctab->Number_of_density_bins;
logHydrogen_density_bins.x0 = -8.0;
logHydrogen_density_bins.dx = 0.1;
/* Helium mass fraction */
struct lt3d_RegularBinning Helium_mass_fraction_bins;
Helium_mass_fraction_bins.name = "Helium mass fraction";
Helium_mass_fraction_bins.n = ctab->Number_of_helium_fractions;
Helium_mass_fraction_bins.x0 = 0.238;
Helium_mass_fraction_bins.dx = 0.01;
status = H5Gclose (group);
/*******************/
/* read the tables */
/*******************/
/* cooling function, metal free */
group = H5Gopen( table , METAL_FREE_GRP, H5P_DEFAULT);
strcpy(ctab->lambda_metal_free.name,"Net_Cooling_metal_free");
//ctab->lambda_metal_free.data = readDatasetAsArray3dDouble(group,"Net_Cooling",ctab->lambda_metal_free.dims);
readDatasetAsArray3dDouble(group,"Net_Cooling",ctab->lambda_metal_free.dims,ctab->lambda_metal_free.data);
ctab->lambda_metal_free.x = Helium_mass_fraction_bins;
ctab->lambda_metal_free.y = logTemperature_bins;
ctab->lambda_metal_free.z = logHydrogen_density_bins;
status = H5Gclose (group);
/* cooling function, metal free */
group = H5Gopen( table , TOTAL_METAL_GRP, H5P_DEFAULT);
strcpy(ctab->lambda_metal.name ,"Net_Cooling");
//ctab->lambda_metal.data = readDatasetAsArray2dDouble(group,"Net_cooling",ctab->lambda_metal.dims);
readDatasetAsArray2dDouble(group,"Net_cooling",ctab->lambda_metal.dims,ctab->lambda_metal.data);
ctab->lambda_metal.x = logTemperature_bins;
ctab->lambda_metal.y = logHydrogen_density_bins;
status = H5Gclose (group);
/* electron density of hydrogen density */
group = H5Gopen( table , METAL_FREE_GRP, H5P_DEFAULT);
strcpy(ctab->ne_n_H.name,"Electron_density_over_n_h");
//ctab->ne_n_H.data = readDatasetAsArray3dDouble(group,"Electron_density_over_n_h",ctab->ne_n_H.dims);
readDatasetAsArray3dDouble(group,"Electron_density_over_n_h",ctab->ne_n_H.dims,ctab->ne_n_H.data);
ctab->ne_n_H.x = Helium_mass_fraction_bins;
ctab->ne_n_H.y = logTemperature_bins;
ctab->ne_n_H.z = logHydrogen_density_bins;
status = H5Gclose (group);
/* electron density of hydrogen density solar */
group = H5Gopen( table , SOLAR_GRP, H5P_DEFAULT);
strcpy(ctab->ne_n_H_solar.name,"Electron_density_over_n_h_solar");
//ctab->ne_n_H_solar.data = readDatasetAsArray2dDouble(group,"Electron_density_over_n_h",ctab->ne_n_H_solar.dims);
readDatasetAsArray2dDouble(group,"Electron_density_over_n_h",ctab->ne_n_H_solar.dims,ctab->ne_n_H_solar.data);
ctab->ne_n_H_solar.x = logTemperature_bins;
ctab->ne_n_H_solar.y = logHydrogen_density_bins;
status = H5Gclose (group);
status = H5Fclose(table);
return 0;
}
int printCoolingTable(struct CoolingTableStruct ctab)
{
//int i;
//char *space="\t\t";
printf("\n");
printf("table name: %s\n\n",ctab.filename);
printf("\t Redshift = %g\n\n",ctab.Redshift);
printf("\t Number_of_temperature_bins = %d\n\n",ctab.Number_of_temperature_bins);
//for(i=0;i<3;i++)
// printf("%s\t [%d] %g\n",space,i,ctab.Temperature_bins[i]);
//printf("%s\t ...\n",space);
//for(i=ctab.Number_of_temperature_bins-3;i<ctab.Number_of_temperature_bins;i++)
// printf("%s\t [%d] %g\n",space,i,ctab.Temperature_bins[i]);
//printf("\n");
printf("\t Number_of_density_bins = %d\n\n",ctab.Number_of_density_bins );
//for(i=0;i<3;i++)
// printf("%s\t [%d] %g\n",space,i,ctab.Hydrogen_density_bins[i]);
//printf("%s\t ...\n",space);
//for(i=ctab.Number_of_density_bins-3;i<ctab.Number_of_density_bins;i++)
// printf("%s\t [%d] %g\n",space,i,ctab.Hydrogen_density_bins[i]);
//printf("\n");
printf("\t Number_of_helium_fractions = %d\n\n",ctab.Number_of_helium_fractions);
//for(i=0;i<3;i++)
// printf("%s\t [%d] %g\n",space,i,ctab.Helium_mass_fraction_bins[i]);
//printf("%s\t ...\n",space);
//for(i=ctab.Number_of_helium_fractions-3;i<ctab.Number_of_helium_fractions;i++)
// printf("%s\t [%d] %g\n",space,i,ctab.Helium_mass_fraction_bins[i]);
//printf("\n");
printf("\t Number_of_species = %d\n",ctab.Number_of_species );
printf("\n");
return 0;
}
int setTablesFromRedshift(double RedShift)
{
int idx0,idx1;
currentRedshift=RedShift; /* make it global */
if (locThisTask==0)
{
findRedshiftij(RedShift,&idx0,&idx1);
//Redshift2ij(RedShift,&idx0,&idx1),
if(locThisTask==0)
printf("SeTables : z=%8.5f z0=%8.5f(%d) z1=%8.5f(%d) \n",RedShift, RedshiftFiles[idx0].Redshift,idx0, RedshiftFiles[idx1].Redshift,idx1 );
if (cTable0.index!=idx0)
{
readCoolingTable(RedshiftFiles[idx0].filename,&cTable0);
cTable0.index = idx0;
}
if (RedshiftInterpolation)
{
if (cTable1.index!=idx1)
{
readCoolingTable(RedshiftFiles[idx1].filename,&cTable1);
cTable1.index = idx1;
}
}
//printCoolingTable(cTable0);
//printCoolingTable(cTable1);
}
/* at this point, we needs to broadcast the tables */
MPI_Bcast(&cTable0, sizeof(struct CoolingTableStruct), MPI_BYTE, 0, MPI_COMM_WORLD);
MPI_Bcast(&cTable1, sizeof(struct CoolingTableStruct), MPI_BYTE, 0, MPI_COMM_WORLD);
return 0;
}
double compute_LambdaMetalFree(double T,double nH,double HeMfrac)
{
double lognH;
double logT;
double L,L0;
lognH = log10(nH);
logT = log10(T);
L0 = lt3d_getInterpolationFromRegularTable3D(&cTable0.lambda_metal_free,HeMfrac,logT,lognH);
if (RedshiftInterpolation)
{
if (cTable0.index==cTable1.index)
L = L0;
else
{
double L1;
double xd;
L1 = lt3d_getInterpolationFromRegularTable3D(&cTable1.lambda_metal_free,HeMfrac,logT,lognH);
xd = (currentRedshift-cTable0.Redshift) / (cTable1.Redshift-cTable0.Redshift);
L = L0*(1-xd) + L1*xd;
}
}
else
L = L0;
return L;
}
double compute_LambdaMetal(double T,double nH)
{
double lognH;
double logT;
double L,L0;
lognH = log10(nH);
logT = log10(T);
L0 = lt3d_getInterpolationFromRegularTable2D(&cTable0.lambda_metal,logT,lognH);
if (RedshiftInterpolation)
{
if (cTable0.index==cTable1.index)
L = L0;
else
{
double L1;
double xd;
L1 = lt3d_getInterpolationFromRegularTable2D(&cTable1.lambda_metal,logT,lognH);
xd = (currentRedshift-cTable0.Redshift) / (cTable1.Redshift-cTable0.Redshift);
L = L0*(1-xd) + L1*xd;
}
}
else
L = L0;
return L;
}
double compute_ElectronHydrogenRatio(double T,double nH,double HeMfrac)
{
double lognH;
double logT;
double L,L0;
lognH = log10(nH);
logT = log10(T);
L0 = lt3d_getInterpolationFromRegularTable3D(&cTable0.ne_n_H,HeMfrac,logT,lognH);
if (RedshiftInterpolation)
{
if (cTable0.index==cTable1.index)
L = L0;
else
{
double L1;
double xd;
L1 = lt3d_getInterpolationFromRegularTable3D(&cTable1.ne_n_H,HeMfrac,logT,lognH);
xd = (currentRedshift-cTable0.Redshift) / (cTable1.Redshift-cTable0.Redshift);
L = L0*(1-xd) + L1*xd;
}
}
else
L = L0;
return L;
}
double compute_ElectronHydrogenRatioSolar(double T,double nH)
{
double lognH;
double logT;
double L,L0;
lognH = log10(nH);
logT = log10(T);
L0 = lt3d_getInterpolationFromRegularTable2D(&cTable0.ne_n_H_solar,logT,lognH);
if (RedshiftInterpolation)
{
if (cTable0.index==cTable1.index)
L = L0;
else
{
double L1;
double xd;
L1 = lt3d_getInterpolationFromRegularTable2D(&cTable1.ne_n_H_solar,logT,lognH);
xd = (currentRedshift-cTable0.Redshift) / (cTable1.Redshift-cTable0.Redshift);
L = L0*(1-xd) + L1*xd;
}
}
else
L = L0;
return L;
}
double compute_LambdaTotal(double T,double nH,double HeMfrac,double Z, double Zsol)
{
double Lf;
double Lm;
double nEnH;
double nEnHsol;
double L;
Lf = compute_LambdaMetalFree(T,nH,HeMfrac);
Lm = compute_LambdaMetal(T,nH);
nEnH = compute_ElectronHydrogenRatio(T,nH,HeMfrac);
nEnHsol = compute_ElectronHydrogenRatioSolar(T,nH);
L= Lf + Lm * (nEnH/nEnHsol)*(Z/Zsol);
return L;
}
/*****************************************************************/
/* python interface */
/*****************************************************************/
#ifdef PY_INTERFACE
static PyObject * CoolingTables_Init(PyObject *self, PyObject *args, PyObject *kwds)
{
//static char *kwlist[] = {"directory", NULL};
PyObject *directory=PyString_FromString(TABLE_DIR);
/* this fails with python2.6, I do not know why ??? */
/*
if (! PyArg_ParseTupleAndKeywords(args, kwds, "|O",kwlist,&directory))
{
PyErr_SetString(PyExc_ValueError,"Init, error in parsing arguments.");
return NULL;
}
*/
if (! PyArg_ParseTuple(args, "O",&directory))
{
PyErr_SetString(PyExc_ValueError,"Init, error in parsing arguments.");
return NULL;
}
if (!PyString_Check(directory))
{
PyErr_SetString(PyExc_ValueError,"Argument must be a string.");
return NULL;
}
InitWiersmaCooling(PyString_AsString(directory));
return Py_BuildValue("i",1);
}
static PyObject * CoolingTables_Free(PyObject *self, PyObject *args, PyObject *kwds)
{
Free();
return Py_BuildValue("i",1);
}
static PyObject * CoolingTables_setTablesFromRedshift(self, args)
PyObject *self;
PyObject *args;
{
double Redshift;
if (! PyArg_ParseTuple(args, "d",&Redshift))
{
PyErr_SetString(PyExc_ValueError,"Init, error in parsing arguments.");
return NULL;
}
setTablesFromRedshift(Redshift);
return Py_BuildValue("i",1);
}
static PyObject * CoolingTables_SetRedshiftInterpolationOn(PyObject *self, PyObject *args, PyObject *kwds)
{
SetRedshiftInterpolationOn();
return Py_BuildValue("i",1);
}
static PyObject * CoolingTables_SetRedshiftInterpolationOff(PyObject *self, PyObject *args, PyObject *kwds)
{
SetRedshiftInterpolationOff();
return Py_BuildValue("i",1);
}
static PyObject * CoolingTables_compute_LambdaMetalFree(PyObject *self, PyObject *args, PyObject *kwds)
{
double nH,HeMfrac,T;
double L;
if (! PyArg_ParseTuple(args, "ddd",&T,&nH,&HeMfrac))
{
PyErr_SetString(PyExc_ValueError,"Init, error in parsing arguments.");
return NULL;
}
L = compute_LambdaMetalFree(T,nH,HeMfrac);
return Py_BuildValue("d",L);
}
static PyObject * CoolingTables_compute_LambdaMetal(PyObject *self, PyObject *args, PyObject *kwds)
{
double nH,T;
double L;
if (! PyArg_ParseTuple(args, "dd",&T,&nH))
{
PyErr_SetString(PyExc_ValueError,"Init, error in parsing arguments.");
return NULL;
}
L = compute_LambdaMetal(T,nH);
return Py_BuildValue("d",L);
}
static PyObject * CoolingTables_compute_LambdaTotal(PyObject *self, PyObject *args, PyObject *kwds)
{
double nH,HeMfrac,T;
double L;
//double Lf,Lm;
//double nEnH,nEnHsol;
double Z,Zsol;
if (! PyArg_ParseTuple(args, "ddddd",&T,&nH,&HeMfrac,&Z,&Zsol))
{
PyErr_SetString(PyExc_ValueError,"Init, error in parsing arguments.");
return NULL;
}
/*
Lf = compute_LambdaMetalFree(T,nH,HeMfrac);
Lm = compute_LambdaMetal(T,nH);
nEnH = compute_ElectronHydrogenRatio(T,nH,HeMfrac);
nEnHsol = compute_ElectronHydrogenRatioSolar(T,nH);
L= Lf + Lm * (nEnH/nEnHsol)*(Z/Zsol);
*/
L = compute_LambdaTotal(T,nH,HeMfrac,Z,Zsol);
return Py_BuildValue("d",L);
}
static PyObject * CoolingTables_compute_ElectronHydrogenRatio(PyObject *self, PyObject *args, PyObject *kwds)
{
double nH,T,HeMfrac;
double nEnH;
if (! PyArg_ParseTuple(args, "ddd",&T,&nH,&HeMfrac))
{
PyErr_SetString(PyExc_ValueError,"Init, error in parsing arguments.");
return NULL;
}
nEnH = compute_ElectronHydrogenRatio(T,nH,HeMfrac);
return Py_BuildValue("d",nEnH);
}
static PyObject * CoolingTables_compute_ElectronHydrogenRatioSolar(PyObject *self, PyObject *args, PyObject *kwds)
{
double nH,T;
double nEnH;
if (! PyArg_ParseTuple(args, "dd",&T,&nH))
{
PyErr_SetString(PyExc_ValueError,"Init, error in parsing arguments.");
return NULL;
}
nEnH = compute_ElectronHydrogenRatioSolar(T,nH);
return Py_BuildValue("d",nEnH);
}
/*****************************************************/
/* Method table */
/*****************************************************/
static PyMethodDef CoolingTablesMethods[] = {
{"Init", CoolingTables_Init, METH_VARARGS,
"Init Coolings tables."},
{"Free", CoolingTables_Free, METH_VARARGS,
"Free memory."},
{"setTablesFromRedshift", CoolingTables_setTablesFromRedshift, METH_VARARGS,
"Set table from redshift."},
{"setRedshiftInterpolationOn", CoolingTables_SetRedshiftInterpolationOn, METH_VARARGS,
"Switch on the redshift interpolation."},
{"setRedshiftInterpolationOff", CoolingTables_SetRedshiftInterpolationOff, METH_VARARGS,
"Switch off the redshift interpolation."},
{"compute_LambdaMetalFree", CoolingTables_compute_LambdaMetalFree, METH_VARARGS,
"compute LambdaMetalFree from tables."},
{"compute_LambdaMetal", CoolingTables_compute_LambdaMetal, METH_VARARGS,
"compute LambdaMetal from tables."},
{"compute_LambdaTotal", CoolingTables_compute_LambdaTotal, METH_VARARGS,
"compute Lambda Total."},
{"compute_ElectronHydrogenRatio", CoolingTables_compute_ElectronHydrogenRatio, METH_VARARGS,
"Density of electron over density of hydrogen."},
{"compute_ElectronHydrogenRatioSolar", CoolingTables_compute_ElectronHydrogenRatioSolar, METH_VARARGS,
"Density of electron over density of hydrogen at solar metallicities."},
{NULL, NULL, 0, NULL} /* Sentinel */
};
/*****************************************************/
/* initialization */
/*****************************************************/
PyMODINIT_FUNC initCoolingWiersma(void) {
(void) Py_InitModule("CoolingWiersma", CoolingTablesMethods);
}
#endif //PY_INTERFACE
#endif //COOLING_WIERSMA
#endif //COOLING
c4science · Help