"""
tkcrystalbase.py

概要: 結晶構造計算と描画に関連するユーティリティ関数を提供するモジュール。
詳細説明:
    このモジュールは、結晶の格子定数から実空間および逆空間の格子ベクトル、
    メトリックテンソル、単位格子の体積などを計算する機能を提供します。
    また、原子の分数座標とデカルト座標間の変換、原子間距離や角度の計算、
    そしてmatplotlibを使用した3D描画ヘルパー関数を含んでいます。
    物理定数や一般的な格子パラメータ、サイト情報もここで定義されています。

関連リンク:
    :doc:`tkcrystalbase_usage`
"""

from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, sqrt
import numpy as np
from numpy import linalg as la
from pprint import pprint

pi          = 3.14159265358979323846
pi2         = pi + pi
torad       = 0.01745329251944 # rad/deg";
todeg       = 57.29577951472   # deg/rad";
basee       = 2.71828183

h           = 6.6260755e-34    # Js";
h_bar       = 1.05459e-34      # "Js";
hbar        = h_bar
c           = 2.99792458e8     # m/s";
e           = 1.60218e-19      # C";
me          = 9.1093897e-31    # kg";
mp          = 1.6726231e-27    # kg";
mn          = 1.67495e-27      # kg";
u0          = 4.0 * 3.14*1e-7; # . "Ns<sup>2</sup>C<sup>-2</sup>";
e0          = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
e2_4pie0    = 2.30711e-28      # Nm<sup>2</sup>";
a0          = 5.29177e-11      # m";
kB          = 1.380658e-23     # JK<sup>-1</sup>";
NA          = 6.0221367e23     # mol<sup>-1</sup>";
R           = 8.31451          # J/K/mol";
F           = 96485.3          # C/mol";
g           = 9.81             # m/s2";


# Lattice parameters (angstrom and degree)
#lattice_parameters = [ 5.62, 5.62, 5.62, 60.0, 60.0, 60.0]
lattice_parameters = [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0]

# Site information (atom name, site label, atomic number, atomic mass, charge, radius, color, position)
sites = [
         ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.0, 0.0, 0.0])]
        ,['Na', 'Na2', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.0, 0.5, 0.5])]
        ,['Na', 'Na3', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.5, 0.0, 0.5])]
        ,['Na', 'Na4', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.5, 0.5, 0.0])]
        ,['Cl', 'Cl1', 17, 35.4527,  +1.0, 1.4, 'blue', np.array([0.5, 0.0, 0.0])]
        ,['Cl', 'Cl2', 17, 35.4527,  +1.0, 1.4, 'blue', np.array([0.5, 0.5, 0.5])]
        ,['Cl', 'Cl3', 17, 35.4527,  +1.0, 1.4, 'blue', np.array([0.0, 0.0, 0.5])]
        ,['Cl', 'Cl4', 17, 35.4527,  +1.0, 1.4, 'blue', np.array([0.0, 0.5, 0.0])]
        ]


def reduce01(x):
    """
    概要: 数値を0以上1未満の範囲に変換する。
    詳細説明:
        入力された数値の小数部分を返します。
        例えば、3.75は0.75に、-2.3は0.7に変換されます。
    :param x: float: 変換する数値。
    :returns: float: 0以上1未満の範囲に変換された値。
    """
    return x - int(x)

    return x

def round01(x):
    """
    概要: 数値を特定の閾値に基づいて0または1に丸めるか、0以上1未満の範囲に変換する。
    詳細説明:
        `x` が1.0に近い（絶対差が0.0002未満）場合は1.0を返します。
        `x` が0.0に近い（絶対差が0.0002未満）場合は0.0を返します。
        それ以外の場合、`x > 1.0` なら `x - int(x)` を返します。
        `x < 1.0` なら `x - int(x) + 1.0` を返します。
        この関数は、負の数や1を超える数値を0以上1未満の範囲に正規化するために使用されることがあります。
    :param x: float: 丸める、または変換する数値。
    :returns: float: 丸められた、または0以上1未満の範囲に変換された値。
    """
    if abs(x - 1.0) < 0.0002:
        return 1.0
    if abs(x) < 0.0002:
        return 0.0
    if x > 1.0:
        return x - int(x)
    if x < 1.0:
        return x - int(x) + 1.0
    return x

def round_parameter(x, tol):
    """
    概要: 数値を指定された許容範囲 (tolerance) の整数倍に丸める。
    詳細説明:
        入力値 `x` を `tol` の整数倍に丸めます。
        丸め処理は、`(x + 0.1 * tol) / tol` を整数に変換し、`tol` を掛けることで行われます。
        これにより、例えば `tol=0.5` で `x=1.2` なら `1.0` に、`x=1.3` なら `1.5` に丸められます。
    :param x: float: 丸める数値。
    :param tol: float: 丸めの基準となる許容範囲（単位）。
    :returns: float: `tol` の整数倍に丸められた値。
    """
    val = tol * int( (x+0.1*tol) / tol )
    return val 


def cal_lattice_vectors(latt):
    """
    概要: 単位格子の格子定数から実空間の格子ベクトルを計算する。
    詳細説明:
        格子定数 (a, b, c, α, β, γ) を受け取り、
        直交座標系における3つの格子ベクトル a_vec, b_vec, c_vec を含む3x3行列を計算します。
        格子ベクトルは、x軸に沿って a_vec が配置されるように構成されます。
        角度は度数で指定されますが、内部ではラジアンに変換して計算されます。
    :param latt: list of float: 格子定数のリスト `[a, b, c, alpha, beta, gamma]`。
                                長さはオングストローム、角度は度で与えられます。
    :returns: numpy.ndarray: 3x3 の格子ベクトル行列。
                            `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
    """
    cosa = cos(torad * latt[3])
    cosb = cos(torad * latt[4])
    cosg = cos(torad * latt[5])
    sing = sin(torad * latt[5])

    aij = np.empty([3, 3], dtype = float)
    aij[0][0] = latt[0]
    aij[0][1] = 0.0;
    aij[0][2] = 0.0;
    aij[1][0] = latt[1] * cosg
    aij[1][1] = latt[1] * sing
    aij[1][2] = 0.0;
    aij[2][0] = latt[2] * cosb
    aij[2][1] = latt[2] * (cosa - cosb * cosg) / sing
    if abs(aij[2][1]) < 1.0e-8:
        aij[2][1] = 0.0
    else:
        aij[2][1] = aij[2][1] / sing
    aij[2][2] = sqrt(latt[2] * latt[2] - aij[2][0] * aij[2][0] - aij[2][1] * aij[2][1])

    return aij

def cal_metrics(latt):
    """
    概要: 格子定数から格子ベクトルとメトリックテンソルを計算する。
    詳細説明:
        単位格子の格子定数から実空間の格子ベクトル `aij` を計算し、
        それらの内積からメトリックテンソル `gij` を構築します。
        結果は辞書形式で返されます。
    :param latt: list of float: 格子定数のリスト `[a, b, c, alpha, beta, gamma]`。
                                長さはオングストローム、角度は度で与えられます。
    :returns: dict: 計算された格子ベクトル行列 (`'aij'`) とメトリックテンソル (`'gij'`)
                    をキーに持つ辞書。
                    `{'aij': numpy.ndarray, 'gij': numpy.ndarray}` の形式。
    """
    inf = {}

    cosa = cos(torad * latt[3])
    cosb = cos(torad * latt[4])
    cosg = cos(torad * latt[5])

    aij = cal_lattice_vectors(latt)
    inf['aij'] = aij

    gij = np.empty([3, 3], dtype = float)
    for i in range(3):
        for j in range(i, 3):
            gij[i][j] = np.dot(aij[i], aij[j])
            gij[j][i] = gij[i][j]
    inf['gij'] = gij
    
    return inf

def cal_volume(aij):
    """
    概要: 格子ベクトルから単位格子の体積を計算する。
    詳細説明:
        3つの格子ベクトル (a_vec, b_vec, c_vec) のスカラー三重積 (a_vec・(b_vec × c_vec)) を用いて、
        単位格子の体積を計算します。
    :param aij: numpy.ndarray: 3x3 の格子ベクトル行列。
                               `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
    :returns: float: 単位格子の体積。
    """
    axb = np.cross(aij[0], aij[1])      # Outner product
    vol = np.dot(axb, aij[2])           # Inner product
    return vol

def cal_reciprocal_lattice_vectors(aij):
    """
    概要: 実空間格子ベクトルから逆空間格子ベクトルを計算する。
    詳細説明:
        実空間の格子ベクトル `aij` (a_vec, b_vec, c_vec) を用いて、
        逆空間の格子ベクトル (a*_vec, b*_vec, c*_vec) を計算します。
        逆空間ベクトルは、単位格子の体積 V と実空間ベクトルの外積によって定義されます。
    :param aij: numpy.ndarray: 3x3 の実空間格子ベクトル行列。
                               `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
    :returns: list of numpy.ndarray: 逆空間の3つの格子ベクトル `[Ra_vec, Rb_vec, Rc_vec]` のリスト。
    """
    V = cal_volume(aij)
    Ra = np.cross(aij[1], aij[2]) / V
    Rb = np.cross(aij[2], aij[0]) / V
    Rc = np.cross(aij[0], aij[1]) / V

    return [Ra, Rb, Rc]

def cal_reciprocal_lattice_parameters(Raij):
    """
    概要: 逆空間格子ベクトルから逆空間の格子定数を計算する。
    詳細説明:
        逆空間の格子ベクトル `Raij` (a*_vec, b*_vec, c*_vec) を用いて、
        逆空間の格子定数 (a*, b*, c*, α*, β*, γ*) を計算します。
        長さはベクトルのノルムから、角度はベクトルの内積から導出されます。
    :param Raij: list of numpy.ndarray: 逆空間の3つの格子ベクトル `[Ra_vec, Rb_vec, Rc_vec]` のリスト。
    :returns: list of float: 逆空間の格子定数 `[Ra, Rb, Rc, Ralpha, Rbeta, Rgamma]` のリスト。
    """
    Ra = la.norm(Raij[0])
    Rb = la.norm(Raij[1])
    Rc = la.norm(Raij[2])
    Ralpha = todeg * arccos(np.dot(Raij[1], Raij[2]) / Rb / Rc)
    Rbeta  = todeg * arccos(np.dot(Raij[2], Raij[0]) / Rc / Ra)
    Rgamma = todeg * arccos(np.dot(Raij[0], Raij[1]) / Ra / Rb)

    return [Ra, Rb, Rc, Ralpha, Rbeta, Rgamma]

def fractional_to_cartesian(pos, aij):
    """
    概要: 分数座標からデカルト座標に変換する。
    詳細説明:
        単位格子内の原子の分数座標 (x_frac, y_frac, z_frac) と格子ベクトル `aij` を用いて、
        原子のデカルト座標 (x_cart, y_cart, z_cart) を計算します。
    :param pos: numpy.ndarray or list of float: 変換する分数座標 `[x_frac, y_frac, z_frac]`。
    :param aij: numpy.ndarray: 3x3 の格子ベクトル行列。
                               `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
    :returns: tuple of float: 計算されたデカルト座標 `(x_cart, y_cart, z_cart)`。
    """
    x = pos[0] * aij[0][0] + pos[1] * aij[1][0] + pos[2] * aij[2][0]
    y = pos[0] * aij[0][1] + pos[1] * aij[1][1] + pos[2] * aij[2][1]
    z = pos[0] * aij[0][2] + pos[1] * aij[1][2] + pos[2] * aij[2][2]

    return x, y, z

def distance2(pos0, pos1, gij):
    """
    概要: 2つの分数座標間の距離の2乗を計算する。
    詳細説明:
        2つの原子の分数座標 `pos0` と `pos1` の間の距離の2乗を、
        メトリックテンソル `gij` を使用して計算します。
        この関数は周期境界条件を考慮しません。
    :param pos0: numpy.ndarray or list of float: 1点目の分数座標 `[x0, y0, z0]`。
    :param pos1: numpy.ndarray or list of float: 2点目の分数座標 `[x1, y1, z1]`。
    :param gij: numpy.ndarray: 3x3 のメトリックテンソル。
    :returns: float: 2点間の距離の2乗。
    """
    dx = pos1 - pos0
#    dx = [pos1[0] - pos0[0], pos1[1] - pos0[1], pos1[2] - pos0[2]]
    r2 = gij[0][0] * dx[0]*dx[0] + gij[1][1] * dx[1]*dx[1] + gij[2][2] * dx[2]*dx[2] \
       + 2.0 * (gij[0][1] * dx[0]*dx[1] + gij[0][2] * dx[0]*dx[2] + gij[1][2] * dx[1]*dx[2])

    return r2

def distance(pos0, pos1, gij):
    """
    概要: 2つの分数座標間の距離を計算する。
    詳細説明:
        2つの原子の分数座標 `pos0` と `pos1` の間の距離を、
        メトリックテンソル `gij` を使用して計算します。
        この関数は周期境界条件を考慮しません。
    :param pos0: numpy.ndarray or list of float: 1点目の分数座標 `[x0, y0, z0]`。
    :param pos1: numpy.ndarray or list of float: 2点目の分数座標 `[x1, y1, z1]`。
    :param gij: numpy.ndarray: 3x3 のメトリックテンソル。
    :returns: float: 2点間の距離。
    """
    dx = pos1 - pos0
#    dx = [pos1[0] - pos0[0], pos1[1] - pos0[1], pos1[2] - pos0[2]]
    r2 = gij[0][0] * dx[0]*dx[0] + gij[1][1] * dx[1]*dx[1] + gij[2][2] * dx[2]*dx[2] \
       + 2.0 * (gij[0][1] * dx[0]*dx[1] + gij[0][2] * dx[0]*dx[2] + gij[1][2] * dx[1]*dx[2])
    r = sqrt(r2)

    return r

def angle(pos0, pos1, pos2, gij):
    """
    概要: 3つの分数座標で定義される角度を計算する。
    詳細説明:
        `pos0` を頂点とし、`pos0-pos1` と `pos0-pos2` の2つのベクトルがなす角度を計算します。
        角度の計算にはメトリックテンソル `gij` が使用されます。
        結果は度数で返され、180度を超える場合は `360.0 - angle` として調整されます。
    :param pos0: numpy.ndarray or list of float: 角度の頂点となる分数座標 `[x0, y0, z0]`。
    :param pos1: numpy.ndarray or list of float: 1つ目のベクトルを定義する分数座標 `[x1, y1, z1]`。
    :param pos2: numpy.ndarray or list of float: 2つ目のベクトルを定義する分数座標 `[x2, y2, z2]`。
    :param gij: numpy.ndarray: 3x3 のメトリックテンソル。
    :returns: float: 3点間で形成される角度（度）。
                     いずれかの距離が0の場合は0.0を返します。
    """
    dis01 = distance(pos0, pos1, gij)
    if dis01 == 0.0:
        return 0.0
    dis02 = distance(pos0, pos2, gij)
    if dis02 == 0.0:
        return 0.0

    dx01 = pos1 - pos0
    dx02 = pos2 - pos0
#    dx01 = [pos1[0] - pos0[0], pos1[1] - pos0[1], pos1[2] - pos0[2]]
#    dx02 = [pos2[0] - pos0[0], pos2[1] - pos0[1], pos2[2] - pos0[2]]
    ip = gij[0][0] * dx01[0]*dx02[0] + gij[1][1] * dx01[1]*dx02[1] + gij[2][2] * dx01[2]*dx02[2] \
       + 2.0 * (gij[0][1] * dx01[0]*dx02[1] + gij[0][2] * dx01[0]*dx02[2] + gij[1][2] * dx01[1]*dx02[2])

    cosa = ip / dis01 / dis02
    angle = todeg * arccos(cosa)
    if angle > 180.0:
        angle = 360.0 - angle 

    return angle;

def configure_axis_structure(ax, xrange, yrange, zrange, fontsize = 12, legend_fontsize = 12):
    """
    概要: Matplotlibの3D軸オブジェクトの表示設定を行う。
    詳細説明:
        MatplotlibのAxes3Dオブジェクト (`ax`) に対して、軸ラベルのフォントサイズ、
        軸の線、目盛り、グリッド、ペインの表示設定をカスタマイズします。
        デフォルトでは、軸の線、目盛り、ラベルは非表示に設定され、
        箱型のアスペクト比で視覚的に整えられます。
    :param ax: matplotlib.axes.Axes: 設定を行うMatplotlibの3D軸オブジェクト。
    :param xrange: tuple or list of float: x軸の表示範囲 `(xmin, xmax)`。
    :param yrange: tuple or list of float: y軸の表示範囲 `(ymin, ymax)`。
    :param zrange: tuple or list of float: z軸の表示範囲 `(zmin, zmax)`。
    :param fontsize: int, optional: 軸ラベルおよび目盛りのフォントサイズ。
                                    0以下の場合は軸ラベルと目盛りを非表示にします。デフォルトは12。
    :param legend_fontsize: int, optional: 凡例のフォントサイズ（この関数では直接使用されません）。デフォルトは12。
    :returns: None
    """
    ax.tick_params(axis = 'both', which = 'major', labelsize = fontsize)
    if fontsize > 0:
        ax.set_xlabel(f'$x$', fontsize = fontsize, labelpad = -5)
        ax.set_ylabel(f'$y$', fontsize = fontsize, labelpad = -5)
        ax.set_zlabel(f'$z$', fontsize = fontsize, labelpad = 0)

# x軸、y軸、z軸の線を非表示にする
    ax.xaxis.line.set_visible(False)
    ax.yaxis.line.set_visible(False)
    ax.zaxis.line.set_visible(False)
    
    ax.grid(False)
    ax.xaxis.pane.fill = False
    ax.yaxis.pane.fill = False
    ax.zaxis.pane.fill = False
    ax.xaxis.pane.set_edgecolor('w')
    ax.yaxis.pane.set_edgecolor('w')
    ax.zaxis.pane.set_edgecolor('w')
    ax.xaxis.pane.set_alpha(0)
    ax.yaxis.pane.set_alpha(0)
    ax.zaxis.pane.set_alpha(0)
#        ax.w_xaxis.line.set_color('none')
#        ax.w_yaxis.line.set_color('none')
#        ax.w_zaxis.line.set_color('none')
    ax.set_xticks([])
    ax.set_yticks([])
    ax.set_zticks([])
    ax.set_xticklabels([]) 
    ax.set_yticklabels([]) 
    ax.set_zticklabels([]) 

    ax.set_xlim(xrange)
    ax.set_ylim(yrange)
    ax.set_zlim(zrange)
    ax.set_aspect('equal','box')
#    ax.set_box_aspect([xrange[1] - xrange[0], yrange[1] - yrange[0], zrange[1] - zrange[0]]) 

def draw_box(ax, aij, nrange, color = 'black'):
    """
    概要: 3D軸上に単位格子のボックスを描画する。
    詳細説明:
        格子ベクトル `aij` で定義される単位格子の12本の辺を、Matplotlibの3D軸 `ax` に描画します。
        この関数は、単一の単位格子を描画するために設計されており、`nrange` は現在の実装では使用されません。
    :param ax: matplotlib.axes.Axes: 描画対象のMatplotlib 3D軸オブジェクト。
    :param aij: numpy.ndarray: 3x3 の格子ベクトル行列。
                               `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
    :param nrange: list of list of int: 描画するセルの範囲 `[[xmin, xmax], [ymin, ymax], [zmin, zmax]]`。
                                        この関数では使用されません。
    :param color: str, optional: ボックスの線の色。デフォルトは 'black'。
    :returns: None
    """
# (0,0,0) -> ax
    ax.plot([0.0, aij[0][0]], 
            [0.0, aij[0][1]], 
            [0.0, aij[0][2]], color = color)
# (0,0,0) -> ay
    ax.plot([0.0, aij[1][0]], 
            [0.0, aij[1][1]], 
            [0.0, aij[1][2]], color = color)
# (0,0,0) -> az
    ax.plot([0.0, aij[2][0]], 
            [0.0, aij[2][1]], 
            [0.0, aij[2][2]], color = color)

# ax -> ax + ay
    ax.plot([aij[0][0], aij[0][0] + aij[1][0]], 
            [aij[0][1], aij[0][1] + aij[1][1]], 
            [aij[0][2], aij[0][2] + aij[1][2]], color = color)
# ax -> ax + az
    ax.plot([aij[0][0], aij[0][0] + aij[2][0]], 
            [aij[0][1], aij[0][1] + aij[2][1]], 
            [aij[0][2], aij[0][2] + aij[2][2]], color = color)

# ay -> ay + ax
    ax.plot([aij[1][0], aij[1][0] + aij[0][0]], 
            [aij[1][1], aij[1][1] + aij[0][1]], 
            [aij[1][2], aij[1][2] + aij[0][2]], color = color)
# ay -> ay + az
    ax.plot([aij[1][0], aij[1][0] + aij[2][0]], 
            [aij[1][1], aij[1][1] + aij[2][1]], 
            [aij[1][2], aij[1][2] + aij[2][2]], color = color)

# az -> az + ax
    ax.plot([aij[2][0], aij[2][0] + aij[0][0]], 
            [aij[2][1], aij[2][1] + aij[0][1]], 
            [aij[2][2], aij[2][2] + aij[0][2]], color = color)
# az -> ax + ay
    ax.plot([aij[2][0], aij[2][0] + aij[1][0]], 
            [aij[2][1], aij[2][1] + aij[1][1]], 
            [aij[2][2], aij[2][2] + aij[1][2]], color = color)

# ax + ay -> ax + ay + az
    ax.plot([aij[0][0] + aij[1][0], aij[0][0] + aij[1][0] + aij[2][0]], 
            [aij[0][1] + aij[1][1], aij[0][1] + aij[1][1] + aij[2][1]], 
            [aij[0][2] + aij[1][2], aij[0][2] + aij[1][2] + aij[2][2]], color = color)

# ax + az -> ax + ay + az
    ax.plot([aij[0][0] + aij[2][0], aij[0][0] + aij[1][0] + aij[2][0]], 
            [aij[0][1] + aij[2][1], aij[0][1] + aij[1][1] + aij[2][1]], 
            [aij[0][2] + aij[2][2], aij[0][2] + aij[1][2] + aij[2][2]], color = color)

# ay + az -> ax + ay + az
    ax.plot([aij[1][0] + aij[2][0], aij[0][0] + aij[1][0] + aij[2][0]], 
            [aij[1][1] + aij[2][1], aij[0][1] + aij[1][1] + aij[2][1]], 
            [aij[1][2] + aij[2][2], aij[0][2] + aij[1][2] + aij[2][2]], color = color)

def draw_unitcell(ax, sites, aij, nrange, color = 'black', facecolor = 'black', edgecolor = 'white', alpha = 0.7, kr = 1.0):
    """
    概要: 3D軸上に単位格子とその内部のサイト原子を描画する。
    詳細説明:
        まず `draw_box` 関数を呼び出して単位格子の境界線を描画します。
        次に、`sites` リストに含まれる各原子について、その分数座標をデカルト座標に変換し、
        指定された `nrange` の範囲内の周期的な位置に原子を散布図として描画します。
        原子の描画スタイル（色、透明度、サイズ）は引数でカスタマイズ可能です。
    :param ax: matplotlib.axes.Axes: 描画対象のMatplotlib 3D軸オブジェクト。
    :param sites: list of list: サイト情報（原子名、ラベル、原子番号、原子質量、電荷、半径、色、位置）のリスト。
                                各要素は `[name, label, z, M, q, r, color, pos]` の形式。
                                `None` の場合、原子は描画されません。
    :param aij: numpy.ndarray: 3x3 の格子ベクトル行列。
                               `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
    :param nrange: list of list of int: 描画するセルの範囲 `[[xmin, xmax], [ymin, ymax], [zmin, zmax]]`。
                                        原子の繰り返し構造を描画するために使用されます。
    :param color: str, optional: 単位格子ボックスの線の色。デフォルトは 'black'。
    :param facecolor: str, optional: サイト原子の面の描画色。デフォルトは 'black'。
    :param edgecolor: str, optional: サイト原子のエッジの描画色。デフォルトは 'white'。
    :param alpha: float, optional: サイト原子の透明度（0.0から1.0）。デフォルトは 0.7。
    :param kr: float, optional: サイト原子の描画サイズに対するスケール因子。デフォルトは 1.0。
    :returns: None
    """
    draw_box(ax, aij, nrange, color)

    if sites is None: return

    for isite, site in enumerate(sites):
        name, label, z, M, q, r, color, pos = site
        pos01 = [reduce01(pos[0]), reduce01(pos[1]), reduce01(pos[2])]
        for iz in range(int(nrange[2][0]) - 1, int(nrange[2][1]) + 1):
         for iy in range(int(nrange[1][0]) - 1, int(nrange[1][1]) + 1):
          for ix in range(int(nrange[0][0]) - 1, int(nrange[0][1]) + 1):
            posn = [pos01[0] + ix, pos01[1] + iy, pos01[2] + iz]
            if    posn[0] < nrange[0][0] or nrange[0][1] < posn[0]  \
               or posn[1] < nrange[1][0] or nrange[1][1] < posn[1]  \
               or posn[2] < nrange[2][0] or nrange[2][1] < posn[2]:
                  continue

            x, y, z = fractional_to_cartesian(posn, aij)
            ax.scatter([x], [y], [z], marker = 'o', c = facecolor, edgecolors = edgecolor, alpha = alpha, s = kr *r)


def main():
    """
    概要: モジュールの主要な計算と結果の表示を実行する。
    詳細説明:
        定義された `lattice_parameters` を使用して、実空間および逆空間の
        格子ベクトル、メトリックテンソル、単位格子の体積を計算します。
        さらに、サンプルとして原子間距離と角度の計算も行い、
        その結果を標準出力に表示します。
        最後にプログラムを終了します。
    :returns: None
    """
    print("")
    print("Lattice parameters:", lattice_parameters)
    aij = cal_lattice_vectors(lattice_parameters)
    print("Lattice vectors:")
    print("  ax: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[0][0], aij[0][1], aij[0][2]))
    print("  ay: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[1][0], aij[1][1], aij[1][2]))
    print("  az: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[2][0], aij[2][1], aij[2][2]))
    inf = cal_metrics(lattice_parameters)
    gij = inf['gij']
    print("Metric tensor:")
    print("  gij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[0]))
    print("       ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[1]))
    print("       ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[2]))
    print("")
    volume = cal_volume(aij)
    print("Volume: {:12.4g} A^3".format(volume))

    print("")
    print("Unit cell volume: {:12.4g} A^3".format(volume))
    Raij  = cal_reciprocal_lattice_vectors(aij)
    Rlatt = cal_reciprocal_lattice_parameters(Raij)
    Rinf  = cal_metrics(Rlatt)
    Rgij  = Rinf['gij']
    print("Reciprocal lattice parameters:", Rlatt)
    print("Reciprocal lattice vectors:")
    print("  Rax: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[0]))
    print("  Ray: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[1]))
    print("  Raz: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[2]))
    print("Reciprocal lattice metric tensor:")
    print("  Rgij: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[0]))
    print("        ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[1]))
    print("        ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[2]))
    Rvolume = cal_volume(Raij)
    print("Reciprocal unit cell volume: {:12.4g} A^-3".format(Rvolume))

    print("")
    print("dis=", distance(np.array([0,0,0]), np.array([1,1,1]), gij))
    print("angle=", angle (np.array([0,0,0]), np.array([1,1,1]), np.array([1,0,0]), gij))
    print("angle=", angle (np.array([0,0,0]), np.array([1,0,0]), np.array([0,1,0]), gij))
    
    print("")
    exit()


if __name__ == '__main__':
    main()