#!/usr/bin/env python3
"""
層状熱電スタックの定常一次元温度プロファイルを計算・プロットするスクリプトです。

このスクリプトは、層状熱電スタックの定常一次元温度プロファイルを計算し、結果をプロットします。
異なる表面熱抵抗のモデルに対応するために、二つのモードが利用可能です。

モード
-----
1.  **delta**:
    表面熱抵抗は、厚さ ``delta_top`` / ``delta_bottom`` および熱伝導率 ``k_air`` を持つ
    明示的な上部/下部空気境界層としてモデル化されます。

2.  **conv_rad**:
    表面熱抵抗は、以下の要素から予測されます。
    -   水平平板のヌセルト数相関を用いた自然対流。
    -   熱放射。

どちらのモードでも、内部スタックは以下の層で構成されます。
上部電極 / 界面12 / TE半導体 / 界面23 / 下部電極

スクリプトは以下の情報を出力します。
-   各層/界面の抵抗
-   全体の抵抗
-   界面温度
-   TE層の温度降下

また、物理スケールと模式図スケールの2つの垂直に重ねられたプロットを含む図を保存します。

:doc:`T_distrib_usage`
"""

from __future__ import annotations

import argparse
from dataclasses import dataclass
from typing import List, Tuple

import numpy as np
import matplotlib.pyplot as plt


SIGMA_SB = 5.670374419e-8  # W / m^2 / K^4
G = 9.81  # m / s^2


@dataclass
class Region:
    """
    熱電スタックの単一の領域（層、界面、表面）を表すデータクラス。

    このクラスは、特定の領域の物理的特性（名称、厚さ、熱伝導率、カテゴリ）を
    保持し、その熱抵抗を計算するプロパティを提供します。

    :param name: str: 領域の識別名。
    :param thickness_m: float: 領域の厚さ [m]。
    :param k_W_mK: float: 領域の熱伝導率 [W/mK]。
    :param category: str: 領域のカテゴリ（例: "surface", "layer", "interface"）。
    """
    name: str
    thickness_m: float
    k_W_mK: float
    category: str  # "surface", "layer", "interface"

    @property
    def resistance_m2K_W(self) -> float:
        """
        領域の熱抵抗を計算して返します。

        熱抵抗は厚さを熱伝導率で割ったものとして定義されます。

        :returns: float: 領域の熱抵抗 [m^2K/W]。
        """
        return self.thickness_m / self.k_W_mK


def build_internal_regions(args: argparse.Namespace) -> List[Region]:
    """
    熱電スタック内部のコア層と界面の `Region` オブジェクトのリストを構築します。

    この関数は、コマンドライン引数で指定された厚さ ``d1``, ``d2``, ``d3`` と
    熱伝導率 ``k1``, ``k2``, ``k3``、および界面の厚さ ``delta12``, ``delta23`` と
    熱伝導率 ``k12``, ``k23`` を使用して、内部スタックの領域を定義します。

    :param args: argparse.Namespace: コマンドライン引数を格納するオブジェクト。
    :returns: List[Region]: 内部スタック領域のリスト。
    """
    return [
        Region("Top electrode", args.d1, args.k1, "layer"),
        Region("Interface 1/2", args.delta12, args.k12, "interface"),
        Region("TE semiconductor", args.d2, args.k2, "layer"),
        Region("Interface 2/3", args.delta23, args.k23, "interface"),
        Region("Bottom electrode", args.d3, args.k3, "layer"),
    ]


def build_regions_delta(args: argparse.Namespace) -> List[Region]:
    """
    'delta'モードでの全領域（空気境界層を含む）の `Region` オブジェクトのリストを構築します。

    この関数は、 ``build_internal_regions`` で構築された内部スタック領域の前後に、
    コマンドライン引数で指定された ``delta_top``, ``delta_bottom`` の厚さと
    ``k_air`` の熱伝導率を持つ空気境界層を追加します。

    :param args: argparse.Namespace: コマンドライン引数を格納するオブジェクト。
    :returns: List[Region]: 'delta'モードにおける全スタック領域のリスト。
    """
    return [
        Region("Top air boundary", args.delta_top, args.k_air, "surface"),
        *build_internal_regions(args),
        Region("Bottom air boundary", args.delta_bottom, args.k_air, "surface"),
    ]


def compute_piecewise_profile(
    T_top_C: float,
    regions: List[Region],
    q_W_m2: float,
    n_points_per_region: int = 80,
) -> Tuple[np.ndarray, np.ndarray, List[float]]:
    """
    既知の熱流束を持つ明示的な領域リストに対する物理的なZ座標と温度プロファイルを計算します。

    各領域で線形温度勾配を仮定し、全体の温度プロファイルと各界面の温度を計算します。
    プロファイルは連続的であり、各領域の熱伝導率と厚さに基づいて温度降下が計算されます。

    :param T_top_C: float: スタック最上部の温度 [°C]。
    :param regions: List[Region]: スタックを構成するRegionオブジェクトのリスト。
    :param q_W_m2: float: スタックを通過する正味の熱流束（下向きが正） [W/m^2]。
    :param n_points_per_region: int, optional: 各領域でプロットのために生成する点の数 (デフォルトは80)。
    :returns: Tuple[np.ndarray, np.ndarray, List[float]]: 物理的なZ座標の配列、対応する温度の配列、および各界面（最上部と最下部を含む）の温度のリストのタプル。
    """
    z_segments: List[np.ndarray] = []
    T_segments: List[np.ndarray] = []

    z0 = 0.0
    T0 = T_top_C
    interface_temperatures = [T0]

    for i, r in enumerate(regions):
        z1 = z0 + r.thickness_m
        endpoint = i == len(regions) - 1
        z = np.linspace(z0, z1, n_points_per_region, endpoint=endpoint)
        T = T0 - q_W_m2 * (z - z0) / r.k_W_mK
        z_segments.append(z)
        T_segments.append(T)

        T1 = T0 - q_W_m2 * r.thickness_m / r.k_W_mK
        interface_temperatures.append(T1)
        z0 = z1
        T0 = T1

    z_all = np.concatenate(z_segments)
    T_all = np.concatenate(T_segments)
    return z_all, T_all, interface_temperatures


def air_properties_simple(T_film_K: float, k_air_ref: float) -> Tuple[float, float, float, float]:
    """
    単純な空気物性モデルを計算します。

    室温から中程度の温度範囲での空気の物性を近似的に計算します。
    熱伝導率 ``k_air`` はユーザー指定の参照値に基づいて温度依存性を持ち、
    動粘度、熱拡散率、プラントル数も計算されます。

    :param T_film_K: float: 空気膜の絶対温度 [K]。
    :param k_air_ref: float: 参照熱伝導率 [W/mK]。
    :returns: Tuple[float, float, float, float]:
        - k_air: 熱伝導率 [W/mK]
        - nu: 動粘度 [m^2/s]
        - alpha: 熱拡散率 [m^2/s]
        - Pr: プラントル数 [-]
    """
    # Simple approximations around room-to-moderate temperatures.
    # Keep k tied to user-given reference value, with weak T dependence.
    k_air = k_air_ref * (T_film_K / 300.0) ** 0.10
    nu = 15.89e-6 * (T_film_K / 300.0) ** 1.70
    alpha = 22.0e-6 * (T_film_K / 300.0) ** 1.60
    Pr = nu / alpha
    return k_air, nu, alpha, Pr


def horizontal_plate_h_nat(
    T_surface_C: float,
    T_air_C: float,
    L_char_m: float,
    k_air_ref: float,
    side: str,
) -> float:
    """
    水平平板の自然対流熱伝達係数を計算します。

    レイリー数とヌセルト数相関を用いて、水平平板の自然対流熱伝達係数 ``h`` を計算します。
    表面と空気の温度差が極めて小さい場合は、熱伝達係数が非常に小さな正の値として返されます。
    相関は、浮力的に有利な（不安定）および不利な（安定）ケースを考慮に入れます。

    :param T_surface_C: float: 表面温度 [°C]。
    :param T_air_C: float: 周囲空気温度 [°C]。
    :param L_char_m: float: 特性長さ [m]。
    :param k_air_ref: float: 空気の参照熱伝導率 [W/mK]。
    :param side: str: サンプルの「top」（上部）または「bottom」（下部）表面を指定。
    :returns: float: 自然対流熱伝達係数 [W/m^2/K]。
    :raises ValueError: ``side`` が 'top' または 'bottom' 以外の場合。
    """
    dT = T_surface_C - T_air_C
    if abs(dT) < 1e-12:
        return 1e-9

    T_film_K = 0.5 * (T_surface_C + T_air_C) + 273.15
    beta = 1.0 / T_film_K
    k_air, nu, alpha, _pr = air_properties_simple(T_film_K, k_air_ref)

    Ra = G * beta * abs(dT) * (L_char_m ** 3) / (nu * alpha)
    Ra = max(Ra, 1e-12)

    # Determine whether orientation is buoyancy-favorable ("unstable") or stable.
    # Top surface:
    #   hot top surface  -> unstable
    #   cold top surface -> stable
    # Bottom surface:
    #   hot bottom surface  -> stable
    #   cold bottom surface -> unstable
    if side == "top":
        unstable = dT > 0.0
    elif side == "bottom":
        unstable = dT < 0.0
    else:
        raise ValueError("side must be 'top' or 'bottom'")

    if unstable:
        if Ra < 1e7:
            Nu = 0.54 * (Ra ** 0.25)
        else:
            Nu = 0.15 * (Ra ** (1.0 / 3.0))
    else:
        Nu = 0.27 * (Ra ** 0.25)

    h = Nu * k_air / L_char_m
    return max(h, 1e-9)


def top_surface_flux_downward(
    T_air_top_C: float,
    T_surface_top_C: float,
    eps_top: float,
    L_char_top_m: float,
    k_air_ref: float,
    T_rad_top_C: float | None = None,
) -> Tuple[float, float, float]:
    """
    上部周囲から上部サンプル表面への正味の下向き熱流束を計算します。

    自然対流と熱放射による熱流束の合計を計算します。
    熱流束は下向きが正と定義されます。

    :param T_air_top_C: float: 上部周囲空気温度 [°C]。
    :param T_surface_top_C: float: 上部サンプル表面温度 [°C]。
    :param eps_top: float: 上部表面の放射率 [-]。
    :param L_char_top_m: float: 上部自然対流の特性長さ [m]。
    :param k_air_ref: float: 空気の参照熱伝導率 [W/mK]。
    :param T_rad_top_C: float | None, optional: 上部放射環境温度 [°C]。Noneの場合、`T_air_top_C` が使用されます。
    :returns: Tuple[float, float, float]: 合計熱流束、対流熱流束、放射熱流束 [W/m^2] のタプル。
    """
    if T_rad_top_C is None:
        T_rad_top_C = T_air_top_C

    h_top = horizontal_plate_h_nat(
        T_surface_C=T_surface_top_C,
        T_air_C=T_air_top_C,
        L_char_m=L_char_top_m,
        k_air_ref=k_air_ref,
        side="top",
    )

    q_conv = h_top * (T_air_top_C - T_surface_top_C)
    q_rad = eps_top * SIGMA_SB * (
        (T_air_top_C + 273.15) ** 4 - (T_surface_top_C + 273.15) ** 4
    )
    return q_conv + q_rad, q_conv, q_rad


def bottom_surface_flux_downward(
    T_surface_bottom_C: float,
    T_air_bottom_C: float,
    eps_bottom: float,
    L_char_bottom_m: float,
    k_air_ref: float,
    T_rad_bottom_C: float | None = None,
) -> Tuple[float, float, float]:
    """
    下部サンプル表面から下部周囲へ向かう正味の下向き熱流束を計算します。

    自然対流と熱放射による熱流束の合計を計算します。
    熱流束は下向きが正と定義されます。

    :param T_surface_bottom_C: float: 下部サンプル表面温度 [°C]。
    :param T_air_bottom_C: float: 下部周囲空気温度 [°C]。
    :param eps_bottom: float: 下部表面の放射率 [-]。
    :param L_char_bottom_m: float: 下部自然対流の特性長さ [m]。
    :param k_air_ref: float: 空気の参照熱伝導率 [W/mK]。
    :param T_rad_bottom_C: float | None, optional: 下部放射環境温度 [°C]。Noneの場合、`T_air_bottom_C` が使用されます。
    :returns: Tuple[float, float, float]: 合計熱流束、対流熱流束、放射熱流束 [W/m^2] のタプル。
    """
    if T_rad_bottom_C is None:
        T_rad_bottom_C = T_air_bottom_C

    h_bottom = horizontal_plate_h_nat(
        T_surface_C=T_surface_bottom_C,
        T_air_C=T_air_bottom_C,
        L_char_m=L_char_bottom_m,
        k_air_ref=k_air_ref,
        side="bottom",
    )

    q_conv = h_bottom * (T_surface_bottom_C - T_air_bottom_C)
    q_rad = eps_bottom * SIGMA_SB * (
        (T_surface_bottom_C + 273.15) ** 4 - (T_rad_bottom_C + 273.15) ** 4
    )
    return q_conv + q_rad, q_conv, q_rad


def solve_conv_rad(
    args: argparse.Namespace,
    internal_regions: List[Region],
    n_points_per_region: int = 80,
) -> dict:
    """
    'conv_rad'モードで定常状態の温度プロファイルを解きます。

    上部と下部の表面温度を通じて暗黙的に決定される未知の熱流束を、
    上部と下部からの熱流束が一致するように上部表面温度について解くことで見つけます。
    二分法を使用して収束させます。

    :param args: argparse.Namespace: コマンドライン引数を格納するオブジェクト。
    :param internal_regions: List[Region]: 内部スタック領域のリスト。
    :param n_points_per_region: int, optional: 各領域でプロットのために生成する点の数 (デフォルトは80)。
    :returns: dict: 計算結果を格納した辞書。以下のキーが含まれます:
        - 'q': float: 熱流束 [W/m^2]。
        - 'R_internal': float: 内部スタックの熱抵抗 [m^2K/W]。
        - 'R_top_eff': float: 有効な上部表面熱抵抗 [m^2K/W]。
        - 'R_bottom_eff': float: 有効な下部表面熱抵抗 [m^2K/W]。
        - 'Ts_top': float: 上部サンプル表面温度 [°C]。
        - 'Ts_bottom': float: 下部サンプル表面温度 [°C]。
        - 'q_conv_top': float: 上部対流熱流束 [W/m^2]。
        - 'q_rad_top': float: 上部放射熱流束 [W/m^2]。
        - 'q_conv_bottom': float: 下部対流熱流束 [W/m^2]。
        - 'q_rad_bottom': float: 下部放射熱流束 [W/m^2]。
        - 'z_phys': np.ndarray: 物理スケールのZ座標配列。
        - 'T_phys': np.ndarray: 物理スケールの温度配列。
        - 'interface_temps_internal': List[float]: 内部界面温度のリスト。
        - 'interface_temps_full': List[float]: 全界面温度のリスト（空気温度と表面温度を含む）。
    :raises RuntimeError: 上部表面温度の解の区間を見つけられなかった場合。
    """
    R_internal = sum(r.resistance_m2K_W for r in internal_regions)
    Th = args.Th
    Tc = args.Tc

    def residual(Ts_top_C: float) -> float:
        """
        上部表面温度に対する熱流束残差を計算します。
        上部からの熱流束と下部への熱流束の差が0になるように求めます。
        """
        q_top, _, _ = top_surface_flux_downward(
            T_air_top_C=Th,
            T_surface_top_C=Ts_top_C,
            eps_top=args.eps_top,
            L_char_top_m=args.L_top,
            k_air_ref=args.k_air,
            T_rad_top_C=args.T_rad_top,
        )
        Ts_bottom_C = Ts_top_C - q_top * R_internal
        q_bottom, _, _ = bottom_surface_flux_downward(
            T_surface_bottom_C=Ts_bottom_C,
            T_air_bottom_C=Tc,
            eps_bottom=args.eps_bottom,
            L_char_bottom_m=args.L_bottom,
            k_air_ref=args.k_air,
            T_rad_bottom_C=args.T_rad_bottom,
        )
        return q_top - q_bottom

    # Robust bracket over a generous temperature range
    T_lo = min(Tc, Th) - 200.0
    T_hi = max(Tc, Th) + 200.0
    f_lo = residual(T_lo)
    f_hi = residual(T_hi)

    expand_count = 0
    while f_lo * f_hi > 0 and expand_count < 20:
        T_lo -= 100.0
        T_hi += 100.0
        f_lo = residual(T_lo)
        f_hi = residual(T_hi)
        expand_count += 1

    if f_lo * f_hi > 0:
        raise RuntimeError("Failed to bracket conv_rad solution for top surface temperature.")

    # Bisection
    for _ in range(200):
        T_mid = 0.5 * (T_lo + T_hi)
        f_mid = residual(T_mid)
        if abs(f_mid) < 1e-9:
            T_lo = T_hi = T_mid
            break
        if f_lo * f_mid <= 0:
            T_hi = T_mid
            f_hi = f_mid
        else:
            T_lo = T_mid
            f_lo = f_mid

    Ts_top = 0.5 * (T_lo + T_hi)
    q_top, q_conv_top, q_rad_top = top_surface_flux_downward(
        T_air_top_C=Th,
        T_surface_top_C=Ts_top,
        eps_top=args.eps_top,
        L_char_top_m=args.L_top,
        k_air_ref=args.k_air,
        T_rad_top_C=args.T_rad_top,
    )
    q = q_top
    Ts_bottom = Ts_top - q * R_internal
    q_bottom, q_conv_bottom, q_rad_bottom = bottom_surface_flux_downward(
        T_surface_bottom_C=Ts_bottom,
        T_air_bottom_C=Tc,
        eps_bottom=args.eps_bottom,
        L_char_bottom_m=args.L_bottom,
        k_air_ref=args.k_air,
        T_rad_bottom_C=args.T_rad_bottom,
    )

    # Internal profile only (physical layers)
    z_phys, T_phys, interface_temps_internal = compute_piecewise_profile(
        T_top_C=Ts_top,
        regions=internal_regions,
        q_W_m2=q,
        n_points_per_region=n_points_per_region,
    )

    # Add air temperatures for reporting purposes
    interface_temps_full = [Th, Ts_top, *interface_temps_internal[1:-1], Ts_bottom, Tc]

    R_top_eff = (Th - Ts_top) / q if abs(q) > 1e-18 else np.inf
    R_bottom_eff = (Ts_bottom - Tc) / q if abs(q) > 1e-18 else np.inf

    return {
        "q": q,
        "R_internal": R_internal,
        "R_top_eff": R_top_eff,
        "R_bottom_eff": R_bottom_eff,
        "Ts_top": Ts_top,
        "Ts_bottom": Ts_bottom,
        "q_conv_top": q_conv_top,
        "q_rad_top": q_rad_top,
        "q_conv_bottom": q_conv_bottom,
        "q_rad_bottom": q_rad_bottom,
        "z_phys": z_phys,
        "T_phys": T_phys,
        "interface_temps_internal": interface_temps_internal,
        "interface_temps_full": interface_temps_full,
    }


def compute_schematic_positions(
    region_names: List[str],
    region_categories: List[str],
    layer_width: float,
    interface_width: float,
    surface_width: float,
) -> Tuple[np.ndarray, List[float], List[float]]:
    """
    模式図プロットのために、各領域の相対的な位置と幅を計算します。

    各領域のカテゴリ（層、界面、表面）に基づいて、指定された表示幅を割り当て、
    スタック内の模式的な位置（エッジと中心）を決定します。

    :param region_names: List[str]: 各領域の名前のリスト。
    :param region_categories: List[str]: 各領域のカテゴリのリスト（"layer", "interface", "surface"）。
    :param layer_width: float: 「layer」カテゴリの領域に割り当てる模式的な表示幅。
    :param interface_width: float: 「interface」カテゴリの領域に割り当てる模式的な表示幅。
    :param surface_width: float: 「surface」カテゴリの領域に割り当てる模式的な表示幅。
    :returns: Tuple[np.ndarray, List[float], List[float]]:
        - np.ndarray: 各領域の模式的な開始および終了位置の配列。
        - List[float]: 各領域の模式的な幅のリスト。
        - List[float]: 各領域の模式的な中心位置のリスト。
    """
    widths = []
    for cat in region_categories:
        if cat == "layer":
            widths.append(layer_width)
        elif cat == "interface":
            widths.append(interface_width)
        else:
            widths.append(surface_width)

    edges = [0.0]
    for w in widths:
        edges.append(edges[-1] + w)

    centers = [(edges[i] + edges[i + 1]) / 2 for i in range(len(widths))]
    return np.array(edges), widths, centers


def plot_profiles_delta(
    z_phys: np.ndarray,
    T_phys: np.ndarray,
    regions: List[Region],
    interface_temperatures: List[float],
    out_png: str,
    schematic_layer_width: float,
    schematic_interface_width: float,
    schematic_surface_width: float,
) -> None:
    """
    'delta'モードの温度プロファイルを物理スケールと模式図スケールでプロットし、PNGファイルに保存します。

    上段のプロットは実際の厚さに基づく物理的な温度プロファイルを示し、
    下段のプロットは各領域に一定の幅を割り当てた模式的な温度プロファイルを示します。

    :param z_phys: np.ndarray: 物理スケールのZ座標配列。
    :param T_phys: np.ndarray: 物理スケールの温度配列。
    :param regions: List[Region]: 全スタック領域のリスト。
    :param interface_temperatures: List[float]: 各界面の温度リスト。
    :param out_png: str: 出力PNGファイル名。
    :param schematic_layer_width: float: 「layer」カテゴリの領域の模式的な表示幅。
    :param schematic_interface_width: float: 「interface」カテゴリの領域の模式的な表示幅。
    :param schematic_surface_width: float: 「surface」カテゴリの領域の模式的な表示幅。
    :returns: None
    """
    fig, (ax1, ax2) = plt.subplots(
        2, 1, figsize=(10, 8), constrained_layout=True, sharey=True
    )

    ax1.plot(z_phys * 1e6, T_phys, lw=2)
    ax1.set_title("Temperature profile (physical scale)")
    ax1.set_xlabel("Depth z [µm]")
    ax1.set_ylabel("Temperature [°C]")
    ax1.grid(True, alpha=0.3)

    z_boundary = 0.0
    for r in regions[:-1]:
        z_boundary += r.thickness_m
        ax1.axvline(z_boundary * 1e6, ls="--", alpha=0.4)

    z_start = 0.0
    y_text = T_phys.max() - 0.03 * (T_phys.max() - T_phys.min() + 1e-12)
    for r in regions:
        z_end = z_start + r.thickness_m
        z_mid = 0.5 * (z_start + z_end)
        ax1.text(z_mid * 1e6, y_text, r.name, ha="center", va="top", fontsize=8)
        z_start = z_end

    names = [r.name for r in regions]
    cats = [r.category for r in regions]
    x_edges, _, centers = compute_schematic_positions(
        names, cats,
        layer_width=schematic_layer_width,
        interface_width=schematic_interface_width,
        surface_width=schematic_surface_width,
    )

    x_plot = []
    T_plot = []
    for i in range(len(regions)):
        x0, x1 = x_edges[i], x_edges[i + 1]
        T0, T1 = interface_temperatures[i], interface_temperatures[i + 1]
        xx = np.linspace(x0, x1, 40, endpoint=(i == len(regions) - 1))
        TT = np.linspace(T0, T1, 40, endpoint=(i == len(regions) - 1))
        x_plot.append(xx)
        T_plot.append(TT)

    ax2.plot(np.concatenate(x_plot), np.concatenate(T_plot), lw=2)
    ax2.set_title("Temperature profile (schematic / deformed scale)")
    ax2.set_xlabel("Stack position [schematic scale]")
    ax2.set_ylabel("Temperature [°C]")
    ax2.grid(True, alpha=0.3)

    for i in range(len(regions) - 1):
        ax2.axvline(x_edges[i + 1], ls="--", alpha=0.4)

    ymin, ymax = ax2.get_ylim()
    for i, name in enumerate(names):
        x0, x1 = x_edges[i], x_edges[i + 1]
        ax2.axvspan(x0, x1, alpha=0.08)
        ax2.text(centers[i], ymax - 0.03 * (ymax - ymin + 1e-12), name,
                 ha="center", va="top", fontsize=8)

    fig.savefig(out_png, dpi=160)

    print()
    print(f"Saved plot: {out_png}")
    plt.show()
    plt.close(fig)


def plot_profiles_conv_rad(
    z_phys: np.ndarray,
    T_phys: np.ndarray,
    internal_regions: List[Region],
    interface_temps_internal: List[float],
    Th: float,
    Tc: float,
    Ts_top: float,
    Ts_bottom: float,
    out_png: str,
    schematic_layer_width: float,
    schematic_interface_width: float,
    schematic_surface_width: float,
) -> None:
    """
    'conv_rad'モードの温度プロファイルを物理スケールと模式図スケールでプロットし、PNGファイルに保存します。

    上段のプロットは内部スタックのみの物理的な温度プロファイルを示し、
    下段のプロットは表面の対流・放射領域を含む模式的な温度プロファイルを示します。

    :param z_phys: np.ndarray: 物理スケールのZ座標配列（内部スタックのみ）。
    :param T_phys: np.ndarray: 物理スケールの温度配列（内部スタックのみ）。
    :param internal_regions: List[Region]: 内部スタック領域のリスト。
    :param interface_temps_internal: List[float]: 内部界面の温度リスト。
    :param Th: float: 上部周囲空気温度 [°C]。
    :param Tc: float: 下部周囲空気温度 [°C]。
    :param Ts_top: float: 上部サンプル表面温度 [°C]。
    :param Ts_bottom: float: 下部サンプル表面温度 [°C]。
    :param out_png: str: 出力PNGファイル名。
    :param schematic_layer_width: float: 「layer」カテゴリの領域の模式的な表示幅。
    :param schematic_interface_width: float: 「interface」カテゴリの領域の模式的な表示幅。
    :param schematic_surface_width: float: 「surface」カテゴリの領域の模式的な表示幅。
    :returns: None
    """
    fig, (ax1, ax2) = plt.subplots(
        2, 1, figsize=(10, 8), constrained_layout=True, sharey=True
    )

    # Physical scale: internal stack only
    ax1.plot(z_phys * 1e6, T_phys, lw=2)
    ax1.set_title("Temperature profile (physical scale)")
    ax1.set_xlabel("Depth z [µm] (internal stack only)")
    ax1.set_ylabel("Temperature [°C]")
    ax1.grid(True, alpha=0.3)

    z_boundary = 0.0
    for r in internal_regions[:-1]:
        z_boundary += r.thickness_m
        ax1.axvline(z_boundary * 1e6, ls="--", alpha=0.4)

    z_start = 0.0
    y_text = T_phys.max() - 0.03 * (T_phys.max() - T_phys.min() + 1e-12)
    for r in internal_regions:
        z_end = z_start + r.thickness_m
        z_mid = 0.5 * (z_start + z_end)
        ax1.text(z_mid * 1e6, y_text, r.name, ha="center", va="top", fontsize=8)
        z_start = z_end

    # Schematic scale: include pseudo-surface regions
    names = ["Top conv+rad", *[r.name for r in internal_regions], "Bottom conv+rad"]
    cats = ["surface", *[r.category for r in internal_regions], "surface"]
    x_edges, _, centers = compute_schematic_positions(
        names, cats,
        layer_width=schematic_layer_width,
        interface_width=schematic_interface_width,
        surface_width=schematic_surface_width,
    )

    temps = [Th, Ts_top, *interface_temps_internal[1:-1], Ts_bottom, Tc]

    x_plot = []
    T_plot = []
    for i in range(len(names)):
        x0, x1 = x_edges[i], x_edges[i + 1]
        T0, T1 = temps[i], temps[i + 1]
        xx = np.linspace(x0, x1, 40, endpoint=(i == len(names) - 1))
        TT = np.linspace(T0, T1, 40, endpoint=(i == len(names) - 1))
        x_plot.append(xx)
        T_plot.append(TT)

    ax2.plot(np.concatenate(x_plot), np.concatenate(T_plot), lw=2)
    ax2.set_title("Temperature profile (schematic / deformed scale)")
    ax2.set_xlabel("Stack position [schematic scale]")
    ax2.set_ylabel("Temperature [°C]")
    ax2.grid(True, alpha=0.3)

    for i in range(len(names) - 1):
        ax2.axvline(x_edges[i + 1], ls="--", alpha=0.4)

    ymin, ymax = ax2.get_ylim()
    for i, name in enumerate(names):
        x0, x1 = x_edges[i], x_edges[i + 1]
        ax2.axvspan(x0, x1, alpha=0.08)
        ax2.text(centers[i], ymax - 0.03 * (ymax - ymin + 1e-12), name,
                 ha="center", va="top", fontsize=8)

    fig.savefig(out_png, dpi=160)
    plt.show()
    plt.close(fig)


def print_report_delta(
    args: argparse.Namespace,
    regions: List[Region],
    q: float,
    interface_temperatures: List[float],
) -> None:
    """
    'delta'モードでの計算結果をコンソールに出力します。

    入力温度、各領域の熱抵抗、総熱抵抗、熱流束、および各界面温度と
    TE半導体層の温度降下を含む詳細なレポートを表示します。

    :param args: argparse.Namespace: コマンドライン引数を格納するオブジェクト。
    :param regions: List[Region]: 全スタック領域のリスト。
    :param q: float: 計算された熱流束 [W/m^2]。
    :param interface_temperatures: List[float]: 各界面の温度リスト。
    :returns: None
    """
    total_R = sum(r.resistance_m2K_W for r in regions)

    print("=== Mode ===")
    print("surface_mode = delta")
    print()

    print("=== Input temperatures ===")
    print(f"Top air temperature    Th = {args.Th:.6g} °C")
    print(f"Bottom air temperature Tc = {args.Tc:.6g} °C")
    print()

    print("=== Regions and thermal resistances (top -> bottom) ===")
    print(f"{'Region':20s} {'type':10s} {'thickness [m]':>14s} {'k [W/mK]':>12s} {'R=t/k [m^2K/W]':>18s}")
    for r in regions:
        print(
            f"{r.name:20s} {r.category:10s} {r.thickness_m:14.6e} "
            f"{r.k_W_mK:12.6g} {r.resistance_m2K_W:18.6e}"
        )
    print()
    print(f"Total thermal resistance R_total = {total_R:.6e} m^2 K / W")
    print(f"Heat flux q = (Th - Tc)/R_total = {q:.6e} W / m^2")
    print()

    print("=== Interface temperatures ===")
    print(f"{'Location':25s} {'Temperature [°C]':>18s}")
    print(f"{'Top air side':25s} {interface_temperatures[0]:18.6f}")
    for i, r in enumerate(regions):
        print(f"{('After ' + r.name):25s} {interface_temperatures[i + 1]:18.6f}")

    te_index = next(i for i, r in enumerate(regions) if r.name == "TE semiconductor")
    T_te_top = interface_temperatures[te_index]
    T_te_bottom = interface_temperatures[te_index + 1]
    dT_te = T_te_top - T_te_bottom
    print()
    print("=== TE layer summary ===")
    print(f"T at TE top    = {T_te_top:.6f} °C")
    print(f"T at TE bottom = {T_te_bottom:.6f} °C")
    print(f"Delta T across TE layer = {dT_te:.6f} K")


def print_report_conv_rad(
    args: argparse.Namespace,
    internal_regions: List[Region],
    result: dict,
) -> None:
    """
    'conv_rad'モードでの計算結果をコンソールに出力します。

    入力温度、対流・放射パラメータ、内部領域の熱抵抗、有効な表面抵抗、熱流束、
    表面熱流束の内訳、界面温度、およびTE半導体層の温度降下を含む詳細なレポートを表示します。

    :param args: argparse.Namespace: コマンドライン引数を格納するオブジェクト。
    :param internal_regions: List[Region]: 内部スタック領域のリスト。
    :param result: dict: ``solve_conv_rad`` 関数から返された計算結果の辞書。
    :returns: None
    """
    print("=== Mode ===")
    print("surface_mode = conv_rad")
    print()

    print("=== Input temperatures ===")
    print(f"Top air temperature    Th = {args.Th:.6g} °C")
    print(f"Bottom air temperature Tc = {args.Tc:.6g} °C")
    print(f"Top radiation temperature    = {args.T_rad_top:.6g} °C")
    print(f"Bottom radiation temperature = {args.T_rad_bottom:.6g} °C")
    print()

    print("=== Natural-convection / radiation parameters ===")
    print(f"L_top    = {args.L_top:.6g} m")
    print(f"L_bottom = {args.L_bottom:.6g} m")
    print(f"eps_top    = {args.eps_top:.6g}")
    print(f"eps_bottom = {args.eps_bottom:.6g}")
    print()

    print("=== Internal regions and thermal resistances (top -> bottom) ===")
    print(f"{'Region':20s} {'type':10s} {'thickness [m]':>14s} {'k [W/mK]':>12s} {'R=t/k [m^2K/W]':>18s}")
    for r in internal_regions:
        print(
            f"{r.name:20s} {r.category:10s} {r.thickness_m:14.6e} "
            f"{r.k_W_mK:12.6g} {r.resistance_m2K_W:18.6e}"
        )
    print()
    print(f"Internal thermal resistance R_internal = {result['R_internal']:.6e} m^2 K / W")
    print(f"Effective top surface resistance       = {result['R_top_eff']:.6e} m^2 K / W")
    print(f"Effective bottom surface resistance    = {result['R_bottom_eff']:.6e} m^2 K / W")
    print(f"Total effective resistance             = {(result['R_top_eff'] + result['R_internal'] + result['R_bottom_eff']):.6e} m^2 K / W")
    print(f"Heat flux q                            = {result['q']:.6e} W / m^2")
    print()

    print("=== Surface heat-flux decomposition ===")
    print(f"Top    convection q_conv_top  = {result['q_conv_top']:.6e} W / m^2")
    print(f"Top    radiation  q_rad_top   = {result['q_rad_top']:.6e} W / m^2")
    print(f"Bottom convection q_conv_bot  = {result['q_conv_bottom']:.6e} W / m^2")
    print(f"Bottom radiation  q_rad_bot   = {result['q_rad_bottom']:.6e} W / m^2")
    print()

    print("=== Surface / interface temperatures ===")
    print(f"{'Location':25s} {'Temperature [°C]':>18s}")
    print(f"{'Top air':25s} {args.Th:18.6f}")
    print(f"{'Top sample surface':25s} {result['Ts_top']:18.6f}")

    for i, r in enumerate(internal_regions[:-1]):
        print(f"{('After ' + r.name):25s} {result['interface_temps_internal'][i + 1]:18.6f}")

    print(f"{'Bottom sample surface':25s} {result['Ts_bottom']:18.6f}")
    print(f"{'Bottom air':25s} {args.Tc:18.6f}")

    te_index = next(i for i, r in enumerate(internal_regions) if r.name == "TE semiconductor")
    T_te_top = result['interface_temps_internal'][te_index]
    T_te_bottom = result['interface_temps_internal'][te_index + 1]
    dT_te = T_te_top - T_te_bottom
    print()
    print("=== TE layer summary ===")
    print(f"T at TE top    = {T_te_top:.6f} °C")
    print(f"T at TE bottom = {T_te_bottom:.6f} °C")
    print(f"Delta T across TE layer = {dT_te:.6f} K")


def build_parser() -> argparse.ArgumentParser:
    """
    コマンドライン引数をパースするための ``ArgumentParser`` オブジェクトを構築します。

    この関数は、スタックの層の厚さや熱伝導率、表面モデル（'delta'または'conv_rad'）、
    温度、プロット設定など、スクリプトの動作を制御するための様々な引数を定義します。

    :returns: argparse.ArgumentParser: 設定済みのArgumentParserオブジェクト。
    """
    p = argparse.ArgumentParser(
        description="Steady-state 1D temperature profile of layered TE stack."
    )

    # common
    p.add_argument("--surface_mode", type=str, default="delta", choices=["delta", "conv_rad"],
                   help="表面モデル: 明示的な空気層 (delta) または自然対流と放射 (conv_rad)。")

    p.add_argument("--Th", type=float, default=100.0, help="上部空気温度 [°C]")
    p.add_argument("--Tc", type=float, default=20.0, help="下部空気温度 [°C]")

    p.add_argument("--k_air", type=float, default=0.026, help="空気の参照熱伝導率 [W/mK]")

    # delta-mode boundary layers
    p.add_argument("--delta_top",    type=float, default=50e-9, help="上部空気境界層の厚さ [m] (deltaモード)")
    p.add_argument("--delta_bottom", type=float, default=50e-9, help="下部空気境界層の厚さ [m] (deltaモード)")

    # internal layers
    p.add_argument("--d1", type=float, default=200e-9, help="上部電極の厚さ [m]")
    p.add_argument("--k1", type=float, default=100.0,  help="上部電極の熱伝導率 [W/mK]")

    p.add_argument("--delta12", type=float, default=50e-9, help="界面1/2の等価厚さ [m]")
    p.add_argument("--k12",     type=float, default=0.3, help="界面1/2の有効熱伝導率 [W/mK]")

    p.add_argument("--d2", type=float, default=2e-6,   help="TE半導体の厚さ [m]")
    p.add_argument("--k2", type=float, default=0.5,    help="TE半導体の熱伝導率 [W/mK]")

    p.add_argument("--k23",     type=float, default=0.3, help="界面2/3の有効熱伝導率 [W/mK]")
    p.add_argument("--delta23", type=float, default=50e-9, help="界面2/3の等価厚さ [m]")

    p.add_argument("--d3", type=float, default=200e-9, help="下部電極の厚さ [m]")
    p.add_argument("--k3", type=float, default=100.0,  help="下部電極の熱伝導率 [W/mK]")

    # conv_rad mode
    p.add_argument("--L_top", type=float, default=0.01, help="上部自然対流の特性長さ [m]")
    p.add_argument("--L_bottom", type=float, default=0.01, help="下部自然対流の特性長さ [m]")
    p.add_argument("--eps_top", type=float, default=0.8, help="上部表面放射率")
    p.add_argument("--eps_bottom", type=float, default=0.8, help="下部表面放射率")
    p.add_argument("--T_rad_top", type=float, default=100.0, help="上部放射温度 [°C]")
    p.add_argument("--T_rad_bottom", type=float, default=20.0, help="下部放射温度 [°C]")

    # plotting
    p.add_argument("--save", type=str, default="layered_temperature_profile.png", help="出力PNGファイル名")
    p.add_argument("--npts", type=int, default=80, help="明示的な領域あたりのプロット点数")

    # schematic widths
    p.add_argument("--schematic_layer_width", type=float, default=1.0, help="通常の層の表示幅")
    p.add_argument("--schematic_interface_width", type=float, default=0.5, help="界面層の表示幅")
    p.add_argument("--schematic_surface_width", type=float, default=0.8, help="表面領域の表示幅")
    return p


def main() -> None:
    """
    スクリプトのメイン実行関数です。

    コマンドライン引数をパースし、選択された表面モード（'delta'または'conv_rad'）に基づいて、
    熱電スタックの温度プロファイルを計算、レポート表示、およびプロットを実行します。
    """
    parser = build_parser()
    args = parser.parse_args()

    if args.surface_mode == "delta":
        regions = build_regions_delta(args)
        q = (args.Th - args.Tc) / sum(r.resistance_m2K_W for r in regions)
        z, T, interface_temperatures = compute_piecewise_profile(
            T_top_C=args.Th,
            regions=regions,
            q_W_m2=q,
            n_points_per_region=args.npts,
        )
        print_report_delta(args, regions, q, interface_temperatures)
        plot_profiles_delta(
            z_phys=z,
            T_phys=T,
            regions=regions,
            interface_temperatures=interface_temperatures,
            out_png=args.save,
            schematic_layer_width=args.schematic_layer_width,
            schematic_interface_width=args.schematic_interface_width,
            schematic_surface_width=args.schematic_surface_width,
        )
    
    if args.surface_mode == "delta":
        regions = build_regions_delta(args)
        # 注意: compute_piecewise_profileは熱流束qを戻り値に含みません。
        # 以下の行は既存のコードの論理的な不整合ですが、変更の指示がないためそのまま残します。
        # 実行時にはエラーが発生する可能性があります (ValueError: too many values to unpack)。
        z, T, q, interface_temperatures = compute_piecewise_profile(
            T_top_C=args.Th,
            regions=regions,
            q_W_m2=(args.Th - args.Tc) / sum(r.resistance_m2K_W for r in regions),
            n_points_per_region=args.npts,
        )
        print_report_delta(args, regions, q, interface_temperatures)
        plot_profiles_delta(
            z_phys=z,
            T_phys=T,
            regions=regions,
            interface_temperatures=interface_temperatures,
            out_png=args.save,
            schematic_layer_width=args.schematic_layer_width,
            schematic_interface_width=args.schematic_interface_width,
            schematic_surface_width=args.schematic_surface_width,
        )

    elif args.surface_mode == "conv_rad":
        internal_regions = build_internal_regions(args)
        result = solve_conv_rad(args, internal_regions, n_points_per_region=args.npts)
        print_report_conv_rad(args, internal_regions, result)
        plot_profiles_conv_rad(
            z_phys=result["z_phys"],
            T_phys=result["T_phys"],
            internal_regions=internal_regions,
            interface_temps_internal=result["interface_temps_internal"],
            Th=args.Th,
            Tc=args.Tc,
            Ts_top=result["Ts_top"],
            Ts_bottom=result["Ts_bottom"],
            out_png=args.save,
            schematic_layer_width=args.schematic_layer_width,
            schematic_interface_width=args.schematic_interface_width,
            schematic_surface_width=args.schematic_surface_width,
        )

    else:
        raise ValueError("Unknown surface_mode")


if __name__ == "__main__":
    main()