Page MenuHomec4science
Files F72682450

cosmictime.c
No OneTemporary
Actions

File Metadata

Created
Tue, Jul 16, 13:23

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
/*! \file driftfac.c
* \brief compute loop-up tables for prefactors in cosmological integration
*/
static double logTimeBegin;
static double logTimeMax;
/*! 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);
/*
check
*/
//double a0,a1;
//int ti0,ti1;
//
//a0 = 0.7;
//a1 = 1.0;
//
//ti0 = log(a0 / All.TimeBegin) / All.Timebase_interval;
//ti1 = log(a1 / All.TimeBegin) / All.Timebase_interval;
//
//printf("- - - %g ti0=%d ti1=%d\n",get_cosmictime_difference(ti0,ti1),ti0,ti1);
//endrun(-7);
/*
end of check check
*/
}
/*! 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;
}
#endif

Event Timeline

c4science · Help