lattice_plane.py ダウンロード/コピー

lattice_plane.py をダウンロード

lattice_plane.py
lattice_plane.py
  1"""
  2格子面 (hkl) を3Dで可視化するスクリプト
  3
  4このスクリプトは、指定された格子定数とミラー指数 (hkl) に基づいて、
  5結晶の単位格子とその中の (hkl) 面を3Dで描画します。
  6コマンドライン引数としてミラー指数h, k, lを受け入れ、それらの値に応じて面を可視化します。
  7描画されるのは単位格子に最も近い (hkl) 面の一部です。
  8
  9:doc:`lattice_plane_usage`
 10"""
 11import sys
 12import numpy as np
 13import matplotlib.pyplot as plt
 14from mpl_toolkits.mplot3d.art3d import Poly3DCollection
 15
 16# ─── 格子定数とミラー指数 ───────────────────────
 17a = 1.0  # a 軸長さ
 18b = 1.2  # b 軸長さ
 19c = 0.8  # c 軸長さ
 20alpha, beta, gamma = 90, 100, 110  # α, β, γ （度）
 21h, k, l = 1, 2, 3  # ミラー指数
 22# ───────────────────────────────────────────
 23
 24plane_color = 'cyan'
 25edge_color = 'blue'
 26plane_alpha = 0.1
 27
 28# フォントを MS Gothic に設定
 29plt.rcParams['font.family'] = 'MS Gothic'
 30
 31
 32def init_frac_cell():
 33    """
 34    分数座標系における単位格子の頂点と辺を初期化します。
 35
 36    この関数は、立方晶を想定した分数座標の8つの頂点と、
 37    それらの頂点を結ぶ12本の辺のインデックスのリストを生成します。
 38    これらの情報は、単位格子を表現するために使用されます。
 39
 40    :returns: (fverts, edges)
 41    :rtype: (numpy.ndarray, list[tuple[int, int]])
 42    :returns fverts: 分数座標系における単位格子の8つの頂点の配列。
 43    :returns edges: 単位格子の辺を定義する頂点インデックスのペアのリスト。
 44    """
 45    fverts = np.array([
 46        [0,0,0], [1,0,0], [1,1,0], [0,1,0],
 47        [0,0,1], [1,0,1], [1,1,1], [0,1,1]
 48    ])
 49    edges = [
 50        (0,1),(1,2),(2,3),(3,0),
 51        (4,5),(5,6),(6,7),(7,4),
 52        (0,4),(1,5),(2,6),(3,7)
 53    ]
 54    return fverts, edges
 55
 56
 57def compute_frac_polygon(h: int, k: int, l: int, fverts: np.ndarray, edges: list[tuple[int, int]]) -> np.ndarray | None:
 58    """
 59    ミラー指数 (hkl) によって定義される平面と単位格子の辺との交点を計算し、ポリゴンとして返します。
 60
 61    この関数は、平面の方程式 hx + ky + lz = 1 を使用して、
 62    単位格子を構成する各辺との交点を求めます。
 63    交点が辺の内部 (0 <= t <= 1) にある場合のみ採用し、
 64    それらの交点を多角形の頂点として収集します。
 65    最後に、これらの頂点を平面上の角度に基づいてソートし、
 66    凸多角形を形成するように並べ替えます。
 67
 68    :param h: ミラー指数 h。
 69    :type h: int
 70    :param k: ミラー指数 k。
 71    :type k: int
 72    :param l: ミラー指数 l。
 73    :type l: int
 74    :param fverts: 分数座標系における単位格子の頂点配列。
 75    :type fverts: numpy.ndarray
 76    :param edges: 単位格子の辺を定義する頂点インデックスのペアのリスト。
 77    :type edges: list[tuple[int, int]]
 78    :returns: 計算された交点を角度順に並べたポリゴン頂点の配列。交点が存在しない場合はNone。
 79    :rtype: numpy.ndarray or None
 80    """
 81    pts = []
 82    for i0, i1 in edges:
 83        v0, v1 = fverts[i0], fverts[i1]
 84        d = v1 - v0
 85        denom = h*d[0] + k*d[1] + l*d[2]
 86        if abs(denom) < 1e-8: # 平面と辺が平行な場合
 87            continue
 88        t = (1 - (h*v0[0] + k*v0[1] + l*v0[2])) / denom
 89        if 0 <= t <= 1:
 90            pts.append(v0 + t*d)
 91    if not pts:
 92        return None
 93    
 94    # 浮動小数点数の誤差により重複する点を考慮し、ユニークな点を抽出
 95    # rtolとatolを調整して重複判定の閾値を設定することも可能だが、ここではデフォルトを使用
 96    pts = np.unique(np.array(pts), axis=0)
 97
 98    # 交点をソートしてポリゴンを形成
 99    # 平面内の任意の中心点と、平面内の直交する2つのベクトルu, vを計算
100    cen = pts.mean(axis=0)
101    n = np.array([h, k, l], dtype=float)
102    n /= np.linalg.norm(n) # 法線ベクトル
103
104    # nに直交するベクトルuを計算
105    # nが[1,0,0]に平行でない限り、[1,0,0]との外積が有効
106    # nが[1,0,0]に平行な場合は[0,1,0]との外積を使用
107    if abs(np.dot(n, [1,0,0])) < 0.99: # ほぼ平行でない場合
108        u = np.cross(n, [1,0,0])
109    else: # ほぼ平行な場合
110        u = np.cross(n, [0,1,0])
111    u /= np.linalg.norm(u) # 単位ベクトル化
112    v = np.cross(n, u) # uとnに直交するvを計算
113    
114    # 各点を中心からのベクトルとu, vの内積を使って角度を計算しソート
115    ang = [np.arctan2(np.dot(p-cen, v), np.dot(p-cen, u)) for p in pts]
116    return pts[np.argsort(ang)]
117
118
119def main():
120    """
121    スクリプトのメインエントリポイント。単位格子と指定されたミラー指数面を描画します。
122
123    この関数は、コマンドライン引数からミラー指数 (h, k, l) を読み込み、
124    結晶の格子定数と角度に基づいて実空間の格子ベクトルを計算します。
125    その後、`init_frac_cell` と `compute_frac_polygon` を呼び出して、
126    分数座標系の単位格子頂点と (hkl) 面のポリゴンを取得します。
127    これらの分数座標を実空間座標に変換し、Matplotlibの3Dプロットを用いて、
128    単位格子、(hkl) 面、格子軸ラベル、およびミラーインターセプトの注釈を描画します。
129    最終的に、生成されたプロットを画像ファイルとして保存し、画面に表示します。
130
131    :returns: なし
132    :rtype: None
133    """
134    global h, k, l
135
136    # コマンドライン引数処理
137    argv = sys.argv
138    nargs = len(argv)
139    if nargs > 1: h = int(argv[1])
140    if nargs > 2: k = int(argv[2])
141    if nargs > 3: l = int(argv[3])
142
143    # ラジアン変換
144    alpha_r, beta_r, gamma_r = np.deg2rad(alpha), np.deg2rad(beta), np.deg2rad(gamma)
145    cos_alpha, cos_beta, cos_gamma = np.cos(alpha_r), np.cos(beta_r), np.cos(gamma_r)
146    sin_gamma = np.sin(gamma_r)
147
148    # 格子ベクトル (実座標系)
149    # 結晶学の慣例に従い、a1をx軸に沿わせ、a2をxy平面に配置、a3は空間に配置
150    a1 = np.array([a, 0.0, 0.0])
151    a2 = np.array([b * cos_gamma, b * sin_gamma, 0.0])
152    
153    # a3の成分計算
154    a3_x = c * cos_beta
155    a3_y = c * (cos_alpha - cos_beta * cos_gamma) / sin_gamma
156    # a3_zは三つのベクトルが作る体積が正になるように計算
157    a3_z_squared = 1 - cos_beta**2 - ((cos_alpha - cos_beta * cos_gamma)/sin_gamma)**2
158    a3_z = c * np.sqrt(max(0, a3_z_squared)) # 浮動小数点誤差で負になる可能性があるのでmax(0, ...)
159    a3 = np.array([a3_x, a3_y, a3_z])
160    
161    # 分数座標から実座標への変換行列M
162    M = np.column_stack((a1, a2, a3))
163
164    # 分数格子
165    fverts, edges = init_frac_cell()
166
167    # 平面ポリゴン
168    frac_poly = compute_frac_polygon(h, k, l, fverts, edges)
169    # 分数座標ポリゴンを実座標に変換
170    cart_poly = np.dot(M, frac_poly.T).T if frac_poly is not None else None
171    # 分数座標単位格子頂点を実座標に変換
172    cverts = np.dot(M, fverts.T).T
173
174    # プロット設定
175    fig = plt.figure(figsize=(6,6))
176    ax = fig.add_subplot(projection='3d')
177    ax.set_box_aspect((1,1,1)) # 3Dプロットのアスペクト比を揃える
178    ax.grid(False) # グリッドを非表示に
179    ax.set_axis_off() # 軸を非表示に
180    # 軸ペインの背景色を白に設定
181    ax.xaxis.set_pane_color((1,1,1,1))
182    ax.yaxis.set_pane_color((1,1,1,1))
183    ax.zaxis.set_pane_color((1,1,1,1))
184
185    # 単位格子描画
186    for i0, i1 in edges:
187        xs, ys, zs = zip(cverts[i0], cverts[i1])
188        ax.plot(xs, ys, zs, 'k-', lw=1)
189
190    # (hkl) 面描画
191    if cart_poly is not None and len(cart_poly) >= 3: # ポリゴンが有効かつ3点以上で構成されている場合
192        ax.add_collection3d(
193            Poly3DCollection(
194                [cart_poly], facecolors=plane_color,
195                edgecolors=edge_color, alpha=plane_alpha
196            )
197        )
198
199    # 格子軸ベクトルラベル（少し外側にオフセット）
200    offset_factor = 1.05
201    label_pos_a = a1 * offset_factor
202    label_pos_b = a2 * offset_factor
203    label_pos_c = a3 * offset_factor
204
205    def label_vec(vec: np.ndarray, label: str):
206        """格子軸ラベルをプロットに配置するヘルパー関数。"""
207        ax.text(*vec, label, fontsize=14,
208                fontfamily='MS Gothic', ha='center', va='center')
209
210    label_vec(label_pos_a, 'a')
211    label_vec(label_pos_b, 'b')
212    label_vec(label_pos_c, 'c')
213
214    # ミラーインターセプト注記（オフセット適用）
215    offset_alabel = -0.1
216    offset_blabel = -0.2
217    offset_clabel = -0.05
218
219    if h != 0: # hが0でない場合、a軸との交点が存在
220        pos_h = np.dot(M, np.array([1/h, 0, 0])) # 分数座標[1/h,0,0]を実座標に変換
221        pos_h_offset = [pos_h[0], pos_h[1] + offset_alabel, pos_h[2]]
222        ax.text(*pos_h_offset, f'$1/{{{h}}}$', fontsize=12,
223                fontfamily='MS Gothic', ha='center', va='bottom')
224
225    if k != 0: # kが0でない場合、b軸との交点が存在
226        pos_k = np.dot(M, np.array([0, 1/k, 0])) # 分数座標[0,1/k,0]を実座標に変換
227        pos_k_offset = [pos_k[0] + offset_blabel, pos_k[1], pos_k[2]]
228        ax.text(*pos_k_offset, f'$1/{{{k}}}$', fontsize=12,
229                fontfamily='MS Gothic', ha='right', va='center')
230
231    if l != 0: # lが0でない場合、c軸との交点が存在
232        pos_l = np.dot(M, np.array([0, 0, 1/l])) # 分数座標[0,0,1/l]を実座標に変換
233        pos_l_offset = [pos_l[0], pos_l[1] + offset_clabel, pos_l[2]]
234        ax.text(*pos_l_offset, f'$1/{{{l}}}$', fontsize=12,
235                fontfamily='MS Gothic', ha='right', va='bottom')
236
237    # 視点調整
238    ax.view_init(elev=20, azim=30) # 仰角と方位角を設定
239    ax.dist = 12 # カメラとオブジェクト間の距離を設定
240
241    plt.tight_layout() # レイアウトを自動調整
242    plt.savefig(f"lattice_plane{h}{k}{l}.png",
243                dpi=300, bbox_inches='tight', transparent=True) # 画像として保存
244    plt.show() # プロットを表示
245
246
247# ───────────────────────────────────────────
248# ここから実行
249# ───────────────────────────────────────────
250if __name__ == "__main__":
251    main()