draw_rhombohedral_cells.py のソースコードを解析し、Sphinx(MyST)でビルド可能なMarkdownドキュメントを作成します。


draw_rhombohedral_cells.py の解析とドキュメント

概要

このスクリプトは、六方格子の基本ベクトル a_1, a_2, a_3 を定義し、変換行列 T を用いて菱面体格子の基本ベクトル r_1, r_2, r_3 へ変換します。変換後の格子ベクトルおよび単位格子を 3D 空間上に描画し、両者の関係を視覚化します。

非標準ライブラリ

このプログラムは以下の非標準ライブラリを使用しています。

  • numpy

  • matplotlib

環境構築

このプログラムを実行するには、上記の非標準ライブラリをインストールする必要があります。以下のコマンドでインストールできます。

pip install numpy matplotlib

ファイル構成

このドキュメントの対象となるファイルは以下の通りです。

  • draw_rhombohedral_cells.py

スクリプトの実行方法

スクリプトはコマンドラインから直接実行することを想定しています。特別な引数やオプションはありません。

python draw_rhombohedral_cells.py

実行後、Matplotlib の 3D プロットウィンドウが表示されます。

プログラムの仕様

概要

このモジュールは、六方格子(Hexagonal)と菱面体格子(Rhombohedral)の幾何学的関係を可視化します。 六方格子の基本ベクトルを定義し、変換行列を適用して菱面体格子の基本ベクトルを導出します。 両方の格子ベクトルと単位格子を 3D 空間に描画することで、その関係性を視覚的に示します。

詳細説明

このスクリプトは、以下の手順で六方格子と菱面体格子の関係を可視化します。

  1. 六方格子の軸長 a_Hc_H (それぞれ 1.01.6)に基づき、直交座標系における六方格子の基本ベクトル a1, a2, a3 を生成します。

  2. 六方格子から菱面体格子への変換行列 T を定義します。 変換式は以下で表されます。

    R_{rhombo} = T \cdot A_{hex}
    

    ここで、R_{rhombo} は菱面体格子ベクトル、A_{hex} は六方格子ベクトルを表します。

  3. この変換行列 T を六方格子ベクトルに適用し、菱面体格子の基本ベクトル r1, r2, r3 を計算します。

  4. Arrow3D クラスを使用して、六方格子と菱面体格子の基本ベクトルを 3D 空間に矢印として描画します。

  5. draw_unit_cell() 関数を使用して、それぞれの格子の単位格子(平行六面体)を描画します。

  6. 計算された六方格子点と菱面体格子点もプロットされます。

  7. 視覚的な参照として、六方格子の外郭を示す六角柱も描画されます。

  8. 描画された 3D プロットは Matplotlib ウィンドウに表示されます。軸は非表示に設定され、等方的なボックスアスペクト比で表示されます。

グローバル設定

matplotlib.rcParams['font.family']'MS Gothic' に設定されており、日本語フォントの表示に対応しています。

Arrow3D クラス

matplotlib.patches.FancyArrowPatch を継承したクラスで、3D 空間に矢印を描画するための機能を提供します。 主なメソッドは以下の通りです。

  • __init__(self, xs, ys, zs, *args, **kwargs): 矢印の開始点と終了点の 3D 座標を受け取ります。

  • draw(self, renderer): 矢印をレンダリングします。3D 座標を 2D スクリーン座標に変換してから描画します。

  • do_3d_projection(self, renderer=None): 3D 投影変換を実行し、深度情報(zs の最小値)を返します。

draw_unit_cell() 関数

指定された 3 つのベクトル v1, v2, v3 と原点 origin に基づいて、平行六面体(単位格子)を描画します。

  • ax: Matplotlib の 3D Axes オブジェクト。

  • origin: 原点の 3D 座標 (numpy.array)。

  • v1, v2, v3: 格子の基底ベクトル (それぞれ numpy.array)。

  • color: 線の色(デフォルトは 'black')。

  • lw: 線の太さ(デフォルトは 0.5)。

draw_vector() 関数

3D 空間にベクトル(矢印)とそのラベルを描画します。

  • ax: Matplotlib の 3D Axes オブジェクト。

  • vec: 描画するベクトル (numpy.array)。

  • color: 矢印とラベルの色。

  • label: ベクトルのラベル文字列。

  • fontsize: ラベルのフォントサイズ。

main() 関数

格子変換の計算と描画を実行するメインルーチンです。

  1. 六方格子の軸長 a_Hc_H を定義します。

  2. 六方格子ベクトル a1, a2, a3 を直交座標系で計算します。

  3. 六方格子点(0 または 1 の係数で a1, a2, a3 を組み合わせた点)を計算します。

  4. 六方→菱面体変換行列 T を定義します。

  5. TA_hex (六方格子ベクトルを行ベクトルとして並べたもの)に適用し、菱面体格子ベクトル r1, r2, r3 を計算します。

  6. 菱面体格子点(0 または 1 の係数で r1, r2, r3 を組み合わせた点)を計算します。

  7. Matplotlib を使用して 3D プロットを設定します。

  8. draw_vector() 関数で六方格子の基本ベクトル(青)と菱面体格子の基本ベクトル(赤)を描画します。

  9. draw_unit_cell() 関数で六方格子と菱面体格子の単位格子を描画します。

  10. 計算された六方格子点と菱面体格子点をそれぞれ青と赤でプロットします。

  11. 六方格子の外郭を示す灰色の六角柱を描画します。

  12. グラフの軸を非表示にし、ボックスアスペクト比を設定し、カメラの視点(elev=20, azim=30)を設定して、プロットを表示します。

入出力仕様

  • 入力: このプログラムは外部からの入力を受け付けません。必要な軸長や変換行列などのすべての定数は、スクリプトのコード内で直接定義されています。

  • 出力: 実行すると、六方格子と菱面体格子の関係を視覚化した 3D プロットを含む Matplotlib のグラフィックウィンドウが表示されます。このプログラムは画像ファイルとして出力を保存する機能は含まれていません。

主要な変数

  • a_H: 六方格子の水平方向の軸長(コードでは 1.0)。

  • c_H: 六方格子の垂直方向の軸長(コードでは 1.6)。

  • a1, a2, a3: 六方格子の基本ベクトル。numpy.array で表現されます。

  • A_hex: 六方格子の基本ベクトルを行ベクトルとして並べた numpy.array

  • T: 六方格子から菱面体格子への変換行列。numpy.array で定義されます。

  • R_rhombo: 菱面体格子の基本ベクトルを計算した結果。T @ A_hex で得られます。

  • r1, r2, r3: 菱面体格子の基本ベクトル。R_rhombo から抽出されます。

  • hex_points: 六方格子の単位格子を構成する頂点。

  • rhombo_points: 菱面体格子の単位格子を構成する頂点。

コード全体

"""
六方格子(Hexagonal)と菱面体格子(Rhombohedral)の幾何学的関係を可視化するモジュール。

概要:
    このスクリプトは、六方格子の基本ベクトル $a_1, a_2, a_3$ を定義し、
    変換行列 $T$ を用いて菱面体格子の基本ベクトル $r_1, r_2, r_3$ へ変換します。
    変換後の格子ベクトルおよび単位格子を 3D 空間上に描画し、両者の関係を視覚化します。

詳細説明:
    1. 六方格子の軸長 $a_H, c_H$ に基づき、直交座標系での格子ベクトルを生成します。
    2. 六方→菱面体の変換行列 $T$ を適用し、菱面体格子の基底を計算します。
    3. 3D 矢印(Arrow3D クラス)を用いて格子ベクトルを描画します。
    4. 平行六面体(単位格子)を描画し、空間的な重なりを示します。

変換式:
    菱面体格子ベクトルを $R_{rhombo}$、六方格子ベクトルを $A_{hex}$ とすると、
    ```math
    R_{rhombo} = T \cdot A_{hex}
    ```
    で表されます。
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import Axes3D, proj3d

#===================
# グローバル設定
#===================
# フォント設定(日本語化対策)
rcParams['font.family'] = 'MS Gothic'

class Arrow3D(FancyArrowPatch):
    """
    3D空間に矢印を描画するための matplotlib 拡張クラス。
    """
    def __init__(self, xs, ys, zs, *args, **kwargs):
        super().__init__((0, 0), (0, 0), *args, **kwargs)
        self._verts3d = xs, ys, zs

    def draw(self, renderer):
        xs3d, ys3d, zs3d = self._verts3d
        proj = self.axes.get_proj()
        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, proj)
        self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
        super().draw(renderer)

    def do_3d_projection(self, renderer=None):
        xs3d, ys3d, zs3d = self._verts3d
        proj = self.axes.get_proj()
        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, proj)
        self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
        return np.min(zs)

def draw_unit_cell(ax, origin, v1, v2, v3, color='black', lw=0.5):
    """
    指定された3つのベクトルに基づいて平行六面体(単位格子)を描画します。
    """
    O = origin
    A = O + v1
    B = O + v2
    C = O + v3
    D = O + v1 + v2
    E = O + v2 + v3
    F = O + v3 + v1
    G = O + v1 + v2 + v3

    edges = [
        (O, A), (O, B), (O, C),
        (A, D), (A, F),
        (B, D), (B, E),
        (C, E), (C, F),
        (D, G), (E, G), (F, G)
    ]

    for start, end in edges:
        ax.plot([start[0], end[0]],
                [start[1], end[1]],
                [start[2], end[2]],
                color=color, lw=lw)

def draw_vector(ax, vec, color, label, fontsize):
    """
    3D空間にベクトル(矢印)とそのラベルを描画します。
    """
    arrow = Arrow3D([0, vec[0]], [0, vec[1]], [0, vec[2]],
                    mutation_scale=20, lw=2, arrowstyle="-|>", color=color)
    ax.add_artist(arrow)
    ax.text(vec[0] * 1.05, vec[1] * 1.05, vec[2] * 1.05, label, fontsize=fontsize, color=color)

def main():
    """
    格子変換の計算と描画を実行するメインルーチン。
    """
    axis_label_font_size = 24

    # 六方格子の軸長
    a_H = 1.0
    c_H = 1.6

    # 六方格子ベクトル
    a1 = np.array([a_H, 0, 0])
    a2 = np.array([-a_H/2, a_H*np.sqrt(3)/2, 0])
    a3 = np.array([0, 0, c_H])

    # 六方格子点
    hex_points = []
    for i in range(2):
        for j in range(2):
            for k in range(2):
                pt = i*a1 + j*a2 + k*a3
                hex_points.append(pt)
    hex_points = np.array(hex_points)

    # 六方格子ベクトルを行ベクトルとして並べる
    A_hex = np.array([a1, a2, a3])

    # 六方→菱面体変換行列 T
    T = np.array([
        [ 2/3,  1/3,  1/3],
        [-1/3,  1/3,  1/3],
        [-1/3, -2/3,  1/3]
    ])

    # 菱面体格子ベクトルを計算(T × A_hex)
    R_rhombo = T @ A_hex
    r1, r2, r3 = R_rhombo

    # 菱面体格子点
    rhombo_points = []
    for i in range(2):
        for j in range(2):
            for k in range(2):
                pt = i*r1 + j*r2 + k*r3
                rhombo_points.append(pt)
    rhombo_points = np.array(rhombo_points)

    # --- 描画開始 ---
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')

    # 格子ベクトル描画
    draw_vector(ax, a1, 'blue', 'a(H)', axis_label_font_size)
    draw_vector(ax, a2, 'blue', 'b(H)', axis_label_font_size)
    draw_vector(ax, a3, 'blue', 'c(H)', axis_label_font_size)
    draw_vector(ax, r1, 'red',  'a(R)', axis_label_font_size)
    draw_vector(ax, r2, 'red',  'b(R)', axis_label_font_size)
    draw_vector(ax, r3, 'red',  'c(R)', axis_label_font_size)

    # 単位格子描画
    draw_unit_cell(ax, origin=np.array([0,0,0]), v1=a1, v2=a2, v3=a3, color='blue', lw=0.5)
    draw_unit_cell(ax, origin=np.array([0,0,0]), v1=r1, v2=r2, v3=r3, color='red', lw=0.5)

    # 格子点のプロット
    ax.scatter(hex_points[:,0], hex_points[:,1], hex_points[:,2], color='blue', alpha=0.6)
    ax.scatter(rhombo_points[:,0], rhombo_points[:,1], rhombo_points[:,2], color='red', alpha=0.6)

    # 六角柱描画(外郭)
    theta = np.linspace(0, 2*np.pi, 7)
    r = a_H
    z_vals = [0, c_H]
    for z in z_vals:
        x = r * np.cos(theta)
        y = r * np.sin(theta)
        ax.plot(x, y, z, color='gray', alpha=0.5)
    for i in range(6):
        x = [r*np.cos(theta[i]), r*np.cos(theta[i])]
        y = [r*np.sin(theta[i]), r*np.sin(theta[i])]
        z = [0, c_H]
        ax.plot(x, y, z, color='gray', alpha=0.5)

    # グラフ設定
    ax.set_axis_off()
    ax.grid(False)
    ax.set_box_aspect([1,1,1])
    ax.view_init(elev=20, azim=30)
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    main()