"""
等温・灰色・非散乱ガス層(両壁黒体)の P1 解析解と厳密解の比較
  - G(x) の分布
  - 壁面正味放射熱流束 q_w の光学厚さ依存性
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import expn

SIGMA = 5.670374419e-8  # Stefan-Boltzmann定数 [W/m^2/K^4]

# ---- 条件 ----
kappa = 1.0      # 吸収係数 [1/m]
Tg = 1000.0      # ガス温度 [K]
Tw = 300.0       # 壁温 [K]
Ls = [0.1, 1.0, 10.0]   # 層厚さ [m] (tau = 0.1, 1, 10)

Delta = Tg**4 - Tw**4

# ---- OpenFOAM 結果の設定 ----
# 層厚さ L [m] -> 2列(x, G)のテキストファイル(空白/カンマ区切り、# 行はコメント扱い)
FOAM_FILES = {
    0.1: "foam_G_L0.1.dat",    # tau = 0.1
    1.0: "foam_G_L1.dat",      # tau = 1
    10.0: "foam_G_L10.dat",    # tau = 10
}
FOAM_X_ORIGIN = "wall"     # "wall": x が壁(0)から L まで / "center": x が中心原点(-L/2〜L/2)

# OpenFOAM の壁面正味放射熱流束 [W/m^2](qr の大きさ。符号は無視して絶対値を使う)。
# None のままなら、G の表の壁面端点の値から q_w = (G_w - 4σTw^4)/2 で算出する
# (表の端点 x=0, L が壁面の値である場合のみ有効)。
FOAM_QW = {0.1: 1.020314e+04, 1.0: 5.025475e+04, 10.0: 6.017844e+04 }


def load_foam(path, L, origin="wall"):
    """(x, G) の表を読み込み、層中心原点の x [m] と G [W/m^2] を返す。"""
    with open(path) as f:
        lines = [ln.replace(",", " ").replace("(", " ").replace(")", " ")
                 for ln in f if ln.strip() and not ln.lstrip().startswith("#")]
    data = np.array([[float(v) for v in ln.split()[:2]] for ln in lines])
    x, G = data[:, 0], data[:, 1]
    if origin == "wall":
        x = x - L / 2.0
    return x, G


foam = {}   # L -> (x, G)
for L_key, path in FOAM_FILES.items():
    try:
        foam[L_key] = load_foam(path, L_key, FOAM_X_ORIGIN)
    except OSError:
        print(f"[注意] {path} が見つからないため、L = {L_key:g} m の OpenFOAM 結果は表示しません。")

# 各ケースの OpenFOAM 壁面熱流束 q_w [W/m^2]
q_foam = {}
for L_key, qr_val in FOAM_QW.items():
    if qr_val is not None:
        q_foam[L_key] = abs(qr_val)
    elif L_key in foam:
        xf, Gf = foam[L_key]
        G_wall = 0.5 * (Gf[np.argmin(xf)] + Gf[np.argmax(xf)])   # 両壁の平均
        q_foam[L_key] = (G_wall - 4.0 * SIGMA * Tw**4) / 2.0


def G_p1(x, L, kappa):
    """P1 解析解 G(x)。x は層中心を原点とする。"""
    tau = kappa * L
    lam = np.sqrt(3.0) * kappa
    a = np.sqrt(3.0) * tau / 2.0
    D = np.sqrt(3.0) * np.sinh(a) + 1.5 * np.cosh(a)
    C = -6.0 * SIGMA * Delta / D
    return 4.0 * SIGMA * Tg**4 + C * np.cosh(lam * x)


def q_p1(tau):
    """P1 による壁面正味熱流束 [W/m^2]"""
    a = np.sqrt(3.0) * tau / 2.0
    D = np.sqrt(3.0) * np.sinh(a) + 1.5 * np.cosh(a)
    return SIGMA * Delta * (2.0 - 3.0 * np.cosh(a) / D)


def G_exact(x, L, kappa):
    """厳密解 G(x)(指数積分関数 E2)"""
    e = 0.5 * (expn(2, kappa * (L / 2 - x)) + expn(2, kappa * (L / 2 + x)))
    return 4.0 * SIGMA * Tw**4 + 4.0 * SIGMA * Delta * (1.0 - e)


def q_exact(tau):
    """厳密解による壁面正味熱流束 [W/m^2]"""
    return SIGMA * Delta * (1.0 - 2.0 * expn(3, tau))


# ---- G(x) の比較プロット ----
fig, axes = plt.subplots(1, 3, figsize=(14, 4), sharey=False)
for ax, L in zip(axes, Ls):
    x = np.linspace(-L / 2, L / 2, 400)
    ax.plot(x / L, G_exact(x, L, kappa) / (4 * SIGMA * Tg**4), "k-", label="Exact")
    ax.plot(x / L, G_p1(x, L, kappa) / (4 * SIGMA * Tg**4), "r--", label="P1")
    if L in foam:
        xf, Gf = foam[L]
        ax.plot(xf / L, Gf / (4 * SIGMA * Tg**4), "bo", ms=3,
                label="OpenFOAM (P1)")
    ax.set_title(f"tau = {kappa * L:g}")
    ax.set_xlabel("x / L")
    ax.set_ylabel(r"G / (4$\sigma T_g^4$)")
    ax.grid(alpha=0.3)
    ax.legend()
fig.tight_layout()

# ---- q_w の光学厚さ依存性 ----
tau = np.logspace(-2, 2, 200)
fig2, ax2 = plt.subplots(figsize=(6, 4.5))
ax2.semilogx(tau, q_exact(tau) / (SIGMA * Delta), "k-", label="Exact: 1 - 2E3(tau)")
ax2.semilogx(tau, q_p1(tau) / (SIGMA * Delta), "r--", label="P1")
if q_foam:
    ax2.semilogx([kappa * L for L in q_foam],
                 [q / (SIGMA * Delta) for q in q_foam.values()],
                 "bo", ms=6, label="OpenFOAM (P1)")
ax2.set_xlabel(r"optical thickness $\tau = \kappa L$")
ax2.set_ylabel(r"$q_w / (\sigma (T_g^4 - T_w^4))$")
ax2.grid(alpha=0.3, which="both")
ax2.legend()
fig2.tight_layout()

# ---- 数値の表示 ----
print(f"{'tau':>6} {'q_exact [W/m2]':>16} {'q_P1 [W/m2]':>14} {'誤差[%]':>9}")
for L in Ls:
    t = kappa * L
    qe, qp = q_exact(t), q_p1(t)
    print(f"{t:6.2f} {qe:16.4e} {qp:14.4e} {100 * (qp - qe) / qe:9.2f}")

for L_key, (xf, Gf) in foam.items():
    ref = G_p1(xf, L_key, kappa)
    rel = np.abs(Gf - ref) / ref
    print(f"OpenFOAM vs P1解析解 (L={L_key:g} m): "
          f"最大相対誤差 {100 * rel.max():.3f} %, RMS {100 * np.sqrt(np.mean(rel**2)):.3f} %")

if q_foam:
    print(f"\n{'tau':>6} {'q_foam [W/m2]':>15} {'q_P1 [W/m2]':>14} {'q_exact [W/m2]':>16} "
          f"{'foam/P1-1[%]':>13} {'foam/exact-1[%]':>16}")
    for L_key, qf in q_foam.items():
        t = kappa * L_key
        qp, qe = q_p1(t), q_exact(t)
        print(f"{t:6.2f} {qf:15.4e} {qp:14.4e} {qe:16.4e} "
              f"{100 * (qf / qp - 1):13.2f} {100 * (qf / qe - 1):16.2f}")

plt.show()