Skip to content

Commit 058f456

Browse files
rajeshmondal18rajeshmondal18
authoredJan 15, 2019
Add files via upload
0 parents  commit 058f456

10 files changed

+2406
-0
lines changed
 

‎allotarrays.c

+51
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
#include<stdlib.h>
2+
#include<fftw3.h>
3+
4+
float ***allocate_fftwf_3d(long N1,long N2,long N3)
5+
{
6+
long ii,jj;
7+
long asize,index;
8+
float ***phia, *phi;
9+
10+
phia=(float ***) fftwf_malloc (N1 * sizeof(float **));
11+
12+
13+
for(ii=0;ii<N1;++ii)
14+
phia[ii]=(float **) fftwf_malloc (N2 * sizeof(float *));
15+
16+
asize = N1*N2;
17+
asize = asize*N3;
18+
19+
if(!(phi = (float *) calloc(asize,sizeof(float))))
20+
{
21+
printf("error in allocate_fftwf_3d");
22+
exit(0);
23+
}
24+
25+
for(ii=0;ii<N1;++ii)
26+
for(jj=0;jj<N2;++jj)
27+
{
28+
index = N2*N3;
29+
index = index*ii + N3*jj;
30+
phia[ii][jj]=phi+ index;
31+
}
32+
return(phia);
33+
}
34+
35+
float **allocate_float_2d(long N1,int N2) /* Typically N1 is number of particles and N2 is number of dimensions namely 3 */
36+
{
37+
float **xxa, *xx;
38+
long ii;
39+
40+
xxa=(float**)malloc(N1 * sizeof(float*));
41+
if(!(xx = (float *) calloc((size_t)(N1*N2),sizeof(float))))
42+
{
43+
printf("error in allocate_float_2d");
44+
exit(0);
45+
}
46+
47+
for(ii=0;ii<N1;++ii)
48+
xxa[ii]=xx + N2*ii ;
49+
50+
return(xxa);
51+
}

‎funcs.c

+200
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,200 @@
1+
/* ----------------------------------------------------
2+
This has various functions - needed in cosmology
3+
4+
Df(aa) - calculates normalized D(a) (growing mode) at scale factor
5+
aa=1/(1+z) D(1) = 1
6+
ff(aa) - calculates f(a) (dlnD/dlna) at scale factor aa=1/(1+z)
7+
8+
ONLY works for FRW universe with ordinary matter (vomegam), cosmological
9+
constant(vomegalam) and curvature(1.0-vomegam-vomegalam).
10+
11+
Also calculates a random number (routine copied from NRC) ran1(long *)
12+
13+
/*---------------------------------------------------------------------------*/
14+
15+
#include<stdio.h>
16+
#include<math.h>
17+
#define FUNC(x) ((*func)(x))
18+
#define EPS1 1.0e-4 /* Fractional accuracy */
19+
#define JMAX 20
20+
21+
extern float vhh,vomegam,vomegalam;
22+
/* Hubble parameter,Omega-Matter,Cosmological Const */
23+
24+
25+
/*---------------------------------------------------------------------------*/
26+
/* Does integration using trapezoidal rule */
27+
/*---------------------------------------------------------------------------*/
28+
29+
float trapzd(float (*func)(float), float a, float b, int n)
30+
{
31+
float x,tnm,sum,del;
32+
static float s;
33+
int it,j;
34+
35+
if (n == 1) {
36+
return (s=0.5*(b-a)*(FUNC(a)+FUNC(b)));
37+
} else {
38+
for (it=1,j=1;j<n-1;j++) it <<= 1;
39+
tnm=it;
40+
del=(b-a)/tnm;
41+
x=a+0.5*del;
42+
for (sum=0.0,j=1;j<=it;j++,x+=del) sum += FUNC(x);
43+
s=0.5*(s+(b-a)*sum/tnm);
44+
return s;
45+
}
46+
}
47+
48+
49+
float qtrap(float (*func)(float), float a, float b)
50+
{
51+
float trapzd(float (*func)(float), float a, float b, int n);
52+
int j;
53+
float s,olds;
54+
55+
float dummy;
56+
57+
olds = -1.0e30;
58+
for (j=1;j<=JMAX;j++) {
59+
s=trapzd(func,a,b,j);
60+
if (fabs(s-olds) < EPS1*fabs(olds)) return s;
61+
olds=s;
62+
}
63+
printf("\nToo many steps in routine qtrap...exiting routine...\n");
64+
return 0.0; /* Never gets here normally */
65+
}
66+
/*---------------------------------------------------------------------------*/
67+
/*Trapezoidal Rule done */
68+
/*---------------------------------------------------------------------------*/
69+
70+
71+
float Hf(float aa)
72+
{
73+
return(sqrt(vomegam*pow(aa,-3.) + (1.0-vomegam-vomegalam)*pow(aa,-2.)+ vomegalam));
74+
}
75+
76+
float qf(float aa)
77+
{
78+
return(0.5*(vomegam*pow(aa,-3.0) - 2.0*vomegalam)/(Hf(aa)*Hf(aa)));
79+
}
80+
81+
float Integrandf(float aa)
82+
{
83+
return(pow(aa*Hf(aa),-3.0));
84+
}
85+
86+
float Integralf(float aa)
87+
{
88+
/* The Integrand blows up at 0.0 so we take a small step 0.0001 */
89+
return( qtrap(Integrandf,0.00001,aa) );
90+
}
91+
92+
/*---------------------------------------------------------------------------*/
93+
/*Returns the growing mode of denisty perturbation normalized to unity at present */
94+
/*---------------------------------------------------------------------------*/
95+
96+
float Df(float aa)
97+
{
98+
return( Hf(aa)*Integralf(aa)/( Hf(1.0)*Integralf(1.0) ) );
99+
}
100+
/*---------------------------------------------------------------------------*/
101+
102+
103+
104+
/*---------------------------------------------------------------------------*/
105+
/*Returns the logarithmic derivative of the growing mode of denisty perturbatio */
106+
/*---------------------------------------------------------------------------*/
107+
108+
float ff(float aa)
109+
{
110+
return( (1.0/( aa*aa*pow(Hf(aa),3.0)*Integralf(aa) ) ) - (1.0 + qf(aa) ) );
111+
}
112+
/*---------------------------------------------------------------------------*/
113+
114+
115+
/*---------------------------------------------------------------------------*/
116+
/* function ran1 copied from Numerical Recipes*/
117+
/*---------------------------------------------------------------------------*/
118+
119+
#define IA 16807
120+
#define IM 2147483647
121+
#define AM (1.0/IM)
122+
#define IQ 127773
123+
#define IR 2836
124+
#define NTAB 32
125+
#define NDIV (1+(IM-1)/NTAB)
126+
#define EPS 1.2e-7
127+
#define RNMX (1.0-EPS)
128+
129+
float ran1(long *idum)
130+
{
131+
int j;
132+
long k;
133+
static long iy=0;
134+
static long iv[NTAB];
135+
float temp;
136+
137+
if (*idum <= 0 || !iy) {
138+
if (-(*idum) < 1) *idum=1;
139+
else *idum = -(*idum);
140+
for (j=NTAB+7;j>=0;j--) {
141+
k=(*idum)/IQ;
142+
*idum=IA*(*idum-k*IQ)-IR*k;
143+
if (*idum < 0) *idum += IM;
144+
if (j < NTAB) iv[j] = *idum;
145+
}
146+
iy=iv[0];
147+
}
148+
k=(*idum)/IQ;
149+
*idum=IA*(*idum-k*IQ)-IR*k;
150+
if (*idum < 0) *idum += IM;
151+
j=iy/NDIV;
152+
iy=iv[j];
153+
iv[j] = *idum;
154+
if ((temp=AM*iy) > RNMX) return RNMX;
155+
else return temp;
156+
}
157+
#undef IA
158+
#undef IM
159+
#undef AM
160+
#undef IQ
161+
#undef IR
162+
#undef NTAB
163+
#undef NDIV
164+
#undef EPS
165+
#undef RNMX
166+
#undef EPS1
167+
#undef FUNC
168+
#undef JMAX
169+
/*---------------------------------------------------------------------------*/
170+
/* ran1 done*/
171+
/*---------------------------------------------------------------------------*/
172+
173+
/*---------------------------------------------------------------------------*/
174+
/* function gasdev copied from Numerical Recipes*/
175+
/*---------------------------------------------------------------------------*/
176+
float gasdev(long *idum)
177+
{
178+
float ran1(long *idum);
179+
static int iset=0;
180+
static float gset;
181+
float fac,rsq,v1,v2;
182+
183+
if (iset == 0) {
184+
do {
185+
v1=2.0*ran1(idum)-1.0;
186+
v2=2.0*ran1(idum)-1.0;
187+
rsq=v1*v1+v2*v2;
188+
} while (rsq >= 1.0 || rsq == 0.0);
189+
fac=sqrt(-2.0*log(rsq)/rsq);
190+
gset=v1*fac;
191+
iset=1;
192+
return v2*fac;
193+
} else {
194+
iset=0;
195+
return gset;
196+
}
197+
}
198+
/*---------------------------------------------------------------------------*/
199+
/* gasdev done*/
200+
/*---------------------------------------------------------------------------*/

‎input.nbody_comp

+53
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,53 @@
1+
-100012 10
2+
0.6704 0.3183 0.6817 0.9619
3+
0.04902 0.8347
4+
1024 1024 1024 134217728 0.07
5+
0.008 2 2 0
6+
0.01
7+
15
8+
14.000 10.804 9.821 9.210 8.745 8.369 8.076 8.000 7.856 7.676 7.483 7.226 6.853 6.314 6.000
9+
10+
11+
12+
13+
8 8 8 512 0.14
14+
16 16 16 4096 0.14
15+
32 32 32 32768 0.14
16+
128 128 128 2097152 0.14
17+
192 192 192 7077888 0.14
18+
256 256 256 16777216 0.14
19+
384 384 384 56623104 0.14
20+
512 512 512 134217728 0.14
21+
640 640 640 262144000 0.14
22+
768 768 768 452984832 0.14
23+
1280 1280 1280 2097152000 0.14
24+
25+
16 16 16 512 0.07
26+
32 32 32 4096 0.07
27+
256 256 256 2097152 0.07
28+
512 512 512 16777216 0.07
29+
768 768 768 56623104 0.07
30+
1024 1024 1024 134217728 0.07
31+
1280 1280 1280 262144000 0.07
32+
1344 1344 1344 303464448 0.07
33+
1408 1408 1408 348913664 0.07
34+
1536 1536 1536 452984832 0.07
35+
2048 2048 2048 1073741824 0.07
36+
2080 2080 2080 1124864000 0.07
37+
2112 2112 2112 1177583616 0.07
38+
2144 2144 2144 1231925248 0.07
39+
40+
256 256 256 2097152 0.033
41+
512 512 512 16777216 0.033
42+
2144 2144 2144 1231925248 0.033
43+
44+
45+
2144 2144 2144 1231925248 0.017
46+
47+
48+
seed Nbin
49+
hh omega_m omega_l spectral_index
50+
omega_baryon sigma_8
51+
N1 N2 N3 MM LL
52+
aa_initial fraction_fill output_flag
53+
delta_aa aa_final

‎makefile

+60
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
LINKLIB= -L/usr/local/lib/ -fopenmp -lfftw3f_omp -lfftw3f -lm
2+
3+
#LINKLIB= -lm -lc -L/usr/local/lib/ -lfftw3
4+
5+
INCLUDE=-I/usr/local/include/
6+
7+
8+
CFLAGS=-g
9+
CC=gcc
10+
11+
12+
nbody_comp: nbody_comp.o nbody_funcs.o allotarrays.o funcs.o powerspec.o tf_fit.o
13+
$(CC) $(CFLAGS) -o nbody_comp nbody_comp.o nbody_funcs.o funcs.o allotarrays.o powerspec.o tf_fit.o $(LINKLIB)
14+
rm -rf *.o
15+
rm -rf *~
16+
17+
nbody_comp.o: nbody_comp.c
18+
$(CC) -c $(CFLAGS) $(INCLUDE) nbody_comp.c $(LINKLIB)
19+
20+
nbody_funcs.o: nbody_funcs.c
21+
$(CC) -c $(CFLAGS) $(INCLUDE) nbody_funcs.c $(LINKLIB)
22+
23+
allotarrays.o: allotarrays.c
24+
$(CC) -c $(CFLAGS) $(INCLUDE) allotarrays.c $(LINKLIB)
25+
26+
powerspec.o: powerspec.c
27+
$(CC) -c $(CFLAGS) $(INCLUDE) powerspec.c $(LINKLIB)
28+
29+
funcs.o: funcs.c
30+
$(CC) -c $(CFLAGS) $(INCLUDE) funcs.c $(LINKLIB)
31+
32+
tf_fit.o: tf_fit.c
33+
$(CC) -c $(CFLAGS) $(INCLUDE) tf_fit.c $(LINKLIB)
34+
35+
36+
fof_main: fof_main.o nbody_funcs.o allotarrays.o funcs.o powerspec.o tf_fit.o
37+
$(CC) $(CFLAGS) -o fof_main fof_main.o nbody_funcs.o funcs.o allotarrays.o powerspec.o tf_fit.o $(LINKLIB)
38+
rm -rf *.o
39+
rm -rf *~
40+
41+
fof_main.o: fof_main.c
42+
$(CC) -c $(CFLAGS) $(INCLUDE) fof_main.c $(LINKLIB)
43+
44+
45+
46+
mass_func: mass_func.o nbody_funcs.o allotarrays.o tf_fit.o powerspec.o funcs.o
47+
$(CC) $(CFLAGS) -o mass_func mass_func.o nbody_funcs.o allotarrays.o powerspec.o tf_fit.o funcs.o $(LINKLIB)
48+
rm -rf *.o
49+
rm -rf *~
50+
51+
mass_func.o: mass_func.c
52+
$(CC) -c $(CFLAGS) $(INCLUDE) mass_func.c $(LINKLIB)
53+
54+
funcs_fit.o: funcs_fit.c
55+
$(CC) -c $(CFLAGS) $(INCLUDE) funcs_fit.c $(LINKLIB)
56+
57+
58+
clean:
59+
rm -rf *.o
60+
rm -rf *~

‎nbody.h

+52
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,52 @@
1+
/* in allotarrays.c */
2+
float **allocate_float_2d(long N1,int N2);
3+
float ***allocate_fftwf_3d(long N1,long N2,long N3);
4+
float Hf(float aa),Df(float aa),ff(float aa);
5+
void Setting_Up_Memory_For_Ro(float av);
6+
void cic(float** rra);
7+
void Get_Phi(int i);
8+
void Update_v(float aa,float delta_aa,float **rra,float **vva);
9+
void Update_x(float aa,float delta_aa,float **rra,float **vva);
10+
void calpow(int f_flag,int Nbin,double* power, double* powerk, double* kmode,long *no);
11+
void calpow_k(int f_flag,float kmin, float kmax, int Nbin,double* power, double* powerk, double* kmode,long *no);
12+
13+
void Zel_move_gradphi(float av,float **rra,float **vva);
14+
void grad_phi(int ix);
15+
int write_output(char *fname,long int seed,int output_flag,float **rra,float **vva,float vaa);
16+
int write_multiout(char *fname,long int seed,int output_flag,float **rra,float **vva,float vaa);
17+
int read_output(char *fname, int read_flag,long int *seed,int *output_flag,int *in_flag,float **rra,float **vva,float *aa);
18+
19+
//Functions required to generate the power spectrum by Hu and Eisentine
20+
void TFset_parameters(float omega0hh, float f_baryon, float Tcmb);
21+
float TFfit_onek(float k, float *tf_baryon, float *tf_cdm);
22+
float Pk(float kk); // power spectrum
23+
float sigma_func(float kk); // for calculating sigma 8
24+
float simp(float (*func)(float),float a,float b,int N); // for integration
25+
void delta_fill(long*); // fills value of phases for Fourier modes
26+
27+
typedef struct
28+
{
29+
long npart[6];
30+
double mass[6];
31+
double time;
32+
double redshift;
33+
int flag_sfr;
34+
int flag_feedback;
35+
long npartTotal[6];
36+
int flag_cooling;
37+
int num_files;
38+
double BoxSize;
39+
double Omega0;
40+
double OmegaLambda;
41+
double HubbleParam;
42+
double Omegab;
43+
double sigma_8_present;
44+
long Nx;
45+
long Ny;
46+
long Nz;
47+
float LL;
48+
int output_flag;
49+
int in_flag;
50+
long int seed;
51+
char fill[256- 6*4- 6*8- 2*8- 2*4- 6*4- 2*4 - 6*8 -6*4 -8]; /* fills to 256 Bytes */
52+
} io_header;

‎nbody_comp.c

+276
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,276 @@
1+
#include<stdio.h>
2+
#include<math.h>
3+
#include<stdlib.h>
4+
#include<string.h>
5+
#include<fftw3.h>
6+
#include"nbody.h"
7+
#include<omp.h>
8+
9+
/* GLOBAL VARIABLES */
10+
11+
// cosmological parameters read from input file "input.nbody_com"
12+
13+
float vhh, // Hubble parameter in units of 100 km/s/Mpc
14+
vomegam, // Omega_matter; total matter density (baryons+CDM) parameter
15+
vomegalam, // Cosmological Constant
16+
vomegab, //Omega_baryon
17+
sigma_8_present ,// Last updated value of sigma_8 (Presently WMAP)
18+
vnn; // Spectral index of primordial Power spectrum
19+
20+
long N1,N2,N3;// box dimension (grid)
21+
int NF, // Fill every NF grid point
22+
Nbin; // Number of bins to calculate final P(k) (output)
23+
24+
float LL; // grid spacing in Mpc
25+
26+
long MM; // Number of particles
27+
// global variables (calculated )
28+
29+
int zel_flag=1, // memory allocation for zel is 3 times that for nbody
30+
fourier_flag; //for fourier transfrom
31+
32+
float DM_m, // Darm matter mass of simulation particle in 10^10 M_sun h^-1 unit
33+
norm, // normalize Pk
34+
pi=M_PI;
35+
36+
io_header header1;
37+
38+
// arrays for storing data
39+
float ***ro; // for density/potential
40+
fftwf_plan p_ro; // for FFT
41+
fftwf_plan q_ro; // for FFT
42+
43+
44+
//end of declaration of global variables
45+
46+
main()
47+
{
48+
FILE *inp, *outpp;
49+
int ii,jj;
50+
int Nflag, Noutput;
51+
int oflag, // desired output format
52+
pk_flag; //for power spectrum calculation
53+
long seed,*no; // seed for random phases
54+
float **rra, **vva; // particle position and velocity
55+
float aa,delta_aa,delta_aap,afin,aa_i,vpk;
56+
float DD,vaa,*nz;
57+
char fname[20];
58+
59+
char file[100],num[8];
60+
61+
double t,T=omp_get_wtime(), // for timing
62+
*power, *powerk, *kmode;
63+
64+
float rho_c=2.7755*1.e11; //rho_c in units of h^2 M_sun/Mpc^3
65+
66+
/*---------------------------------------------------------------------------*/
67+
/* Read input parameters for the simulation from the file "input.nbody_comp" */
68+
/*---------------------------------------------------------------------------*/
69+
inp=fopen("input.nbody_comp","r");
70+
fscanf(inp,"%ld%d",&seed,&Nbin);
71+
fscanf(inp,"%f%f%f%f",&vhh,&vomegam,&vomegalam,&vnn);
72+
fscanf(inp,"%f%f",&vomegab,&sigma_8_present);
73+
fscanf(inp,"%ld%ld%ld%ld%f",&N1,&N2,&N3,&MM,&LL);
74+
fscanf(inp,"%f",&vaa);
75+
fscanf(inp,"%d",&NF);
76+
fscanf(inp,"%d%d",&oflag,&pk_flag);
77+
fscanf(inp,"%f",&delta_aa); /* time step, final scale factor*/
78+
fscanf(inp,"%d",&Noutput);
79+
80+
nz=(float*)calloc(Noutput,sizeof(float)); // array to store Noutput
81+
82+
for(ii=0;ii<Noutput;ii++)
83+
fscanf(inp,"%f",&nz[ii]);
84+
85+
fclose(inp);
86+
87+
DD=Df(vaa); // growing mode
88+
89+
DM_m=vomegam*rho_c*pow(vhh,3.)*(1.0*N1)*(1.0*N2)*(1.0*N3)*powf(LL,3.)/(MM*1.e10); //mass per particle in 10^10 M_sun h^-1 unit
90+
printf("DM_m= %e 10^10 M_sun \t L_box=%e Mpc\n", DM_m/vhh, N1*LL);
91+
92+
/*---------------------------done inputting parameters-----------------*/
93+
94+
/*---------------------------------------------------------------------*/
95+
/* initialize power spectrum */
96+
/*---------------------------------------------------------------------*/
97+
98+
float tcmb=2.728; // CMBR temperature
99+
TFset_parameters(vomegam*vhh*vhh,vomegab/vomegam,tcmb);
100+
101+
/*---------------------------------------------------------------------*/
102+
/* done intitializing power spectrum */
103+
/*---------------------------------------------------------------------*/
104+
/* normalizing the power spectrum using sigma_8 */
105+
/*---------------------------------------------------------------------*/
106+
norm=1.;
107+
norm=simp(sigma_func,0.00001,3.5,100000);
108+
norm=pow(sigma_8_present,2.)/norm; //normalization factor for Pk(k)
109+
/*---------------------------------------------------------------------*/
110+
/* Normalization of powerspectrum done */
111+
/*---------------------------------------------------------------------*/
112+
/*---------------------------------------------------------------------*/
113+
/* allocate memory for particle positions */
114+
/*---------------------------------------------------------------------*/
115+
116+
rra= allocate_float_2d(MM,3);
117+
vva= allocate_float_2d(MM,3);
118+
119+
/*---------------------------------------------------------------------*/
120+
/*----------allocate memory for power spectrum and k modes-------------*/
121+
power = calloc((size_t)Nbin,sizeof(double));
122+
powerk = calloc((size_t)Nbin,sizeof(double));
123+
kmode = calloc((size_t)Nbin,sizeof(double));
124+
no = calloc((size_t)Nbin,sizeof(long));
125+
/*----------------------------------------------------------------*/
126+
127+
Setting_Up_Memory_For_Ro(vaa);
128+
129+
/*----------------------------------------------------------------*/
130+
131+
/* get values of delta in fourier Space */
132+
133+
t=omp_get_wtime();
134+
135+
delta_fill(&seed); //random phases for Fourier modes
136+
137+
printf("ok delta_fill time = %e\n",omp_get_wtime()-t);
138+
139+
//*------------------------------------------------------------------------*//
140+
141+
if(pk_flag==1)
142+
{
143+
t=omp_get_wtime();
144+
145+
outpp=fopen("pk.inp","w");
146+
147+
calpow(0, Nbin, power, powerk, kmode, no); // calculates power spectrum of (Delta(k))
148+
149+
for(ii=0;ii<Nbin;++ii)
150+
fprintf(outpp,"%e %e %e %ld\n",kmode[ii],power[ii],DD*DD*2.*pi*pi*powerk[ii],no[ii]);
151+
152+
fclose(outpp);
153+
154+
printf("ok cal_pow time = %e\n",omp_get_wtime()-t);
155+
}
156+
157+
//---------------------------------------------------------------------
158+
159+
t=omp_get_wtime();
160+
161+
Zel_move_gradphi(vaa,rra,vva); // move particles using ZA
162+
163+
printf("ok Zel_move time = %e\n",omp_get_wtime()-t);
164+
165+
/*--------------------------ZA complete--------------------------------*/
166+
167+
if(pk_flag==1)
168+
{
169+
t=omp_get_wtime();
170+
cic(rra);
171+
printf("ok cic = %e\n",omp_get_wtime()-t);
172+
173+
outpp=fopen("pk.zel","w");
174+
calpow(1, Nbin, power, powerk, kmode, no); // calculates power spectrum of (Delta rho)
175+
176+
for(ii=0;ii<Nbin;++ii)
177+
fprintf(outpp,"%e %e %e %ld\n", kmode[ii], power[ii], DD*DD*2.*pi*pi*powerk[ii], no[ii]);
178+
179+
fclose(outpp);
180+
}
181+
/*----------------------------------------------------------------------*/
182+
183+
aa=vaa; //initial scale factor
184+
185+
Update_v(aa, -1.*delta_aa, rra, vva); // Initialize v after ZA
186+
187+
Nflag=0;
188+
189+
//-------------------------Starting Nbody Block---------------------------//
190+
191+
for(jj=0;jj<Noutput;jj++)
192+
{
193+
afin=1.0/(nz[jj]+1.0);
194+
195+
t=omp_get_wtime();
196+
197+
/*---------------------------first step ---------------------------------*/
198+
199+
while((afin-aa)<delta_aa)
200+
{
201+
delta_aa=delta_aa*0.5;
202+
}
203+
204+
Update_v(aa, 0.5*delta_aa, rra, vva); // half step for v
205+
Update_x(aa+0.5*delta_aa, delta_aa, rra, vva); // full step for x
206+
207+
/*----------------------first step done---------------------------------*/
208+
209+
aa=aa+delta_aa;
210+
211+
while (aa + delta_aa <=afin)
212+
{
213+
Update_v( aa, delta_aa, rra, vva);
214+
Update_x( aa+0.5*delta_aa, delta_aa, rra, vva);
215+
aa=aa+delta_aa;
216+
Nflag++;
217+
} /* NBody loop end */
218+
219+
/*-----------------------last step--------------------------------------*/
220+
221+
delta_aap=afin-aa;
222+
printf("delta_aap=%e\n",delta_aap);
223+
224+
if(delta_aap>0.)
225+
{
226+
Update_v(aa, 0.5*delta_aa, rra, vva);
227+
Update_x(aa, delta_aap, rra, vva);
228+
229+
aa=aa+delta_aap;
230+
Update_v(aa, delta_aap, rra, vva);
231+
}
232+
233+
printf("(z=%3.1f) time for Nbody loop = %e\t",nz[jj], omp_get_wtime()-t);
234+
235+
/*------------------- done last step--------------------------------------*/
236+
237+
/*----------------------print final power spectrum -----------------------------------*/
238+
239+
if(pk_flag==1)
240+
{
241+
strcpy(file,"pk.nbody");
242+
sprintf(num,"%3.1f",nz[jj]);
243+
strcat(file,num);
244+
outpp=fopen(file,"w");
245+
246+
cic(rra);
247+
calpow(1, Nbin, power, powerk, kmode, no);
248+
DD=Df(aa); // growing mode
249+
250+
for(ii=0;ii<Nbin;++ii)
251+
fprintf(outpp,"%e %e %e %ld\n", kmode[ii], power[ii], DD*DD*2.*pi*pi*powerk[ii], no[ii]);
252+
253+
fclose(outpp);
254+
}
255+
/*----------------------print nbody output -----------------------------------*/
256+
257+
t=omp_get_wtime();
258+
259+
strcpy(file,"output.nbody_");
260+
sprintf(num,"%3.3f",nz[jj]);
261+
strcat(file,num);
262+
263+
write_multiout(file,seed,oflag,rra,vva,aa);
264+
printf("time for write output = %e\n",omp_get_wtime()-t);
265+
266+
aa=afin;
267+
}
268+
269+
printf("number of step = %d\n",Nflag);
270+
271+
free(rra);
272+
free(vva);
273+
free(ro);
274+
275+
printf("done. Total time taken = %dhr %dmin %dsec\n",(int)((omp_get_wtime()-T)/3600), (int)((omp_get_wtime()-T)/60)%60, (int)(omp_get_wtime()-T)%60);
276+
}

‎nbody_funcs.c

+1,305
Large diffs are not rendered by default.

‎power1.h

+8
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
float TFfit_onek(float k, float *tf_baryon, float *tf_cdm);
2+
void TFset_parameters(float omega0hh, float f_baryon, float Tcmb);
3+
float TFzerobaryon(float omega0, float hubble, float Tcmb, float k_hmpc);
4+
/*
5+
int TFmdm_set_cosm(float omega_matter, float omega_baryon, float omega_hdm,int degen_hdm, float omega_lambda, float hubble, float redshift, float spectral);
6+
float TFmdm_onek_mpc(float kk);
7+
float TFmdm_onek_hmpc(float kk);
8+
*/

‎powerspec.c

+62
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
#include<math.h>
2+
#include<stdio.h>
3+
#include <stdlib.h>
4+
#include "power1.h"
5+
6+
/*
7+
PK(k) is the power spectrum divided by (2 \pi^2) : Pk = A k^n T(k) /(2 pi^2)
8+
kk is wavenumber (in Mpc^{-1}) and power spectrum is in units (Mpc^3).
9+
Notice there are no 'h's in the definitions.
10+
*/
11+
extern float norm,vhh;
12+
13+
float Pk(float kk)
14+
{
15+
float y,ffb,ffcdm;
16+
17+
/* Modified on 15.02.2008, using WMAP 7 year data if required please change the value of sigma_8 in the line below. It is presently using sigma_8 = 0.809 */
18+
19+
//norm=pow(0.809,2.)/(1.424307e-06);
20+
21+
if (kk>1.e-6)
22+
{
23+
y=norm*kk*pow(TFfit_onek(kk,&ffb,&ffcdm),2.);
24+
25+
}
26+
else
27+
{
28+
y=0.;
29+
}
30+
31+
return(y);
32+
}
33+
34+
float sigma_func(float kk)
35+
{
36+
float y,R;
37+
38+
R=8./vhh;
39+
y=3.*(sin(R*kk)-(R*kk)*cos(R*kk))*pow(R*kk,-3);
40+
y=(kk*kk*Pk(kk)*y*y);
41+
return(y);
42+
}
43+
44+
45+
float simp(float (*func)(float),float a,float b,int N)
46+
{
47+
float x,sum1=0.0,sum2=0.0;
48+
int i;
49+
float h=(float)((b-a)/N),z;
50+
51+
52+
for(x=a+fabs(h);x<=b;x=x+(2*fabs(h)))
53+
{
54+
sum1=sum1+(4.0*(*func)(x));
55+
}
56+
for(x=a+(2*fabs(h));x<=b;x=x+(2*fabs(h)))
57+
{
58+
sum2=sum2+(2.0*(*func)(x));
59+
}
60+
z=((sum1+sum2+(*func)(a)+(*func)(b))*fabs(h)/3);
61+
return(z);
62+
}

‎tf_fit.c

+339
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,339 @@
1+
/* The following routines implement all of the fitting formulae in
2+
Eisenstein \& Hu (1997) */
3+
4+
/* There are two sets of routines here. The first set,
5+
6+
TFfit_hmpc(), TFset_parameters(), and TFfit_onek(),
7+
8+
calculate the transfer function for an arbitrary CDM+baryon universe using
9+
the fitting formula in Section 3 of the paper. The second set,
10+
11+
TFsound_horizon_fit(), TFk_peak(), TFnowiggles(), and TFzerobaryon(),
12+
13+
calculate other quantities given in Section 4 of the paper. */
14+
15+
#include <math.h>
16+
#include <stdio.h>
17+
#include<stdlib.h>
18+
19+
void TFset_parameters(float omega0hh, float f_baryon, float Tcmb);
20+
float TFfit_onek(float k, float *tf_baryon, float *tf_cdm);
21+
22+
void TFfit_hmpc(float omega0, float f_baryon, float hubble, float Tcmb,
23+
int numk, float *k, float *tf_full, float *tf_baryon, float *tf_cdm);
24+
25+
float TFsound_horizon_fit(float omega0, float f_baryon, float hubble);
26+
float TFk_peak(float omega0, float f_baryon, float hubble);
27+
float TFnowiggles(float omega0, float f_baryon, float hubble,
28+
float Tcmb, float k_hmpc);
29+
float TFzerobaryon(float omega0, float hubble, float Tcmb, float k_hmpc);
30+
31+
/* ------------------------ DRIVER ROUTINE --------------------------- */
32+
/* The following is an example of a driver routine you might use. */
33+
/* Basically, the driver routine needs to call TFset_parameters() to
34+
set all the scalar parameters, and then call TFfit_onek() for each
35+
wavenumber k you desire. */
36+
37+
/* While the routines use Mpc^-1 units internally, this driver has been
38+
written to take an array of wavenumbers in units of h Mpc^-1. On the
39+
other hand, if you want to use Mpc^-1 externally, you can do this by
40+
altering the variables you pass to the driver:
41+
omega0 -> omega0*hubble*hubble, hubble -> 1.0 */
42+
43+
/* INPUT: omega0 -- the matter density (baryons+CDM) in units of critical
44+
f_baryon -- the ratio of baryon density to matter density
45+
hubble -- the Hubble constant, in units of 100 km/s/Mpc
46+
Tcmb -- the CMB temperature in Kelvin. T<=0 uses the COBE value 2.728.
47+
numk -- the length of the following zero-offset array
48+
k[] -- the array of wavevectors k[0..numk-1] */
49+
50+
/* INPUT/OUTPUT: There are three output arrays of transfer functions.
51+
All are zero-offset and, if used, must have storage [0..numk-1] declared
52+
in the calling program. However, if you substitute the NULL pointer for
53+
one or more of the arrays, then that particular transfer function won't
54+
be outputted. The transfer functions are:
55+
56+
tf_full[] -- The full fitting formula, eq. (16), for the matter
57+
transfer function.
58+
tf_baryon[] -- The baryonic piece of the full fitting formula, eq. 21.
59+
tf_cdm[] -- The CDM piece of the full fitting formula, eq. 17. */
60+
61+
/* Again, you can set these pointers to NULL in the function call if
62+
you don't want a particular output. */
63+
64+
/* Various intermediate scalar quantities are stored in global variables,
65+
so that you might more easily access them. However, this also means that
66+
you would be better off not simply #include'ing this file in your programs,
67+
but rather compiling it separately, calling only the driver, and using
68+
extern declarations to access the intermediate quantities. */
69+
70+
/* ----------------------------- DRIVER ------------------------------- */
71+
72+
void TFfit_hmpc(float omega0, float f_baryon, float hubble, float Tcmb,
73+
int numk, float *k, float *tf_full, float *tf_baryon, float *tf_cdm)
74+
/* Remember: k[0..numk-1] is in units of h Mpc^-1. */
75+
{
76+
int j;
77+
float tf_thisk, baryon_piece, cdm_piece;
78+
79+
TFset_parameters(omega0*hubble*hubble, f_baryon, Tcmb);
80+
81+
for (j=0;j<numk;j++) {
82+
tf_thisk = TFfit_onek(k[j]*hubble, &baryon_piece, &cdm_piece);
83+
if (tf_full!=NULL) tf_full[j] = tf_thisk;
84+
if (tf_baryon!=NULL) tf_baryon[j] = baryon_piece;
85+
if (tf_cdm!=NULL) tf_cdm[j] = cdm_piece;
86+
}
87+
return;
88+
}
89+
90+
/* ------------------------ FITTING FORMULAE ROUTINES ----------------- */
91+
92+
/* There are two routines here. TFset_parameters() sets all the scalar
93+
parameters, while TFfit_onek() calculates the transfer function for a
94+
given wavenumber k. TFfit_onek() may be called many times after a single
95+
call to TFset_parameters() */
96+
97+
/* Global variables -- We've left many of the intermediate results as
98+
global variables in case you wish to access them, e.g. by declaring
99+
them as extern variables in your main program. */
100+
/* Note that all internal scales are in Mpc, without any Hubble constants! */
101+
102+
float omhh, /* Omega_matter*h^2 */
103+
obhh, /* Omega_baryon*h^2 */
104+
theta_cmb, /* Tcmb in units of 2.7 K */
105+
z_equality, /* Redshift of matter-radiation equality, really 1+z */
106+
k_equality, /* Scale of equality, in Mpc^-1 */
107+
z_drag, /* Redshift of drag epoch */
108+
R_drag, /* Photon-baryon ratio at drag epoch */
109+
R_equality, /* Photon-baryon ratio at equality epoch */
110+
sound_horizon, /* Sound horizon at drag epoch, in Mpc */
111+
k_silk, /* Silk damping scale, in Mpc^-1 */
112+
alpha_c, /* CDM suppression */
113+
beta_c, /* CDM log shift */
114+
alpha_b, /* Baryon suppression */
115+
beta_b, /* Baryon envelope shift */
116+
beta_node, /* Sound horizon shift */
117+
k_peak, /* Fit to wavenumber of first peak, in Mpc^-1 */
118+
sound_horizon_fit, /* Fit to sound horizon, in Mpc */
119+
alpha_gamma; /* Gamma suppression in approximate TF */
120+
121+
/* Convenience from Numerical Recipes in C, 2nd edition */
122+
static float sqrarg;
123+
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
124+
static float cubearg;
125+
#define CUBE(a) ((cubearg=(a)) == 0.0 ? 0.0 : cubearg*cubearg*cubearg)
126+
static float pow4arg;
127+
#define POW4(a) ((pow4arg=(a)) == 0.0 ? 0.0 : pow4arg*pow4arg*pow4arg*pow4arg)
128+
/* Yes, I know the last one isn't optimal; it doesn't appear much */
129+
130+
void TFset_parameters(float omega0hh, float f_baryon, float Tcmb)
131+
/* Set all the scalars quantities for Eisenstein & Hu 1997 fitting formula */
132+
/* Input: omega0hh -- The density of CDM and baryons, in units of critical dens,
133+
multiplied by the square of the Hubble constant, in units
134+
of 100 km/s/Mpc */
135+
/* f_baryon -- The fraction of baryons to CDM */
136+
/* Tcmb -- The temperature of the CMB in Kelvin. Tcmb<=0 forces use
137+
of the COBE value of 2.728 K. */
138+
/* Output: Nothing, but set many global variables used in TFfit_onek().
139+
You can access them yourself, if you want. */
140+
/* Note: Units are always Mpc, never h^-1 Mpc. */
141+
{
142+
float z_drag_b1, z_drag_b2;
143+
float alpha_c_a1, alpha_c_a2, beta_c_b1, beta_c_b2, alpha_b_G, y;
144+
145+
if (f_baryon<=0.0 || omega0hh<=0.0) {
146+
fprintf(stderr, "TFset_parameters(): Illegal input.\n");
147+
exit(1);
148+
}
149+
omhh = omega0hh;
150+
obhh = omhh*f_baryon;
151+
if (Tcmb<=0.0) Tcmb=2.728; /* COBE FIRAS */
152+
theta_cmb = Tcmb/2.7;
153+
154+
z_equality = 2.50e4*omhh/POW4(theta_cmb); /* Really 1+z */
155+
k_equality = 0.0746*omhh/SQR(theta_cmb);
156+
157+
z_drag_b1 = 0.313*pow(omhh,-0.419)*(1+0.607*pow(omhh,0.674));
158+
z_drag_b2 = 0.238*pow(omhh,0.223);
159+
z_drag = 1291*pow(omhh,0.251)/(1+0.659*pow(omhh,0.828))*
160+
(1+z_drag_b1*pow(obhh,z_drag_b2));
161+
162+
R_drag = 31.5*obhh/POW4(theta_cmb)*(1000/(1+z_drag));
163+
R_equality = 31.5*obhh/POW4(theta_cmb)*(1000/z_equality);
164+
165+
sound_horizon = 2./3./k_equality*sqrt(6./R_equality)*
166+
log((sqrt(1+R_drag)+sqrt(R_drag+R_equality))/(1+sqrt(R_equality)));
167+
168+
k_silk = 1.6*pow(obhh,0.52)*pow(omhh,0.73)*(1+pow(10.4*omhh,-0.95));
169+
170+
alpha_c_a1 = pow(46.9*omhh,0.670)*(1+pow(32.1*omhh,-0.532));
171+
alpha_c_a2 = pow(12.0*omhh,0.424)*(1+pow(45.0*omhh,-0.582));
172+
alpha_c = pow(alpha_c_a1,-f_baryon)*
173+
pow(alpha_c_a2,-CUBE(f_baryon));
174+
175+
beta_c_b1 = 0.944/(1+pow(458*omhh,-0.708));
176+
beta_c_b2 = pow(0.395*omhh, -0.0266);
177+
beta_c = 1.0/(1+beta_c_b1*(pow(1-f_baryon, beta_c_b2)-1));
178+
179+
y = z_equality/(1+z_drag);
180+
alpha_b_G = y*(-6.*sqrt(1+y)+(2.+3.*y)*log((sqrt(1+y)+1)/(sqrt(1+y)-1)));
181+
alpha_b = 2.07*k_equality*sound_horizon*pow(1+R_drag,-0.75)*alpha_b_G;
182+
183+
beta_node = 8.41*pow(omhh, 0.435);
184+
beta_b = 0.5+f_baryon+(3.-2.*f_baryon)*sqrt(pow(17.2*omhh,2.0)+1);
185+
186+
k_peak = 2.5*3.14159*(1+0.217*omhh)/sound_horizon;
187+
sound_horizon_fit = 44.5*log(9.83/omhh)/sqrt(1+10.0*pow(obhh,0.75));
188+
189+
alpha_gamma = 1-0.328*log(431.0*omhh)*f_baryon + 0.38*log(22.3*omhh)*
190+
SQR(f_baryon);
191+
192+
return;
193+
}
194+
195+
float TFfit_onek(float k, float *tf_baryon, float *tf_cdm)
196+
/* Input: k -- Wavenumber at which to calculate transfer function, in Mpc^-1.
197+
*tf_baryon, *tf_cdm -- Input value not used; replaced on output if
198+
the input was not NULL. */
199+
/* Output: Returns the value of the full transfer function fitting formula.
200+
This is the form given in Section 3 of Eisenstein & Hu (1997).
201+
*tf_baryon -- The baryonic contribution to the full fit.
202+
*tf_cdm -- The CDM contribution to the full fit. */
203+
/* Notes: Units are Mpc, not h^-1 Mpc. */
204+
{
205+
float T_c_ln_beta, T_c_ln_nobeta, T_c_C_alpha, T_c_C_noalpha;
206+
float q, xx, xx_tilde, q_eff;
207+
float T_c_f, T_c, s_tilde, T_b_T0, T_b, f_baryon, T_full;
208+
float T_0_L0, T_0_C0, T_0, gamma_eff;
209+
float T_nowiggles_L0, T_nowiggles_C0, T_nowiggles;
210+
211+
k = fabs(k); /* Just define negative k as positive */
212+
if (k==0.0) {
213+
if (tf_baryon!=NULL) *tf_baryon = 1.0;
214+
if (tf_cdm!=NULL) *tf_cdm = 1.0;
215+
return 1.0;
216+
}
217+
218+
q = k/13.41/k_equality;
219+
xx = k*sound_horizon;
220+
221+
T_c_ln_beta = log(2.718282+1.8*beta_c*q);
222+
T_c_ln_nobeta = log(2.718282+1.8*q);
223+
T_c_C_alpha = 14.2/alpha_c + 386.0/(1+69.9*pow(q,1.08));
224+
T_c_C_noalpha = 14.2 + 386.0/(1+69.9*pow(q,1.08));
225+
226+
T_c_f = 1.0/(1.0+POW4(xx/5.4));
227+
T_c = T_c_f*T_c_ln_beta/(T_c_ln_beta+T_c_C_noalpha*SQR(q)) +
228+
(1-T_c_f)*T_c_ln_beta/(T_c_ln_beta+T_c_C_alpha*SQR(q));
229+
230+
s_tilde = sound_horizon*pow(1+CUBE(beta_node/xx),-1./3.);
231+
xx_tilde = k*s_tilde;
232+
233+
T_b_T0 = T_c_ln_nobeta/(T_c_ln_nobeta+T_c_C_noalpha*SQR(q));
234+
T_b = sin(xx_tilde)/(xx_tilde)*(T_b_T0/(1+SQR(xx/5.2))+
235+
alpha_b/(1+CUBE(beta_b/xx))*exp(-pow(k/k_silk,1.4)));
236+
237+
f_baryon = obhh/omhh;
238+
T_full = f_baryon*T_b + (1-f_baryon)*T_c;
239+
240+
/* Now to store these transfer functions */
241+
if (tf_baryon!=NULL) *tf_baryon = T_b;
242+
if (tf_cdm!=NULL) *tf_cdm = T_c;
243+
return T_full;
244+
}
245+
246+
/* ======================= Approximate forms =========================== */
247+
248+
float TFsound_horizon_fit(float omega0, float f_baryon, float hubble)
249+
/* Input: omega0 -- CDM density, in units of critical density
250+
f_baryon -- Baryon fraction, the ratio of baryon to CDM density.
251+
hubble -- Hubble constant, in units of 100 km/s/Mpc
252+
/* Output: The approximate value of the sound horizon, in h^-1 Mpc. */
253+
/* Note: If you prefer to have the answer in units of Mpc, use hubble -> 1
254+
and omega0 -> omega0*hubble^2. */
255+
{
256+
float omhh, sound_horizon_fit_mpc;
257+
omhh = omega0*hubble*hubble;
258+
sound_horizon_fit_mpc =
259+
44.5*log(9.83/omhh)/sqrt(1+10.0*pow(omhh*f_baryon,0.75));
260+
return sound_horizon_fit_mpc*hubble;
261+
}
262+
263+
float TFk_peak(float omega0, float f_baryon, float hubble)
264+
/* Input: omega0 -- CDM density, in units of critical density
265+
f_baryon -- Baryon fraction, the ratio of baryon to CDM density.
266+
hubble -- Hubble constant, in units of 100 km/s/Mpc
267+
/* Output: The approximate location of the first baryonic peak, in h Mpc^-1 */
268+
/* Note: If you prefer to have the answer in units of Mpc^-1, use hubble -> 1
269+
and omega0 -> omega0*hubble^2. */
270+
{
271+
float omhh, k_peak_mpc;
272+
omhh = omega0*hubble*hubble;
273+
k_peak_mpc = 2.5*3.14159*(1+0.217*omhh)/TFsound_horizon_fit(omhh,f_baryon,1.0);
274+
return k_peak_mpc/hubble;
275+
}
276+
277+
float TFnowiggles(float omega0, float f_baryon, float hubble,
278+
float Tcmb, float k_hmpc)
279+
/* Input: omega0 -- CDM density, in units of critical density
280+
f_baryon -- Baryon fraction, the ratio of baryon to CDM density.
281+
hubble -- Hubble constant, in units of 100 km/s/Mpc
282+
Tcmb -- Temperature of the CMB in Kelvin; Tcmb<=0 forces use of
283+
COBE FIRAS value of 2.728 K
284+
k_hmpc -- Wavenumber in units of (h Mpc^-1). */
285+
/* Output: The value of an approximate transfer function that captures the
286+
non-oscillatory part of a partial baryon transfer function. In other words,
287+
the baryon oscillations are left out, but the suppression of power below
288+
the sound horizon is included. See equations (30) and (31). */
289+
/* Note: If you prefer to use wavenumbers in units of Mpc^-1, use hubble -> 1
290+
and omega0 -> omega0*hubble^2. */
291+
{
292+
float k, omhh, theta_cmb, k_equality, q, xx, alpha_gamma, gamma_eff;
293+
float q_eff, T_nowiggles_L0, T_nowiggles_C0;
294+
295+
k = k_hmpc*hubble; /* Convert to Mpc^-1 */
296+
omhh = omega0*hubble*hubble;
297+
if (Tcmb<=0.0) Tcmb=2.728; /* COBE FIRAS */
298+
theta_cmb = Tcmb/2.7;
299+
300+
k_equality = 0.0746*omhh/SQR(theta_cmb);
301+
q = k/13.41/k_equality;
302+
xx = k*TFsound_horizon_fit(omhh, f_baryon, 1.0);
303+
304+
alpha_gamma = 1-0.328*log(431.0*omhh)*f_baryon + 0.38*log(22.3*omhh)*
305+
SQR(f_baryon);
306+
gamma_eff = omhh*(alpha_gamma+(1-alpha_gamma)/(1+POW4(0.43*xx)));
307+
q_eff = q*omhh/gamma_eff;
308+
309+
T_nowiggles_L0 = log(2.0*2.718282+1.8*q_eff);
310+
T_nowiggles_C0 = 14.2 + 731.0/(1+62.5*q_eff);
311+
return T_nowiggles_L0/(T_nowiggles_L0+T_nowiggles_C0*SQR(q_eff));
312+
}
313+
314+
/* ======================= Zero Baryon Formula =========================== */
315+
316+
float TFzerobaryon(float omega0, float hubble, float Tcmb, float k_hmpc)
317+
/* Input: omega0 -- CDM density, in units of critical density
318+
hubble -- Hubble constant, in units of 100 km/s/Mpc
319+
Tcmb -- Temperature of the CMB in Kelvin; Tcmb<=0 forces use of
320+
COBE FIRAS value of 2.728 K
321+
k_hmpc -- Wavenumber in units of (h Mpc^-1). */
322+
/* Output: The value of the transfer function for a zero-baryon universe. */
323+
/* Note: If you prefer to use wavenumbers in units of Mpc^-1, use hubble -> 1
324+
and omega0 -> omega0*hubble^2. */
325+
{
326+
float k, omhh, theta_cmb, k_equality, q, T_0_L0, T_0_C0;
327+
328+
k = k_hmpc*hubble; /* Convert to Mpc^-1 */
329+
omhh = omega0*hubble*hubble;
330+
if (Tcmb<=0.0) Tcmb=2.728; /* COBE FIRAS */
331+
theta_cmb = Tcmb/2.7;
332+
333+
k_equality = 0.0746*omhh/SQR(theta_cmb);
334+
q = k/13.41/k_equality;
335+
336+
T_0_L0 = log(2.0*2.718282+1.8*q);
337+
T_0_C0 = 14.2 + 731.0/(1+62.5*q);
338+
return T_0_L0/(T_0_L0+T_0_C0*q*q);
339+
}

0 commit comments

Comments
 (0)
Please sign in to comment.