Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85237054
cosmictime.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
Fri, Sep 27, 17:31
Size
20 KB
Mime Type
text/x-c
Expires
Sun, Sep 29, 17:31 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
21144699
Attached To
rGEAR Gear
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;
}
}
/*! 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
Event Timeline
Log In to Comment