#!/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 ***")