!mkdir -p ./bh_gas

import pandas as pd
import numpy as np
import minemine

h = minemine.h

snaplist = []
masklist = []
starlistmedian = []
bhlistmedian = []
xx = [1,2,3,4,5,6,7,8,9,10]
for x in xx:
    for j, snap in enumerate(range(20, 5, -1)):
        for i in [1,2,3,4]:
            star = np.hstack(pd.read_feather(f"../VirialRadius/catalog_history_mask{i}/fiducial/mask{i}_catalog.fiducial.subhalo.{snap}.feather")[f"GasSurfaceDensityIn{x}kpc"].values)
            bh =  np.hstack(pd.read_feather(f"../VirialRadius/catalog_history_mask{i}/fiducial/mask{i}_catalog.fiducial.subhalo.{snap}.feather")["SubhaloFirstBlackHoleMass"].values*1e10/h)

            snaplist.append([snap for _ in range(len(star))])
            masklist.append([i for _ in range(len(star))])
            starlistmedian.append(star)
            bhlistmedian.append(bh)


    all_data = pd.DataFrame(
            {
                    "Snap": np.hstack(snaplist),
                    "Mask": np.hstack(masklist),
                    "GasMass": np.hstack(starlistmedian),
                    "BlackHoleMass": np.hstack(bhlistmedian),
            }
    )

    all_data.to_feather(f"./bh_gas/bh_gas_mass.nonstatic2.in{x}.feather")

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os

plt.style.use("default")
plt.rcParams.update({
    "figure.dpi": 300,
    "axes.grid": True,
    "grid.linestyle": "--",
    "grid.alpha": 0.3,
    "axes.spines.top": False,
    "axes.spines.right": False,
    "legend.frameon": False,
})

# 列ごとの (上段, 下段) 半径[kpc]
# 例: [[2,3],[4,5],[6,7],[8,9],[10, None]]
x_pairs = [(1,2), (3,4), (5,6), (7,8), (9, 10)]

mm = {
    1: r"$12 \leq \log (M_{\mathrm{Halo}}/M_\odot) \leq 12.5$",
    2: r"$12.5 \leq \log (M_{\mathrm{Halo}}/M_\odot) \leq 13$",
    3: r"$13 \leq \log (M_{\mathrm{Halo}}/M_\odot) \leq 13.5$",
    4: r"$13.5 \leq \log (M_{\mathrm{Halo}}/M_\odot) \leq 14$",
}
palette = ["tab:blue", "tab:orange", "tab:green", "tab:red"]
EPS = 1e-12

def safe_log10(a):
    a = np.asarray(a, dtype=float)
    a = np.where(a > 0, a, np.nan)
    return np.log10(a + EPS)

def median_p16_p84(x):
    x = np.asarray(x, dtype=float)
    x = x[~np.isnan(x)]
    if x.size == 0:
        return np.nan, np.nan, np.nan
    med = np.nanmedian(x)
    p16 = np.nanpercentile(x, 16)
    p84 = np.nanpercentile(x, 84)
    return med, p16, p84

def compute_track(sub_df, snaps):
    x_med, x_p16, x_p84 = [], [], []
    y_med, y_p16, y_p84 = [], [], []
    for s in snaps[::2]:
        chunk = sub_df[sub_df["Snap"] == s]
        if len(chunk) == 0:
            continue
        xm, x16, x84 = median_p16_p84(safe_log10(chunk["GasMass"].values))
        ym, y16, y84 = median_p16_p84(safe_log10(chunk["BlackHoleMass"].values))
        if np.isnan(xm) or np.isnan(ym):
            continue
        x_med.append(xm); x_p16.append(x16); x_p84.append(x84)
        y_med.append(ym); y_p16.append(y16); y_p84.append(y84)
    return (np.array(x_med), np.array(x_p16), np.array(x_p84),
            np.array(y_med), np.array(y_p16), np.array(y_p84))

def plot_panel(ax, x_radius, grab_legend=False):
    if x_radius is None:
        ax.axis("off")
        return [], []
    path = f"./bh_gas/bh_gas_mass.nonstatic2.in{x_radius}.feather"
    if not os.path.exists(path):
        ax.text(0.5, 0.5, f"File not found:\n{x_radius} kpc",
                ha="center", va="center", fontsize=9)
        ax.set_title(fr"$r < {x_radius}\,\mathrm{{kpc}}$")
        return [], []

    df = pd.read_feather(path)
    snaps = sorted(df["Snap"].unique(), reverse=True)
    masks = sorted(df["Mask"].unique())

    for i, m in enumerate(masks):
        sub = df[df["Mask"] == m]
        (x_med, x_p16, x_p84,
         y_med, y_p16, y_p84) = compute_track(sub, snaps)
        if x_med.size == 0:
            continue
        xerr = np.vstack([x_med - x_p16, x_p84 - x_med])
        yerr = np.vstack([y_med - y_p16, y_p84 - y_med])
        c = palette[i % len(palette)]
        ax.errorbar(x_med, y_med, xerr=xerr, yerr=yerr,
                    fmt="none", ecolor=c, alpha=0.3, elinewidth=1, capsize=2)
        ax.plot(x_med, y_med, marker="o", markersize=3.5, markevery=5,
                color=c, linewidth=1.8, label=mm.get(m, f"Mask {m}"))

    ax.set_title(fr"$r < {x_radius}\,\mathrm{{kpc}}$", fontsize=11)
    return ax.get_legend_handles_labels() if grab_legend else ([], [])

# ---------- 2×5 図 ----------
fig, axes = plt.subplots(
    2, 5, figsize=(17, 6.2), sharex=True, sharey=True, constrained_layout=True
)
global_h, global_l = [], []

for col, (xtop, xbot) in enumerate(x_pairs):
    h, l = plot_panel(axes[0, col], xtop, grab_legend=(col == 0))
    if h and not global_l:
        global_h, global_l = h, l
    plot_panel(axes[1, col], xbot)

# 軸ラベル（共通）
fig.supxlabel(r"$\log_{10}\, (\Sigma_{\rm Gas}(<r)\,[M_{\odot}/{\rm kpc}^2])$")
fig.supylabel(r"$\log_{10}\, M_{\mathrm{BH}}$")
fig.suptitle("Gas Surface Density vs. Black Hole Mass (median ± [16,84] percentiles)")

# 目盛は外側のみ表示してスッキリ
for ax in axes.flat:
    ax.label_outer()

# 凡例：図の右外だが幅を取りすぎない位置
if global_l:
    fig.legend(global_h, global_l, loc="center left",
               bbox_to_anchor=(1.005, 0.5), fontsize=9)

plt.show()

gas¶