-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathstacking_big_box.py
362 lines (310 loc) · 11.3 KB
/
stacking_big_box.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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
import numpy as np
import matplotlib.pyplot as plt
import splashback as sp
from scipy.optimize import curve_fit
plt.style.use("mnras.mplstyle")
axes_labels = {
"mass": "$M_{\\rm{200m}} / \\rm{M_{\odot}}$",
"accretion": "$\Gamma$",
"energy": "$X_{\\rm{E}}$"}
def bin_profiles(data, accretion_bins, mass_bins, energy_bins,
bootstrap=False):
"""
Takes a given run object and bins the density profiles according to
3 given bin arrays.
Parameters
----------
data : obj
Simulation data.
accretion_bins : array, float
Bin edges for stacking according to the accretion rate
mass_bins : array, float
Bin edges for stacking following the cluster mass
energy_bins : array. float
bin edges for stacking following the cluster's energy ratio
Returns
-------
None.
"""
sp.stack_and_find_3D(data, "accretion", accretion_bins, bootstrap=bootstrap)
sp.stack_and_find_3D(data, "mass", mass_bins, bootstrap=bootstrap)
sp.stack_and_find_3D(data, "energy", energy_bins, bootstrap=bootstrap)
def second_caustic(split):
"""
Determines which minima is the splashback and which is the second
caustic and saves the values accordingly
Parameters
----------
split : string
Name of the criteria previously used to stack the profiles.
Returns
-------
None.
"""
second_all = getattr(flm, "2_DM_"+split)
R_sp = getattr(flm, "R_DM_"+split)
second_mask = np.where(np.isfinite(second_all))[0]
for i in range(len(second_mask)):
index = second_mask[i]
if R_sp[index] < second_all[index]:
larger = second_all[index]
smaller = R_sp[index]
R_sp[index] = larger
second_all[index] = smaller
setattr(flm, "R_DM_"+split, R_sp)
setattr(flm, "2_DM_"+split, second_all)
def plot_profiles_3D_all(flm, accretion_bins, mass_bins, energy_bins):
"""
Plots stacked profiles for one run to compare the effect of using different
stacking criteria to define the bins used to stack.
Parameters
----------
flm : obj
Run data to plot
accretion_bins : array (N_bins-1)
Edge values of accretion rate bins used for labelling.
mass_bins : array (N_bins-1)
Edge values of mass bins used for labelling.
energy_bins : array (N_bins-1)
Edge values of energy ratio bins used for labelling.
Returns
-------
None.
"""
bins = np.vstack((accretion_bins, mass_bins, energy_bins))
bin_type = np.array(["accretion", "mass", "energy"])
labels = np.array(["$\Gamma$", "$\log M_{\\rm{200m}}$", "$X_{\\rm{E}}$"])
N_bins = len(accretion_bins) - 1
ylim = (-4.1,1.5)
fig, ax = plt.subplots(nrows=3, ncols=2,
figsize=(3.3,4.5),
sharey=True,
gridspec_kw={'hspace' : 0, 'wspace' : 0})
cm1 = plt.cm.autumn(np.linspace(0,0.95,N_bins))
cm2 = plt.cm.winter(np.linspace(0,1,N_bins))
cm3 = plt.cm.copper(np.linspace(0,1,N_bins))
lw = 0.8
for i in range(N_bins):
for j in range(3):
if i == 0:
label = labels[j] + "$<$" + str(np.round(bins[j,i+1],2))
elif i == N_bins-1:
label = labels[j] + "$>$" + str(np.round(bins[j,i],2))
else:
label = str(np.round(bins[j,i],2)) \
+ "$<$" + labels[j] + "$<$" \
+ str(np.round(bins[j,i+1],2))
if j == 0:
cm = cm1
elif j==1:
cm = cm2
else:
cm = cm3
if i <= 2:
label_DM = label
label_gas = ""
else:
label_DM = ""
label_gas = label
ax[j,0].semilogx(flm.rad_mid, getattr(flm, bin_type[j] + "_log_DM")[i,:],
color=cm[i], linewidth=lw,
label=label_DM)
ax[j,1].semilogx(flm.rad_mid, getattr(flm, bin_type[j] + "_log_gas")[i,:],
color=cm[i], linewidth=lw,
label=label_gas)
ax[0,0].set_ylim(ylim)
ax[0,0].legend(loc="upper right")
ax[1,0].legend(loc="upper right")
ax[2,0].legend(loc="upper right")
ax[0,1].legend(loc="upper left")
ax[1,1].legend(loc="upper left")
ax[2,1].legend(loc="upper left")
ax[1,0].set_ylabel(r"$d \log \rho / d \log r$")
plt.subplots_adjust(bottom=0.1)
ax[0,0].set_xticklabels([])
ax[1,0].set_xticklabels([])
fig.text(0.48, 0.05, "$r/R_{\\rm{200m}}$", transform=fig.transFigure)
ax[0,0].text(0.05, 0.07, "$\\rho=\\rho_{\\rm{DM}}$", transform=ax[0,0].transAxes)
ax[0,1].text(0.05, 0.07, "$\\rho=\\rho_{\\rm{gas}}$", transform=ax[0,1].transAxes)
filename = "splashback_data/flamingo/plots/HF_compare_bins_all.png"
plt.savefig(filename, dpi=300)
plt.show()
def fit_model(accretion, a, b, c, d):
"""
Model to get splashback radius from accretion rate.
Parameters
----------
accretion : float
Accretion rate values to use to estimate the splashback radius.
a : float
Fitting parameter.
b : float
Fitting parameter.
c : float
Fitting parameter.
d : float
Fitting parameter.
Returns
-------
Rsp : float
Predicted splashabck radius values.
"""
omega_m = 0.307
Rsp = a * (1+b*omega_m) * (1+ c*np.exp(-1*accretion/d))
return Rsp
def Rsp_model(accretion, model="More", ydata=[]):
"""
Sorts which model paramters to use/fits the model to get the parameters
and then gets an estimate of the splashback radius.
Parameters
----------
accretion : float
Accretion rate values to use to estimate the splashback radius
model : str, optional
Name of the model parameters to use in the model that predictes the
splashback radius. The default is "More".
ydata : array, optional
If model=="mine", this is necessary to pass onto the fitting function.
The default is [].
Returns
-------
Rsp : float
Estimated splashback radius values.
"""
more = [0.54, 0.53, 1.36, 3.04]
oneil = [0.8, 0.26, 1.14, 1.25]
if model == "O'Neil":
params = oneil
elif model == "mine":
popt, pcov = curve_fit(fit_model, accretion, ydata,
p0=oneil)
params = popt
else:
params = more
Rsp = fit_model(accretion, *params)
return Rsp
def plot_param_correlations(split, ax, plot="R"):
"""
Plots one panel of a plot with relationship between a named stacking
criteria and the obtained splashback radius.
Parameters
----------
split : str
Name of the stacking criteria used to create the profile to plot.
ax : obj
Pyplot axes object.
plot : str, optional
Which value to plot on the y-axis. Either "R" for the splashback radius
or "depth" to plot the depth of the feature. The default is "R".
Returns
-------
None.
"""
plot_DM = getattr(flm, plot+"_DM_"+split)
plot_gas = getattr(flm, plot+"_gas_"+split)
plot_pressure = getattr(flm, plot+"_P_"+split)
mids = getattr(flm, split+"_mid")
errors_DM = getattr(flm, "error_"+plot+"_DM_"+split)
errors_gas = getattr(flm, "error_"+plot+"_gas_"+split)
errors_pressure = getattr(flm, "error_"+plot+"_P_"+split)
label_DM = "Dark matter density"
label_gas = "Gas density"
label_P = "Gas pressure"
fmt = ""
# plt.figure()
if split == "mass":
ax.set_xscale('log')
elif split == "accretion" and plot == "R":
Rsp_more = Rsp_model(mids)
Rsp_oneil = Rsp_model(mids, model="O'Neil")
Rsp_mine = Rsp_model(mids, model="mine", ydata=plot_DM)
ax.plot(mids, Rsp_more,
color="darkmagenta", label="More 2015", linestyle="--")
ax.plot(mids, Rsp_oneil,
color="darkmagenta", label="O'Neil 2021", linestyle=":")
ax.plot(mids, Rsp_mine,
color="darkmagenta", label="Our model",)
label_DM = "Our data"
label_gas = ""
label_P = ""
fmt = "."
ax.errorbar(mids, plot_DM, yerr=3*errors_DM, fmt=fmt,
color="darkmagenta", label=label_DM, capsize=2)
ax.errorbar(mids, plot_gas, yerr=3*errors_gas,
color="gold", label=label_gas, capsize=2)
ax.errorbar(mids, plot_pressure, yerr=3*errors_pressure,
color="c", label=label_P, capsize=2)
ax.set_xlabel(axes_labels[split])
def stack_for_profiles(flm):
"""
Defines bins and then plots 3D density profiles.
Parameters
----------
flm : obj
Simulation data.
Returns
-------
None.
"""
N_bins = 4
mass_bins = np.linspace(14, 15, N_bins+1)
mass_bins = np.append(mass_bins, 16)
accretion_bins = np.linspace(0, 4, N_bins+1)
accretion_bins = np.append(accretion_bins, 20)
energy_bins = np.linspace(0.05, 0.3, N_bins+1)
energy_bins = np.append(energy_bins, 1)
bin_profiles(flm, accretion_bins, mass_bins, energy_bins)
plot_profiles_3D_all(flm, accretion_bins, mass_bins, energy_bins)
def stack_for_params():
"""
Defines bins and then plots relationship between the different stacking
criteria used to stack the profiles and either the radius or depth of the
splashback feature.
Returns
-------
None.
"""
N_bins = 10
mass_bins = np.linspace(14, 15.6, N_bins+1)
accretion_bins = np.linspace(0, 4.2, N_bins+1)
energy_bins = np.linspace(0.05, 0.35, N_bins+1)
bin_profiles(flm, accretion_bins, mass_bins, energy_bins,
bootstrap=True)
second_caustic("accretion")
flm.mass_mid = 10**((mass_bins[:-1] + mass_bins[1:])/2)
flm.accretion_mid = (accretion_bins[:-1] + accretion_bins[1:])/2
flm.energy_mid = (energy_bins[:-1] + energy_bins[1:])/2
fig, axes = plt.subplots(nrows=1, ncols=3,
sharey=True,
figsize=(7,2),
gridspec_kw={'hspace' : 0.1, 'wspace' : 0})
plot_param_correlations("mass", axes[1])
plot_param_correlations("accretion", axes[0])
plot_param_correlations("energy", axes[2])
axes[0].set_ylabel("$R_{\\rm{SP}} / R_{\\rm{200m}}$")
axes[0].legend()
axes[2].legend()
plt.subplots_adjust(bottom=0.18)
filename = "splashback_data/flamingo/plots/parameter_dependence_R.png"
plt.savefig(filename, dpi=300)
plt.show()
fig, axes = plt.subplots(nrows=1, ncols=3,
sharey=True,
figsize=(7,2),
gridspec_kw={'hspace' : 0.1, 'wspace' : 0})
plot_param_correlations("mass", axes[1], plot="depth")
plot_param_correlations("accretion", axes[0], plot="depth")
plot_param_correlations("energy", axes[2], plot="depth")
axes[0].set_ylabel("$\gamma_{\\rm{SP}}$")
axes[0].legend()
plt.subplots_adjust(bottom=0.18)
filename = "splashback_data/flamingo/plots/parameter_dependence_gamma.png"
plt.savefig(filename, dpi=300)
plt.show()
if __name__ == "__main__":
box = "L2800N5040"
flm = sp.flamingo(box, "HF")
flm.read_properties()
flm.read_pressure()
# stack_for_profiles(flm)
stack_for_params()