#!/usr/bin/env python3
"""
第一ブリルアンゾーンまたは実空間のウィグナー・ザイツセル描画スクリプト。

このスクリプトは、指定された結晶学的パラメータと格子タイプに基づいて、
第一ブリルアンゾーン（逆空間）またはウィグナー・ザイツセル（実空間）を3Dで可視化します。
さらに、周囲の周期的なセルをタイリング表示する機能、高対称点のパス表示（逆空間の場合）、
およびインタラクティブなパラメータ調整機能（スライダーとApplyボタン）を提供します。

:doc:`bz_draw_usage`
"""
import argparse
import numpy as np
import sys, math, functools
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from scipy.spatial import Voronoi
import warnings
warnings.filterwarnings(
  "ignore",
  category=DeprecationWarning,
  module=r"spglib\.spglib",
  message=r"dict interface \(SpglibDataset\['.*'\]\) is deprecated"
)

# ================== 幾何ユーティリティ ==================
def build_conventional_basis(a, b, c, alpha, beta, gamma):
    """
    従来の結晶学的パラメータから直交座標系における格子基底ベクトルを構築します。

    最初の基底ベクトルをx軸に沿わせ、2番目の基底ベクトルをxy平面に置き、
    3番目の基底ベクトルを正のz成分を持つように配置します。

    :param a: 格子定数 a (Å)
    :type a: float
    :param b: 格子定数 b (Å)
    :type b: float
    :param c: 格子定数 c (Å)
    :type c: float
    :param alpha: 角度 alpha (度)
    :type alpha: float
    :param beta: 角度 beta (度)
    :type beta: float
    :param gamma: 角度 gamma (度)
    :type gamma: float
    :returns: 構築された3x3の格子基底行列。各列が基底ベクトル (a_vec, b_vec, c_vec) を表します。
    :rtype: numpy.ndarray
    """
    ar, br, gr = np.radians([alpha, beta, gamma])
    a_vec = np.array([a, 0.0, 0.0])
    b_vec = np.array([b*np.cos(gr), b*np.sin(gr), 0.0])
    cx = c*np.cos(br)
    denom = np.sin(gr) if abs(np.sin(gr)) > 1e-12 else 1e-12
    cy = c*(np.cos(ar) - np.cos(br)*np.cos(gr)) / denom
    cz2 = c**2 - cx**2 - cy**2
    cz = math.sqrt(max(cz2, 0.0))
    c_vec = np.array([cx, cy, cz])
    return np.column_stack((a_vec, b_vec, c_vec))  # 列: a,b,c

def primitive_from_centering(A_conv, lattice_tag):
    """
    従来の格子基底と格子タイプから原始格子基底を計算します。

    与えられた従来の格子基底に対し、指定されたセンタリングタイプ（P, I, Fなど）に
    対応する変換行列を適用して原始格子基底を導出します。

    :param A_conv: 従来の3x3格子基底行列（列が基底ベクトル）
    :type A_conv: numpy.ndarray
    :param lattice_tag: 格子タイプを示す文字列 (例: 'P', 'I', 'F', 'SC', 'BCC', 'FCC', 'A', 'B', 'C')
    :type lattice_tag: str
    :returns: 原始3x3格子基底行列（列が原始基底ベクトル）
    :rtype: numpy.ndarray
    :raises ValueError: 未知の格子タイプまたは未サポートの格子タイプが指定された場合
    """
    tag = lattice_tag.strip().upper()
    alias = {'SC':'P','P':'P','BCC':'I','I':'I','FCC':'F','F':'F','A':'A','B':'B','C':'C'}
    if tag not in alias: raise ValueError(f"Unknown lattice '{lattice_tag}'")
    tag = alias[tag]
    if tag == 'P':
        P = np.eye(3)
    elif tag == 'I':
        P = np.array([[-0.5,  0.5,  0.5],
                      [ 0.5, -0.5,  0.5],
                      [ 0.5,  0.5, -0.5]])
    elif tag == 'F':
        P = np.array([[0.0, 0.5, 0.5],
                      [0.5, 0.0, 0.5],
                      [0.5, 0.5, 0.0]])
    elif tag == 'R':
        P = np.array([[ 2/3, -1/3, -1/3],
                      [ 1/3,  1/3, -2/3],
                      [ 1/3,  1/3,  1/3]])
    elif tag == 'A':
        P = np.array([[1.0,  0.0,  0.0],
                      [0.0,  0.5,  0.5],
                      [0.0, -0.5,  0.5]])
    elif tag == 'B':
        P = np.array([[0.5,  0.0,  0.5],
                      [0.0,  1.0,  0.0],
                      [0.5,  0.0, -0.5]])
    elif tag == 'C':
        P = np.array([[ 0.5,  0.5, 0.0],
                      [ 0.5, -0.5, 0.0],
                      [ 0.0,  0.0, 1.0]])
    else:
        raise ValueError(f"Unsupported lattice '{lattice_tag}'")
    return A_conv @ P  # 列が原始基底

def order_polygon_vertices_on_plane(P):
    """
    平面上にある多角形の頂点を反時計回りに並べ替えるインデックスを返します。

    頂点の重心を計算し、SVD（特異値分解）を用いて平面の法線ベクトルを決定します。
    その後、基準ベクトルと法線に垂直なベクトルを生成し、これらの座標系で各頂点の角度を計算してソートします。
    この順序はMatplotlibで多角形を描画する際に重要です。

    :param P: 多角形の頂点群を表すNx3のNumPy配列。
    :type P: numpy.ndarray
    :returns: 頂点を反時計回りに並べ替えるためのインデックスの配列。
    :rtype: numpy.ndarray
    """
    C = P.mean(axis=0); Q = P - C
    if len(P) <= 2: return np.arange(len(P))
    _, _, Vt = np.linalg.svd(Q, full_matrices=False)
    n = Vt[-1]
    ref = Q[0]
    if np.linalg.norm(np.cross(ref, n)) < 1e-12:
        ref = Q[1] if len(Q)>1 else np.array([1.0,0.0,0.0])
    u = ref - np.dot(ref, n)*n
    if np.linalg.norm(u) < 1e-12:
        tmp = np.array([1.0,0.0,0.0]); 
        if abs(np.dot(tmp,n)) > 0.9: tmp = np.array([0.0,1.0,0.0])
        u = tmp - np.dot(tmp, n)*n
    u /= np.linalg.norm(u)
    v = np.cross(n, u)
    ang = np.arctan2(Q @ v, Q @ u)
    return np.argsort(ang)

def reciprocal_points_from_basis(B, n):
    """
    与えられた基底ベクトルとシェル数に基づいて、格子点または逆格子点の配列を生成します。

    各基底ベクトルの整数倍（-nからnまで）の組み合わせで、中心の点を含む広範囲の点を生成します。
    これはVoronoi図の計算に必要な点群を提供します。

    :param B: 3x3の基底ベクトル行列（列が基底ベクトル）
    :type B: numpy.ndarray
    :param n: 各基底ベクトルに沿って考慮するシェル数（例: n=1なら-1,0,1の係数を使用）
    :type n: int
    :returns: 生成された点群のNx3のNumPy配列
    :rtype: numpy.ndarray
    """
    pts = []
    b1,b2,b3 = B[:,0], B[:,1], B[:,2]
    for i in range(-n,n+1):
        for j in range(-n,n+1):
            for k in range(-n,n+1):
                pts.append(i*b1 + j*b2 + k*b3)
    return np.array(pts)

def get_cell_faces_and_neighbors(vor, point_index):
    """
    Voronoi図の中心点に関連する面の頂点インデックスと隣接点のインデックスを取得します。

    Voronoiオブジェクトのridge_pointsとridge_verticesを走査し、
    指定された点に関連する有限な面（-1を含まない、3頂点以上）を抽出します。

    :param vor: scipy.spatial.Voronoiオブジェクト
    :type vor: scipy.spatial.Voronoi
    :param point_index: 中心となる点の、Voronoi計算に用いた点群におけるインデックス
    :type point_index: int
    :returns: (faces, neighbors) のタプル。
              facesは各面の頂点インデックスのリストのリスト、
              neighborsは各面に対応する隣接点のインデックスのリスト。
    :rtype: tuple[list[numpy.ndarray], list[int]]
    """
    faces, neighbors = [], []
    for (p,q), v_idx in zip(vor.ridge_points, vor.ridge_vertices):
        if point_index not in (p,q): continue
        if -1 in v_idx or len(v_idx) < 3: continue
        faces.append(np.array(v_idx))
        neighbors.append(q if p==point_index else p)
    return faces, neighbors

def cell_polygons(vor, point_index):
    """
    Voronoi図の中心点に関連する多角形（セル面）と隣接点のリストを生成します。

    `get_cell_faces_and_neighbors` を利用して面の頂点インデックスを取得し、
    `order_polygon_vertices_on_plane` で頂点を反時計回りに順序付けて多角形オブジェクトを構築します。

    :param vor: scipy.spatial.Voronoiオブジェクト
    :type vor: scipy.spatial.Voronoi
    :param point_index: 中心となる点の、Voronoi計算に用いた点群におけるインデックス
    :type point_index: int
    :returns: (polygons, neighbors) のタプル。
              polygonsは各面の頂点座標のNx3 NumPy配列のリスト、
              neighborsは各面に対応する隣接点のインデックスのリスト。
    :rtype: tuple[list[numpy.ndarray], list[int]]
    """
    polys = []
    faces, neighbors = get_cell_faces_and_neighbors(vor, point_index)
    for f in faces:
        poly = vor.vertices[f]
        poly = poly[order_polygon_vertices_on_plane(poly)]
        polys.append(poly)
    return polys, neighbors

# ================== SeeK-path（キャッシュ付き） ==================
def _key_from_Aprim(A_prim, symprec):
    """
    原始格子基底行列と対称操作の許容誤差からキャッシュキーを生成します。

    `functools.lru_cache` で使用するためのキーです。
    浮動小数点数の比較を安定させるため、行列の要素は指定された精度で丸められます。

    :param A_prim: 原始3x3格子基底行列（列が基底ベクトル）
    :type A_prim: numpy.ndarray
    :param symprec: 対称操作の許容誤差
    :type symprec: float
    :returns: キャッシュキーとして使用されるタプル
    :rtype: tuple
    """
    # 行が a1,a2,a3 のフラット（数値は丸めてキャッシュ安定化）
    Arows = A_prim.T
    return (float(symprec),) + tuple(np.round(Arows.ravel(), 12))

@functools.lru_cache(maxsize=128)
def _seekpath_cached(key):
    """
    seekpath.get_path を呼び出し、結果をキャッシュする内部関数です。

    `functools.lru_cache` デコレータにより、同じ入力（キャッシュキー）に対しては
    実際に `seekpath` を呼び出さずに、キャッシュされた結果を返します。
    これにより、インタラクティブな描画での再計算を高速化します。

    :param key: `_key_from_Aprim` で生成されたキャッシュキー
    :type key: tuple
    :returns: (Bstd, kpts_abs, path) のタプル。
              Bstdは標準逆格子基底行列、
              kpts_absは高対称点の絶対座標辞書、
              pathは高対称点パスのリスト。
    :rtype: tuple[numpy.ndarray, dict[str, numpy.ndarray], list[tuple[str, str]]]
    """
    symprec = key[0]
    Arows = np.array(key[1:]).reshape(3,3)
    import seekpath
    cell_tuple = (Arows, np.array([[0.0,0.0,0.0]]), [1])
    res = seekpath.get_path(cell_tuple, with_time_reversal=True, symprec=symprec)
    B_rows = np.array(res['reciprocal_primitive_lattice'])
    Bstd = B_rows.T  # 列が b1,b2,b3（2π含む）
    kpts_abs = {lbl: Bstd @ np.array(frac) for lbl, frac in res['point_coords'].items()}
    path = res['path']
    return Bstd, kpts_abs, path

def bz_and_kpath_from_seekpath_cached(A_prim, symprec=1e-7):
    """
    `seekpath` を用いて、第一ブリルアンゾーンの標準逆格子基底、高対称点、およびパスを取得します。

    `_seekpath_cached` を通じて `functools.lru_cache` によるキャッシュを利用するため、
    同じ原始格子基底に対しては `seekpath` の計算がスキップされ、高速に結果を返します。

    :param A_prim: 原始3x3格子基底行列（列が基底ベクトル）
    :type A_prim: numpy.ndarray
    :param symprec: 対称操作の許容誤差
    :type symprec: float
    :returns: (Bstd, kpts_abs, path) のタプル。
              Bstdは標準逆格子基底行列、
              kpts_absは高対称点の絶対座標辞書、
              pathは高対称点パスのリスト。
    :rtype: tuple[numpy.ndarray, dict[str, numpy.ndarray], list[tuple[str, str]]]
    """
    key = _key_from_Aprim(A_prim, symprec)
    return _seekpath_cached(key)

def pretty_label(lbl):
    """
    高対称点のラベルを整形します（例: "GAMMA" を "Γ" に変換）。

    :param lbl: 整形前の高対称点ラベル文字列
    :type lbl: str
    :returns: 整形された高対称点ラベル文字列
    :rtype: str
    """
    return lbl.replace("GAMMA", "Γ")

# ================== 着色＆透明度ユーティリティ ==================
def color_checker(i,j,k):
    """
    三次元インデックス (i,j,k) に基づいてチェッカーボード状の色を返します。

    各インデックスの偶奇性から0から7までのインデックスを生成し、
    定義されたパレットから色を選択します。

    :param i: 第1軸方向のインデックス
    :type i: int
    :param j: 第2軸方向のインデックス
    :type j: int
    :param k: 第3軸方向のインデックス
    :type k: int
    :returns: HEX形式の色コード文字列 (例: '#1f77b4')
    :rtype: str
    """
    palette = ['#1f77b4','#ff7f0e','#2ca02c','#d62728','#9467bd','#8c564b','#e377c2','#7f7f7f']
    idx = ((i & 1) << 0) | ((j & 1) << 1) | ((k & 1) << 2)  # 0..7
    return palette[idx]

def color_hash(i,j,k, cmap=plt.cm.tab20):
    """
    三次元インデックス (i,j,k) に基づいてハッシュ値から色を生成します。

    与えられた整数インデックスから一意のハッシュ値を生成し、
    そのハッシュ値を正規化して指定されたカラーマップから色を抽出します。

    :param i: 第1軸方向のインデックス
    :type i: int
    :param j: 第2軸方向のインデックス
    :type j: int
    :param k: 第3軸方向のインデックス
    :type k: int
    :param cmap: 使用するMatplotlibカラーマップ
    :type cmap: matplotlib.colors.Colormap
    :returns: RGBA形式の色タプル
    :rtype: tuple
    """
    h = (i*73856093 ^ j*19349663 ^ k*83492791) & 0xffffffff
    return cmap((h % 256)/255.0)

def color_distance(t, rmax, cmap=plt.cm.viridis):
    """
    原点からの距離に基づいてグラデーション色を返します。

    点の位置ベクトル `t` の原点からの距離を計算し、`rmax` で正規化して
    指定されたカラーマップ `cmap` から色を抽出します。

    :param t: 点の3D座標ベクトル
    :type t: numpy.ndarray
    :param rmax: 色のグラデーションの最大半径。これを超えると色が飽和します。
                 Noneまたは0以下の場合は常に最大値の色を返します。
    :type rmax: float or None
    :param cmap: 使用するMatplotlibカラーマップ
    :type cmap: matplotlib.colors.Colormap
    :returns: RGBA形式の色タプル
    :rtype: tuple
    """
    r = np.linalg.norm(t); u = 0.0 if rmax is None or rmax<=0 else min(1.0, r/rmax)
    return cmap(u)

def alpha_fade(t, rmax, a_min=0.06, a_max=0.22, power=1.0):
    """
    原点からの距離に基づいて透明度をフェードさせる値を返します。

    点の位置ベクトル `t` の原点からの距離が `rmax` に近づくにつれて、
    透明度が `a_max` から `a_min` へとフェードアウトするように計算されます。
    `power` パラメータでフェードの非線形性を調整できます。

    :param t: 点の3D座標ベクトル
    :type t: numpy.ndarray
    :param rmax: 透明度フェードの最大半径。これを超えると透明度は `a_min` に近づきます。
                 Noneまたは0以下の場合は常に `a_max` を返します。
    :type rmax: float or None
    :param a_min: 最小透明度 (0.0=完全透明, 1.0=完全不透明)
    :type a_min: float
    :param a_max: 最大透明度 (0.0=完全透明, 1.0=完全不透明)
    :type a_max: float
    :param power: フェードの指数（例: 1.0=線形, 2.0=二乗）
    :type power: float
    :returns: 計算された透明度の値
    :rtype: float
    """
    if rmax is None or rmax<=0: return a_max
    r = np.linalg.norm(t); u = min(1.0, r/rmax); w = (1.0 - u)**power
    return a_min + (a_max - a_min)*w

# ================== 描画マネージャ（差分更新） ==================
class SceneManager:
    """
    Matplotlibの3Dシーンを管理し、描画オブジェクト（アーティスト）の追加、クリア、更新を行います。

    特に、描画のたびに既存のオブジェクトを効率的に削除し、新しいオブジェクトを追加することで、
    インタラクティブな描画でのちらつきを抑え、高速な更新を可能にします。
    """
    def __init__(self, args):
        """
        シーンマネージャを初期化します。

        MatplotlibのFigureと3D Axesを作成し、以前の描画オブジェクトを保持するためのリストを準備します。

        :param args: コマンドライン引数を格納するargparse.Namespaceオブジェクト
        :type args: argparse.Namespace
        """
        self.args = args
        self.fig = plt.figure(figsize=(10,9))
        self.ax = self.fig.add_subplot(111, projection='3d')
        self.artists = []  # 前回描画したアーティストを保持（差し替え用）

    def clear_artists(self):
        """
        前回描画したすべてのMatplotlibアーティストをシーンから削除します。

        これにより、Axes自体は維持しつつ、内容だけを更新することができます。
        """
        # Axesは維持し、前回のアーティストだけ remove()
        for a in self.artists:
            try: a.remove()
            except Exception: pass
        self.artists = []

    def add_artist(self, *objs):
        """
        新しいMatplotlibアーティストを管理リストに追加します。

        描画されたアーティストは、後で `clear_artists` メソッドで効率的に削除できるように
        このリストに格納されます。リストやタプルの形式で複数のアーティストを渡すことも可能です。

        :param objs: シーンに追加するMatplotlibアーティストオブジェクト（一つまたは複数）
        :type objs: matplotlib.artist.Artist or list or tuple
        """
        for o in objs:
            if isinstance(o, (list, tuple)):
                self.artists.extend(o)
            elif o is not None:
                self.artists.append(o)

    def draw_scene(self, space, A_prim, Bstd, n_shell,
                   tile_nx, tile_ny, tile_nz, tile_rmax,
                   show_outside_points, highlight_face_points,
                   draw_connectors_to_faces, draw_kpath,
                   color_mode='checker', fade_power=1.0,
                   kpts_abs=None, kpath=None,
                   label_fontsize=10):
        """
        第一ブリルアンゾーンまたはウィグナー・ザイツセル、およびそのタイリングされたビューを描画します。

        空間タイプ（実空間または逆空間）に応じてVoronoi図を構築し、
        中心セル、面生成点、周期タイル、高対称点とパス（逆空間の場合のみ）を描画します。
        既存のアーティストをクリアし、新しいアーティストを追加する差分更新方式を採用しています。

        :param space: 描画する空間タイプ ('reciprocal' または 'real')
        :type space: str
        :param A_prim: 原始3x3格子基底行列（実空間の基底）
        :type A_prim: numpy.ndarray
        :param Bstd: 標準逆格子3x3基底行列（逆空間の基底、2π因子を含む）
        :type Bstd: numpy.ndarray or None
        :param n_shell: Voronoi計算のために考慮する格子点/逆格子点のシェル数
        :type n_shell: int
        :param tile_nx: x軸方向のタイルの最大数（0で非表示）
        :type tile_nx: int
        :param tile_ny: y軸方向のタイルの最大数（0で非表示）
        :type tile_ny: int
        :param tile_nz: z軸方向のタイルの最大数（0で非表示）
        :type tile_nz: int
        :param tile_rmax: タイル描画の半径上限（中心から距離 r <= rmax のセルのみ描画）
                          Noneまたは0以下の場合、距離による制限なし
        :type tile_rmax: float or None
        :param show_outside_points: セル外部の格子点/逆格子点を表示するかどうか
        :type show_outside_points: bool
        :param highlight_face_points: セルの面を定義する点を強調表示するかどうか
        :type highlight_face_points: bool
        :param draw_connectors_to_faces: 面を定義する点と原点を結ぶ破線を描画するかどうか
        :type draw_connectors_to_faces: bool
        :param draw_kpath: 高対称点パスを描画するかどうか（reciprocal空間のみ有効）
        :type draw_kpath: bool
        :param color_mode: タイルの色付けモード ('checker', 'byindex', 'distance')
        :type color_mode: str
        :param fade_power: 距離に応じた透明度フェードの指数
        :type fade_power: float
        :param kpts_abs: 高対称点の絶対座標を格納した辞書 (reciprocal空間の場合のみ)
        :type kpts_abs: dict[str, numpy.ndarray] or None
        :param kpath: 高対称点パスのリスト (reciprocal空間の場合のみ)
        :type kpath: list[tuple[str, str]] or None
        :param label_fontsize: 高対称点ラベルのフォントサイズ
        :type label_fontsize: int
        """
        self.clear_artists()
        ax = self.ax

        if space == 'reciprocal':
            G = reciprocal_points_from_basis(Bstd, n_shell)
            vor = Voronoi(G)
            idx_center = np.argmin(np.linalg.norm(G, axis=1))
            polys, neighbor_idx = cell_polygons(vor, idx_center)
            tile_basis = Bstd
            label_axes = (r'$k_x$', r'$k_y$', r'$k_z$')
            outside_pts = G
            face_color_main = 'lightskyblue'
            origin_label = 'Γ (origin)'
        else:
            R = reciprocal_points_from_basis(A_prim, n_shell)
            vor = Voronoi(R)
            idx_center = np.argmin(np.linalg.norm(R, axis=1))
            polys, neighbor_idx = cell_polygons(vor, idx_center)
            tile_basis = A_prim
            label_axes = ('x','y','z')
            outside_pts = R
            face_color_main = 'peachpuff'
            origin_label = 'origin'

        # 外部点群
        if show_outside_points:
            mask = np.ones(len(outside_pts), dtype=bool); mask[idx_center] = False
            sc = ax.scatter(outside_pts[mask,0], outside_pts[mask,1], outside_pts[mask,2],
                            s=10, color='gray', alpha=0.45,
                            label=('Reciprocal points' if space=='reciprocal' else 'Lattice points'))
            self.add_artist(sc)

        # 中心セル（外周のみ）
        for poly in polys:
            col = Poly3DCollection([poly], alpha=0.32, facecolor=face_color_main, edgecolor='none')
            ax.add_collection3d(col)
            self.add_artist(col)
            ring = np.vstack([poly, poly[0]])
            ln, = ax.plot(ring[:,0], ring[:,1], ring[:,2], color='k', linewidth=1.8)
            self.add_artist(ln)

        # 面生成点と点線
        centers = outside_pts[np.array(neighbor_idx)] if len(neighbor_idx)>0 else np.empty((0,3))
        if highlight_face_points and len(centers)>0:
            sc2 = ax.scatter(centers[:,0], centers[:,1], centers[:,2], s=28, color='tab:blue', alpha=0.9, label='Face-defining points')
            self.add_artist(sc2)
        if draw_connectors_to_faces:
            lines = []
            for p in centers:
                ln, = ax.plot([0,p[0]],[0,p[1]],[0,p[2]], linestyle=':', linewidth=1.2, color='k', alpha=0.9)
                lines.append(ln)
            self.add_artist(lines)

        # 周期タイル
        if any(n>0 for n in (tile_nx, tile_ny, tile_nz)):
            b1,b2,b3 = tile_basis[:,0], tile_basis[:,1], tile_basis[:,2]
            for i in range(-tile_nx, tile_nx+1):
                for j in range(-tile_ny, tile_ny+1):
                    for k in range(-tile_nz, tile_nz+1):
                        if i==0 and j==0 and k==0: continue
                        t = i*b1 + j*b2 + k*b3
                        if tile_rmax is not None and np.linalg.norm(t) > tile_rmax + 1e-12:
                            continue
                        # 色
                        if color_mode == 'checker':
                            fc = color_checker(i,j,k)
                        elif color_mode == 'distance':
                            fc = color_distance(t, tile_rmax)
                        else:
                            fc = color_hash(i,j,k)
                        # 透明度
                        a = alpha_fade(t, tile_rmax, a_min=0.06, a_max=0.22, power=fade_power)
                        # セル描画
                        cols, lns = [], []
                        for poly in polys:
                            P = poly + t
                            c = Poly3DCollection([P], alpha=a, facecolor=fc, edgecolor='none')
                            ax.add_collection3d(c); cols.append(c)
                            ring = np.vstack([P, P[0]])
                            ln, = ax.plot(ring[:,0], ring[:,1], ring[:,2], color='#666666', linewidth=0.6); lns.append(ln)
                        self.add_artist(cols); self.add_artist(lns)

        # 高対称点と k-path（逆格子のみ）
        if space=='reciprocal' and kpts_abs:
            texts, pts, segs = [], [], []
            for lbl, k in kpts_abs.items():
                sct = ax.scatter(k[0],k[1],k[2], s=40, color='crimson'); pts.append(sct)
                txt = ax.text(k[0],k[1],k[2], " "+pretty_label(lbl), fontsize=label_fontsize, color='crimson'); texts.append(txt)
            self.add_artist(pts); self.add_artist(texts)
            if draw_kpath and kpath:
                for a_lbl,b_lbl in kpath:
                    ka,kb = kpts_abs[a_lbl], kpts_abs[b_lbl]
                    seg = np.vstack([ka,kb])
                    ln, = ax.plot(seg[:,0], seg[:,1], seg[:,2], color='crimson', linewidth=1.6, alpha=0.85)
                    segs.append(ln)
                self.add_artist(segs)

        # 原点
        origin = ax.scatter([0],[0],[0], color='k', s=50, label=origin_label)
        self.add_artist(origin)

        # 軸や凡例は1回だけ設定（差分更新では再設定しない）
        basis = (Bstd if space=='reciprocal' else A_prim)
        scale = 1.35 * np.linalg.norm(basis[:,0])
        ax.set_xlim(-scale, scale); ax.set_ylim(-scale, scale); ax.set_zlim(-scale, scale)

        if self.args.plot_limit is not None:
            L = self.args.plot_limit
            ax.set_xlim(-L, L)
            ax.set_ylim(-L, L)
            ax.set_zlim(-L, L)
        ax.set_box_aspect([1,1,1])
        ax.set_xlabel(label_axes[0]); ax.set_ylabel(label_axes[1]); ax.set_zlabel(label_axes[2])
        ax.set_title('First Brillouin Zone (tiled)' if space=='reciprocal' else 'Wigner–Seitz cell (tiled)')
        # 凡例は毎回作り直す（重複回避のため一旦消してから）
        leg = ax.legend(loc='upper right'); self.add_artist(leg)

        self.fig.tight_layout()
        self.fig.canvas.draw_idle()

# ================== 実行系 ==================
def plot_once(args):
    """
    指定されたパラメータに基づいて、第一ブリルアンゾーンまたはウィグナー・ザイツセルを一度だけ描画します。

    結晶学的パラメータから格子基底を構築し、`SceneManager` を使用してシーンを描画し、
    結果をMatplotlibウィンドウに表示します。

    :param args: コマンドライン引数を格納するargparse.Namespaceオブジェクト
    :type args: argparse.Namespace
    :raises SystemExit: `seekpath` の実行に失敗した場合
    """
    A_conv = build_conventional_basis(args.a, args.b, args.c, args.alpha, args.beta, args.gamma)
    A_prim = primitive_from_centering(A_conv, args.lattice)

    mgr = SceneManager(args)
    if args.space == 'reciprocal':
        try:
            Bstd, kpts_abs, kpath = bz_and_kpath_from_seekpath_cached(A_prim)
        except RuntimeError as e:
            print(e, file=sys.stderr); sys.exit(1)
        mgr.draw_scene('reciprocal', A_prim, Bstd, args.n_shell,
                       args.tile_nx, args.tile_ny, args.tile_nz, args.tile_rmax,
                       (not args.no_outside_points), (not args.no_highlight_face_points),
                       (not args.no_connectors), (not args.no_kpath),
                       color_mode=args.tile_color_mode, fade_power=args.fade_power,
                       kpts_abs=kpts_abs, kpath=kpath, label_fontsize=args.label_fontsize)
    else:
        mgr.draw_scene('real', A_prim, None, args.n_shell,
                       args.tile_nx, args.tile_ny, args.tile_nz, args.tile_rmax,
                       (not args.no_outside_points), (not args.no_highlight_face_points),
                       (not args.no_connectors), False,
                       color_mode=args.tile_color_mode, fade_power=args.fade_power,
                       kpts_abs=None, kpath=None, label_fontsize=args.label_fontsize)
    plt.show()

def plot_interactive_apply(args):
    """
    スライダーと「Apply」ボタンを備えたインタラクティブな描画インターフェースを提供します。

    `a, b, c` の格子定数をスライダーで調整でき、「Apply」ボタンをクリックしたときにのみ
    シーンが再描画されます。`seekpath` の結果はキャッシュされるため、
    再描画時の計算コストが削減され、高速な更新が可能です。

    :param args: コマンドライン引数を格納するargparse.Namespaceオブジェクト
    :type args: argparse.Namespace
    """
    # 初期値
    A_conv = build_conventional_basis(args.a, args.b, args.c, args.alpha, args.beta, args.gamma)
    A_prim = primitive_from_centering(A_conv, args.lattice)
    mgr = SceneManager(args)

    def do_draw(a, b, c):
        A_conv = build_conventional_basis(a, b, c, args.alpha, args.beta, args.gamma)
        A_prim = primitive_from_centering(A_conv, args.lattice)
        if args.space == 'reciprocal':
            Bstd, kpts_abs, kpath = bz_and_kpath_from_seekpath_cached(A_prim)  # キャッシュで高速化
            mgr.draw_scene('reciprocal', A_prim, Bstd, args.n_shell,
                           args.tile_nx, args.tile_ny, args.tile_nz, args.tile_rmax,
                           (not args.no_outside_points), (not args.no_highlight_face_points),
                           (not args.no_connectors), (not args.no_kpath),
                           color_mode=args.tile_color_mode, fade_power=args.fade_power,
                           kpts_abs=kpts_abs, kpath=kpath, label_fontsize=args.label_fontsize)
        else:
            mgr.draw_scene('real', A_prim, None, args.n_shell,
                           args.tile_nx, args.tile_ny, args.tile_nz, args.tile_rmax,
                           (not args.no_outside_points), (not args.no_highlight_face_points),
                           (not args.no_connectors), False,
                           color_mode=args.tile_color_mode, fade_power=args.fade_power,
                           kpts_abs=None, kpath=None, label_fontsize=args.label_fontsize)

    # 初回描画
    do_draw(args.a, args.b, args.c)

    # UI: スライダ + Apply ボタン（押したときのみ再描画）
    plt.subplots_adjust(bottom=0.20)
    ax_a = plt.axes([0.15, 0.08, 0.65, 0.03])
    ax_b = plt.axes([0.15, 0.05, 0.65, 0.03])
    ax_c = plt.axes([0.15, 0.02, 0.65, 0.03])
    s_a = Slider(ax_a, 'a', 0.5*args.a, 2.0*args.a, valinit=args.a)
    s_b = Slider(ax_b, 'b', 0.5*args.b, 2.0*args.b, valinit=args.b)
    s_c = Slider(ax_c, 'c', 0.5*args.c, 2.0*args.c, valinit=args.c)

    ax_btn = plt.axes([0.83, 0.02, 0.12, 0.09])
    btn = Button(ax_btn, 'Apply')

    busy = {'flag': False}
    def on_apply(event):
        if busy['flag']:  # 処理中は無視
            return
        busy['flag'] = True
        # ボタン表示を Redrawing... に変更
        btn.label.set_text('Redrawing...')
        btn.ax.set_facecolor('#dddddd')
        mgr.fig.canvas.draw_idle()

        try:
            do_draw(s_a.val, s_b.val, s_c.val)
        finally:
            # ボタン表示を戻す
            btn.label.set_text('Apply')
            btn.ax.set_facecolor('#f0f0f0')
            mgr.fig.canvas.draw_idle()
            busy['flag'] = False

    btn.on_clicked(on_apply)
    plt.show()

# ================== CLI ==================
def main():
    """
    プログラムのエントリーポイントです。

    コマンドライン引数を解析し、ユーザーの選択に応じて
    第一ブリルアンゾーンまたはウィグナー・ザイツセルの描画機能
    （一度きりの描画またはインタラクティブな描画）を呼び出します。
    """
    ap = argparse.ArgumentParser(description="Voronoi 多面体（第一BZ / Wigner–Seitz）: キャッシュ・差分更新・Applyボタン対応")
    ap.add_argument('--space', choices=['reciprocal','real'], default='reciprocal',
                    help="reciprocal: 第一BZ（SeeK-pathで高対称点） / real: 実空間 Wigner–Seitz")
    ap.add_argument('--lattice', default='P', help="格子タイプ: P|SC, I|BCC, F|FCC, A, B, C")
    ap.add_argument('--a', type=float, default=1.0); ap.add_argument('--b', type=float, default=1.0); ap.add_argument('--c', type=float, default=1.0)
    ap.add_argument('--alpha', type=float, default=90.0); ap.add_argument('--beta', type=float, default=90.0); ap.add_argument('--gamma', type=float, default=90.0)
    ap.add_argument('--n-shell', type=int, default=3, help='Voronoi 計算用に並べる格子点の範囲')
    ap.add_argument('--tile-nx', type=int, default=0); ap.add_argument('--tile-ny', type=int, default=0); ap.add_argument('--tile-nz', type=int, default=0)
    ap.add_argument('--tile-rmax', type=float, default=None,
                    help='タイル描画の半径上限（中心から距離 r ≤ rmax のセルのみ）')
    ap.add_argument('--plot-limit', type=float, default=None,
                    help='軸の表示範囲を [-L, L] に強制的に制限する（内部座標、L=1.0で[-2pi/a, 2pi/a]に対応）')
    ap.add_argument('--tile-color-mode', choices=['checker','byindex','distance'], default='checker',
                    help="タイル色: checker=チェッカー, byindex=ハッシュ, distance=距離グラデーション")
    ap.add_argument('--fade-power', type=float, default=1.0, help='距離フェードの指数（1=線形, 2=二乗など）')
    ap.add_argument('--no-outside-points', action='store_true')
    ap.add_argument('--no-highlight-face-points', action='store_true')
    ap.add_argument('--no-connectors', action='store_true')
    ap.add_argument('--no-kpath', action='store_true')
    ap.add_argument('--label-fontsize', type=int, default=10)
    ap.add_argument('--interactive-apply', action='store_true', help='スライダ + Applyボタンで再描画')
    ap.add_argument('--save', default=None, help='（一括描画時のみ）画像保存パス')

    args = ap.parse_args()
    if args.interactive_apply:
        plot_interactive_apply(args)
    else:
        plot_once(args)

if __name__ == "__main__":
    main()