-
Notifications
You must be signed in to change notification settings - Fork 1
/
mead_BAHAMAS.py
133 lines (106 loc) · 3.85 KB
/
mead_BAHAMAS.py
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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
# Standard imports
import os
# Third-party imports
import numpy as np
power_dir = '/Users/Mead/Physics/BAHAMAS/power'
def power_file_name(mesh, model, snap, field_pair):
# Measured BAHAMAS power spectra file names
field1 = field_pair[0]
field2 = field_pair[1]
file_name1 = power_dir+'/M' + \
str(mesh)+'/'+model+'_L400N1024_WMAP9_snap' + \
str(snap)+'_'+field1+'_'+field2+'_power.dat'
file_name2 = power_dir+'/M' + \
str(mesh)+'/'+model+'_L400N1024_WMAP9_snap' + \
str(snap)+'_'+field2+'_'+field1+'_power.dat'
if os.path.isfile(file_name1):
return file_name1
else:
return file_name2
def error_file_name(mesh, snap, field_pair):
# Measured BAHAMAS errors between different realisations of the AGN_TUNED_nu0 model
field1 = field_pair[0]
field2 = field_pair[1]
return power_dir+'/M'+str(mesh)+'/L400N1024_WMAP9_snap'+str(snap)+'_'+field1+'_'+field2+'_error.dat'
def HMcode_file_name(mesh, snap):
# Corresponding HMcode power file
return power_dir+'/M'+str(mesh)+'/HMcode/snap'+str(snap)+'_HMcode.dat'
def z_to_snap(z):
# Get the snapshot number corresponding to different BAHAMAS redshifts
if z == 0.000:
snap = 32
elif z == 0.125:
snap = 31
elif z == 0.250:
snap = 30
elif z == 0.375:
snap = 29
elif z == 0.500:
snap = 28
elif z == 0.750:
snap = 27
elif z == 1.000:
snap = 26
elif z == 1.250:
snap = 25
elif z == 1.500:
snap = 24
elif z == 1.750:
snap = 23
elif z == 2.000:
snap = 22
else:
print('Redshift = ', z)
raise ValueError('Snapshot not stored corresponding to this z')
return snap
def get_measured_power(mesh, model, z, field_pair, realisation_errors=False, correct_shot_noise=False):
# Read a BAHAMAS power/error file and output k, power, error
column_k = 0
column_power = 1
column_shot = 2
column_modes = 3
column_error = 4
snap = z_to_snap(z)
infile = power_file_name(mesh, model, snap, field_pair)
data = np.loadtxt(infile)
k = data[:, column_k]
power = data[:, column_power]
modes = data[:, column_modes]
shot = data[:, column_shot]
if realisation_errors:
infile = error_file_name(mesh, snap, field_pair)
data = np.loadtxt(infile)
error = data[:, column_error]
# Subtract shot noise
if correct_shot_noise:
power = power - shot
return k, power, shot, modes, error
def get_measured_response(mesh, model, z, field_pair):
# Get the power from the hydro model
k, power, _, _, _ = get_measured_power(mesh, model, z, field_pair,
realisation_errors=False,
correct_shot_noise=True,
)
# Get the power from the DMONLY model
dmonly_name = 'DMONLY_2fluid_nu0'
_, dmonly, _, _, _ = get_measured_power(mesh, dmonly_name, z, field_pair=('all', 'all'),
realisation_errors=False,
correct_shot_noise=True,
)
# Make the response
response = power/dmonly
return k, response
def get_HMcode_power(mesh, z):
snap = z_to_snap(z)
infile = HMcode_file_name(mesh, snap)
data = np.loadtxt(infile)
k = data[:, 0]
power = data[:, 1]
return k, power
def get_HMcode_corrected_power(mesh, model, z, field_pair, realisation_errors=False, correct_shot_noise=False):
_, _, shot, modes, error = get_measured_power(
mesh, model, z, field_pair, realisation_errors, correct_shot_noise)
_, response = get_measured_response(mesh, model, z, field_pair)
k, power = get_HMcode_power(mesh, z)
power = power*response
return k, power, shot, modes, error