Mercurial > repos > astroteam > sgwb_astro_tool

comparison Phase_transition_parameters.py @ 0:c9fc89ee996e draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit d5be3628a0bdad9747451669fa195f1d2acaf3f2
author astroteam
date Wed, 17 Apr 2024 10:29:23 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:c9fc89ee996e
1 #!/usr/bin/env python
2 # coding: utf-8
3
4 # flake8: noqa
5
6 import json
7 import os
8 import shutil
9
10 from oda_api.json import CustomJSONEncoder
11
12 get_ipython().run_cell_magic( # noqa: F821
13 "bash",
14 "",
15 "wget https://gitlab.renkulab.io/astronomy/mmoda/sgwb/-/raw/master/EPTA.csv\nwget https://gitlab.renkulab.io/astronomy/mmoda/sgwb/-/raw/master/NANOGrav23.csv\nwget https://gitlab.renkulab.io/astronomy/mmoda/sgwb/-/raw/master/PPTA.csv\nwget https://gitlab.renkulab.io/astronomy/mmoda/sgwb/-/raw/master/RoperPol22.csv\nwget https://gitlab.renkulab.io/astronomy/mmoda/sgwb/-/raw/master/MAGIC.csv\nwget https://gitlab.renkulab.io/astronomy/mmoda/sgwb/-/raw/master/Ellis.csv\nwget https://gitlab.renkulab.io/astronomy/mmoda/sgwb/-/raw/master/NANO_alpha_T.csv\nwget https://gitlab.renkulab.io/astronomy/mmoda/sgwb/-/raw/master/NANO_alpha_T1.csv\nwget https://gitlab.renkulab.io/astronomy/mmoda/sgwb/-/raw/master/NANO_bubble.csv\nwget https://gitlab.renkulab.io/astronomy/mmoda/sgwb/-/raw/master/NANO_sound.csv\n",
16 )
17
18 import matplotlib.pyplot as plt
19
20 # Loading necessary packages
21 import numpy as np
22
23 get_ipython().run_line_magic("matplotlib", "inline") # noqa: F821
24 import astropy.units as u
25 from astropy.constants import c
26 from numpy import exp, log, pi, sqrt
27 from oda_api.api import ProgressReporter
28 from oda_api.data_products import PictureProduct
29
30 # Parameters of the phase transition
31
32 # fraction of turbulent energy that goes to gw N.B. arXiv:1004.4187 claims that epsilon_turb=0.05, but checks below show that it is rather 0.01
33 epsilon_turb = 1.0 # http://odahub.io/ontology#Float
34
35 _galaxy_wd = os.getcwd()
36
37 with open("inputs.json", "r") as fd:
38 inp_dic = json.load(fd)
39 if "_data_product" in inp_dic.keys():
40 inp_pdic = inp_dic["_data_product"]
41 else:
42 inp_pdic = inp_dic
43
44 for vn, vv in inp_pdic.items():
45 if vn != "_selector":
46 globals()[vn] = type(globals()[vn])(vv)
47
48 Tmax_yrs = 10.0 # http://odahub.io/ontology#Float
49 Tmin_yrs = 2.0 # http://odahub.io/ontology#Float
50 # terminal velocity of bubbles
51 v_w = 0.999 # http://odahub.io/ontology#Float
52 h = 0.7 # http://odahub.io/ontology#Float
53 # Numbers of relativistic degrees of freedom
54 g_star = 20 # http://odahub.io/ontology#Integer
55
56 c_s = 3 ** (-0.5) # speed of sound
57 F0_GW = 1.64e-5 / h**2 * (100 / g_star) ** (1 / 3.0)
58 pr = ProgressReporter()
59 pr.report_progress(stage="Progress", progress=1.0)
60
61 # Eq.20 of arXiv:1512.06239, corrected according to Appendix A of arXiv:1004.4187
62
63 def ka_v(alpha_v, v):
64 zeta = ((2.0 / 3.0 * alpha_v + alpha_v**2) ** 0.5 + (1 / 3.0) ** 0.5) / (
65 1 + alpha_v
66 )
67 kappa_a = (
68 6.9
69 * v ** (6.0 / 5.0)
70 * alpha_v
71 / (1.36 - 0.037 * alpha_v**0.5 + alpha_v)
72 )
73 kappa_b = alpha_v ** (2.0 / 5.0) / (
74 0.017 + (0.997 + alpha_v) ** (2.0 / 5.0)
75 )
76 kappa_c = alpha_v**0.5 / (0.135 + (0.98 + alpha_v) ** 0.5)
77 kappa_d = alpha_v / (0.73 + 0.083 * alpha_v**0.6 + alpha_v)
78 deltak = -0.9 * log(alpha_v**0.5 / (1 + alpha_v**0.5))
79 if v < c_s:
80 return (
81 c_s ** (11.0 / 5.0)
82 * kappa_a
83 * kappa_b
84 / (
85 (c_s ** (11.0 / 5.0) - v ** (11.0 / 5.0)) * kappa_b
86 + v * c_s ** (6.0 / 5.0) * kappa_a
87 )
88 )
89 elif v > zeta:
90 return (
91 (zeta - 1) ** 3
92 * (zeta / v) ** (5.0 / 2)
93 * kappa_c
94 * kappa_d
95 / (
96 ((zeta - 1) ** 3 - (v - 1) ** 3)
97 * zeta ** (5.0 / 2.0)
98 * kappa_c
99 + (v - 1) ** 3 * kappa_d
100 )
101 )
102 else:
103 return (
104 kappa_b
105 + (v - c_s) * deltak
106 + (v - c_s) ** 3
107 / (zeta - c_s) ** 3
108 * (kappa_c - kappa_b - (zeta - c_s) * deltak)
109 )
110
111 T_star = 0.2
112
113 # Comoving Hubble rate at the phase transition Eq. 11 of arXiv:1512.06239
114 def H_star(T_star):
115 return (
116 16.5e-6
117 * (T_star / (100.0 * u.GeV))
118 * (g_star / 100) ** (1.0 / 6.0)
119 * u.Hz
120 )
121
122 Hstar = H_star(T_star * u.GeV)
123 logHstar = np.log10(Hstar / u.Hz)
124
125 fmin = logHstar.value - 5
126 fmax = logHstar.value + 5
127 ff = np.logspace(fmin, fmax, 101)
128
129 # HL model formula
130 def GW_sound(f, T_star, alpha, beta_H, v_w):
131 Hstar = H_star(T_star)
132 kappa_v = ka_v(alpha, v_w)
133 K = kappa_v * alpha / (1 + alpha)
134 lambda_star = (8 * pi) ** (1 / 3) * max([v_w, c_s]) / (beta_H * Hstar) * c
135 Delta_w = sqrt((v_w - c_s) ** 2) / v_w
136 s2 = Delta_w * lambda_star * f / c
137 # print((c/lambda_star).cgs,(c/(Delta_w*lambda_star)).cgs)
138 s1 = lambda_star * f / c
139 M = (
140 16
141 * (1 + Delta_w ** (-3)) ** (2 / 3.0)
142 * (Delta_w * s1) ** 3
143 / ((1 + s1**3) ** (2.0 / 3.0) * (3 + s2**2) ** 2)
144 )
145 factor = (K * lambda_star * Hstar) ** 2 / (
146 sqrt(K) + lambda_star * Hstar / c
147 )
148 mu = 4.78 - 6.27 * Delta_w + 3.34 * Delta_w**2
149 ff = f / u.Hz
150 dlnf = (ff[1] - ff[0]) / ff[0]
151 mu = sum(M) * dlnf
152 B = 1e-2 / mu
153 return (3 * B * factor / c**2 * F0_GW * M).cgs
154
155 # Eq 1 of the new Overleaf, from Alberto
156
157 def GW_turb_Andrii(f, T_star, alpha, beta_H, v_w, epsilon_turb):
158 Hstar = H_star(T_star)
159
160 # Eq. 1 of 2307.10744
161 lambda_star = (
162 (8 * pi) ** (1 / 3) * max([v_w, c_s]) / (beta_H * Hstar)
163 ) * c # characteristic light-crossing distance scale
164
165 kappa_v = ka_v(alpha, v_w)
166 K = kappa_v * alpha / (1 + alpha)
167
168 # Eq. 10 of 2307.10744
169 Omega_star = epsilon_turb * K
170
171 # Eq. 13 of 2307.10744 and text after it
172 u_star = sqrt(0.75 * Omega_star)
173 dt_fin = 2 * lambda_star / u_star / c
174 s3 = dt_fin * f
175 s1 = f * lambda_star / c
176 # print('u_star=',u_star,'lambda_star=',lambda_star.cgs)
177
178 # Eq. 15 of 2307.10744
179 T_GW = np.log(1 + Hstar * dt_fin / (2 * pi)) * (s3 < 1) + np.log(
180 1 + lambda_star * Hstar / c / (2 * pi * s1)
181 ) * (s3 >= 1)
182 T_GW = T_GW**2
183
184 # Eq. 17 of 2307.10744
185 alpha_pi = 2.15
186 s_pi = 2.2
187 P_pi = (1 + (s1 / s_pi) ** alpha_pi) ** (-11 / (3 * alpha_pi))
188
189 # Eq 18,19,20 of 2307.10744
190 A = 2e-3 * 1.4 * 0.6
191
192 # Eq. 14 of 2307.10744
193 Sturb = (
194 4
195 * pi**2
196 * s1**3
197 * T_GW
198 / (lambda_star * Hstar / c) ** 2
199 * P_pi
200 / 1.4
201 / 0.6
202 )
203
204 # Eq 9 of 2307.10744
205 res = (
206 3 * A * Omega_star**2 * (lambda_star * Hstar / c) ** 2 * F0_GW * Sturb
207 )
208
209 Omega_B = Omega_star / 2.0
210 Omega_gamma = 2 / g_star
211 B = 3e-6 * (Omega_B / Omega_gamma) ** 0.5
212 return res, lambda_star, B
213
214 H0 = 70 * (u.km / u.s) / u.Mpc
215
216 d = np.genfromtxt("NANOGrav23.csv")
217 gammas_nano = d[:, 0]
218 As_nano = 10 ** d[:, 1]
219 ps_nano = 5 - gammas_nano
220
221 d = np.genfromtxt("EPTA.csv")
222 gammas_epta = d[:, 0]
223 As_epta = 10 ** d[:, 1]
224 ps_epta = 5 - gammas_epta
225
226 d = np.genfromtxt("PPTA.csv")
227 gammas_ppta = d[:, 0]
228 As_ppta = 10 ** d[:, 1]
229 ps_ppta = 5 - gammas_ppta
230
231 fref = (1 / u.yr).cgs.value
232 lgfmin = np.log10(fref / Tmax_yrs)
233 lgfmax = np.log10(fref / Tmin_yrs)
234 fff = np.logspace(lgfmin, lgfmax, 10) * u.Hz
235 min_nano = np.ones(len(fff))
236 max_nano = np.zeros(len(fff))
237 for i in range(len(As_nano)):
238 spec = (
239 2
240 * pi**2
241 / 3
242 / H0**2
243 * fff**2
244 * As_nano[i] ** 2
245 * (fff / fref) ** (3 - gammas_nano[i])
246 ).cgs.value
247 min_nano = np.minimum(spec, min_nano)
248 max_nano = np.maximum(spec, max_nano)
249 # plt.plot(ff,spec)
250 plt.fill_between(
251 fff.value, min_nano, max_nano, color="red", alpha=0.5, label="NANOGrav"
252 )
253 min_epta = np.ones(len(fff))
254 max_epta = np.zeros(len(fff))
255 for i in range(len(As_epta)):
256 spec = (
257 2
258 * pi**2
259 / 3
260 / H0**2
261 * fff**2
262 * As_epta[i] ** 2
263 * (fff / fref) ** (3 - gammas_epta[i])
264 ).cgs.value
265 min_epta = np.minimum(spec, min_epta)
266 max_epta = np.maximum(spec, max_epta)
267 # plt.plot(ff,spec)
268 plt.fill_between(
269 fff.value, min_epta, max_epta, color="blue", alpha=0.5, label="EPTA"
270 )
271 min_ppta = np.ones(len(fff))
272 max_ppta = np.zeros(len(fff))
273 for i in range(len(As_ppta)):
274 spec = (
275 2
276 * pi**2
277 / 3
278 / H0**2
279 * fff**2
280 * As_ppta[i] ** 2
281 * (fff / fref) ** (3 - gammas_ppta[i])
282 ).cgs.value
283 min_ppta = np.minimum(spec, min_ppta)
284 max_ppta = np.maximum(spec, max_ppta)
285 # plt.plot(ff,spec)
286 plt.fill_between(
287 fff.value, min_ppta, max_ppta, color="green", alpha=0.5, label="PPTA"
288 )
289
290 # PBH constraints
291 def PBH(alpha):
292 return (5.2 * (1.0 - exp(-1.1 * (abs(alpha - 1.0)) ** (1 / 3))) + 1.0) * (
293 alpha > 1
294 )
295
296 Tgrid = np.logspace(-2.6, 1, 81)
297 agrid = np.logspace(-1, 2, 61)
298 bgrid = np.logspace(0, 2, 41)
299 amin_b_nano_eps1 = 100.0 * np.ones(len(bgrid))
300 amax_b_nano_eps1 = 0.0 * np.ones(len(bgrid))
301 amin_b_epta_eps1 = 100.0 * np.ones(len(bgrid))
302 amax_b_epta_eps1 = 0.0 * np.ones(len(bgrid))
303 amin_b_ppta_eps1 = 100.0 * np.ones(len(bgrid))
304 amax_b_ppta_eps1 = 0.0 * np.ones(len(bgrid))
305 amin_T_nano_eps1 = 100.0 * np.ones(len(Tgrid))
306 amax_T_nano_eps1 = 0.0 * np.ones(len(Tgrid))
307 amin_T_epta_eps1 = 100.0 * np.ones(len(Tgrid))
308 amax_T_epta_eps1 = 0.0 * np.ones(len(Tgrid))
309 amin_T_ppta_eps1 = 100.0 * np.ones(len(Tgrid))
310 amax_T_ppta_eps1 = 0.0 * np.ones(len(Tgrid))
311
312 Tmin_b_nano_eps1 = 100.0 * np.ones(len(bgrid))
313 Tmax_b_nano_eps1 = 0.0 * np.ones(len(bgrid))
314 Tmin_b_epta_eps1 = 100.0 * np.ones(len(bgrid))
315 Tmax_b_epta_eps1 = 0.0 * np.ones(len(bgrid))
316 Tmin_b_ppta_eps1 = 100.0 * np.ones(len(bgrid))
317 Tmax_b_ppta_eps1 = 0.0 * np.ones(len(bgrid))
318 B_nano_eps1 = []
319 lam_nano_eps1 = []
320 B_ppta_eps1 = []
321 lam_ppta_eps1 = []
322 B_epta_eps1 = []
323 lam_epta_eps1 = []
324 for i, T in enumerate(Tgrid):
325 print(T, len(B_nano_eps1), len(B_epta_eps1), len(B_ppta_eps1))
326 pr.report_progress(stage="Progress", progress=1 + 94.0 * (i / len(Tgrid)))
327 for j, a in enumerate(agrid):
328 for k, b in enumerate(bgrid):
329 if b > PBH(a):
330 GW_t = GW_turb_Andrii(
331 ff * u.Hz, T * u.GeV, a, b, v_w, epsilon_turb
332 )
333 lam = GW_t[1].cgs.value
334 B = GW_t[2]
335 GW_t = GW_t[0]
336 GW_s = GW_sound(ff * u.Hz, T * u.GeV, a, b, v_w)
337 GW = GW_s + GW_t
338 GW_interp = np.interp(fff.value, ff, GW)
339 if sum((GW_interp < max_nano) * (GW_interp > min_nano)) == len(
340 fff
341 ):
342 if a < amin_b_nano_eps1[k]:
343 amin_b_nano_eps1[k] = a
344 if a > amax_b_nano_eps1[k]:
345 amax_b_nano_eps1[k] = a
346 if a < amin_T_nano_eps1[i]:
347 amin_T_nano_eps1[i] = a
348 if a > amax_T_nano_eps1[i]:
349 amax_T_nano_eps1[i] = a
350 if T < Tmin_b_nano_eps1[k]:
351 Tmin_b_nano_eps1[k] = T
352 if T > Tmax_b_nano_eps1[k]:
353 Tmax_b_nano_eps1[k] = T
354 B_nano_eps1.append(B)
355 lam_nano_eps1.append(lam)
356 if sum((GW_interp < max_epta) * (GW_interp > min_epta)) == len(
357 fff
358 ):
359 if a < amin_b_epta_eps1[k]:
360 amin_b_epta_eps1[k] = a
361 if a > amax_b_epta_eps1[k]:
362 amax_b_epta_eps1[k] = a
363 if a < amin_T_epta_eps1[i]:
364 amin_T_epta_eps1[i] = a
365 if a > amax_T_epta_eps1[i]:
366 amax_T_epta_eps1[i] = a
367 if T < Tmin_b_epta_eps1[k]:
368 Tmin_b_epta_eps1[k] = T
369 if T > Tmax_b_epta_eps1[k]:
370 Tmax_b_epta_eps1[k] = T
371 B_epta_eps1.append(B)
372 lam_epta_eps1.append(lam)
373 if sum((GW_interp < max_ppta) * (GW_interp > min_ppta)) == len(
374 fff
375 ):
376 if a < amin_b_ppta_eps1[k]:
377 amin_b_ppta_eps1[k] = a
378 if a > amax_b_ppta_eps1[k]:
379 amax_b_ppta_eps1[k] = a
380 if a < amin_T_ppta_eps1[i]:
381 amin_T_ppta_eps1[i] = a
382 if a > amax_T_ppta_eps1[i]:
383 amax_T_ppta_eps1[i] = a
384 if T < Tmin_b_ppta_eps1[k]:
385 Tmin_b_ppta_eps1[k] = T
386 if T > Tmax_b_ppta_eps1[k]:
387 Tmax_b_ppta_eps1[k] = T
388 B_ppta_eps1.append(B)
389 lam_ppta_eps1.append(lam)
390
391 lam_bins = np.logspace(-7, -5, 41)
392 B_bins = np.logspace(-7, -5, 41)
393
394 h = plt.hist2d(
395 np.array(lam_nano_eps1) / 3e24,
396 B_nano_eps1,
397 bins=[lam_bins, B_bins],
398 vmax=1,
399 )
400 plt.colorbar()
401 cnt = plt.contour(
402 sqrt(lam_bins[1:] * lam_bins[:-1]),
403 sqrt(B_bins[1:] * B_bins[:-1]),
404 np.transpose(h[0]),
405 levels=[1],
406 colors="white",
407 )
408 cont = cnt.get_paths()[0].vertices
409 B_nanos_eps1 = cont[:, 1]
410 lam_nanos_eps1 = cont[:, 0]
411
412 plt.xscale("log")
413 plt.yscale("log")
414
415 lam_bins = np.logspace(-7, -5, 41)
416 B_bins = np.logspace(-7, -5, 41)
417
418 h = plt.hist2d(
419 np.array(lam_epta_eps1) / 3e24,
420 B_epta_eps1,
421 bins=[lam_bins, B_bins],
422 vmax=1,
423 )
424 plt.colorbar()
425 cnt = plt.contour(
426 sqrt(lam_bins[1:] * lam_bins[:-1]),
427 sqrt(B_bins[1:] * B_bins[:-1]),
428 np.transpose(h[0]),
429 levels=[1],
430 colors="white",
431 )
432 cont = cnt.get_paths()[0].vertices
433 B_eptas_eps1 = cont[:, 1]
434 lam_eptas_eps1 = cont[:, 0]
435
436 plt.xscale("log")
437 plt.yscale("log")
438
439 lam_bins = np.logspace(-7, -5, 41)
440 B_bins = np.logspace(-7, -5, 41)
441
442 h = plt.hist2d(
443 np.array(lam_ppta_eps1) / 3e24,
444 B_ppta_eps1,
445 bins=[lam_bins, B_bins],
446 vmax=1,
447 )
448 plt.colorbar()
449 cnt = plt.contour(
450 sqrt(lam_bins[1:] * lam_bins[:-1]),
451 sqrt(B_bins[1:] * B_bins[:-1]),
452 np.transpose(h[0]),
453 levels=[1],
454 colors="white",
455 )
456 cont = cnt.get_paths()[0].vertices
457 B_pptas_eps1 = cont[:, 1]
458 lam_pptas_eps1 = cont[:, 0]
459
460 plt.xscale("log")
461 plt.yscale("log")
462
463 import numpy as np
464 from matplotlib import pyplot as plt
465
466 fig, ax = plt.subplots(figsize=(7, 7))
467 x = [1e-8, 1e3]
468 y1 = [3e-6, 3e-6]
469 y2 = [3e-5, 3e-5]
470 plt.fill_between(x, y1, y2, color="grey", alpha=0.5)
471 ax.fill(lam_nanos_eps1, B_nanos_eps1, color="red", alpha=0.3, label="NANOGrav")
472 ax.fill(lam_eptas_eps1, B_eptas_eps1, color="blue", alpha=0.3, label="EPTA")
473 ax.fill(lam_pptas_eps1, B_pptas_eps1, color="green", alpha=0.3, label="PPTA")
474 ax.set_xscale("log")
475 ax.set_yscale("log")
476
477 lam_med = np.median(lam_nanos_eps1)
478 B_med = np.median(B_nanos_eps1)
479 lam_evol = np.logspace(np.log10(lam_med), np.log10(lam_med) + 10, 100)
480 B_evol = B_med * (lam_evol / lam_med) ** (-5 / 4.0)
481 mask = B_evol > 10 ** (-8.5) * lam_evol
482 lam_evol = lam_evol[mask]
483 B_evol = B_evol[mask]
484 ax.annotate(
485 "",
486 xy=(lam_evol[-1], B_evol[-1]),
487 xytext=(lam_evol[0], B_evol[0]),
488 arrowprops={"arrowstyle": "->", "color": "black"},
489 )
490
491 lam_evol = np.logspace(np.log10(lam_med), np.log10(lam_med) + 10, 100)
492 B_evol = B_med * (lam_evol / lam_med) ** (-5 / 2.0)
493 mask = B_evol > 10 ** (-6.8) * lam_evol
494 lam_evol = lam_evol[mask]
495 B_evol = B_evol[mask]
496 ax.annotate(
497 "",
498 xy=(lam_evol[-1], B_evol[-1]),
499 xytext=(lam_evol[0], B_evol[0]),
500 arrowprops={"arrowstyle": "->", "linestyle": "dashed", "color": "black"},
501 )
502
503 x = np.logspace(-8, 3, 10)
504 y = 10 ** (-8.5) * x
505 ax.plot(x, y, color="grey", linewidth=4, linestyle="dashed")
506 y = 10 ** (-6.8) * x
507 ax.plot(x, y, color="grey", linewidth=4)
508
509 d = np.genfromtxt("MAGIC.csv")
510 ax.fill_between(d[:, 0], np.zeros(len(d)), d[:, 1], color="grey", alpha=0.5)
511 ax.set_xlim(1e-8, 1e3)
512 ax.set_ylim(5e-18, 2e-5)
513 ax.set_xlabel(r"$\lambda_B$, Mpc")
514 ax.set_ylabel(r"$B$, G")
515 plt.text(1, 3e-17, "MAGIC '22", color="black")
516 plt.text(
517 1e-6,
518 0.1e-14,
519 "Alfvenic larges processed eddies",
520 color="black",
521 rotation=42,
522 )
523
524 plt.text(
525 0.15e-7,
526 0.01e-13,
527 "Helicity fluctuaitons / reconnection controlled",
528 color="black",
529 rotation=41.5,
530 )
531 ax.legend(loc="lower left")
532
533 y = [1.5e-12 * 1e1, 1.5e-12, 1.5e-12]
534 x = [0.3e-3, 0.03, 1e3]
535 plt.plot(x, y, color="olive", linewidth=4, linestyle="dotted")
536 ax.annotate(
537 "",
538 xytext=(1, 1.0e-13),
539 xy=(1, 1.5e-12),
540 arrowprops={"arrowstyle": "->", "color": "olive", "linewidth": 4},
541 )
542 plt.text(0.5, 2e-12, "CTA", color="olive")
543
544 x = [0.7e-3, 1e3]
545 y = [3e-11, 3e-11]
546 plt.plot(x, y, color="olive", linewidth=4, linestyle="dotted")
547 ax.annotate(
548 "",
549 xytext=(1, 3e-10),
550 xy=(1, 3e-11),
551 arrowprops={"arrowstyle": "->", "color": "olive", "linewidth": 4},
552 )
553 plt.text(0.5, 1.1e-11, "CMB", color="olive")
554
555 plt.savefig("B_lambdaB.png", bbox_inches="tight")
556
557 fig = plt.figure()
558 mask = Tmax_b_nano_eps1 > 0
559 plt.fill_between(
560 bgrid[mask],
561 Tmin_b_nano_eps1[mask],
562 Tmax_b_nano_eps1[mask],
563 alpha=0.3,
564 color="red",
565 label="NANOGrav",
566 )
567 mask = Tmax_b_epta_eps1 > 0
568 plt.fill_between(
569 bgrid[mask],
570 Tmin_b_epta_eps1[mask],
571 Tmax_b_epta_eps1[mask],
572 alpha=0.3,
573 color="blue",
574 label="EPTA",
575 )
576 mask = Tmax_b_ppta_eps1 > 0
577 plt.fill_between(
578 bgrid[mask],
579 Tmin_b_ppta_eps1[mask],
580 Tmax_b_ppta_eps1[mask],
581 alpha=0.3,
582 color="green",
583 label="PPTA",
584 )
585 d = np.genfromtxt("Ellis.csv")
586 lgT = d[:, 0]
587 lgb = d[:, 1]
588 plt.plot(10**lgb, 10**lgT, color="black", label="2308.08546")
589 d = np.genfromtxt("NANO_bubble.csv")
590 lgT = d[:, 0]
591 lgb = -d[:, 1]
592 plt.plot(
593 10**lgb, 10**lgT, color="black", linestyle="dashed", label="2306.162196"
594 )
595 d = np.genfromtxt("NANO_sound.csv")
596 lgT = d[:, 0]
597 lgb = -d[:, 1]
598 plt.plot(10**lgb, 10**lgT, color="black", linestyle="dashed")
599 plt.axhline(0.160, color="black", linestyle="dotted")
600
601 plt.xscale("log")
602 plt.yscale("log")
603 plt.ylabel(r"$T$ [GeV]")
604 plt.xlabel(r"$\beta/H$")
605 plt.xlim(0.8, 80)
606 plt.ylim(0.001, 10)
607 plt.legend(loc="lower left")
608 plt.savefig("T_Beta.png", bbox_inches="tight")
609
610 fig = plt.figure()
611 mask = amax_b_nano_eps1 > 0
612 plt.fill_between(
613 bgrid[mask],
614 amin_b_nano_eps1[mask],
615 amax_b_nano_eps1[mask],
616 alpha=0.3,
617 color="red",
618 label="NANOGrav",
619 )
620 mask = amax_b_epta_eps1 > 0
621 plt.fill_between(
622 bgrid[mask],
623 amin_b_epta_eps1[mask],
624 amax_b_epta_eps1[mask],
625 alpha=0.3,
626 color="blue",
627 label="EPTA",
628 )
629 mask = amax_b_ppta_eps1 > 0
630 plt.fill_between(
631 bgrid[mask],
632 amin_b_ppta_eps1[mask],
633 amax_b_ppta_eps1[mask],
634 alpha=0.3,
635 color="green",
636 label="PPTA",
637 )
638 bpbh = PBH(agrid)
639 x = [2, 3]
640 y = [1, 1]
641 y1 = [2, 2]
642 plt.fill_between(x, y, y1, color="white", linewidth=0)
643 plt.fill(bpbh, agrid, color="white")
644 plt.xscale("log")
645 plt.yscale("log")
646 plt.xlim(0.8, 80)
647 plt.ylim(0.2, 80)
648 plt.xlabel(r"$\beta/H$")
649 plt.ylabel(r"$\alpha$")
650 plt.legend(loc="upper left")
651 plt.savefig("Alpha_Beta.png", bbox_inches="tight")
652
653 fig = plt.figure()
654 mask = amax_T_nano_eps1 > 0
655 plt.fill_between(
656 Tgrid[mask],
657 amin_T_nano_eps1[mask],
658 amax_T_nano_eps1[mask],
659 alpha=0.3,
660 color="red",
661 label="NANOGrav",
662 )
663 mask = amax_T_epta_eps1 > 0
664 plt.fill_between(
665 Tgrid[mask],
666 amin_T_epta_eps1[mask],
667 amax_T_epta_eps1[mask],
668 alpha=0.3,
669 color="blue",
670 label="EPTA",
671 linewidth=0,
672 )
673 mask = amax_T_ppta_eps1 > 0
674 plt.fill_between(
675 Tgrid[mask],
676 amin_T_ppta_eps1[mask],
677 amax_T_ppta_eps1[mask],
678 alpha=0.3,
679 color="green",
680 label="PPTA",
681 linewidth=0,
682 )
683 d = np.genfromtxt("NANO_alpha_T.csv")
684 lgT = d[:, 0]
685 lga = d[:, 1]
686 plt.plot(
687 10**lgT, 10**lga, color="black", label="2306.16219", linestyle="dashed"
688 )
689 d = np.genfromtxt("NANO_alpha_T1.csv")
690 lgT = d[:, 0]
691 lga = d[:, 1]
692 plt.plot(10**lgT, 10**lga, color="black", linestyle="dashed")
693 plt.axvline(0.16, color="black", linestyle="dotted")
694 plt.xscale("log")
695 plt.yscale("log")
696 plt.xlim(0.001, 10)
697 plt.ylim(0.2, 80)
698 plt.xlabel(r"$T$ [GeV]")
699 plt.ylabel(r"$\alpha$")
700 plt.legend(loc="upper left")
701 plt.savefig("Alpha_T.png", bbox_inches="tight")
702
703 bin_image1 = PictureProduct.from_file("B_lambdaB.png")
704 bin_image2 = PictureProduct.from_file("T_Beta.png")
705 bin_image3 = PictureProduct.from_file("Alpha_Beta.png")
706 bin_image4 = PictureProduct.from_file("Alpha_T.png")
707
708 B_lambdaB_png = bin_image1 # http://odahub.io/ontology#ODAPictureProduct
709 T_Beta_png = bin_image2 # http://odahub.io/ontology#ODAPictureProduct
710 Alpha_Beta_png = bin_image3 # http://odahub.io/ontology#ODAPictureProduct
711 Alpha_T_png = bin_image4 # http://odahub.io/ontology#ODAPictureProduct
712
713 # output gathering
714 _galaxy_meta_data = {}
715 _oda_outs = []
716 _oda_outs.append(
717 (
718 "out_Phase_transition_parameters_B_lambdaB_png",
719 "B_lambdaB_png_galaxy.output",
720 B_lambdaB_png,
721 )
722 )
723 _oda_outs.append(
724 (
725 "out_Phase_transition_parameters_T_Beta_png",
726 "T_Beta_png_galaxy.output",
727 T_Beta_png,
728 )
729 )
730 _oda_outs.append(
731 (
732 "out_Phase_transition_parameters_Alpha_Beta_png",
733 "Alpha_Beta_png_galaxy.output",
734 Alpha_Beta_png,
735 )
736 )
737 _oda_outs.append(
738 (
739 "out_Phase_transition_parameters_Alpha_T_png",
740 "Alpha_T_png_galaxy.output",
741 Alpha_T_png,
742 )
743 )
744
745 for _outn, _outfn, _outv in _oda_outs:
746 _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn)
747 if isinstance(_outv, str) and os.path.isfile(_outv):
748 shutil.move(_outv, _galaxy_outfile_name)
749 _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
750 elif getattr(_outv, "write_fits_file", None):
751 _outv.write_fits_file(_galaxy_outfile_name)
752 _galaxy_meta_data[_outn] = {"ext": "fits"}
753 elif getattr(_outv, "write_file", None):
754 _outv.write_file(_galaxy_outfile_name)
755 _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
756 else:
757 with open(_galaxy_outfile_name, "w") as fd:
758 json.dump(_outv, fd, cls=CustomJSONEncoder)
759 _galaxy_meta_data[_outn] = {"ext": "json"}
760
761 with open(os.path.join(_galaxy_wd, "galaxy.json"), "w") as fd:
762 json.dump(_galaxy_meta_data, fd)
763 print("*** Job finished successfully ***")