Mercurial > repos > astroteam > sgwb_astro_tool
changeset 0:c9fc89ee996e draft default tip
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 | |
| files | Model_spectrum.py Phase_transition_parameters.py sgwb_astro_tool.xml |
| diffstat | 3 files changed, 1336 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Model_spectrum.py Wed Apr 17 10:29:23 2024 +0000 @@ -0,0 +1,451 @@ +#!/usr/bin/env python +# coding: utf-8 + +# flake8: noqa + +import json +import os +import shutil + +from oda_api.json import CustomJSONEncoder + +get_ipython().run_cell_magic( # noqa: F821 + "bash", + "", + "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", +) + +import matplotlib.pyplot as plt + +# Loading necessary packages +import numpy as np + +get_ipython().run_line_magic("matplotlib", "inline") # noqa: F821 +import astropy.constants as const +import astropy.units as u +from astropy.constants import c +from numpy import log, pi, sqrt +from oda_api.data_products import ODAAstropyTable, PictureProduct + +# Parameters of the phase transition +# Temperature in GeV +T_star = 0.178 # http://odahub.io/ontology#Energy_GeV + +# Numbers of relativistic degrees of freedom +g_star = 20 # http://odahub.io/ontology#Integer + +# ratio of the energy density deposited in the bubbles to the radiation energy density +alpha = 1.0 # http://odahub.io/ontology#Float + +# beta/H : rate of the phase transition compared to Hubble rate +beta_H = 3.3 + +# 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 +epsilon_turb = 1 # http://odahub.io/ontology#Float + +# terminal velocity of bubbles +v_w = 0.999 # http://odahub.io/ontology#Float +h = 0.7 # http://odahub.io/ontology#Float + +# Parameterisation='RoperPol2023' # http://odahub.io/ontology#String ; oda:allowed_value "Lewicki2022","RoperPol2023" + +_galaxy_wd = os.getcwd() + +with open("inputs.json", "r") as fd: + inp_dic = json.load(fd) +if "_data_product" in inp_dic.keys(): + inp_pdic = inp_dic["_data_product"] +else: + inp_pdic = inp_dic + +for vn, vv in inp_pdic.items(): + if vn != "_selector": + globals()[vn] = type(globals()[vn])(vv) + +c_s = 3 ** (-0.5) # speed of sound +F0_GW = 1.64e-5 / h**2 * (100 / g_star) ** (1 / 3.0) + +# Eq.20 of arXiv:1512.06239, corrected according to Appendix A of arXiv:1004.4187 + +def ka_v(alpha_v, v): + zeta = ((2.0 / 3.0 * alpha_v + alpha_v**2) ** 0.5 + (1 / 3.0) ** 0.5) / ( + 1 + alpha_v + ) + kappa_a = ( + 6.9 + * v ** (6.0 / 5.0) + * alpha_v + / (1.36 - 0.037 * alpha_v**0.5 + alpha_v) + ) + kappa_b = alpha_v ** (2.0 / 5.0) / ( + 0.017 + (0.997 + alpha_v) ** (2.0 / 5.0) + ) + kappa_c = alpha_v**0.5 / (0.135 + (0.98 + alpha_v) ** 0.5) + kappa_d = alpha_v / (0.73 + 0.083 * alpha_v**0.6 + alpha_v) + deltak = -0.9 * log(alpha_v**0.5 / (1 + alpha_v**0.5)) + if v < c_s: + return ( + c_s ** (11.0 / 5.0) + * kappa_a + * kappa_b + / ( + (c_s ** (11.0 / 5.0) - v ** (11.0 / 5.0)) * kappa_b + + v * c_s ** (6.0 / 5.0) * kappa_a + ) + ) + elif v > zeta: + return ( + (zeta - 1) ** 3 + * (zeta / v) ** (5.0 / 2) + * kappa_c + * kappa_d + / ( + ((zeta - 1) ** 3 - (v - 1) ** 3) + * zeta ** (5.0 / 2.0) + * kappa_c + + (v - 1) ** 3 * kappa_d + ) + ) + else: + return ( + kappa_b + + (v - c_s) * deltak + + (v - c_s) ** 3 + / (zeta - c_s) ** 3 + * (kappa_c - kappa_b - (zeta - c_s) * deltak) + ) + +kappa_v = ka_v(alpha, v_w) +kappa_therm = 1 - kappa_v +print( + "Fraction of released energy in bulk motion:", + kappa_v, + ", Thermal energy fraction:", + kappa_therm, +) +Omega_star = (kappa_v) * alpha / (1 + alpha) +print("Omega in bulk motion:", Omega_star) + +# Comoving Hubble rate at the phase transition Eq. 11 of arXiv:1512.06239 +def H_star(T_star): + return ( + 16.5e-6 + * (T_star / (100.0 * u.GeV)) + * (g_star / 100) ** (1.0 / 6.0) + * u.Hz + ) + +Hstar = H_star(T_star * u.GeV) +logHstar = np.log10(Hstar / u.Hz) + +fmin = logHstar.value - 5 +fmax = logHstar.value + 5 +ff = np.logspace(fmin, fmax, 101) + +# HL model formula +def GW_sound(f, T_star, alpha, beta_H, v_w): + kappa_v = ka_v(alpha, v_w) + K = kappa_v * alpha / (1 + alpha) + lambda_star = (8 * pi) ** (1 / 3) * max([v_w, c_s]) / (beta_H * Hstar) * c + Delta_w = sqrt((v_w - c_s) ** 2) / v_w + s2 = Delta_w * lambda_star * f / c + print((c / lambda_star).cgs, (c / (Delta_w * lambda_star)).cgs) + s1 = lambda_star * f / c + M = ( + 16 + * (1 + Delta_w ** (-3)) ** (2 / 3.0) + * (Delta_w * s1) ** 3 + / ((1 + s1**3) ** (2.0 / 3.0) * (3 + s2**2) ** 2) + ) + factor = (K * lambda_star * Hstar) ** 2 / ( + sqrt(K) + lambda_star * Hstar / c + ) + mu = 4.78 - 6.27 * Delta_w + 3.34 * Delta_w**2 + ff = f / u.Hz + dlnf = (ff[1] - ff[0]) / ff[0] + mu = sum(M) * dlnf + B = 1e-2 / mu + return (3 * B * factor / c**2 * F0_GW * M).cgs + +# Eq 1 of the new Overleaf, from Alberto + +def GW_turb_Theo(f, T_star, alpha, beta_H, v_w, epsilon_turb): + (const.c / H_star(T_star)).to(u.Mpc) + kappa_v = ka_v(alpha, v_w) + Omega_star = (kappa_v) * alpha / (1 + alpha) + lambda_star = ( + (8 * pi) ** (1 / 3) * max([v_w, c_s]) / (beta_H * H_star(T_star)) + ) * c # characteristic light-crossing distance scale + + Delta_w = ( + sqrt((v_w - c_s) ** 2) / v_w + ) # <1, the width of the bubble wall in units of R_star + R_peak = ( + Delta_w * lambda_star + ) # the GW spectrum peaks ad the frequency correspoding to the bubble wall width + u_star = sqrt(0.75 * epsilon_turb * Omega_star) # alfven velocity + dt_fin = 2 * lambda_star / u_star / c # eddy processing time + + T_GW = np.log10(1 + (H_star(T_star) * dt_fin / (2 * pi))) * ( + f < 1 / dt_fin + ) + np.log10(1 + (H_star(T_star) / (2 * pi * f))) * (f >= 1 / dt_fin) + T_GW = T_GW**2 + alpha_pi = 2.15 + P_pi = (1 + ((lambda_star * f) / (2.2 * c)) ** alpha_pi) ** ( + -11 / (3 * alpha_pi) + ) + M = ( + (lambda_star * f / c) ** 3 + * (4 * pi**2 * T_GW * P_pi) + / (0.84 * (lambda_star * H_star(T_star) / c) ** 2) + ) + res = ( + 3 + * (1.75 * 10 ** (-3)) + * (Omega_star * epsilon_turb) ** 2 + * (lambda_star * H_star(T_star) / c) ** 2 + * 1.64 + * 10 ** (-5) + * (100 / g_star) ** (1 / 3.0) + / h**2 + ) + return res * M + +def GW_turb_Andrii(f, T_star, alpha, beta_H, v_w, epsilon_turb): + + # Eq. 1 of 2307.10744 + lambda_star = ( + (8 * pi) ** (1 / 3) * max([v_w, c_s]) / (beta_H * Hstar) + ) * c # characteristic light-crossing distance scale + + kappa_v = ka_v(alpha, v_w) + K = kappa_v * alpha / (1 + alpha) + + # Eq. 10 of 2307.10744 + Omega_star = epsilon_turb * K + + # Eq. 13 of 2307.10744 and text after it + u_star = sqrt(0.75 * Omega_star) + dt_fin = 2 * lambda_star / u_star / c + s3 = dt_fin * f + s1 = f * lambda_star / c + print("u_star=", u_star, "lambda_star=", lambda_star.cgs) + + # Eq. 15 of 2307.10744 + T_GW = np.log(1 + Hstar * dt_fin / (2 * pi)) * (s3 < 1) + np.log( + 1 + lambda_star * Hstar / c / (2 * pi * s1) + ) * (s3 >= 1) + T_GW = T_GW**2 + + # Eq. 17 of 2307.10744 + alpha_pi = 2.15 + s_pi = 2.2 + P_pi = (1 + (s1 / s_pi) ** alpha_pi) ** (-11 / (3 * alpha_pi)) + + # Eq 18,19,20 of 2307.10744 + A = 2e-3 * 1.4 * 0.6 + + # Eq. 14 of 2307.10744 + Sturb = ( + 4 + * pi**2 + * s1**3 + * T_GW + / (lambda_star * Hstar / c) ** 2 + * P_pi + / 1.4 + / 0.6 + ) + + # Eq 9 of 2307.10744 + res = ( + 3 * A * Omega_star**2 * (lambda_star * Hstar / c) ** 2 * F0_GW * Sturb + ) + return res + +# Eq 19 of 2308.08546 +def GW_Ellis(f, T_star, alpha, beta_H, v_w, epsilon_turb): + A = 5.1e-2 + a = 2.4 + b = 2.4 + c = 4.0 + Hstar = H_star(T_star).cgs + f_H = Hstar / 2 / pi + fp = 0.7 * beta_H * f_H + lambda_star = (2 * pi * 3e10 * u.cm / u.s / fp).cgs + + SH = 1 / (1 + (f / f_H) ** (a - 3)) + Omega_star = beta_H ** (-2) * alpha / (1 + alpha) * A + Omega_B = epsilon_turb * Omega_star / 2.0 + Omega_gamma = 2 / g_star + B = 3e-6 * (Omega_B / Omega_gamma) ** 0.5 + res = ( + beta_H ** (-2) + * alpha + / (1 + alpha) + * A + * (a + b) ** c + / (b * (f / fp) ** (-a / c) + a * (f / fp) ** (b / c)) ** c + * SH + ) + return ( + 1.6e-5 * (g_star / 100) ** (-1 / 3.0) * res / h**2, + lambda_star.cgs.value, + B, + ) + +d = np.genfromtxt("NANOGrav23.csv") +gammas_nano = d[:, 0] +As_nano = 10 ** d[:, 1] +ps_nano = 5 - gammas_nano + +d = np.genfromtxt("EPTA.csv") +gammas_epta = d[:, 0] +As_epta = 10 ** d[:, 1] +ps_epta = 5 - gammas_epta + +d = np.genfromtxt("PPTA.csv") +gammas_ppta = d[:, 0] +As_ppta = 10 ** d[:, 1] +ps_ppta = 5 - gammas_ppta + +H0 = 70 * (u.km / u.s) / u.Mpc + +GW_t = GW_turb_Theo( + ff * u.Hz, T_star * u.GeV, alpha, beta_H, v_w, epsilon_turb +) +# plt.plot(ff,GW_t,color='magenta') +GW_t = GW_turb_Andrii( + ff * u.Hz, T_star * u.GeV, alpha, beta_H, v_w, epsilon_turb +) +GW_s = GW_sound(ff * u.Hz, T_star * u.GeV, alpha, beta_H, v_w) +GW = GW_s + GW_t + +GW1 = GW_Ellis(ff * u.Hz, T_star * u.GeV, alpha, beta_H, v_w, epsilon_turb)[0] + +# if(Parameterisation=='Lewicki2022'): +# plt.plot(ff,GW1,color='magenta',alpha=0.5,linewidth=4,label='2208.11697') +# else: +plt.plot(ff, GW_t, color="blue", linestyle="dashed", label="turbulence") +plt.plot(ff, GW_s, color="red", linestyle="dotted", label="sound waves") +plt.plot(ff, GW, linewidth=4, color="black", alpha=0.5, label="total") + +fref = (1 / u.yr).cgs.value +lgfmin = np.log10(fref / 10.0) +lgfmax = np.log10(fref / 2.0) +fff = np.logspace(lgfmin, lgfmax, 10) * u.Hz +min_nano = np.ones(len(fff)) +max_nano = np.zeros(len(fff)) +for i in range(len(As_nano)): + spec = ( + 2 + * pi**2 + / 3 + / H0**2 + * fff**2 + * As_nano[i] ** 2 + * (fff / fref) ** (3 - gammas_nano[i]) + ).cgs.value + min_nano = np.minimum(spec, min_nano) + max_nano = np.maximum(spec, max_nano) + # plt.plot(ff,spec) +plt.fill_between( + fff.value, min_nano, max_nano, color="red", alpha=0.5, label="NANOGrav" +) +min_epta = np.ones(len(fff)) +max_epta = np.zeros(len(fff)) +for i in range(len(As_epta)): + spec = ( + 2 + * pi**2 + / 3 + / H0**2 + * fff**2 + * As_epta[i] ** 2 + * (fff / fref) ** (3 - gammas_epta[i]) + ).cgs.value + min_epta = np.minimum(spec, min_epta) + max_epta = np.maximum(spec, max_epta) + # plt.plot(ff,spec) +plt.fill_between( + fff.value, min_epta, max_epta, color="blue", alpha=0.5, label="EPTA" +) +min_ppta = np.ones(len(fff)) +max_ppta = np.zeros(len(fff)) +for i in range(len(As_ppta)): + spec = ( + 2 + * pi**2 + / 3 + / H0**2 + * fff**2 + * As_ppta[i] ** 2 + * (fff / fref) ** (3 - gammas_ppta[i]) + ).cgs.value + min_ppta = np.minimum(spec, min_ppta) + max_ppta = np.maximum(spec, max_ppta) + # plt.plot(ff,spec) +plt.fill_between( + fff.value, min_ppta, max_ppta, color="green", alpha=0.5, label="PPTA" +) + +maxGW = max(GW) +ind = np.argmax(GW) +fmax = ff[ind] +plt.xscale("log") +plt.yscale("log") +plt.ylim(maxGW / 1e5, maxGW * 10) +plt.xlim(fmax / 1e3, fmax * 1e3) +plt.grid() +plt.xscale("log") +plt.yscale("log") +plt.legend(loc="upper left") +plt.xlabel("$f$, Hz") +plt.ylabel("$\Omega_{gw}(f)$") +plt.savefig("Spectrum.png", format="png", bbox_inches="tight") + +bin_image = PictureProduct.from_file("Spectrum.png") +from astropy.table import Table + +# if(Parameterisation=='Lewicki2022'): +# data=[ff,GW1] +# names=('f[Hz]','Omega_gw') +# else: +data = [ff, GW_s, GW_t] +names = ("f[Hz]", "Omega_sound_waves", "Omega_turbulence") +spectrum = ODAAstropyTable(Table(data, names=names)) + +png = bin_image # http://odahub.io/ontology#ODAPictureProduct +astropy_table = spectrum # http://odahub.io/ontology#ODAAstropyTable + +# output gathering +_galaxy_meta_data = {} +_oda_outs = [] +_oda_outs.append(("out_Model_spectrum_png", "png_galaxy.output", png)) +_oda_outs.append( + ( + "out_Model_spectrum_astropy_table", + "astropy_table_galaxy.output", + astropy_table, + ) +) + +for _outn, _outfn, _outv in _oda_outs: + _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn) + if isinstance(_outv, str) and os.path.isfile(_outv): + shutil.move(_outv, _galaxy_outfile_name) + _galaxy_meta_data[_outn] = {"ext": "_sniff_"} + elif getattr(_outv, "write_fits_file", None): + _outv.write_fits_file(_galaxy_outfile_name) + _galaxy_meta_data[_outn] = {"ext": "fits"} + elif getattr(_outv, "write_file", None): + _outv.write_file(_galaxy_outfile_name) + _galaxy_meta_data[_outn] = {"ext": "_sniff_"} + else: + with open(_galaxy_outfile_name, "w") as fd: + json.dump(_outv, fd, cls=CustomJSONEncoder) + _galaxy_meta_data[_outn] = {"ext": "json"} + +with open(os.path.join(_galaxy_wd, "galaxy.json"), "w") as fd: + json.dump(_galaxy_meta_data, fd) +print("*** Job finished successfully ***")
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Phase_transition_parameters.py Wed Apr 17 10:29:23 2024 +0000 @@ -0,0 +1,763 @@ +#!/usr/bin/env python +# coding: utf-8 + +# flake8: noqa + +import json +import os +import shutil + +from oda_api.json import CustomJSONEncoder + +get_ipython().run_cell_magic( # noqa: F821 + "bash", + "", + "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", +) + +import matplotlib.pyplot as plt + +# Loading necessary packages +import numpy as np + +get_ipython().run_line_magic("matplotlib", "inline") # noqa: F821 +import astropy.units as u +from astropy.constants import c +from numpy import exp, log, pi, sqrt +from oda_api.api import ProgressReporter +from oda_api.data_products import PictureProduct + +# Parameters of the phase transition + +# 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 +epsilon_turb = 1.0 # http://odahub.io/ontology#Float + +_galaxy_wd = os.getcwd() + +with open("inputs.json", "r") as fd: + inp_dic = json.load(fd) +if "_data_product" in inp_dic.keys(): + inp_pdic = inp_dic["_data_product"] +else: + inp_pdic = inp_dic + +for vn, vv in inp_pdic.items(): + if vn != "_selector": + globals()[vn] = type(globals()[vn])(vv) + +Tmax_yrs = 10.0 # http://odahub.io/ontology#Float +Tmin_yrs = 2.0 # http://odahub.io/ontology#Float +# terminal velocity of bubbles +v_w = 0.999 # http://odahub.io/ontology#Float +h = 0.7 # http://odahub.io/ontology#Float +# Numbers of relativistic degrees of freedom +g_star = 20 # http://odahub.io/ontology#Integer + +c_s = 3 ** (-0.5) # speed of sound +F0_GW = 1.64e-5 / h**2 * (100 / g_star) ** (1 / 3.0) +pr = ProgressReporter() +pr.report_progress(stage="Progress", progress=1.0) + +# Eq.20 of arXiv:1512.06239, corrected according to Appendix A of arXiv:1004.4187 + +def ka_v(alpha_v, v): + zeta = ((2.0 / 3.0 * alpha_v + alpha_v**2) ** 0.5 + (1 / 3.0) ** 0.5) / ( + 1 + alpha_v + ) + kappa_a = ( + 6.9 + * v ** (6.0 / 5.0) + * alpha_v + / (1.36 - 0.037 * alpha_v**0.5 + alpha_v) + ) + kappa_b = alpha_v ** (2.0 / 5.0) / ( + 0.017 + (0.997 + alpha_v) ** (2.0 / 5.0) + ) + kappa_c = alpha_v**0.5 / (0.135 + (0.98 + alpha_v) ** 0.5) + kappa_d = alpha_v / (0.73 + 0.083 * alpha_v**0.6 + alpha_v) + deltak = -0.9 * log(alpha_v**0.5 / (1 + alpha_v**0.5)) + if v < c_s: + return ( + c_s ** (11.0 / 5.0) + * kappa_a + * kappa_b + / ( + (c_s ** (11.0 / 5.0) - v ** (11.0 / 5.0)) * kappa_b + + v * c_s ** (6.0 / 5.0) * kappa_a + ) + ) + elif v > zeta: + return ( + (zeta - 1) ** 3 + * (zeta / v) ** (5.0 / 2) + * kappa_c + * kappa_d + / ( + ((zeta - 1) ** 3 - (v - 1) ** 3) + * zeta ** (5.0 / 2.0) + * kappa_c + + (v - 1) ** 3 * kappa_d + ) + ) + else: + return ( + kappa_b + + (v - c_s) * deltak + + (v - c_s) ** 3 + / (zeta - c_s) ** 3 + * (kappa_c - kappa_b - (zeta - c_s) * deltak) + ) + +T_star = 0.2 + +# Comoving Hubble rate at the phase transition Eq. 11 of arXiv:1512.06239 +def H_star(T_star): + return ( + 16.5e-6 + * (T_star / (100.0 * u.GeV)) + * (g_star / 100) ** (1.0 / 6.0) + * u.Hz + ) + +Hstar = H_star(T_star * u.GeV) +logHstar = np.log10(Hstar / u.Hz) + +fmin = logHstar.value - 5 +fmax = logHstar.value + 5 +ff = np.logspace(fmin, fmax, 101) + +# HL model formula +def GW_sound(f, T_star, alpha, beta_H, v_w): + Hstar = H_star(T_star) + kappa_v = ka_v(alpha, v_w) + K = kappa_v * alpha / (1 + alpha) + lambda_star = (8 * pi) ** (1 / 3) * max([v_w, c_s]) / (beta_H * Hstar) * c + Delta_w = sqrt((v_w - c_s) ** 2) / v_w + s2 = Delta_w * lambda_star * f / c + # print((c/lambda_star).cgs,(c/(Delta_w*lambda_star)).cgs) + s1 = lambda_star * f / c + M = ( + 16 + * (1 + Delta_w ** (-3)) ** (2 / 3.0) + * (Delta_w * s1) ** 3 + / ((1 + s1**3) ** (2.0 / 3.0) * (3 + s2**2) ** 2) + ) + factor = (K * lambda_star * Hstar) ** 2 / ( + sqrt(K) + lambda_star * Hstar / c + ) + mu = 4.78 - 6.27 * Delta_w + 3.34 * Delta_w**2 + ff = f / u.Hz + dlnf = (ff[1] - ff[0]) / ff[0] + mu = sum(M) * dlnf + B = 1e-2 / mu + return (3 * B * factor / c**2 * F0_GW * M).cgs + +# Eq 1 of the new Overleaf, from Alberto + +def GW_turb_Andrii(f, T_star, alpha, beta_H, v_w, epsilon_turb): + Hstar = H_star(T_star) + + # Eq. 1 of 2307.10744 + lambda_star = ( + (8 * pi) ** (1 / 3) * max([v_w, c_s]) / (beta_H * Hstar) + ) * c # characteristic light-crossing distance scale + + kappa_v = ka_v(alpha, v_w) + K = kappa_v * alpha / (1 + alpha) + + # Eq. 10 of 2307.10744 + Omega_star = epsilon_turb * K + + # Eq. 13 of 2307.10744 and text after it + u_star = sqrt(0.75 * Omega_star) + dt_fin = 2 * lambda_star / u_star / c + s3 = dt_fin * f + s1 = f * lambda_star / c + # print('u_star=',u_star,'lambda_star=',lambda_star.cgs) + + # Eq. 15 of 2307.10744 + T_GW = np.log(1 + Hstar * dt_fin / (2 * pi)) * (s3 < 1) + np.log( + 1 + lambda_star * Hstar / c / (2 * pi * s1) + ) * (s3 >= 1) + T_GW = T_GW**2 + + # Eq. 17 of 2307.10744 + alpha_pi = 2.15 + s_pi = 2.2 + P_pi = (1 + (s1 / s_pi) ** alpha_pi) ** (-11 / (3 * alpha_pi)) + + # Eq 18,19,20 of 2307.10744 + A = 2e-3 * 1.4 * 0.6 + + # Eq. 14 of 2307.10744 + Sturb = ( + 4 + * pi**2 + * s1**3 + * T_GW + / (lambda_star * Hstar / c) ** 2 + * P_pi + / 1.4 + / 0.6 + ) + + # Eq 9 of 2307.10744 + res = ( + 3 * A * Omega_star**2 * (lambda_star * Hstar / c) ** 2 * F0_GW * Sturb + ) + + Omega_B = Omega_star / 2.0 + Omega_gamma = 2 / g_star + B = 3e-6 * (Omega_B / Omega_gamma) ** 0.5 + return res, lambda_star, B + +H0 = 70 * (u.km / u.s) / u.Mpc + +d = np.genfromtxt("NANOGrav23.csv") +gammas_nano = d[:, 0] +As_nano = 10 ** d[:, 1] +ps_nano = 5 - gammas_nano + +d = np.genfromtxt("EPTA.csv") +gammas_epta = d[:, 0] +As_epta = 10 ** d[:, 1] +ps_epta = 5 - gammas_epta + +d = np.genfromtxt("PPTA.csv") +gammas_ppta = d[:, 0] +As_ppta = 10 ** d[:, 1] +ps_ppta = 5 - gammas_ppta + +fref = (1 / u.yr).cgs.value +lgfmin = np.log10(fref / Tmax_yrs) +lgfmax = np.log10(fref / Tmin_yrs) +fff = np.logspace(lgfmin, lgfmax, 10) * u.Hz +min_nano = np.ones(len(fff)) +max_nano = np.zeros(len(fff)) +for i in range(len(As_nano)): + spec = ( + 2 + * pi**2 + / 3 + / H0**2 + * fff**2 + * As_nano[i] ** 2 + * (fff / fref) ** (3 - gammas_nano[i]) + ).cgs.value + min_nano = np.minimum(spec, min_nano) + max_nano = np.maximum(spec, max_nano) + # plt.plot(ff,spec) +plt.fill_between( + fff.value, min_nano, max_nano, color="red", alpha=0.5, label="NANOGrav" +) +min_epta = np.ones(len(fff)) +max_epta = np.zeros(len(fff)) +for i in range(len(As_epta)): + spec = ( + 2 + * pi**2 + / 3 + / H0**2 + * fff**2 + * As_epta[i] ** 2 + * (fff / fref) ** (3 - gammas_epta[i]) + ).cgs.value + min_epta = np.minimum(spec, min_epta) + max_epta = np.maximum(spec, max_epta) + # plt.plot(ff,spec) +plt.fill_between( + fff.value, min_epta, max_epta, color="blue", alpha=0.5, label="EPTA" +) +min_ppta = np.ones(len(fff)) +max_ppta = np.zeros(len(fff)) +for i in range(len(As_ppta)): + spec = ( + 2 + * pi**2 + / 3 + / H0**2 + * fff**2 + * As_ppta[i] ** 2 + * (fff / fref) ** (3 - gammas_ppta[i]) + ).cgs.value + min_ppta = np.minimum(spec, min_ppta) + max_ppta = np.maximum(spec, max_ppta) + # plt.plot(ff,spec) +plt.fill_between( + fff.value, min_ppta, max_ppta, color="green", alpha=0.5, label="PPTA" +) + +# PBH constraints +def PBH(alpha): + return (5.2 * (1.0 - exp(-1.1 * (abs(alpha - 1.0)) ** (1 / 3))) + 1.0) * ( + alpha > 1 + ) + +Tgrid = np.logspace(-2.6, 1, 81) +agrid = np.logspace(-1, 2, 61) +bgrid = np.logspace(0, 2, 41) +amin_b_nano_eps1 = 100.0 * np.ones(len(bgrid)) +amax_b_nano_eps1 = 0.0 * np.ones(len(bgrid)) +amin_b_epta_eps1 = 100.0 * np.ones(len(bgrid)) +amax_b_epta_eps1 = 0.0 * np.ones(len(bgrid)) +amin_b_ppta_eps1 = 100.0 * np.ones(len(bgrid)) +amax_b_ppta_eps1 = 0.0 * np.ones(len(bgrid)) +amin_T_nano_eps1 = 100.0 * np.ones(len(Tgrid)) +amax_T_nano_eps1 = 0.0 * np.ones(len(Tgrid)) +amin_T_epta_eps1 = 100.0 * np.ones(len(Tgrid)) +amax_T_epta_eps1 = 0.0 * np.ones(len(Tgrid)) +amin_T_ppta_eps1 = 100.0 * np.ones(len(Tgrid)) +amax_T_ppta_eps1 = 0.0 * np.ones(len(Tgrid)) + +Tmin_b_nano_eps1 = 100.0 * np.ones(len(bgrid)) +Tmax_b_nano_eps1 = 0.0 * np.ones(len(bgrid)) +Tmin_b_epta_eps1 = 100.0 * np.ones(len(bgrid)) +Tmax_b_epta_eps1 = 0.0 * np.ones(len(bgrid)) +Tmin_b_ppta_eps1 = 100.0 * np.ones(len(bgrid)) +Tmax_b_ppta_eps1 = 0.0 * np.ones(len(bgrid)) +B_nano_eps1 = [] +lam_nano_eps1 = [] +B_ppta_eps1 = [] +lam_ppta_eps1 = [] +B_epta_eps1 = [] +lam_epta_eps1 = [] +for i, T in enumerate(Tgrid): + print(T, len(B_nano_eps1), len(B_epta_eps1), len(B_ppta_eps1)) + pr.report_progress(stage="Progress", progress=1 + 94.0 * (i / len(Tgrid))) + for j, a in enumerate(agrid): + for k, b in enumerate(bgrid): + if b > PBH(a): + GW_t = GW_turb_Andrii( + ff * u.Hz, T * u.GeV, a, b, v_w, epsilon_turb + ) + lam = GW_t[1].cgs.value + B = GW_t[2] + GW_t = GW_t[0] + GW_s = GW_sound(ff * u.Hz, T * u.GeV, a, b, v_w) + GW = GW_s + GW_t + GW_interp = np.interp(fff.value, ff, GW) + if sum((GW_interp < max_nano) * (GW_interp > min_nano)) == len( + fff + ): + if a < amin_b_nano_eps1[k]: + amin_b_nano_eps1[k] = a + if a > amax_b_nano_eps1[k]: + amax_b_nano_eps1[k] = a + if a < amin_T_nano_eps1[i]: + amin_T_nano_eps1[i] = a + if a > amax_T_nano_eps1[i]: + amax_T_nano_eps1[i] = a + if T < Tmin_b_nano_eps1[k]: + Tmin_b_nano_eps1[k] = T + if T > Tmax_b_nano_eps1[k]: + Tmax_b_nano_eps1[k] = T + B_nano_eps1.append(B) + lam_nano_eps1.append(lam) + if sum((GW_interp < max_epta) * (GW_interp > min_epta)) == len( + fff + ): + if a < amin_b_epta_eps1[k]: + amin_b_epta_eps1[k] = a + if a > amax_b_epta_eps1[k]: + amax_b_epta_eps1[k] = a + if a < amin_T_epta_eps1[i]: + amin_T_epta_eps1[i] = a + if a > amax_T_epta_eps1[i]: + amax_T_epta_eps1[i] = a + if T < Tmin_b_epta_eps1[k]: + Tmin_b_epta_eps1[k] = T + if T > Tmax_b_epta_eps1[k]: + Tmax_b_epta_eps1[k] = T + B_epta_eps1.append(B) + lam_epta_eps1.append(lam) + if sum((GW_interp < max_ppta) * (GW_interp > min_ppta)) == len( + fff + ): + if a < amin_b_ppta_eps1[k]: + amin_b_ppta_eps1[k] = a + if a > amax_b_ppta_eps1[k]: + amax_b_ppta_eps1[k] = a + if a < amin_T_ppta_eps1[i]: + amin_T_ppta_eps1[i] = a + if a > amax_T_ppta_eps1[i]: + amax_T_ppta_eps1[i] = a + if T < Tmin_b_ppta_eps1[k]: + Tmin_b_ppta_eps1[k] = T + if T > Tmax_b_ppta_eps1[k]: + Tmax_b_ppta_eps1[k] = T + B_ppta_eps1.append(B) + lam_ppta_eps1.append(lam) + +lam_bins = np.logspace(-7, -5, 41) +B_bins = np.logspace(-7, -5, 41) + +h = plt.hist2d( + np.array(lam_nano_eps1) / 3e24, + B_nano_eps1, + bins=[lam_bins, B_bins], + vmax=1, +) +plt.colorbar() +cnt = plt.contour( + sqrt(lam_bins[1:] * lam_bins[:-1]), + sqrt(B_bins[1:] * B_bins[:-1]), + np.transpose(h[0]), + levels=[1], + colors="white", +) +cont = cnt.get_paths()[0].vertices +B_nanos_eps1 = cont[:, 1] +lam_nanos_eps1 = cont[:, 0] + +plt.xscale("log") +plt.yscale("log") + +lam_bins = np.logspace(-7, -5, 41) +B_bins = np.logspace(-7, -5, 41) + +h = plt.hist2d( + np.array(lam_epta_eps1) / 3e24, + B_epta_eps1, + bins=[lam_bins, B_bins], + vmax=1, +) +plt.colorbar() +cnt = plt.contour( + sqrt(lam_bins[1:] * lam_bins[:-1]), + sqrt(B_bins[1:] * B_bins[:-1]), + np.transpose(h[0]), + levels=[1], + colors="white", +) +cont = cnt.get_paths()[0].vertices +B_eptas_eps1 = cont[:, 1] +lam_eptas_eps1 = cont[:, 0] + +plt.xscale("log") +plt.yscale("log") + +lam_bins = np.logspace(-7, -5, 41) +B_bins = np.logspace(-7, -5, 41) + +h = plt.hist2d( + np.array(lam_ppta_eps1) / 3e24, + B_ppta_eps1, + bins=[lam_bins, B_bins], + vmax=1, +) +plt.colorbar() +cnt = plt.contour( + sqrt(lam_bins[1:] * lam_bins[:-1]), + sqrt(B_bins[1:] * B_bins[:-1]), + np.transpose(h[0]), + levels=[1], + colors="white", +) +cont = cnt.get_paths()[0].vertices +B_pptas_eps1 = cont[:, 1] +lam_pptas_eps1 = cont[:, 0] + +plt.xscale("log") +plt.yscale("log") + +import numpy as np +from matplotlib import pyplot as plt + +fig, ax = plt.subplots(figsize=(7, 7)) +x = [1e-8, 1e3] +y1 = [3e-6, 3e-6] +y2 = [3e-5, 3e-5] +plt.fill_between(x, y1, y2, color="grey", alpha=0.5) +ax.fill(lam_nanos_eps1, B_nanos_eps1, color="red", alpha=0.3, label="NANOGrav") +ax.fill(lam_eptas_eps1, B_eptas_eps1, color="blue", alpha=0.3, label="EPTA") +ax.fill(lam_pptas_eps1, B_pptas_eps1, color="green", alpha=0.3, label="PPTA") +ax.set_xscale("log") +ax.set_yscale("log") + +lam_med = np.median(lam_nanos_eps1) +B_med = np.median(B_nanos_eps1) +lam_evol = np.logspace(np.log10(lam_med), np.log10(lam_med) + 10, 100) +B_evol = B_med * (lam_evol / lam_med) ** (-5 / 4.0) +mask = B_evol > 10 ** (-8.5) * lam_evol +lam_evol = lam_evol[mask] +B_evol = B_evol[mask] +ax.annotate( + "", + xy=(lam_evol[-1], B_evol[-1]), + xytext=(lam_evol[0], B_evol[0]), + arrowprops={"arrowstyle": "->", "color": "black"}, +) + +lam_evol = np.logspace(np.log10(lam_med), np.log10(lam_med) + 10, 100) +B_evol = B_med * (lam_evol / lam_med) ** (-5 / 2.0) +mask = B_evol > 10 ** (-6.8) * lam_evol +lam_evol = lam_evol[mask] +B_evol = B_evol[mask] +ax.annotate( + "", + xy=(lam_evol[-1], B_evol[-1]), + xytext=(lam_evol[0], B_evol[0]), + arrowprops={"arrowstyle": "->", "linestyle": "dashed", "color": "black"}, +) + +x = np.logspace(-8, 3, 10) +y = 10 ** (-8.5) * x +ax.plot(x, y, color="grey", linewidth=4, linestyle="dashed") +y = 10 ** (-6.8) * x +ax.plot(x, y, color="grey", linewidth=4) + +d = np.genfromtxt("MAGIC.csv") +ax.fill_between(d[:, 0], np.zeros(len(d)), d[:, 1], color="grey", alpha=0.5) +ax.set_xlim(1e-8, 1e3) +ax.set_ylim(5e-18, 2e-5) +ax.set_xlabel(r"$\lambda_B$, Mpc") +ax.set_ylabel(r"$B$, G") +plt.text(1, 3e-17, "MAGIC '22", color="black") +plt.text( + 1e-6, + 0.1e-14, + "Alfvenic larges processed eddies", + color="black", + rotation=42, +) + +plt.text( + 0.15e-7, + 0.01e-13, + "Helicity fluctuaitons / reconnection controlled", + color="black", + rotation=41.5, +) +ax.legend(loc="lower left") + +y = [1.5e-12 * 1e1, 1.5e-12, 1.5e-12] +x = [0.3e-3, 0.03, 1e3] +plt.plot(x, y, color="olive", linewidth=4, linestyle="dotted") +ax.annotate( + "", + xytext=(1, 1.0e-13), + xy=(1, 1.5e-12), + arrowprops={"arrowstyle": "->", "color": "olive", "linewidth": 4}, +) +plt.text(0.5, 2e-12, "CTA", color="olive") + +x = [0.7e-3, 1e3] +y = [3e-11, 3e-11] +plt.plot(x, y, color="olive", linewidth=4, linestyle="dotted") +ax.annotate( + "", + xytext=(1, 3e-10), + xy=(1, 3e-11), + arrowprops={"arrowstyle": "->", "color": "olive", "linewidth": 4}, +) +plt.text(0.5, 1.1e-11, "CMB", color="olive") + +plt.savefig("B_lambdaB.png", bbox_inches="tight") + +fig = plt.figure() +mask = Tmax_b_nano_eps1 > 0 +plt.fill_between( + bgrid[mask], + Tmin_b_nano_eps1[mask], + Tmax_b_nano_eps1[mask], + alpha=0.3, + color="red", + label="NANOGrav", +) +mask = Tmax_b_epta_eps1 > 0 +plt.fill_between( + bgrid[mask], + Tmin_b_epta_eps1[mask], + Tmax_b_epta_eps1[mask], + alpha=0.3, + color="blue", + label="EPTA", +) +mask = Tmax_b_ppta_eps1 > 0 +plt.fill_between( + bgrid[mask], + Tmin_b_ppta_eps1[mask], + Tmax_b_ppta_eps1[mask], + alpha=0.3, + color="green", + label="PPTA", +) +d = np.genfromtxt("Ellis.csv") +lgT = d[:, 0] +lgb = d[:, 1] +plt.plot(10**lgb, 10**lgT, color="black", label="2308.08546") +d = np.genfromtxt("NANO_bubble.csv") +lgT = d[:, 0] +lgb = -d[:, 1] +plt.plot( + 10**lgb, 10**lgT, color="black", linestyle="dashed", label="2306.162196" +) +d = np.genfromtxt("NANO_sound.csv") +lgT = d[:, 0] +lgb = -d[:, 1] +plt.plot(10**lgb, 10**lgT, color="black", linestyle="dashed") +plt.axhline(0.160, color="black", linestyle="dotted") + +plt.xscale("log") +plt.yscale("log") +plt.ylabel(r"$T$ [GeV]") +plt.xlabel(r"$\beta/H$") +plt.xlim(0.8, 80) +plt.ylim(0.001, 10) +plt.legend(loc="lower left") +plt.savefig("T_Beta.png", bbox_inches="tight") + +fig = plt.figure() +mask = amax_b_nano_eps1 > 0 +plt.fill_between( + bgrid[mask], + amin_b_nano_eps1[mask], + amax_b_nano_eps1[mask], + alpha=0.3, + color="red", + label="NANOGrav", +) +mask = amax_b_epta_eps1 > 0 +plt.fill_between( + bgrid[mask], + amin_b_epta_eps1[mask], + amax_b_epta_eps1[mask], + alpha=0.3, + color="blue", + label="EPTA", +) +mask = amax_b_ppta_eps1 > 0 +plt.fill_between( + bgrid[mask], + amin_b_ppta_eps1[mask], + amax_b_ppta_eps1[mask], + alpha=0.3, + color="green", + label="PPTA", +) +bpbh = PBH(agrid) +x = [2, 3] +y = [1, 1] +y1 = [2, 2] +plt.fill_between(x, y, y1, color="white", linewidth=0) +plt.fill(bpbh, agrid, color="white") +plt.xscale("log") +plt.yscale("log") +plt.xlim(0.8, 80) +plt.ylim(0.2, 80) +plt.xlabel(r"$\beta/H$") +plt.ylabel(r"$\alpha$") +plt.legend(loc="upper left") +plt.savefig("Alpha_Beta.png", bbox_inches="tight") + +fig = plt.figure() +mask = amax_T_nano_eps1 > 0 +plt.fill_between( + Tgrid[mask], + amin_T_nano_eps1[mask], + amax_T_nano_eps1[mask], + alpha=0.3, + color="red", + label="NANOGrav", +) +mask = amax_T_epta_eps1 > 0 +plt.fill_between( + Tgrid[mask], + amin_T_epta_eps1[mask], + amax_T_epta_eps1[mask], + alpha=0.3, + color="blue", + label="EPTA", + linewidth=0, +) +mask = amax_T_ppta_eps1 > 0 +plt.fill_between( + Tgrid[mask], + amin_T_ppta_eps1[mask], + amax_T_ppta_eps1[mask], + alpha=0.3, + color="green", + label="PPTA", + linewidth=0, +) +d = np.genfromtxt("NANO_alpha_T.csv") +lgT = d[:, 0] +lga = d[:, 1] +plt.plot( + 10**lgT, 10**lga, color="black", label="2306.16219", linestyle="dashed" +) +d = np.genfromtxt("NANO_alpha_T1.csv") +lgT = d[:, 0] +lga = d[:, 1] +plt.plot(10**lgT, 10**lga, color="black", linestyle="dashed") +plt.axvline(0.16, color="black", linestyle="dotted") +plt.xscale("log") +plt.yscale("log") +plt.xlim(0.001, 10) +plt.ylim(0.2, 80) +plt.xlabel(r"$T$ [GeV]") +plt.ylabel(r"$\alpha$") +plt.legend(loc="upper left") +plt.savefig("Alpha_T.png", bbox_inches="tight") + +bin_image1 = PictureProduct.from_file("B_lambdaB.png") +bin_image2 = PictureProduct.from_file("T_Beta.png") +bin_image3 = PictureProduct.from_file("Alpha_Beta.png") +bin_image4 = PictureProduct.from_file("Alpha_T.png") + +B_lambdaB_png = bin_image1 # http://odahub.io/ontology#ODAPictureProduct +T_Beta_png = bin_image2 # http://odahub.io/ontology#ODAPictureProduct +Alpha_Beta_png = bin_image3 # http://odahub.io/ontology#ODAPictureProduct +Alpha_T_png = bin_image4 # http://odahub.io/ontology#ODAPictureProduct + +# output gathering +_galaxy_meta_data = {} +_oda_outs = [] +_oda_outs.append( + ( + "out_Phase_transition_parameters_B_lambdaB_png", + "B_lambdaB_png_galaxy.output", + B_lambdaB_png, + ) +) +_oda_outs.append( + ( + "out_Phase_transition_parameters_T_Beta_png", + "T_Beta_png_galaxy.output", + T_Beta_png, + ) +) +_oda_outs.append( + ( + "out_Phase_transition_parameters_Alpha_Beta_png", + "Alpha_Beta_png_galaxy.output", + Alpha_Beta_png, + ) +) +_oda_outs.append( + ( + "out_Phase_transition_parameters_Alpha_T_png", + "Alpha_T_png_galaxy.output", + Alpha_T_png, + ) +) + +for _outn, _outfn, _outv in _oda_outs: + _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn) + if isinstance(_outv, str) and os.path.isfile(_outv): + shutil.move(_outv, _galaxy_outfile_name) + _galaxy_meta_data[_outn] = {"ext": "_sniff_"} + elif getattr(_outv, "write_fits_file", None): + _outv.write_fits_file(_galaxy_outfile_name) + _galaxy_meta_data[_outn] = {"ext": "fits"} + elif getattr(_outv, "write_file", None): + _outv.write_file(_galaxy_outfile_name) + _galaxy_meta_data[_outn] = {"ext": "_sniff_"} + else: + with open(_galaxy_outfile_name, "w") as fd: + json.dump(_outv, fd, cls=CustomJSONEncoder) + _galaxy_meta_data[_outn] = {"ext": "json"} + +with open(os.path.join(_galaxy_wd, "galaxy.json"), "w") as fd: + json.dump(_galaxy_meta_data, fd) +print("*** Job finished successfully ***")
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sgwb_astro_tool.xml Wed Apr 17 10:29:23 2024 +0000 @@ -0,0 +1,122 @@ +<tool id="sgwb_astro_tool" name="SGWB" version="0.0.1+galaxy0" profile="23.0"> + <requirements> + <requirement type="package" version="8.22.2">ipython</requirement> + <requirement type="package" version="1.26.4">numpy</requirement> + <requirement type="package" version="3.8.4">matplotlib</requirement> + <requirement type="package" version="6.0.1">astropy</requirement> + <requirement type="package" version="2.2.2">pandas</requirement> + <requirement type="package" version="5.19.0">plotly</requirement> + <requirement type="package" version="1.13.0">scipy</requirement> + <requirement type="package" version="1.21.4">wget</requirement> + <requirement type="package" version="1.2.15">oda-api</requirement> + <!--Requirements string 'nb2workflow[cwl,service,rdf,mmoda]>=1.3.30 +' can't be converted automatically. Please add the galaxy/conda requirement manually or modify the requirements file!--> + <requirement type="package" version="7.16.3">nbconvert</requirement> + </requirements> + <command detect_errors="exit_code">ipython '$__tool_directory__/${_data_product._selector}.py'</command> + <configfiles> + <inputs name="inputs" filename="inputs.json" data_style="paths" /> + </configfiles> + <inputs> + <conditional name="_data_product"> + <param name="_selector" type="select" label="Data Product"> + <option value="Model_spectrum" selected="true">Model_spectrum</option> + <option value="Phase_transition_parameters" selected="false">Phase_transition_parameters</option> + </param> + <when value="Model_spectrum"> + <param name="T_star" type="float" value="0.178" label="T_star (unit: GeV)" /> + <param name="g_star" type="integer" value="20" label="g_star" /> + <param name="alpha" type="float" value="1.0" label="alpha" /> + <param name="beta_H" type="float" value="3.3" label="beta_H" /> + <param name="epsilon_turb" type="integer" value="1" label="epsilon_turb" /> + <param name="v_w" type="float" value="0.999" label="v_w" /> + <param name="h" type="float" value="0.7" label="h" /> + </when> + <when value="Phase_transition_parameters"> + <param name="epsilon_turb" type="float" value="1.0" label="epsilon_turb" /> + </when> + </conditional> + </inputs> + <outputs> + <data label="${tool.name} -> Model_spectrum png" name="out_Model_spectrum_png" format="auto" from_work_dir="png_galaxy.output"> + <filter>_data_product['_selector'] == 'Model_spectrum'</filter> + </data> + <data label="${tool.name} -> Model_spectrum astropy_table" name="out_Model_spectrum_astropy_table" format="auto" from_work_dir="astropy_table_galaxy.output"> + <filter>_data_product['_selector'] == 'Model_spectrum'</filter> + </data> + <data label="${tool.name} -> Phase_transition_parameters B_lambdaB_png" name="out_Phase_transition_parameters_B_lambdaB_png" format="auto" from_work_dir="B_lambdaB_png_galaxy.output"> + <filter>_data_product['_selector'] == 'Phase_transition_parameters'</filter> + </data> + <data label="${tool.name} -> Phase_transition_parameters T_Beta_png" name="out_Phase_transition_parameters_T_Beta_png" format="auto" from_work_dir="T_Beta_png_galaxy.output"> + <filter>_data_product['_selector'] == 'Phase_transition_parameters'</filter> + </data> + <data label="${tool.name} -> Phase_transition_parameters Alpha_Beta_png" name="out_Phase_transition_parameters_Alpha_Beta_png" format="auto" from_work_dir="Alpha_Beta_png_galaxy.output"> + <filter>_data_product['_selector'] == 'Phase_transition_parameters'</filter> + </data> + <data label="${tool.name} -> Phase_transition_parameters Alpha_T_png" name="out_Phase_transition_parameters_Alpha_T_png" format="auto" from_work_dir="Alpha_T_png_galaxy.output"> + <filter>_data_product['_selector'] == 'Phase_transition_parameters'</filter> + </data> + </outputs> + <tests> + <test expect_num_outputs="2"> + <conditional name="_data_product"> + <param name="_selector" value="Model_spectrum" /> + <param name="T_star" value="0.178" /> + <param name="g_star" value="20" /> + <param name="alpha" value="1.0" /> + <param name="beta_H" value="3.3" /> + <param name="epsilon_turb" value="1" /> + <param name="v_w" value="0.999" /> + <param name="h" value="0.7" /> + </conditional> + <assert_stdout> + <has_text text="*** Job finished successfully ***" /> + </assert_stdout> + </test> + <test expect_num_outputs="4"> + <conditional name="_data_product"> + <param name="_selector" value="Phase_transition_parameters" /> + <param name="epsilon_turb" value="1.0" /> + </conditional> + <assert_stdout> + <has_text text="*** Job finished successfully ***" /> + </assert_stdout> + </test> + </tests> + <help>This service provides a calculaiton of the power spectrum of Stochastic +Gravitational Wave Backgorund (SGWB) from a first-order cosmological +phase transition based on the parameterisations of `Roper Pol et +al. (2023) <https://arxiv.org/abs/2307.10744>`__. The power spectrum +includes two components: from the sound waves excited by collisions of +bubbles of the new phase and from the turbulence that is induced by +these collisions. + +The cosmological epoch of the phase transition is described by the +temperature, ``T_star`` and by the number(s) of relativistic degrees of +freedom, ``g_star`` that should be specified as parameters. + +The phase transition itself is characterised by phenomenological +parameters, ``alpha``, ``beta_H`` and ``epsilon_turb``, the latent heat, +the ratio of the Hubble radius to the bubble size at percolation and the +fraction of the energy otuput of the phase transition that goes into +turbulence. + +The tool ``Model spectrum`` outputs the power spectrum for fixed values +of these parameters. The tool ``Phase transition parameters`` reproduces +the constraints on the phase transition parameters from the Pulsar +Timing Array gravitational wave detectors, reported by Boyer & Neronov +(2024), including the estimate of the cosmological magnetic field +induced by turbulence. +</help> + <citations> + <citation type="bibtex">@article{RoperPol:2023bqa, + author = {Roper Pol, A. and Neronov, A. and Caprini, C. and Boyer, T. and Semikoz, D.}, + title = {{LISA and $\gamma$-ray telescopes as multi-messenger probes of a first-order cosmological phase transition}}, + eprint = {2307.10744}, + archivePrefix = {arXiv}, + primaryClass = {astro-ph.CO}, + month = {7}, + year = {2023} + }</citation> + </citations> +</tool> \ No newline at end of file