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

draw_rhombohedral_cells.py をダウンロード

draw_rhombohedral_cells.py
draw_rhombohedral_cells.py
  1"""
  2六方格子（Hexagonal）と菱面体格子（Rhombohedral）の幾何学的関係を可視化するモジュール。
  3
  4概要:
  5    このスクリプトは、六方格子の基本ベクトル a1, a2, a3 を定義し、
  6    変換行列 T を用いて菱面体格子の基本ベクトル r1, r2, r3 へ変換します。
  7    変換後の格子ベクトルおよび単位格子を 3D 空間上に描画し、両者の関係を視覚化します。
  8
  9詳細説明:
 10    このスクリプトは以下の手順で処理を実行します。
 11    1. 六方格子の軸長 aH, cH に基づき、直交座標系での格子ベクトルを生成します。
 12    2. 六方→菱面体の変換行列 T を適用し、菱面体格子の基底を計算します。
 13    3. 3D 矢印（Arrow3D クラス）を用いて格子ベクトルを描画します。
 14    4. 平行六面体（単位格子）を描画し、空間的な重なりを示します。
 15
 16変換式:
 17    菱面体格子ベクトルを Rrhombo、六方格子ベクトルを Ahex とすると、
 18    Rrhombo = T @ Ahex で表されます。
 19"""
 20
 21import numpy as np
 22import matplotlib.pyplot as plt
 23from matplotlib import rcParams
 24from matplotlib.patches import FancyArrowPatch
 25from mpl_toolkits.mplot3d import Axes3D, proj3d
 26
 27#===================
 28# グローバル設定
 29#===================
 30# フォント設定（日本語化対策）
 31rcParams['font.family'] = 'MS Gothic'
 32
 33class Arrow3D(FancyArrowPatch):
 34    """
 35    概要:
 36        3D空間に矢印を描画するためのmatplotlib拡張クラスです。
 37
 38    詳細説明:
 39        matplotlib.patches.FancyArrowPatch を継承し、
 40        3D座標を受け取って2D投影された矢印として描画します。
 41        主に draw メソッドと do_3d_projection メソッドが3D描画のためにオーバーライドされています。
 42
 43    引数:
 44        :param xs: 矢印のX座標のリストまたはタプル (開始点X, 終了点X)。
 45        :type xs: list or tuple
 46        :param ys: 矢印のY座標のリストまたはタプル (開始点Y, 終了点Y)。
 47        :type ys: list or tuple
 48        :param zs: 矢印のZ座標のリストまたはタプル (開始点Z, 終了点Z)。
 49        :type zs: list or tuple
 50        :param args: FancyArrowPatch に渡される追加の位置引数。
 51        :param kwargs: FancyArrowPatch に渡される追加のキーワード引数。
 52    """
 53    def __init__(self, xs, ys, zs, *args, **kwargs):
 54        super().__init__((0, 0), (0, 0), *args, **kwargs)
 55        self._verts3d = xs, ys, zs
 56
 57    def draw(self, renderer):
 58        xs3d, ys3d, zs3d = self._verts3d
 59        proj = self.axes.get_proj()
 60        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, proj)
 61        self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
 62        super().draw(renderer)
 63
 64    def do_3d_projection(self, renderer=None):
 65        xs3d, ys3d, zs3d = self._verts3d
 66        proj = self.axes.get_proj()
 67        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, proj)
 68        self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
 69        return np.min(zs)
 70
 71def draw_unit_cell(ax, origin, v1, v2, v3, color='black', lw=0.5):
 72    """
 73    概要:
 74        指定された3つのベクトルに基づいて平行六面体（単位格子）を描画します。
 75
 76    詳細説明:
 77        原点 origin から始まり、3つの基底ベクトル v1, v2, v3 で定義される
 78        平行六面体の全ての辺を描画します。これにより単位格子が可視化されます。
 79
 80    引数:
 81        :param ax: 描画対象の3D Axesオブジェクト。
 82        :type ax: mpl_toolkits.mplot3d.axes3d.Axes3D
 83        :param origin: 格子の原点座標。
 84        :type origin: numpy.ndarray
 85        :param v1: 単位格子を構成する第一のベクトル。
 86        :type v1: numpy.ndarray
 87        :param v2: 単位格子を構成する第二のベクトル。
 88        :type v2: numpy.ndarray
 89        :param v3: 単位格子を構成する第三のベクトル。
 90        :type v3: numpy.ndarray
 91        :param color: 描画する線の色。デフォルトは 'black'。
 92        :type color: str
 93        :param lw: 描画する線の太さ。デフォルトは 0.5。
 94        :type lw: float
 95    """
 96    O = origin
 97    A = O + v1
 98    B = O + v2
 99    C = O + v3
100    D = O + v1 + v2
101    E = O + v2 + v3
102    F = O + v3 + v1
103    G = O + v1 + v2 + v3
104
105    edges = [
106        (O, A), (O, B), (O, C),
107        (A, D), (A, F),
108        (B, D), (B, E),
109        (C, E), (C, F),
110        (D, G), (E, G), (F, G)
111    ]
112
113    for start, end in edges:
114        ax.plot([start[0], end[0]],
115                [start[1], end[1]],
116                [start[2], end[2]],
117                color=color, lw=lw)
118
119def draw_vector(ax, vec, color, label, fontsize):
120    """
121    概要:
122        3D空間にベクトル（矢印）とそのラベルを描画します。
123
124    詳細説明:
125        指定されたベクトルを原点 (0,0,0) から開始する矢印として Axes オブジェクト ax に追加します。
126        矢印の終点付近に、指定されたラベルを配置します。
127
128    引数:
129        :param ax: 描画対象の3D Axesオブジェクト。
130        :type ax: mpl_toolkits.mplot3d.axes3d.Axes3D
131        :param vec: 描画するベクトルの終点座標。原点 (0,0,0) から開始します。
132        :type vec: numpy.ndarray
133        :param color: 矢印とラベルの色。
134        :type color: str
135        :param label: ベクトルに表示するテキストラベル。
136        :type label: str
137        :param fontsize: ラベルのフォントサイズ。
138        :type fontsize: int or float
139    """
140    arrow = Arrow3D([0, vec[0]], [0, vec[1]], [0, vec[2]],
141                    mutation_scale=20, lw=2, arrowstyle="-|>", color=color)
142    ax.add_artist(arrow)
143    ax.text(vec[0] * 1.05, vec[1] * 1.05, vec[2] * 1.05, label, fontsize=fontsize, color=color)
144
145def main():
146    """
147    概要:
148        格子変換の計算と描画を実行するメインルーチンです。
149
150    詳細説明:
151        六方格子の軸長を定義し、六方格子ベクトルと格子点を計算します。
152        その後、六方→菱面体の変換行列を適用して菱面体格子ベクトルと格子点を導出します。
153        最後に、matplotlib を使用してこれらの格子ベクトルと単位格子を3D空間に描画し、
154        六方格子の外郭も表示して、両格子の関係性を視覚的に示します。
155    """
156    axis_label_font_size = 24
157
158    # 六方格子の軸長
159    a_H = 1.0
160    c_H = 1.6
161
162    # 六方格子ベクトル
163    a1 = np.array([a_H, 0, 0])
164    a2 = np.array([-a_H/2, a_H*np.sqrt(3)/2, 0])
165    a3 = np.array([0, 0, c_H])
166
167    # 六方格子点
168    hex_points = []
169    for i in range(2):
170        for j in range(2):
171            for k in range(2):
172                pt = i*a1 + j*a2 + k*a3
173                hex_points.append(pt)
174    hex_points = np.array(hex_points)
175
176    # 六方格子ベクトルを行ベクトルとして並べる
177    A_hex = np.array([a1, a2, a3])
178
179    # 六方→菱面体変換行列 T
180    T = np.array([
181        [ 2/3,  1/3,  1/3],
182        [-1/3,  1/3,  1/3],
183        [-1/3, -2/3,  1/3]
184    ])
185
186    # 菱面体格子ベクトルを計算（T @ A_hex）
187    R_rhombo = T @ A_hex
188    r1, r2, r3 = R_rhombo
189
190    # 菱面体格子点
191    rhombo_points = []
192    for i in range(2):
193        for j in range(2):
194            for k in range(2):
195                pt = i*r1 + j*r2 + k*r3
196                rhombo_points.append(pt)
197    rhombo_points = np.array(rhombo_points)
198
199    # --- 描画開始 ---
200    fig = plt.figure(figsize=(10, 8))
201    ax = fig.add_subplot(111, projection='3d')
202
203    # 格子ベクトル描画
204    draw_vector(ax, a1, 'blue', 'a(H)', axis_label_font_size)
205    draw_vector(ax, a2, 'blue', 'b(H)', axis_label_font_size)
206    draw_vector(ax, a3, 'blue', 'c(H)', axis_label_font_size)
207    draw_vector(ax, r1, 'red',  'a(R)', axis_label_font_size)
208    draw_vector(ax, r2, 'red',  'b(R)', axis_label_font_size)
209    draw_vector(ax, r3, 'red',  'c(R)', axis_label_font_size)
210
211    # 単位格子描画
212    draw_unit_cell(ax, origin=np.array([0,0,0]), v1=a1, v2=a2, v3=a3, color='blue', lw=0.5)
213    draw_unit_cell(ax, origin=np.array([0,0,0]), v1=r1, v2=r2, v3=r3, color='red', lw=0.5)
214
215    # 格子点のプロット
216    ax.scatter(hex_points[:,0], hex_points[:,1], hex_points[:,2], color='blue', alpha=0.6)
217    ax.scatter(rhombo_points[:,0], rhombo_points[:,1], rhombo_points[:,2], color='red', alpha=0.6)
218
219    # 六角柱描画（外郭）
220    theta = np.linspace(0, 2*np.pi, 7)
221    r = a_H
222    z_vals = [0, c_H]
223    for z in z_vals:
224        x = r * np.cos(theta)
225        y = r * np.sin(theta)
226        ax.plot(x, y, z, color='gray', alpha=0.5)
227    for i in range(6):
228        x = [r*np.cos(theta[i]), r*np.cos(theta[i])]
229        y = [r*np.sin(theta[i]), r*np.sin(theta[i])]
230        z = [0, c_H]
231        ax.plot(x, y, z, color='gray', alpha=0.5)
232
233    # グラフ設定
234    ax.set_axis_off()
235    ax.grid(False)
236    ax.set_box_aspect([1,1,1])
237    ax.view_init(elev=20, azim=30)
238    plt.tight_layout()
239    plt.show()
240
241if __name__ == "__main__":
242    main()