Homec4science
Diffusion Gear (master) Browse

Gear/src/cosmictime.c
1c0fbaf83a18master

Last Change
Raw File

cosmictime.c
View Options

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include "allvars.h"
#include "proto.h"
#ifdef COSMICTIME
#define MINSCALINGFACTOR 1e-5
/*! \file driftfac.c
* \brief compute loop-up tables for prefactors in cosmological integration
*/
static double logTimeBegin;
static double logTimeMax;
static double logFullTimeBegin;
static double logFullTimeMax;
/*! Integration kernel for cosmictime factor computation.
*/
double cosmictime_integ(double a, void *param)
{
double h;
h = All.Omega0 / (a * a * a) + (1 - All.Omega0 - All.OmegaLambda) / (a * a) + All.OmegaLambda;
h = All.Hubble * sqrt(h);
return 1 / (h * a);
}
/*! This function computes look-up tables for factors needed in
* cosmological integrations. The (simple) integrations are carried out
* with the GSL library. Separate factors are computed for the "drift",
* and the gravitational and hydrodynamical "kicks". The lookup-table is
* used for reasons of speed.
*/
void init_cosmictime_table(void)
{
#define WORKSIZE 100000
int i;
double result, abserr;
gsl_function F;
gsl_integration_workspace *workspace;
logTimeBegin = log(All.TimeBegin);
logTimeMax = log(All.TimeMax);
workspace = gsl_integration_workspace_alloc(WORKSIZE);
for(i = 0; i < COSMICTIME_TABLE_LENGTH; i++)
{
F.function = &cosmictime_integ;
gsl_integration_qag(&F, exp(logTimeBegin), exp(logTimeBegin + ((logTimeMax - logTimeBegin) / COSMICTIME_TABLE_LENGTH) * (i + 1)), All.Hubble, /* note: absolute error just a dummy */
1.0e-8, WORKSIZE, GSL_INTEG_GAUSS41, workspace, &result, &abserr);
CosmicTimeTable[i] = result;
}
gsl_integration_workspace_free(workspace);
}
/*! Same as init_cosmictime_table, except that
* it uses a full range in a
* This function is used to compute an absolute cosmic time value.
*/
void init_full_cosmictime_table(void)
{
#define WORKSIZE 100000
int i;
double result, abserr;
gsl_function F;
gsl_integration_workspace *workspace;
/* here, we take the full range */
logFullTimeBegin = log(MINSCALINGFACTOR);
logFullTimeMax = log(1.0);
workspace = gsl_integration_workspace_alloc(WORKSIZE);
for(i = 0; i < COSMICTIME_TABLE_LENGTH; i++)
{
F.function = &cosmictime_integ;
//gsl_integration_qag(&F, exp(logFullTimeBegin), exp(logFullTimeBegin + ((logFullTimeMax - logFullTimeBegin) / COSMICTIME_TABLE_LENGTH) * (i + 1)), All.Hubble, /* note: absolute error just a dummy */
// 1.0e-8, WORKSIZE, GSL_INTEG_GAUSS41, workspace, &result, &abserr);
gsl_integration_qag(&F, exp(logFullTimeBegin + ((logFullTimeMax - logFullTimeBegin) / COSMICTIME_TABLE_LENGTH) * (i + 1)) , exp(logFullTimeMax) , All.Hubble, /* note: absolute error just a dummy */
1.0e-8, WORKSIZE, GSL_INTEG_GAUSS41, workspace, &result, &abserr);
FullCosmicTimeTable[i] = -result;
}
/* normalize put CosmicTime=0 -> MINSCALINGFACTOR */
result = FullCosmicTimeTable[0];
for(i = 0; i < COSMICTIME_TABLE_LENGTH; i++)
{
FullCosmicTimeTable[i] = FullCosmicTimeTable[i] - result;
}
gsl_integration_workspace_free(workspace);
double Time;
double a0,a1,t0,t1;
int j;
for(i = 0; i < COSMICTIME_TABLE_LENGTH; i++)
{
Time = FullCosmicTimeTable[0] + ((FullCosmicTimeTable[COSMICTIME_TABLE_LENGTH-1] - FullCosmicTimeTable[0]) / COSMICTIME_TABLE_LENGTH) * (i + 1) ;
j=0;
while ( FullCosmicTimeTable[j] < Time )
{
if (j==(COSMICTIME_TABLE_LENGTH-1))
break;
j++;
}
a0 = exp(logFullTimeBegin + ((logFullTimeMax - logFullTimeBegin) / COSMICTIME_TABLE_LENGTH) * (j + 0));
a1 = exp(logFullTimeBegin + ((logFullTimeMax - logFullTimeBegin) / COSMICTIME_TABLE_LENGTH) * (j + 1));
t0 = FullCosmicTimeTable[j-1];
t1 = FullCosmicTimeTable[j];
FullCosmicTimeTableInv[i]=(Time-t0)/(t1-t0)*(a1-a0) + a0;
}
}
/*! Get da from a given dt, at scaling factor a.
*/
double get_da_from_dt(double a, double dt)
{
double hubble_a;
double da;
hubble_a = All.Omega0 / (All.Time * All.Time * All.Time)
+ (1 - All.Omega0 - All.OmegaLambda) / (All.Time * All.Time) + All.OmegaLambda;
da = dt * a * hubble_a;
return da;
}
/*! Get dt from a given da, at scaling factor a.
*/
double get_dt_from_da(double a, double da)
{
double hubble_a;
double dt;
hubble_a = All.Omega0 / (All.Time * All.Time * All.Time)
+ (1 - All.Omega0 - All.OmegaLambda) / (All.Time * All.Time) + All.OmegaLambda;
dt = da / ( a * hubble_a);
return dt;
}
/*! This function return cosmic time difference for time in the time line
*/
double get_cosmictime_difference(int time0, int time1)
{
double a1, a2, df1, df2, u1, u2;
int i1, i2;
/* note: will only be called for cosmological integration */
a1 = logTimeBegin + time0 * All.Timebase_interval;
a2 = logTimeBegin + time1 * All.Timebase_interval;
u1 = (a1 - logTimeBegin) / (logTimeMax - logTimeBegin) * COSMICTIME_TABLE_LENGTH;
i1 = (int) u1;
if(i1 >= COSMICTIME_TABLE_LENGTH)
i1 = COSMICTIME_TABLE_LENGTH - 1;
if(i1 <= 1)
df1 = u1 * CosmicTimeTable[0];
else
df1 = CosmicTimeTable[i1 - 1] + (CosmicTimeTable[i1] - CosmicTimeTable[i1 - 1]) * (u1 - i1);
u2 = (a2 - logTimeBegin) / (logTimeMax - logTimeBegin) * COSMICTIME_TABLE_LENGTH;
i2 = (int) u2;
if(i2 >= COSMICTIME_TABLE_LENGTH)
i2 = COSMICTIME_TABLE_LENGTH - 1;
if(i2 <= 1)
df2 = u2 * CosmicTimeTable[0];
else
df2 = CosmicTimeTable[i2 - 1] + (CosmicTimeTable[i2] - CosmicTimeTable[i2 - 1]) * (u2 - i2);
return df2 - df1;
}
/*! This function return cosmic time difference for time in scaling factor
*/
double get_cosmictime_difference_from_a(double a1, double a2)
{
double df1, df2, u1, u2;
int i1, i2;
/* note: will only be called for cosmological integration */
//a1 = logTimeBegin + time0 * All.Timebase_interval;
//a2 = logTimeBegin + time1 * All.Timebase_interval;
u1 = (a1 - logTimeBegin) / (logTimeMax - logTimeBegin) * COSMICTIME_TABLE_LENGTH;
i1 = (int) u1;
if(i1 >= COSMICTIME_TABLE_LENGTH)
i1 = COSMICTIME_TABLE_LENGTH - 1;
if(i1 <= 1)
df1 = u1 * CosmicTimeTable[0];
else
df1 = CosmicTimeTable[i1 - 1] + (CosmicTimeTable[i1] - CosmicTimeTable[i1 - 1]) * (u1 - i1);
u2 = (a2 - logTimeBegin) / (logTimeMax - logTimeBegin) * COSMICTIME_TABLE_LENGTH;
i2 = (int) u2;
if(i2 >= COSMICTIME_TABLE_LENGTH)
i2 = COSMICTIME_TABLE_LENGTH - 1;
if(i2 <= 1)
df2 = u2 * CosmicTimeTable[0];
else
df2 = CosmicTimeTable[i2 - 1] + (CosmicTimeTable[i2] - CosmicTimeTable[i2 - 1]) * (u2 - i2);
return df2 - df1;
}
/*! This function return cosmic time from a given scaling factor
*/
double get_CosmicTime_from_a(double a1)
{
double df1, df2, u1, u2;
double a2=1.0;
int i1, i2;
a1 = dmax(a1,MINSCALINGFACTOR);
u1 = (log(a1) - logFullTimeBegin) / (logFullTimeMax - logFullTimeBegin) * COSMICTIME_TABLE_LENGTH;
i1 = (int) u1;
if(i1 >= COSMICTIME_TABLE_LENGTH)
i1 = COSMICTIME_TABLE_LENGTH - 1;
if(i1 < 1)
df1 = FullCosmicTimeTable[0];
else
df1 = FullCosmicTimeTable[i1 - 1] + (FullCosmicTimeTable[i1] - FullCosmicTimeTable[i1 - 1]) * (u1 - i1);
return df1;
}
/*! This function return the scaling factor from a given cosmic time
*/
double get_a_from_CosmicTime(double t)
{
double a1;
double u1;
int i1;
u1 = (t - FullCosmicTimeTable[0]) / (FullCosmicTimeTable[COSMICTIME_TABLE_LENGTH-1] - FullCosmicTimeTable[0]) * COSMICTIME_TABLE_LENGTH;
i1 = (int) u1;
if(i1 < 1)
a1 = FullCosmicTimeTableInv[0];
else
a1 = FullCosmicTimeTableInv[i1 - 1] + (FullCosmicTimeTableInv[i1] - FullCosmicTimeTableInv[i1 - 1]) * (u1 - i1);
return a1;
}
/*! This function return cosmic time from a given scaling factor
*/
double get_Redshift_from_a(double a)
{
return 1./a - 1;
}
/*! This function return the scaling factor from a given cosmic time
*/
double get_a_from_Redshift(double z)
{
return 1/(z+1);
}
/****************************************************************************************/
/*
/*
/*
/* PYTHON INTERFACE
/*
/*
/*
/****************************************************************************************/
#ifdef PY_INTERFACE
#include <Python.h>
#include <numpy/arrayobject.h>
static PyObject * cosmic_InitDefaultParameters(void)
{
/* list of Gadget parameters */
/* System of units */
All.UnitLength_in_cm = 3.085e+21; /* 1.0 kpc */
All.UnitMass_in_g = 1.989e+43; /* 1.0e10 solar masses */
All.UnitVelocity_in_cm_per_s = 20725573.785998672; /* 207 km/sec */
All.ComovingIntegrationOn=1;
All.TimeBegin=0.1;
All.TimeMax=1;
All.Omega0 = 0;
All.OmegaLambda = 0;
All.OmegaBaryon = 0;
All.HubbleParam = 0;
All.BoxSize = 0;
return Py_BuildValue("i",1);
}
static PyObject * SetParameters(PyObject *dict)
{
PyObject *key;
PyObject *value;
int ivalue;
float fvalue;
double dvalue;
/* check that it is a PyDictObject */
if(!PyDict_Check(dict))
{
PyErr_SetString(PyExc_AttributeError, "argument is not a dictionary.");
return NULL;
}
if (PyDict_Size(dict)==0)
return Py_BuildValue("i",0);
Py_ssize_t pos=0;
while(PyDict_Next(dict,&pos,&key,&value))
{
if(PyString_Check(key))
{
/* System of units */
if(strcmp(PyString_AsString(key), "UnitLength_in_cm")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.UnitLength_in_cm = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "UnitMass_in_g")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.UnitMass_in_g = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "UnitVelocity_in_cm_per_s")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.UnitVelocity_in_cm_per_s = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "ComovingIntegrationOn")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.ComovingIntegrationOn = PyInt_AsLong(value);
}
/* Caracteristics of run */
if(strcmp(PyString_AsString(key), "TimeBegin")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TimeBegin = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "TimeMax")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TimeMax = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "Omega0")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.Omega0 = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "OmegaLambda")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.OmegaLambda = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "OmegaBaryon")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.OmegaBaryon = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "HubbleParam")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.HubbleParam = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "BoxSize")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.BoxSize = PyFloat_AsDouble(value);
}
}
}
return Py_BuildValue("i",1);
}
static PyObject * cosmic_SetParameters(PyObject *self, PyObject *args)
{
PyObject *dict;
/* here, we can have either arguments or dict directly */
if(PyDict_Check(args))
{
dict = args;
}
else
{
if (! PyArg_ParseTuple(args, "O",&dict))
return NULL;
}
SetParameters(dict);
return Py_BuildValue("i",1);
}
static PyObject * cosmic_GetParameters()
{
PyObject *dict;
PyObject *key;
PyObject *value;
dict = PyDict_New();
/* System of units */
key = PyString_FromString("UnitLength_in_cm");
value = PyFloat_FromDouble(All.UnitLength_in_cm);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("UnitMass_in_g");
value = PyFloat_FromDouble(All.UnitMass_in_g);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("UnitVelocity_in_cm_per_s");
value = PyFloat_FromDouble(All.UnitVelocity_in_cm_per_s);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("ComovingIntegrationOn");
value = PyInt_FromLong(All.ComovingIntegrationOn);
PyDict_SetItem(dict,key,value);
/* Caracteristics of run */
key = PyString_FromString("TimeBegin");
value = PyFloat_FromDouble(All.TimeBegin);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("TimeMax");
value = PyFloat_FromDouble(All.TimeMax);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("Omega0");
value = PyFloat_FromDouble(All.Omega0);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("OmegaLambda");
value = PyFloat_FromDouble(All.OmegaLambda);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("OmegaBaryon");
value = PyFloat_FromDouble(All.OmegaBaryon);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("HubbleParam");
value = PyFloat_FromDouble(All.HubbleParam);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("BoxSize");
value = PyFloat_FromDouble(All.BoxSize);
PyDict_SetItem(dict,key,value);
return Py_BuildValue("O",dict);
}
static PyObject *
cosmic_Init(PyObject *self, PyObject *args, PyObject *kwds)
{
if(All.ComovingIntegrationOn)
All.Timebase_interval = (log(All.TimeMax) - log(All.TimeBegin)) / TIMEBASE;
else
All.Timebase_interval = (All.TimeMax - All.TimeBegin) / TIMEBASE;
All.UnitTime_in_s = All.UnitLength_in_cm / All.UnitVelocity_in_cm_per_s;
All.Hubble = HUBBLE * All.UnitTime_in_s;
if(All.ComovingIntegrationOn)
init_cosmictime_table();
init_full_cosmictime_table();
return Py_BuildValue("O",Py_None);
}
static PyObject *
cosmic_get_cosmictime_difference(PyObject *self, PyObject *args, PyObject *kwds)
{
double a1,a2,t;
if (!PyArg_ParseTuple(args, "dd",&a1,&a2))
return NULL;
t = get_cosmictime_difference_from_a(a1,a2);
return Py_BuildValue("d",t);
}
static PyObject *
cosmic_get_CosmicTime_from_a(PyObject *self, PyObject *args, PyObject *kwds)
{
double a,t;
if (!PyArg_ParseTuple(args, "d",&a))
return NULL;
t = get_CosmicTime_from_a(a);
return Py_BuildValue("d",t);
}
static PyObject *
cosmic_get_a_from_CosmicTime(PyObject *self, PyObject *args, PyObject *kwds)
{
double a,t;
if (!PyArg_ParseTuple(args, "d",&t))
return NULL;
a = get_a_from_CosmicTime(t);
return Py_BuildValue("d",a);
}
static PyObject *
cosmic_get_Redshift_from_a(PyObject *self, PyObject *args, PyObject *kwds)
{
double a,z;
if (!PyArg_ParseTuple(args, "d",&a))
return NULL;
z = get_Redshift_from_a(a);
return Py_BuildValue("d",z);
}
static PyObject *
cosmic_get_a_from_Redshift(PyObject *self, PyObject *args, PyObject *kwds)
{
double a,z;
if (!PyArg_ParseTuple(args, "d",&z))
return NULL;
a = get_a_from_Redshift(z);
return Py_BuildValue("d",a);
}
static PyObject *
cosmic_get_CosmicTimeTable(PyObject *self, PyObject *args, PyObject *kwds)
{
int i;
double a;
double t;
PyArrayObject *va;
PyArrayObject *vt;
npy_intp ld[1];
ld[0] = COSMICTIME_TABLE_LENGTH;
va = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
vt = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
for(i = 0; i < COSMICTIME_TABLE_LENGTH; i++)
{
a = exp(logTimeBegin + ((logTimeMax - logTimeBegin) / COSMICTIME_TABLE_LENGTH) * (i + 1));
t = CosmicTimeTable[i];
*(double *) (va->data + i*(va->strides[0])) = a;
*(double *) (vt->data + i*(vt->strides[0])) = t;
}
return Py_BuildValue("(OO)",va,vt);
}
static PyObject *
cosmic_get_FullCosmicTimeTable(PyObject *self, PyObject *args, PyObject *kwds)
{
int i;
double a;
double t;
PyArrayObject *va;
PyArrayObject *vt;
npy_intp ld[1];
ld[0] = COSMICTIME_TABLE_LENGTH;
va = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
vt = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
for(i = 0; i < COSMICTIME_TABLE_LENGTH; i++)
{
a = exp(logFullTimeBegin + ((logFullTimeMax - logFullTimeBegin) / COSMICTIME_TABLE_LENGTH) * (i + 1));
t = FullCosmicTimeTable[i];
*(double *) (va->data + i*(va->strides[0])) = a;
*(double *) (vt->data + i*(vt->strides[0])) = t;
}
return Py_BuildValue("(OO)",va,vt);
}
static PyObject *
cosmic_get_FullCosmicTimeTableInv(PyObject *self, PyObject *args, PyObject *kwds)
{
int i;
double a;
double t;
PyArrayObject *va;
PyArrayObject *vt;
npy_intp ld[1];
ld[0] = COSMICTIME_TABLE_LENGTH;
va = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
vt = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
for(i = 0; i < COSMICTIME_TABLE_LENGTH; i++)
{
t = FullCosmicTimeTable[0] + ((FullCosmicTimeTable[COSMICTIME_TABLE_LENGTH-1] - FullCosmicTimeTable[0]) / COSMICTIME_TABLE_LENGTH) * (i + 1) ;
a = FullCosmicTimeTableInv[i];
*(double *) (va->data + i*(va->strides[0])) = a;
*(double *) (vt->data + i*(vt->strides[0])) = t;
}
return Py_BuildValue("(OO)",va,vt);
}
static PyMethodDef cosmicMethods[] = {
{"Init", (PyCFunction)cosmic_Init, METH_VARARGS,
"Initialise the cosmic module"},
{"InitDefaultParameters", (PyCFunction)cosmic_InitDefaultParameters, METH_VARARGS,
"Init default parameters"},
{"SetParameters", (PyCFunction)cosmic_SetParameters, METH_VARARGS,
"Set gadget parameters"},
{"GetParameters", (PyCFunction)cosmic_GetParameters, METH_VARARGS,
"get some gadget parameters"},
{"get_cosmictime_difference", (PyCFunction)cosmic_get_cosmictime_difference, METH_VARARGS,
"get_cosmictime_difference, input in scaling factor"},
{"get_CosmicTime_from_a", (PyCFunction)cosmic_get_CosmicTime_from_a, METH_VARARGS,
"get the cosmic time value from a given scaling factor"},
{"get_a_from_CosmicTime", (PyCFunction)cosmic_get_a_from_CosmicTime, METH_VARARGS,
"get the scaling factor value from a given cosmic time"},
{"get_Redshift_from_a", (PyCFunction)cosmic_get_Redshift_from_a, METH_VARARGS,
"get the redshift value from a given scaling factor"},
{"get_a_from_Redshift", (PyCFunction)cosmic_get_a_from_Redshift, METH_VARARGS,
"get the scaling factor value from a given redshift"},
{"get_CosmicTimeTable", (PyCFunction)cosmic_get_CosmicTimeTable, METH_VARARGS,
"get the cosmic time table and the corresponding scaling factors"},
{"get_FullCosmicTimeTable", (PyCFunction)cosmic_get_FullCosmicTimeTable, METH_VARARGS,
"get the full cosmic time table and the corresponding scaling factors"},
{"get_FullCosmicTimeTableInv", (PyCFunction)cosmic_get_FullCosmicTimeTableInv, METH_VARARGS,
"get the cosmic time table (inv) and the corresponding time"},
{NULL, NULL, 0, NULL} /* Sentinel */
};
void initcosmic(void)
{
(void) Py_InitModule("cosmic", cosmicMethods);
import_array();
}
#endif /* PY_INTERFACE */
#endif
c4science · Help