1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98 | #include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include "power.h"
#include "msg.h"
PowerSpectrum::PowerSpectrum(const char filename[])
{
// Read k P from filename
FILE* fp= fopen(filename, "r");
if(fp == 0) {
msg_printf(msg_error, "Error: Unable to open input power spectrum file: %s\n",filename);
throw PowerFileError();
}
int nalloc= 1000;
double* buf= (double*) malloc(sizeof(double)*2*nalloc);
char line[128];
int nlines= 0;
double kk, PP;
// Read lines and push to buf as k1,P1,k2,P2, ...
// Doubles the length of buf when the length of the array is not enough
while(fgets(line, 127, fp)) {
if(nlines == nalloc) {
msg_printf(msg_debug, "Reallocating power spectrum table %d -> %d\n",
nalloc, 2*nalloc);
nalloc *= 2;
buf= (double*) realloc(buf, sizeof(double)*2*nalloc); assert(buf);
}
if(line[0] == '#')
continue;
else if(sscanf(line, "%lg %lg", &kk, &PP) == 2) {
buf[2*nlines ]= kk;
buf[2*nlines + 1]= PP;
nlines++;
}
else {
msg_printf(msg_warn,
"Warning: Unable to understand a line in the power spectrum file;"
"following data are ignored: %s", line);
break;
}
}
int ret= fclose(fp); assert(ret == 0);
msg_printf(msg_verbose, "Read %d pairs of k P(k) from %s\n", nlines, filename);
k = (double*) malloc(2*nlines*sizeof(double)); assert(this->k);
P = k + nlines;
for(int j=0; j<nlines; j++) {
this->k[j] = buf[2*j];
this->P[j] = buf[2*j + 1];
}
free(buf);
this->n= nlines;
}
PowerSpectrum::~PowerSpectrum()
{
free(k);
}
double PowerSpectrum::compute_sigma(const double R) const
{
// Computes sigma (rms amplituede) smoothed on scale R
// R: smoothing length [/h Mpc] (8 for sigma_8)
// 1/(2*pi^2) \int P(k) W(k*R) k^2 dk
const double fac= 1.0/(2.0*M_PI*M_PI);
double k0= k[0];
double f0= 0.0;
double sigma2= 0.0;
for(int i=0; i<n; i++) {
double kk= k[i];
double x= kk*R;
double w= 3.0*(sin(x)-x*cos(x))/(x*x*x);
double f1= P[i]*kk*kk*w*w;
sigma2 += 0.5*(f0 + f1)*(kk - k0);
k0= kk;
f0= f1;
}
return sqrt(fac*sigma2);
}
|