draw_unit_cell.py ダウンロード/コピー
draw_unit_cell.py
draw_unit_cell.py
1"""
2単位格子を描画するスクリプト。
3
4概要:
5 結晶学における単位格子を3Dで可視化する機能を提供します。
6
7詳細説明:
8 任意の格子定数と角度から単位格子の形状を計算し、Matplotlibを使用して描画します。
9 基本ベクトル（a, b, c）の描画には、パースペクティブを考慮したカスタム3D矢印クラスを使用します。
10
11関連リンク:
12 draw_unit_cell_usage
13"""
14
15import numpy as np
16import matplotlib.pyplot as plt
17from mpl_toolkits.mplot3d import Axes3D
18from matplotlib.patches import FancyArrowPatch
19from mpl_toolkits.mplot3d import Axes3D, proj3d
20
21class Arrow3D(FancyArrowPatch):
22 """
23 概要:
24 Matplotlibの3Dプロットでカスタムの3D矢印を描画するためのクラス。
25
26 詳細説明:
27 FancyArrowPatch を継承し、3D空間内の座標を2Dスクリーン座標に変換して描画することで、
28 常に適切なパースペクティブで矢印が表示されるようにします。
29 """
30 def __init__(self, xs, ys, zs, *args, **kwargs):
31 """
32 概要:
33 Arrow3Dクラスのコンストラクタ。
34
35 引数:
36 :param xs: 矢印のX座標のリストまたはタプル。
37 :type xs: list[float]またはtuple[float]
38 :param ys: 矢印のY座標のリストまたはタプル。
39 :type ys: list[float]またはtuple[float]
40 :param zs: 矢印のZ座標のリストまたはタプル。
41 :type zs: list[float]またはtuple[float]
42 :param args: FancyArrowPatch に渡される追加の引数。
43 :type args: tuple
44 :param kwargs: FancyArrowPatch に渡される追加のキーワード引数。
45 :type kwargs: dict
46 """
47 super().__init__((0, 0), (0, 0), *args, **kwargs)
48 self._verts3d = xs, ys, zs
49
50 def draw(self, renderer):
51 """
52 概要:
53 矢印をレンダリングします。
54
55 詳細説明:
56 3D座標を2Dスクリーン座標に変換し、変換された座標を使用して矢印を描画します。
57 引数:
58 :param renderer: Matplotlibのレンダラーオブジェクト。
59 :type renderer: matplotlib.backend_bases.RendererBase
60 戻り値:
61 :returns: なし
62 :rtype: None
63 """
64 xs3d, ys3d, zs3d = self._verts3d
65 proj = self.axes.get_proj()
66 xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, proj)
67 self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
68 super().draw(renderer)
69
70 def do_3d_projection(self, renderer=None):
71 """
72 概要:
73 3Dプロジェクションを実行し、矢印の深度を計算します。
74
75 詳細説明:
76 3D座標を2Dスクリーン座標に変換し、Z軸方向の最小値（深度）を返すことで、
77 3Dビューにおけるオブジェクトの描画順序を適切に管理します。
78 引数:
79 :param renderer: Matplotlibのレンダラーオブジェクト（オプション）。
80 :type renderer: matplotlib.backend_bases.RendererBase
81 戻り値:
82 :returns: 矢印のZ軸方向の最小値（深度）。
83 :rtype: numpy.float
84 """
85 xs3d, ys3d, zs3d = self._verts3d
86 proj = self.axes.get_proj()
87 xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, proj)
88 self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
89 return np.min(zs)
90
91def draw_vector(ax, vec, color, label, fontsize):
92 """
93 概要:
94 3D空間にベクトルと対応するラベルを描画します。
95
96 詳細説明:
97 Arrow3D クラスを利用して、原点から指定されたベクトルまでの矢印を描画し、
98 ベクトルの先端にラベルを配置します。
99 引数:
100 :param ax: 3DプロットのAxesオブジェクト。
101 :type ax: matplotlib.axes.Axes
102 :param vec: 描画するベクトルの3D座標（[x, y, z]）。
103 :type vec: numpy.ndarray
104 :param color: 矢印とラベルの色。
105 :type color: str
106 :param label: ベクトルに付けるラベル。
107 :type label: str
108 :param fontsize: ラベルのフォントサイズ。
109 :type fontsize: int
110 戻り値:
111 :returns: なし
112 :rtype: None
113 """
114 arrow = Arrow3D([0, vec[0]], [0, vec[1]], [0, vec[2]],
115 mutation_scale=20, lw=2, arrowstyle="-|>", color=color)
116 ax.add_artist(arrow)
117 ax.text(vec[0] * 1.05, vec[1] * 1.05, vec[2] * 1.05, label, fontsize = fontsize, color = color)
118
119def lattice_vectors(a, b, c, alpha, beta, gamma):
120 """
121 概要:
122 格子定数と角度から基本格子ベクトルを計算します。
123
124 詳細説明:
125 結晶学で用いられる格子定数 (a, b, c) と軸間角 (alpha, beta, gamma) を用いて、
126 直交座標系における3つの基本格子ベクトル (va, vb, vc) を計算します。
127 引数:
128 :param a: 格子定数aの長さ。
129 :type a: float
130 :param b: 格子定数bの長さ。
131 :type b: float
132 :param c: 格子定数cの長さ。
133 :type c: float
134 :param alpha: b軸とc軸の間の角度（度数）。
135 :type alpha: float
136 :param beta: a軸とc軸の間の角度（度数）。
137 :type beta: float
138 :param gamma: a軸とb軸の間の角度（度数）。
139 :type gamma: float
140 戻り値:
141 :returns: 3つの基本格子ベクトル (va, vb, vc)。各ベクトルは [x, y, z] の形式の numpy.ndarray。
142 :rtype: tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]
143 """
144 alpha_r, beta_r, gamma_r = np.radians([alpha, beta, gamma])
145 va = np.array([a, 0, 0])
146 vb = np.array([b * np.cos(gamma_r), b * np.sin(gamma_r), 0])
147 cx = c * np.cos(beta_r)
148 cy = c * (np.cos(alpha_r) - np.cos(beta_r) * np.cos(gamma_r)) / np.sin(gamma_r)
149 cz = np.sqrt(max(0, c**2 - cx**2 - cy**2)) # 負の平方根を避けるための安全策
150 vc = np.array([cx, cy, cz])
151 return va, vb, vc
152
153def set_equal_aspect(ax, points):
154 """
155 概要:
156 3Dプロットのアスペクト比を均等に設定し、描画範囲を調整します。
157
158 詳細説明:
159 描画される点の最大・最小座標に基づいて、X, Y, Z軸の表示範囲を同じにします。
160 引数:
161 :param ax: 3DプロットのAxesオブジェクト。
162 :type ax: matplotlib.axes.Axes
163 :param points: 描画される全ての点のNumpy配列。形状は (N, 3)。
164 :type points: numpy.ndarray
165 戻り値:
166 :returns: なし
167 :rtype: None
168 """
169 xlim = [np.min(points[:,0]), np.max(points[:,0])]
170 ylim = [np.min(points[:,1]), np.max(points[:,1])]
171 zlim = [np.min(points[:,2]), np.max(points[:,2])]
172 max_range = max(
173 xlim[1] - xlim[0],
174 ylim[1] - ylim[0],
175 zlim[1] - zlim[0]
176 ) / 2
177
178 mid_x = np.mean(xlim)
179 mid_y = np.mean(ylim)
180 mid_z = np.mean(zlim)
181
182 ax.set_xlim(mid_x - max_range, mid_x + max_range)
183 ax.set_ylim(mid_y - max_range, mid_y + max_range)
184 ax.set_zlim(mid_z - max_range, mid_z + max_range)
185 ax.set_box_aspect([1,1,1])
186
187def draw_unit_cell_plot(a, b, c, alpha, beta, gamma, elev=5.90, azim=-67.55, fontsize=24, color='blue'):
188 """
189 概要:
190 指定された格子定数と角度を持つ単位格子を3Dで描画します。
191
192 詳細説明:
193 lattice_vectors 関数を使用して基本格子ベクトルを計算し、それらから単位格子の8つの頂点を生成します。
194 その後、Matplotlibを用いて単位格子の辺、頂点、および基本格子ベクトルを3D空間に描画します。
195 引数:
196 :param a: 格子定数aの長さ。
197 :type a: float
198 :param b: 格子定数bの長さ。
199 :type b: float
200 :param c: 格子定数cの長さ。
201 :type c: float
202 :param alpha: b軸とc軸の間の角度（度数）。
203 :type alpha: float
204 :param beta: a軸とc軸の間の角度（度数）。
205 :type beta: float
206 :param gamma: a軸とb軸の間の角度（度数）。
207 :type gamma: float
208 :param elev: 3Dビューの仰角（度数）。デフォルトは5.90。
209 :type elev: float
210 :param azim: 3Dビューの方位角（度数）。デフォルトは-67.55。
211 :type azim: float
212 :param fontsize: ラベルのフォントサイズ。デフォルトは24。
213 :type fontsize: int
214 :param color: ベクトルの色。デフォルトは 'blue'。
215 :type color: str
216 戻り値:
217 :returns: なし（グラフを表示し、PNGとして保存します）
218 :rtype: None
219 """
220 va, vb, vc = lattice_vectors(a, b, c, alpha, beta, gamma)
221 origin = np.array([0, 0, 0])
222 points = [
223 origin,
224 va,
225 vb,
226 vc,
227 va + vb,
228 va + vc,
229 vb + vc,
230 va + vb + vc
231 ]
232 points = np.array(points)
233
234 edges = [
235 (0,1), (0,2), (0,3),
236 (1,4), (1,5),
237 (2,4), (2,6),
238 (3,5), (3,6),
239 (4,7), (5,7), (6,7)
240 ]
241
242 fig = plt.figure(figsize=(8,6))
243 ax = fig.add_subplot(111, projection='3d')
244 ax.view_init(elev=elev, azim=azim)
245
246 # インタラクティブな角度確認用のハンドラ（直接実行時のみ意味を持つ）
247 def on_draw(event):
248 curr_elev = ax.elev
249 curr_azim = ax.azim
250 print(f"elev: {curr_elev:.2f}, azim: {curr_azim:.2f}")
251 fig.canvas.mpl_connect('draw_event', on_draw)
252
253 # 単位格子の辺（細い黒線）
254 for i,j in edges:
255 ax.plot(*zip(points[i], points[j]), color='black', linewidth=0.5)
256
257 # 頂点（黒い小さな●）
258 ax.scatter(points[:,0], points[:,1], points[:,2], color='black', s=20)
259
260 # 基本ベクトル（太い矢印）とラベル
261 for vec, label in zip([va, vb, vc], ['a', 'b', 'c']):
262 draw_vector(ax, vec, color, label, fontsize)
263
264 # グラフの見た目を整える
265 ax.set_axis_off()
266 ax.grid(False)
267 ax.xaxis.pane.set_visible(False)
268 ax.yaxis.pane.set_visible(False)
269 ax.zaxis.pane.set_visible(False)
270 ax.set_facecolor((1,1,1,0)) # 背景透明
271
272 set_equal_aspect(ax, points)
273 plt.tight_layout()
274 plt.savefig("unit_cell.png", dpi=300, bbox_inches='tight', transparent=True)
275 plt.show()
276
277def main():
278 """
279 概要:
280 メインの実行ルーチン。特定の格子定数で単位格子を描画します。
281 """
282 # 始点の初期値やフォントサイズなどの設定
283 params = {
284 'a': 5.0, 'b': 5.5, 'c': 4.5,
285 'alpha': 80, 'beta': 70, 'gamma': 100,
286 'elev': 5.90, 'azim': -67.55,
287 'fontsize': 24, 'color': 'blue'
288 }
289
290 print(f"Drawing unit cell with a={params['a']}, b={params['b']}, c={params['c']}")
291 print(f"Angles: alpha={params['alpha']}, beta={params['beta']}, gamma={params['gamma']}")
292
293 draw_unit_cell_plot(**params)
294
295if __name__ == "__main__":
296 main()