Mercurial > repos > astroteam > sgwb_astro_tool

comparison Model_spectrum.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\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.constants as const
25 import astropy.units as u
26 from astropy.constants import c
27 from numpy import log, pi, sqrt
28 from oda_api.data_products import ODAAstropyTable, PictureProduct
29
30 # Parameters of the phase transition
31 # Temperature in GeV
32 T_star = 0.178 # http://odahub.io/ontology#Energy_GeV
33
34 # Numbers of relativistic degrees of freedom
35 g_star = 20 # http://odahub.io/ontology#Integer
36
37 # ratio of the energy density deposited in the bubbles to the radiation energy density
38 alpha = 1.0 # http://odahub.io/ontology#Float
39
40 # beta/H : rate of the phase transition compared to Hubble rate
41 beta_H = 3.3
42
43 # 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
44 epsilon_turb = 1 # http://odahub.io/ontology#Float
45
46 # terminal velocity of bubbles
47 v_w = 0.999 # http://odahub.io/ontology#Float
48 h = 0.7 # http://odahub.io/ontology#Float
49
50 # Parameterisation='RoperPol2023' # http://odahub.io/ontology#String ; oda:allowed_value "Lewicki2022","RoperPol2023"
51
52 _galaxy_wd = os.getcwd()
53
54 with open("inputs.json", "r") as fd:
55 inp_dic = json.load(fd)
56 if "_data_product" in inp_dic.keys():
57 inp_pdic = inp_dic["_data_product"]
58 else:
59 inp_pdic = inp_dic
60
61 for vn, vv in inp_pdic.items():
62 if vn != "_selector":
63 globals()[vn] = type(globals()[vn])(vv)
64
65 c_s = 3 ** (-0.5) # speed of sound
66 F0_GW = 1.64e-5 / h**2 * (100 / g_star) ** (1 / 3.0)
67
68 # Eq.20 of arXiv:1512.06239, corrected according to Appendix A of arXiv:1004.4187
69
70 def ka_v(alpha_v, v):
71 zeta = ((2.0 / 3.0 * alpha_v + alpha_v**2) ** 0.5 + (1 / 3.0) ** 0.5) / (
72 1 + alpha_v
73 )
74 kappa_a = (
75 6.9
76 * v ** (6.0 / 5.0)
77 * alpha_v
78 / (1.36 - 0.037 * alpha_v**0.5 + alpha_v)
79 )
80 kappa_b = alpha_v ** (2.0 / 5.0) / (
81 0.017 + (0.997 + alpha_v) ** (2.0 / 5.0)
82 )
83 kappa_c = alpha_v**0.5 / (0.135 + (0.98 + alpha_v) ** 0.5)
84 kappa_d = alpha_v / (0.73 + 0.083 * alpha_v**0.6 + alpha_v)
85 deltak = -0.9 * log(alpha_v**0.5 / (1 + alpha_v**0.5))
86 if v < c_s:
87 return (
88 c_s ** (11.0 / 5.0)
89 * kappa_a
90 * kappa_b
91 / (
92 (c_s ** (11.0 / 5.0) - v ** (11.0 / 5.0)) * kappa_b
93 + v * c_s ** (6.0 / 5.0) * kappa_a
94 )
95 )
96 elif v > zeta:
97 return (
98 (zeta - 1) ** 3
99 * (zeta / v) ** (5.0 / 2)
100 * kappa_c
101 * kappa_d
102 / (
103 ((zeta - 1) ** 3 - (v - 1) ** 3)
104 * zeta ** (5.0 / 2.0)
105 * kappa_c
106 + (v - 1) ** 3 * kappa_d
107 )
108 )
109 else:
110 return (
111 kappa_b
112 + (v - c_s) * deltak
113 + (v - c_s) ** 3
114 / (zeta - c_s) ** 3
115 * (kappa_c - kappa_b - (zeta - c_s) * deltak)
116 )
117
118 kappa_v = ka_v(alpha, v_w)
119 kappa_therm = 1 - kappa_v
120 print(
121 "Fraction of released energy in bulk motion:",
122 kappa_v,
123 ", Thermal energy fraction:",
124 kappa_therm,
125 )
126 Omega_star = (kappa_v) * alpha / (1 + alpha)
127 print("Omega in bulk motion:", Omega_star)
128
129 # Comoving Hubble rate at the phase transition Eq. 11 of arXiv:1512.06239
130 def H_star(T_star):
131 return (
132 16.5e-6
133 * (T_star / (100.0 * u.GeV))
134 * (g_star / 100) ** (1.0 / 6.0)
135 * u.Hz
136 )
137
138 Hstar = H_star(T_star * u.GeV)
139 logHstar = np.log10(Hstar / u.Hz)
140
141 fmin = logHstar.value - 5
142 fmax = logHstar.value + 5
143 ff = np.logspace(fmin, fmax, 101)
144
145 # HL model formula
146 def GW_sound(f, T_star, alpha, beta_H, v_w):
147 kappa_v = ka_v(alpha, v_w)
148 K = kappa_v * alpha / (1 + alpha)
149 lambda_star = (8 * pi) ** (1 / 3) * max([v_w, c_s]) / (beta_H * Hstar) * c
150 Delta_w = sqrt((v_w - c_s) ** 2) / v_w
151 s2 = Delta_w * lambda_star * f / c
152 print((c / lambda_star).cgs, (c / (Delta_w * lambda_star)).cgs)
153 s1 = lambda_star * f / c
154 M = (
155 16
156 * (1 + Delta_w ** (-3)) ** (2 / 3.0)
157 * (Delta_w * s1) ** 3
158 / ((1 + s1**3) ** (2.0 / 3.0) * (3 + s2**2) ** 2)
159 )
160 factor = (K * lambda_star * Hstar) ** 2 / (
161 sqrt(K) + lambda_star * Hstar / c
162 )
163 mu = 4.78 - 6.27 * Delta_w + 3.34 * Delta_w**2
164 ff = f / u.Hz
165 dlnf = (ff[1] - ff[0]) / ff[0]
166 mu = sum(M) * dlnf
167 B = 1e-2 / mu
168 return (3 * B * factor / c**2 * F0_GW * M).cgs
169
170 # Eq 1 of the new Overleaf, from Alberto
171
172 def GW_turb_Theo(f, T_star, alpha, beta_H, v_w, epsilon_turb):
173 (const.c / H_star(T_star)).to(u.Mpc)
174 kappa_v = ka_v(alpha, v_w)
175 Omega_star = (kappa_v) * alpha / (1 + alpha)
176 lambda_star = (
177 (8 * pi) ** (1 / 3) * max([v_w, c_s]) / (beta_H * H_star(T_star))
178 ) * c # characteristic light-crossing distance scale
179
180 Delta_w = (
181 sqrt((v_w - c_s) ** 2) / v_w
182 ) # <1, the width of the bubble wall in units of R_star
183 R_peak = (
184 Delta_w * lambda_star
185 ) # the GW spectrum peaks ad the frequency correspoding to the bubble wall width
186 u_star = sqrt(0.75 * epsilon_turb * Omega_star) # alfven velocity
187 dt_fin = 2 * lambda_star / u_star / c # eddy processing time
188
189 T_GW = np.log10(1 + (H_star(T_star) * dt_fin / (2 * pi))) * (
190 f < 1 / dt_fin
191 ) + np.log10(1 + (H_star(T_star) / (2 * pi * f))) * (f >= 1 / dt_fin)
192 T_GW = T_GW**2
193 alpha_pi = 2.15
194 P_pi = (1 + ((lambda_star * f) / (2.2 * c)) ** alpha_pi) ** (
195 -11 / (3 * alpha_pi)
196 )
197 M = (
198 (lambda_star * f / c) ** 3
199 * (4 * pi**2 * T_GW * P_pi)
200 / (0.84 * (lambda_star * H_star(T_star) / c) ** 2)
201 )
202 res = (
203 3
204 * (1.75 * 10 ** (-3))
205 * (Omega_star * epsilon_turb) ** 2
206 * (lambda_star * H_star(T_star) / c) ** 2
207 * 1.64
208 * 10 ** (-5)
209 * (100 / g_star) ** (1 / 3.0)
210 / h**2
211 )
212 return res * M
213
214 def GW_turb_Andrii(f, T_star, alpha, beta_H, v_w, epsilon_turb):
215
216 # Eq. 1 of 2307.10744
217 lambda_star = (
218 (8 * pi) ** (1 / 3) * max([v_w, c_s]) / (beta_H * Hstar)
219 ) * c # characteristic light-crossing distance scale
220
221 kappa_v = ka_v(alpha, v_w)
222 K = kappa_v * alpha / (1 + alpha)
223
224 # Eq. 10 of 2307.10744
225 Omega_star = epsilon_turb * K
226
227 # Eq. 13 of 2307.10744 and text after it
228 u_star = sqrt(0.75 * Omega_star)
229 dt_fin = 2 * lambda_star / u_star / c
230 s3 = dt_fin * f
231 s1 = f * lambda_star / c
232 print("u_star=", u_star, "lambda_star=", lambda_star.cgs)
233
234 # Eq. 15 of 2307.10744
235 T_GW = np.log(1 + Hstar * dt_fin / (2 * pi)) * (s3 < 1) + np.log(
236 1 + lambda_star * Hstar / c / (2 * pi * s1)
237 ) * (s3 >= 1)
238 T_GW = T_GW**2
239
240 # Eq. 17 of 2307.10744
241 alpha_pi = 2.15
242 s_pi = 2.2
243 P_pi = (1 + (s1 / s_pi) ** alpha_pi) ** (-11 / (3 * alpha_pi))
244
245 # Eq 18,19,20 of 2307.10744
246 A = 2e-3 * 1.4 * 0.6
247
248 # Eq. 14 of 2307.10744
249 Sturb = (
250 4
251 * pi**2
252 * s1**3
253 * T_GW
254 / (lambda_star * Hstar / c) ** 2
255 * P_pi
256 / 1.4
257 / 0.6
258 )
259
260 # Eq 9 of 2307.10744
261 res = (
262 3 * A * Omega_star**2 * (lambda_star * Hstar / c) ** 2 * F0_GW * Sturb
263 )
264 return res
265
266 # Eq 19 of 2308.08546
267 def GW_Ellis(f, T_star, alpha, beta_H, v_w, epsilon_turb):
268 A = 5.1e-2
269 a = 2.4
270 b = 2.4
271 c = 4.0
272 Hstar = H_star(T_star).cgs
273 f_H = Hstar / 2 / pi
274 fp = 0.7 * beta_H * f_H
275 lambda_star = (2 * pi * 3e10 * u.cm / u.s / fp).cgs
276
277 SH = 1 / (1 + (f / f_H) ** (a - 3))
278 Omega_star = beta_H ** (-2) * alpha / (1 + alpha) * A
279 Omega_B = epsilon_turb * Omega_star / 2.0
280 Omega_gamma = 2 / g_star
281 B = 3e-6 * (Omega_B / Omega_gamma) ** 0.5
282 res = (
283 beta_H ** (-2)
284 * alpha
285 / (1 + alpha)
286 * A
287 * (a + b) ** c
288 / (b * (f / fp) ** (-a / c) + a * (f / fp) ** (b / c)) ** c
289 * SH
290 )
291 return (
292 1.6e-5 * (g_star / 100) ** (-1 / 3.0) * res / h**2,
293 lambda_star.cgs.value,
294 B,
295 )
296
297 d = np.genfromtxt("NANOGrav23.csv")
298 gammas_nano = d[:, 0]
299 As_nano = 10 ** d[:, 1]
300 ps_nano = 5 - gammas_nano
301
302 d = np.genfromtxt("EPTA.csv")
303 gammas_epta = d[:, 0]
304 As_epta = 10 ** d[:, 1]
305 ps_epta = 5 - gammas_epta
306
307 d = np.genfromtxt("PPTA.csv")
308 gammas_ppta = d[:, 0]
309 As_ppta = 10 ** d[:, 1]
310 ps_ppta = 5 - gammas_ppta
311
312 H0 = 70 * (u.km / u.s) / u.Mpc
313
314 GW_t = GW_turb_Theo(
315 ff * u.Hz, T_star * u.GeV, alpha, beta_H, v_w, epsilon_turb
316 )
317 # plt.plot(ff,GW_t,color='magenta')
318 GW_t = GW_turb_Andrii(
319 ff * u.Hz, T_star * u.GeV, alpha, beta_H, v_w, epsilon_turb
320 )
321 GW_s = GW_sound(ff * u.Hz, T_star * u.GeV, alpha, beta_H, v_w)
322 GW = GW_s + GW_t
323
324 GW1 = GW_Ellis(ff * u.Hz, T_star * u.GeV, alpha, beta_H, v_w, epsilon_turb)[0]
325
326 # if(Parameterisation=='Lewicki2022'):
327 # plt.plot(ff,GW1,color='magenta',alpha=0.5,linewidth=4,label='2208.11697')
328 # else:
329 plt.plot(ff, GW_t, color="blue", linestyle="dashed", label="turbulence")
330 plt.plot(ff, GW_s, color="red", linestyle="dotted", label="sound waves")
331 plt.plot(ff, GW, linewidth=4, color="black", alpha=0.5, label="total")
332
333 fref = (1 / u.yr).cgs.value
334 lgfmin = np.log10(fref / 10.0)
335 lgfmax = np.log10(fref / 2.0)
336 fff = np.logspace(lgfmin, lgfmax, 10) * u.Hz
337 min_nano = np.ones(len(fff))
338 max_nano = np.zeros(len(fff))
339 for i in range(len(As_nano)):
340 spec = (
341 2
342 * pi**2
343 / 3
344 / H0**2
345 * fff**2
346 * As_nano[i] ** 2
347 * (fff / fref) ** (3 - gammas_nano[i])
348 ).cgs.value
349 min_nano = np.minimum(spec, min_nano)
350 max_nano = np.maximum(spec, max_nano)
351 # plt.plot(ff,spec)
352 plt.fill_between(
353 fff.value, min_nano, max_nano, color="red", alpha=0.5, label="NANOGrav"
354 )
355 min_epta = np.ones(len(fff))
356 max_epta = np.zeros(len(fff))
357 for i in range(len(As_epta)):
358 spec = (
359 2
360 * pi**2
361 / 3
362 / H0**2
363 * fff**2
364 * As_epta[i] ** 2
365 * (fff / fref) ** (3 - gammas_epta[i])
366 ).cgs.value
367 min_epta = np.minimum(spec, min_epta)
368 max_epta = np.maximum(spec, max_epta)
369 # plt.plot(ff,spec)
370 plt.fill_between(
371 fff.value, min_epta, max_epta, color="blue", alpha=0.5, label="EPTA"
372 )
373 min_ppta = np.ones(len(fff))
374 max_ppta = np.zeros(len(fff))
375 for i in range(len(As_ppta)):
376 spec = (
377 2
378 * pi**2
379 / 3
380 / H0**2
381 * fff**2
382 * As_ppta[i] ** 2
383 * (fff / fref) ** (3 - gammas_ppta[i])
384 ).cgs.value
385 min_ppta = np.minimum(spec, min_ppta)
386 max_ppta = np.maximum(spec, max_ppta)
387 # plt.plot(ff,spec)
388 plt.fill_between(
389 fff.value, min_ppta, max_ppta, color="green", alpha=0.5, label="PPTA"
390 )
391
392 maxGW = max(GW)
393 ind = np.argmax(GW)
394 fmax = ff[ind]
395 plt.xscale("log")
396 plt.yscale("log")
397 plt.ylim(maxGW / 1e5, maxGW * 10)
398 plt.xlim(fmax / 1e3, fmax * 1e3)
399 plt.grid()
400 plt.xscale("log")
401 plt.yscale("log")
402 plt.legend(loc="upper left")
403 plt.xlabel("$f$, Hz")
404 plt.ylabel("$\Omega_{gw}(f)$")
405 plt.savefig("Spectrum.png", format="png", bbox_inches="tight")
406
407 bin_image = PictureProduct.from_file("Spectrum.png")
408 from astropy.table import Table
409
410 # if(Parameterisation=='Lewicki2022'):
411 # data=[ff,GW1]
412 # names=('f[Hz]','Omega_gw')
413 # else:
414 data = [ff, GW_s, GW_t]
415 names = ("f[Hz]", "Omega_sound_waves", "Omega_turbulence")
416 spectrum = ODAAstropyTable(Table(data, names=names))
417
418 png = bin_image # http://odahub.io/ontology#ODAPictureProduct
419 astropy_table = spectrum # http://odahub.io/ontology#ODAAstropyTable
420
421 # output gathering
422 _galaxy_meta_data = {}
423 _oda_outs = []
424 _oda_outs.append(("out_Model_spectrum_png", "png_galaxy.output", png))
425 _oda_outs.append(
426 (
427 "out_Model_spectrum_astropy_table",
428 "astropy_table_galaxy.output",
429 astropy_table,
430 )
431 )
432
433 for _outn, _outfn, _outv in _oda_outs:
434 _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn)
435 if isinstance(_outv, str) and os.path.isfile(_outv):
436 shutil.move(_outv, _galaxy_outfile_name)
437 _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
438 elif getattr(_outv, "write_fits_file", None):
439 _outv.write_fits_file(_galaxy_outfile_name)
440 _galaxy_meta_data[_outn] = {"ext": "fits"}
441 elif getattr(_outv, "write_file", None):
442 _outv.write_file(_galaxy_outfile_name)
443 _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
444 else:
445 with open(_galaxy_outfile_name, "w") as fd:
446 json.dump(_outv, fd, cls=CustomJSONEncoder)
447 _galaxy_meta_data[_outn] = {"ext": "json"}
448
449 with open(os.path.join(_galaxy_wd, "galaxy.json"), "w") as fd:
450 json.dump(_galaxy_meta_data, fd)
451 print("*** Job finished successfully ***")