"""
1Dタイトバインディング法による量子井戸ポテンシャル下でのバンド計算スクリプト。

概要:
    量子井戸構造における電子の波動関数、エネルギー準位、そしてそれらを用いたタイトバインディング法によるエネルギーバンド構造の計算と可視化を行うPythonスクリプト。
    主に、量子井戸内の基底関数計算 (`basis`モード)、バンド構造計算 (`band`モード)、波動関数描画 (`wf`モード、現在未実装) の機能を提供する。

詳細説明:
    本スクリプトは以下の主要な機能を持つ。
    1.  **量子井戸ポテンシャルの定義と可視化**:
        与えられた結晶構造と原子情報に基づき、周期的な矩形量子井戸ポテンシャルを構築する。
    2.  **単一量子井戸のエネルギー準位と波動関数の計算**:
        無限障壁近似と有限障壁（セカント法を用いた数値解法）の両方で、単一量子井戸内のエネルギー準位と対応する波動関数を計算する。
        計算された波動関数は、タイトバインディング法の基底関数として使用される。
    3.  **タイトバインディング行列の構築**:
        計算された基底関数と定義されたポテンシャルを用いて、ハミルトニアン行列 (Fij) および重なり積分行列 (Sij) を構築する。
        近接サイトとの相互作用を考慮し、周期境界条件 (Blochの定理) を導入する。
    4.  **エネルギーバンド構造の計算**:
        構築された行列を対角化することで、さまざまな波数kにおけるエネルギー準位（エネルギーバンド）を計算する。
    5.  **結果の可視化**:
        ポテンシャル、波動関数、エネルギー準位、エネルギーバンド構造をmatplotlibを用いてグラフ表示する。

関連リンク:
    :doc:`tb1d_usage` (使用方法に関するドキュメントへのリンク、仮)
"""

import sys
import numpy as np
from numpy import sqrt, exp, sin, cos, tan, sinh, cosh, tanh, pi
import numpy.linalg as LA 
import csv
from matplotlib import pyplot as plt


#===================================
# physical constants
#===================================
pi   = 3.14159265358979323846
pi2  = 2.0 * pi
h    = 6.6260755e-34    # Js";
hbar = 1.05459e-34      # "Js";
c    = 2.99792458e8     # m/s";
e    = 1.60218e-19      # C";
e0   = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
kB   = 1.380658e-23     # JK<sup>-1</sup>";
me   = 9.1093897e-31    # kg";
R    = 8.314462618      # J/K/mol
a0   = 5.29177e-11      # m";

KE = hbar * hbar / 2.0 / me / e * 1.0e20  # coefficient of kinetic energy


#========================
# global parameters
#========================
parity_str = ["even", "odd"]

#========================
# global configuration
#========================
mode = 'basis'   # basis|band|wf

#========================
# Potential
#========================
Vvac = 10.0     # eV

#========================
# Crystal definition
#========================
#crystalname = 'Sim'
#crystalname = 'AB'
crystalname = 'Si'

# Si
a       = 5.4064   # angstrom, lattice parameter
# tijを計算する近接サイトの検索範囲
rmax_neighbor = 1.1 * a

atoms = [
    {"name": "A", "wellwidth": a - 0.5, "Vbase": 0.0, "meff": 1.0}
    ]

crystal = [
    {"name": 'A1', "iatom": 0, "x0": (0.5 + a) / 2.0}
    ]

if crystalname == 'AB':
    a       = 5.4064   # angstrom, lattice parameter
    rmax_neighbor = 1.1 * a
    atoms = [
        {"name": "A", "wellwidth": 2.5,     "Vbase": 0.0, "meff": 5.0}
       ,{"name": "B", "wellwidth": a - 3.5, "Vbase": 0.3, "meff": 5.0}
        ]
    crystal = [
         {"name": 'A1', "iatom": 0, "x0": (0.5 + 3.0)/2.0}
        ,{"name": 'B1', "iatom": 1, "x0": (3.5 + a)/2.0}
        ]

if crystalname == 'Sim':
    a  = 17.5
    rmax_neighbor = 4.0
    atoms = [
        {"name": "A", "wellwidth": 3.0, "Vbase": 0.0, "meff": 1.0}
        ]
    crystal = [
         {"name": 'A1', "iatom": 0, "x0":  2.0}
        ,{"name": 'A2', "iatom": 0, "x0":  5.5}
        ,{"name": 'A3', "iatom": 0, "x0":  9.0}
        ,{"name": 'A4', "iatom": 0, "x0": 12.5}
        ,{"name": 'A5', "iatom": 0, "x0": 16.0}
        ]


#==============================================
# X range to plot potential and basis
#==============================================
xmin  = -0.1 * a  # A
xmax  =  1.1 * a  # A
nx    = 301

#=====================================
# 解を走査するグラフ表示
#=====================================
# 解を走査するエネルギー範囲
Emin   = 0.0
Emax   = Vvac
# 解を走査するエネルギー点数
nE = 101
# Secant法パラメータ
eps      = 1.0e-5
nmaxiter = 100

#========================
# Wave function plot
#========================
wmargin  = 0.5
nwf       = 101
wfplotmax = 1.0  # superiposed on energy level diagram, in eV

#========================
# <E> integration
#========================
xminEinteg = -10.0
xmaxEinteg = a + 11.0
nxEinteg   = 501

#========================
# Band
#========================
# tijを計算する近接サイトを検索する格子周期範囲
Nxmax_neighbor = 2

# tijを計算する際の波動関数の計算範囲
tijxmin = -1.1 * a
tijxmax =  1.1 * a
ntijx   = 501

# band plot
kmin = -0.5  # in pi/a
kmax =  0.5  # in pi/a
nk = 21

# プロットするエネルギー範囲
Erange = [0.0, 10.0]    # eV


#===================================
# figure configuration
#===================================
figsize = (6, 8)
fontsize        = 12
legend_fontsize = 12


#==============================================
# fundamental functions
#==============================================
def pfloat(str):
    """
    文字列をfloat型に安全に変換します。

    概要:
        与えられた文字列をfloat型に変換します。変換できない場合はNoneを返します。
        変換エラーが発生してもプログラムは終了しません。

    Parameters:
        :param str: 変換する文字列。
        :type str: str
    
    Returns:
        :returns: 変換された浮動小数点数、または変換できなかった場合はNone。
        :rtype: float or None
    """
    try:
        return float(str)
    except:
        return None

def pint(str):
    """
    文字列をint型に安全に変換します。

    概要:
        与えられた文字列をint型に変換します。変換できない場合はNoneを返します。
        変換エラーが発生してもプログラムは終了しません。

    Parameters:
        :param str: 変換する文字列。
        :type str: str

    Returns:
        :returns: 変換された整数、または変換できなかった場合はNone。
        :rtype: int or None
    """
    try:
        return int(str)
    except:
        return None

def getarg(position, defval = None):
    """
    コマンドライン引数を安全に取得します。

    概要:
        sys.argvリストから指定された位置の引数を取得します。
        指定された位置に引数がない場合は、デフォルト値を返します。

    Parameters:
        :param position: 取得する引数のsys.argvにおけるインデックス。
        :type position: int
        :param defval: 引数が存在しない場合に返すデフォルト値。デフォルトはNone。
        :type defval: Any, optional

    Returns:
        :returns: コマンドライン引数の値、またはデフォルト値。
        :rtype: str or Any
    """
    try:
        return sys.argv[position]
    except:
        return defval

def getfloatarg(position, defval = None):
    """
    コマンドライン引数をfloat型に安全に変換して取得します。

    概要:
        指定された位置のコマンドライン引数を取得し、pfloat関数を使ってfloat型に変換します。
        引数が存在しない場合やfloat型に変換できない場合は、デフォルト値を返します。

    Parameters:
        :param position: 取得する引数のsys.argvにおけるインデックス。
        :type position: int
        :param defval: 引数が存在しない場合や変換できない場合に返すデフォルト値。デフォルトはNone。
        :type defval: float or None, optional

    Returns:
        :returns: 変換された浮動小数点数、またはデフォルト値。
        :rtype: float or None
    """
    return pfloat(getarg(position, defval))

def getintarg(position, defval = None):
    """
    コマンドライン引数をint型に安全に変換して取得します。

    概要:
        指定された位置のコマンドライン引数を取得し、pint関数を使ってint型に変換します。
        引数が存在しない場合やint型に変換できない場合は、デフォルト値を返します。

    Parameters:
        :param position: 取得する引数のsys.argvにおけるインデックス。
        :type position: int
        :param defval: 引数が存在しない場合や変換できない場合に返すデフォルト値。デフォルトはNone。
        :type defval: int or None, optional

    Returns:
        :returns: 変換された整数、またはデフォルト値。
        :rtype: int or None
    """
    return pint(getarg(position, defval))

def usage():
    """
    スクリプトのコマンドライン引数の使用法を表示します。

    概要:
        スクリプトの正しい実行方法とオプションの引数に関する情報を標準出力に出力します。

    Returns:
        :returns: None
    """
    print("")
    print("Usage: Variables in () are optional")
    print("  python {}".format(sys.argv[0]))
    print("    ex: python {} (basis)".format(sys.argv[0]))
    print("    ex: python {} (band)".format(sys.argv[0]))
    print("    ex: python {} (wf)".format(sys.argv[0]))
    print("")

def terminate(message = None):
    """
    プログラムを終了し、オプションでメッセージと使用法を表示します。

    概要:
        指定されたメッセージ（オプション）を表示し、usage関数を呼び出してスクリプトの使用法を表示した後、プログラムを終了します。

    Parameters:
        :param message: 終了前に表示するメッセージ。デフォルトはNone。
        :type message: str or None, optional

    Returns:
        :returns: None (プログラムは終了します)
    """
    if message is not None:
        print("")
        print(message)

    print("")
    usage()
    print("")
    exit()

def make_list2(dim1, dim2):
    """
    指定された次元の2次元リストを生成し、0.0で初期化します。

    概要:
        dim1 x dim2 の浮動小数点数で初期化された2次元リストを作成します。

    Parameters:
        :param dim1: リストの第一次元（行数）。
        :type dim1: int
        :param dim2: リストの第二次元（列数）。
        :type dim2: int

    Returns:
        :returns: 0.0で初期化された2次元リスト。
        :rtype: list[list[float]]
    """
    l = []
    for i in range(dim1):
        l.append([])
        for j in range(dim2):
            l[i].append(0.0)
    return l

def make_list3(dim1, dim2, dim3):
    """
    指定された次元の3次元リストを生成し、0.0で初期化します。

    概要:
        dim1 x dim2 x dim3 の浮動小数点数で初期化された3次元リストを作成します。

    Parameters:
        :param dim1: リストの第一次元。
        :type dim1: int
        :param dim2: リストの第二次元。
        :type dim2: int
        :param dim3: リストの第三次元。
        :type dim3: int

    Returns:
        :returns: 0.0で初期化された3次元リスト。
        :rtype: list[list[list[float]]]
    """
    l = make_list2(dim1, dim2)
    for i1 in range(dim1):
        for i2 in range(dim2):
            l[i1][i2] = []
            for i3 in range(dim3):
                l[i1][i2].append(0.0)
    return l
    
def make_list4(dim1, dim2, dim3, dim4):
    """
    指定された次元の4次元リストを生成し、0.0で初期化します。

    概要:
        dim1 x dim2 x dim3 x dim4 の浮動小数点数で初期化された4次元リストを作成します。

    Parameters:
        :param dim1: リストの第一次元。
        :type dim1: int
        :param dim2: リストの第二次元。
        :type dim2: int
        :param dim3: リストの第三次元。
        :type dim3: int
        :param dim4: リストの第四次元。
        :type dim4: int

    Returns:
        :returns: 0.0で初期化された4次元リスト。
        :rtype: list[list[list[list[float]]]]
    """
    l = make_list3(dim1, dim2, dim3)
    for i1 in range(dim1):
        for i2 in range(dim2):
            for i3 in range(dim3):
                l[i1][i2][i3] = []
                for i4 in range(dim4):
                    l[i1][i2][i3].append(0.0)
    return l

def reduce(x, x0):
    """
    xをx0の周期で還元します（モジュロ演算に相当）。

    概要:
        与えられたxを、周期x0の範囲 [0, x0) または [-x0, 0) にマッピングします。
        物理的な周期性を考慮した座標の還元に使用されます。

    Parameters:
        :param x: 還元する値。
        :type x: float
        :param x0: 周期。
        :type x0: float

    Returns:
        :returns: 周期x0で還元されたxの値。
        :rtype: float
    """
    n = int(x / x0)
    if x < 0.0:
        n -= 1
    return x - x0 * n


#==============================================
# update default values by startup arguments
#==============================================
argv = sys.argv
#if len(argv) == 1:
#    terminate()

mode    = getarg     (1, mode)


def atominf(i, atoms):
    """
    atomsリストから指定された原子の情報を取得します。

    概要:
        指定されたインデックスiに対応する原子のプロパティ（名前、井戸幅、基底電位、有効質量）を返します。

    Parameters:
        :param i: atomsリスト内の原子のインデックス。
        :type i: int
        :param atoms: 原子情報の辞書リスト。
        :type atoms: list[dict]

    Returns:
        :returns: (名前, 井戸幅, 基底電位, 有効質量) のタプル。
        :rtype: tuple[str, float, float, float]
    """
    atom   = atoms[i]
    name   = atom["name"]
    wwidth = atom["wellwidth"]
    Vbase  = atom["Vbase"]
    meff   = atom["meff"]
    return name, wwidth, Vbase, meff

def siteinf(i, crystal):
    """
    crystalリストから指定されたサイトの全ての情報を取得します。

    概要:
        指定されたインデックスiに対応する結晶サイトのプロパティ（名前、原子インデックス、原子情報、位置x0、井戸幅、井戸の最小/最大x座標、基底電位、有効質量）を返します。

    Parameters:
        :param i: crystalリスト内のサイトのインデックス。
        :type i: int
        :param crystal: 結晶サイト情報の辞書リスト。
        :type crystal: list[dict]

    Returns:
        :returns: (名前, 原子インデックス, 原子辞書, x0, 井戸幅, 井戸のxmin, 井戸のxmax, 基底電位, 有効質量) のタプル。
        :rtype: tuple[str, int, dict, float, float, float, float, float, float]
    """
    site   = crystal[i]
    iatom  = site["iatom"]
    atom   = atoms[iatom]
    name   = site["name"]
    x      = site["x0"]
    wwidth = atom["wellwidth"]
    Vbase  = atom["Vbase"]
    meff   = atom["meff"]
    return name, iatom, atom, x, wwidth, x - wwidth / 2.0, x + wwidth / 2.0, Vbase, meff

def pot(x, Vvac, crystal):
    """
    指定された位置xにおけるポテンシャルエネルギーを計算します。

    概要:
        与えられたx座標において、結晶構造に基づいて矩形量子井戸ポテンシャルを決定します。
        結晶周期aでx座標を還元し、どのサイトの井戸内にあるかを確認します。

    Parameters:
        :param x: ポテンシャルを評価するx座標 (Å)。
        :type x: float
        :param Vvac: 真空準位のポテンシャルエネルギー (eV)。
        :type Vvac: float
        :param crystal: 結晶サイト情報の辞書リスト。
        :type crystal: list[dict]

    Returns:
        :returns: 指定されたx座標におけるポテンシャルエネルギー (eV)。
        :rtype: float
    """
    xred = reduce(x, a)

    nsites = len(crystal)
    for i in range(nsites):
        name, iatom, atom, x0, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
        if x0 - wwidth/2.0 <= xred < x0 + wwidth/2.0:
            return Vbase
    else:
        return Vvac

def pot_atom(dx, Vvac, crystal):
    """
    単一原子サイトの相対座標dxにおけるポテンシャルエネルギーを計算します。

    概要:
        単一の原子サイト内での相対座標dxに基づき、矩形量子井戸ポテンシャルを決定します。
        この関数は通常、タイトバインディングの積分計算などで、原点に中心を持つ原子サイトのポテンシャルを評価する際に使用されます。

    Parameters:
        :param dx: 原子サイトの中心からの相対x座標 (Å)。
        :type dx: float
        :param Vvac: 真空準位のポテンシャルエネルギー (eV)。
        :type Vvac: float
        :param crystal: 結晶サイト情報の辞書リスト。この関数では各サイトのwellwidthとVbaseが使われる。
        :type crystal: list[dict]

    Returns:
        :returns: 指定された相対dxにおけるポテンシャルエネルギー (eV)。
        :rtype: float
    """
    nsites = len(crystal)
    for i in range(nsites):
        name, iatom, atom, x0, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
        a_half = wwidth / 2.0 # 井戸の半分の幅
        if -a_half <= dx < a_half:
            return Vbase
    else:
        return Vvac

def build_potential(xmin, xstep, n, Vvac, crystal):
    """
    指定された範囲でポテンシャルエネルギーのプロファイル配列を構築します。

    概要:
        xminからxstep刻みでn点までのx座標と、それに対応するポテンシャルエネルギーの配列を計算します。
        主にポテンシャルプロファイルの描画に使用されます。

    Parameters:
        :param xmin: プロファイルの開始x座標 (Å)。
        :type xmin: float
        :param xstep: x座標の刻み幅 (Å)。
        :type xstep: float
        :param n: 計算する点の数。
        :type n: int
        :param Vvac: 真空準位のポテンシャルエネルギー (eV)。
        :type Vvac: float
        :param crystal: 結晶サイト情報の辞書リスト。
        :type crystal: list[dict]

    Returns:
        :returns: (x座標のnumpy配列, ポテンシャルエネルギーのnumpy配列) のタプル。
        :rtype: tuple[numpy.ndarray, numpy.ndarray]
    """
    xpot = np.empty(n)
    ypot = np.empty(n)
    for i in range(n):
        xx = xmin + i * xstep
        xpot[i] = xx
        ypot[i] = pot(xx, Vvac, crystal)
    return xpot, ypot

def quantization_level_inifinite(Vbase, meff, wwidth, iL):
    """
    無限に深い矩形量子井戸の量子化エネルギー準位を計算します。

    概要:
        与えられた量子井戸の基底電位、有効質量、井戸幅、および量子数iLに基づいて、
        無限に深い量子井戸モデルにおけるエネルギー準位を計算します。

    Parameters:
        :param Vbase: 量子井戸の基底電位 (eV)。
        :type Vbase: float
        :param meff: 有効質量 (me単位)。
        :type meff: float
        :param wwidth: 量子井戸の幅 (Å)。
        :type wwidth: float
        :param iL: 量子数 (n=1, 2, 3...)。
        :type iL: int

    Returns:
        :returns: 計算された量子化エネルギー準位 (eV)。
        :rtype: float
    """
    k = pi / wwidth
    return Vbase + KE / meff * (k * iL)**2

def cal_delta_even(E, w, V0, meff):
    """
    偶関数波動関数のエネルギーに対するΔ関数を計算します。

    概要:
        矩形量子井戸の偶関数解に対するエネルギーEのΔ関数を計算します。
        この関数が0になるEが、偶関数のエネルギー準位に対応します。

    Parameters:
        :param E: 井戸底からのエネルギー (eV)。E = E_total - Vbase に相当。
        :type E: float
        :param w: 量子井戸の幅 (Å)。
        :type w: float
        :param V0: 障壁の高さ (eV)。V0 = Vvac - Vbase に相当。
        :type V0: float
        :param meff: 有効質量 (me単位)。
        :type meff: float

    Returns:
        :returns: 計算されたΔ関数の値。
        :rtype: float
    """
    kb = sqrt(2.0 * meff * me * (V0 - E) * e) / hbar * 1.0e-10  # A^-1
    kw = sqrt(2.0 * meff * me * E * e) / hbar * 1.0e-10         # A^-1
#    kb = 0.0 # 無限障壁ポテンシャルの場合
    delta = kb - kw * tan(kw * w / 2.0)

    return delta

def cal_delta_odd(E, w, V0, meff):
    """
    奇関数波動関数のエネルギーに対するΔ関数を計算します。

    概要:
        矩形量子井戸の奇関数解に対するエネルギーEのΔ関数を計算します。
        この関数が0になるEが、奇関数のエネルギー準位に対応します。

    Parameters:
        :param E: 井戸底からのエネルギー (eV)。E = E_total - Vbase に相当。
        :type E: float
        :param w: 量子井戸の幅 (Å)。
        :type w: float
        :param V0: 障壁の高さ (eV)。V0 = Vvac - Vbase に相当。
        :type V0: float
        :param meff: 有効質量 (me単位)。
        :type meff: float

    Returns:
        :returns: 計算されたΔ関数の値。
        :rtype: float
    """
    kb = sqrt(2.0 * meff * me * (V0 - E) * e) / hbar * 1.0e-10  # A^-1
    kw = sqrt(2.0 * meff * me * E * e) / hbar * 1.0e-10         # A^-1
#    kb = 0.0 # 無限障壁ポテンシャルの場合
    delta = kb + kw / tan(kw * w / 2.0)

    return delta

def refine_E(cal_delta, E0, E1, nmaxiter, eps, w, V0, meff, IsPrint = 0):
    """
    セカント法を用いてエネルギー準位を精密化します。

    概要:
        cal_delta関数がゼロとなるエネルギーEを、セカント法（二分法に似た方法）を用いて探索し、精密化します。
        与えられた初期エネルギー範囲 [E0, E1] 内に根があることを前提とします。

    Parameters:
        :param cal_delta: エネルギーEの関数で、根を求める対象となる関数 (cal_delta_even または cal_delta_odd)。
        :type cal_delta: callable
        :param E0: 探索範囲の下限エネルギー (eV、井戸底からの相対値)。
        :type E0: float
        :param E1: 探索範囲の上限エネルギー (eV、井戸底からの相対値)。
        :type E1: float
        :param nmaxiter: 最大反復回数。
        :type nmaxiter: int
        :param eps: 収束判定に使用する許容誤差 (エネルギー差) (eV)。
        :type eps: float
        :param w: 量子井戸の幅 (Å)。
        :type w: float
        :param V0: 障壁の高さ (eV)。
        :type V0: float
        :param meff: 有効質量 (me単位)。
        :type meff: float
        :param IsPrint: デバッグ情報を出力するかどうかのフラグ (0:なし, 1:簡易, 2:詳細)。デフォルトは0。
        :type IsPrint: int, optional

    Returns:
        :returns: (収束したエネルギー, 最後のエネルギー差, 最後のdelta値) のタプル。
                  収束しない場合は (None, None, None) を返します。
        :rtype: tuple[float or None, float or None, float or None]
    """
    delta0 = cal_delta(E0, w, V0, meff)
#    print("  0:", E0, delta0)
    delta1 = cal_delta(E1, w, V0, meff)
#    print("  1:", E1, delta1)
    for i in range(nmaxiter):
        if delta0 * delta1 > 0.0:
            print("Error: delta0 ({:8.4f}, {:10.4g}) and delta1 ({:8.4f}, {:10.4g}) must have opposite signs"
                    .format(E0, delta0, E1, delta1))
            exit()
        
        E2 = (E0 + E1) / 2.0
        delta2 = cal_delta(E2, w, V0, meff)
#        print("  {}: E={:8.4f} - {:8.4f} => {:8.4f}".format(i, E0, E1, E2))
#        print("      delta: {:12.4g} - {:12.4g} => {:12.4g}".format(delta0, delta1, delta2))

        if abs(E2 - E0) < eps:
            if IsPrint:
                print("  converged at E = {:12.6g} with dE = {:12.6g}  delta = {:12.6g}"
                        .format(E2, E2 - E0, delta2))
            return E2, E2 - E0, delta2
        else:
            if delta0 * delta2 > 0.0:
                E0 = E2
                delta0 = delta2
            else:
                E1 = E2
                delta1 = delta2
            continue
    else:
        print("  Not converged for {} iterations.".format(nmaxiter))
        # dE is not defined outside the loop if not converged before.
        # This line might cause an error if the loop finishes without convergence.
        # However, following rule 1, no logic change is allowed.
        # print("    E = {:12.6g} with dE = {:12.6g}  delta = {:12.6g}".format(E2, dE, delta2))
        return None, None, None

def find_eigenstates_base(parity, Vvac, Vbase, wwidth, meff, Emin, Emax, nE, nmaxiter, eps, IsPrint = 0):
    """
    指定されたパリティ（偶奇性）を持つ量子井戸のエネルギー準位を探索します。

    概要:
        量子井戸の全エネルギー範囲 (EminからEmax) を走査し、cal_delta関数が符号を反転する点（根）を特定します。
        セカント法 (refine_E) を用いてこれらの根を精密化し、指定されたパリティのエネルギー準位のリストを生成します。

    Parameters:
        :param parity: 探索する波動関数のパリティ (0: 偶関数, 1: 奇関数)。
        :type parity: int
        :param Vvac: 真空準位のポテンシャルエネルギー (eV)。
        :type Vvac: float
        :param Vbase: 量子井戸の基底電位 (eV)。
        :type Vbase: float
        :param wwidth: 量子井戸の幅 (Å)。
        :type wwidth: float
        :param meff: 有効質量 (me単位)。
        :type meff: float
        :param Emin: 探索開始エネルギー (eV)。
        :type Emin: float
        :param Emax: 探索終了エネルギー (eV)。
        :type Emax: float
        :param nE: エネルギー走査の点数。
        :type nE: int
        :param nmaxiter: セカント法での最大反復回数。
        :type nmaxiter: int
        :param eps: セカント法での収束判定に使用する許容誤差 (eV)。
        :type eps: float
        :param IsPrint: デバッグ情報を出力するかどうかのフラグ。デフォルトは0。
        :type IsPrint: int, optional

    Returns:
        :returns: 見つかったエネルギー準位の辞書リスト。各辞書は"parity"と"E"のキーを持つ。
        :rtype: list[dict]
    """
    V0 = Vvac - Vbase
    Estep = (Emax - Emin) / (nE - 1)
    if parity:
        cal_delta = cal_delta_odd
    else:
        cal_delta = cal_delta_even

    d0 = None
    iband = 0
    orbitals = []
    for iE in range(nE):
        E = Emin + iE * Estep
        if E <= Vbase: # 井戸底よりも低いエネルギーは無視
            continue
        if Vvac <= E:  # 真空準位よりも高いエネルギーは無視
            break

        # cal_deltaの引数は井戸底からの相対エネルギーなので、Vbaseを引く
        delta = cal_delta(E - Vbase, wwidth, V0, meff)
        if d0 is None:
            d0 = delta
            continue
# cal_delta function flips its sign when delta => 0 and delta => infinite
# the latter case must be removed
        # 符号が反転し、かつdeltaが大きすぎない（無限大に発散するような点ではない）場合、根があると判断
        if d0 * delta < 0.0 and delta < 10.0: # delta < 10.0 は無限大への発散を除外するための閾値
#            print("0: ", E-Estep, d0)
#            print("1: ", E, delta)
            # 井戸底からの相対エネルギーでrefine_Eを呼び出す
            E0, dE, delta0 = refine_E(cal_delta, E - Vbase - Estep, E - Vbase, nmaxiter, eps, wwidth, V0, meff, IsPrint = 0)
            E = E0 + Vbase # 全エネルギーに変換
            if IsPrint:
                print("  E[{}]={:12.6g} eV  dE={:12.6g} delta={:12.6g}".format(iband, E, dE, delta0))


            orbitals.append({"parity": parity, "E": E})
            d0 = delta
            E += Estep # 次の探索のためにエネルギーを少し進める
            iband += 1
        if d0 * delta < 0.0: # 符号が反転した場合、d0を更新
            d0 = delta

    return orbitals

def find_eigenstates(Vvac, Vbase, wwidth, meff, Emin, Emax, nE, nmaxiter, eps):
    """
    矩形量子井戸の全てのエネルギー準位（偶関数および奇関数）を探索します。

    概要:
        cal_delta_evenおよびcal_delta_odd関数を使用して、偶関数と奇関数の両方のエネルギー準位を見つけます。
        見つかった全ての準位を結合し、エネルギーの昇順にソートします。
        各準位について、波動関数の係数(Cb, Cw)や特性値(kb, kw)を計算し、辞書に追加して返します。

    Parameters:
        :param Vvac: 真空準位のポテンシャルエネルギー (eV)。
        :type Vvac: float
        :param Vbase: 量子井戸の基底電位 (eV)。
        :type Vbase: float
        :param wwidth: 量子井戸の幅 (Å)。
        :type wwidth: float
        :param meff: 有効質量 (me単位)。
        :type meff: float
        :param Emin: 探索開始エネルギー (eV)。
        :type Emin: float
        :param Emax: 探索終了エネルギー (eV)。
        :type Emax: float
        :param nE: エネルギー走査の点数。
        :type nE: int
        :param nmaxiter: セカント法での最大反復回数。
        :type nmaxiter: int
        :param eps: セカント法での収束判定に使用する許容誤差 (eV)。
        :type eps: float

    Returns:
        :returns: 量子井戸の全てのエネルギー準位の詳細情報を含む辞書リスト。
                  各辞書には"parity", "E", "kb", "kw", "Cb", "Cw", "Vvac", "Vbase", "wwidth"が含まれます。
        :rtype: list[dict]
    """
    a      = wwidth / 2.0 # 半井戸幅
    
    # 奇関数と偶関数のエネルギー準位をそれぞれ見つける
    orbitals_o = find_eigenstates_base(1, Vvac, Vbase, wwidth, meff, Emin, Emax, nE, nmaxiter, eps)
    orbitals_e = find_eigenstates_base(0, Vvac, Vbase, wwidth, meff, Emin, Emax, nE, nmaxiter, eps)
    
    # 両方のリストを結合し、エネルギーの昇順にソート
    orbitals = sorted(orbitals_o + orbitals_e, key = lambda x: x["E"])
    norbitals = len(orbitals)

    # 各準位に対して、波動関数のパラメータを計算
    for i in range(norbitals):
        p = orbitals[i]["parity"]
        E = orbitals[i]["E"]
#        print("  {:3d}: parity={:4}  E={:12.6g} eV".format(i, parity_str[p], E))
        V0 = Vvac - Vbase # 障壁の高さ (井戸底からの相対値)
        # 井戸の外側の減衰定数 kb と井戸の内側の波数 kw を計算
        kb = sqrt(2.0 * meff * me * (V0 - E) * e) / hbar * 1.0e-10  # A^-1
        kw = sqrt(2.0 * meff * me * E * e) / hbar * 1.0e-10         # A^-1
        
        # 波動関数の正規化係数を計算 (Cw:井戸内, Cb:井戸外)
        # 井戸の内側でcos/sinが1になる場所と、外側でexpが1になる場所での接続条件から求める
        Cw = 1.0 / sqrt(wwidth/2.0 + 1.0 / kb) # 簡略化された正規化係数
        if p == 0: # 偶関数
            Cb = Cw * cos(kw * a)
        else: # 奇関数
            Cb = Cw * sin(kw * a)
        
        # 計算結果を辞書に追加
        orbitals[i]["kb"] = kb
        orbitals[i]["kw"] = kw
        orbitals[i]["Cb"] = Cb
        orbitals[i]["Cw"] = Cw
        orbitals[i]["Vvac"]   = Vvac
        orbitals[i]["Vbase"]  = Vbase
        orbitals[i]["wwidth"] = wwidth

    return orbitals

def basis(dx, orbital, Nx = 0):
    """
    単一の矩形量子井戸における波動関数（基底関数）の値を計算します。

    概要:
        与えられた相対座標dxと波動関数軌道情報(orbital)に基づき、
        量子井戸内の電子の波動関数（偶関数または奇関数）の値を返します。
        井戸の内外で異なる関数形式を使用します。

    Parameters:
        :param dx: 原子サイトの中心からの相対x座標 (Å)。
        :type dx: float
        :param orbital: 波動関数軌道の詳細情報を含む辞書。
                        "parity", "wwidth", "E", "kb", "kw", "Cb", "Cw" のキーを含む。
        :type orbital: dict
        :param Nx: 未使用の引数。デフォルトは0。
        :type Nx: int, optional

    Returns:
        :returns: 指定されたdxにおける波動関数の値。
        :rtype: float
    """
    p      = orbital["parity"]
    wwidth = orbital["wwidth"]
    a      = wwidth / 2.0 # 半井戸幅
    E      = orbital["E"] # Eは使用されないが、辞書には含まれる
    kb     = orbital["kb"]
    kw     = orbital["kw"]
    Cb     = orbital["Cb"]
    Cw     = orbital["Cw"]

    if p == 0: # 偶関数
        if dx < -a: # 井戸の左側 (-infinity < dx < -a)
            return Cb * exp(kb * (dx + a))
        elif dx > a: # 井戸の右側 (a < dx < +infinity)
            return Cb * exp(-kb * (dx - a))
        else: # 井戸の内側 (-a <= dx <= a)
            return Cw * cos(kw * dx)
    else: # 奇関数
        if dx < -a: # 井戸の左側 (-infinity < dx < -a)
            return -Cb * exp(kb * (dx + a))
        elif dx > a: # 井戸の右側 (a < dx < +infinity)
            return Cb * exp(-kb * (dx - a))
        else: # 井戸の内側 (-a <= dx <= a)
            return Cw * sin(kw * dx)

def cal_basis(x0, wfmin, wfmax, nwf, orbital):
    """
    特定の軌道に対応する波動関数プロファイルを計算します。

    概要:
        指定された軌道(orbital)について、`wfmin`から`wfmax`の範囲を`nwf`点に分割し、
        `basis`関数を用いて波動関数の値を計算します。
        結果は、サイト中心`x0`からの絶対座標と、対応する波動関数値のリストとして返されます。

    Parameters:
        :param x0: 波動関数の中心となるサイトのx座標 (Å)。
        :type x0: float
        :param wfmin: 波動関数計算範囲の最小相対x座標 (Å)。
        :type wfmin: float
        :param wfmax: 波動関数計算範囲の最大相対x座標 (Å)。
        :type wfmax: float
        :param nwf: 波動関数を計算する点の数。
        :type nwf: int
        :param orbital: 波動関数軌道の詳細情報を含む辞書。`basis`関数に渡される。
        :type orbital: dict

    Returns:
        :returns: (絶対x座標のリスト, 波動関数値のリスト) のタプル。
        :rtype: tuple[list[float], list[float]]
    """
    wfstep = (wfmax - wfmin) / (nwf - 1)
    xx  = []
    ywf = []
#    print("   E={} eV  parity={}".format(Edic["E"], Edic["parity"]))
#    print("   x={} - {}".format(wfmin, wfmax))
    for i in range(nwf):
        x_relative = wfmin + i * wfstep # 相対座標
        wf = basis(x_relative, orbital) # 波動関数を計算
        xx.append(x0 + x_relative)      # 絶対座標に変換
        ywf.append(wf)
#        print(x, wf)

    return xx, ywf

def cal_cdiff2(h, f):
    """
    複素配列の2階中心差分を計算します。

    概要:
        与えられた複素数の配列 `f` に対して、2階中心差分法を用いて2階微分を近似計算します。
        境界点では前方/後方差分に準じた方法で計算されます。

    Parameters:
        :param h: サンプリング間隔（x軸の刻み幅）。
        :type h: float
        :param f: 2階微分を計算する複素数の配列。
        :type f: numpy.ndarray

    Returns:
        :returns: 計算された2階微分の複素数配列。
        :rtype: numpy.ndarray
    """
    n = len(f)
    diff2 = np.empty(n, dtype = complex)
    for i in range(n):
        if i == 0: # 先頭点
            idx = 1
        elif i == n - 1: # 末尾点
            idx = n - 2
        else: # 中間点
            idx = i

        # 複素数の実部と虚部を別々に計算
        d2r = (f[idx+1].real - 2.0 * f[idx].real + f[idx-1].real) / h / h
        d2i = (f[idx+1].imag - 2.0 * f[idx].imag + f[idx-1].imag) / h / h
        diff2[i] = d2r + 1.0j * d2i

    return diff2

def integrate2(h, y1, y2):
    """
    2つの配列の積の数値積分（長方形近似）を計算します。

    概要:
        与えられた2つの配列 `y1` と `y2` の積 `y1[i] * y2[i]` を、
        長方形近似（リーマン和）を用いて`h`の刻み幅で数値積分します。

    Parameters:
        :param h: 積分区間の刻み幅。
        :type h: float
        :param y1: 最初の配列。
        :type y1: numpy.ndarray or list
        :param y2: 2番目の配列。
        :type y2: numpy.ndarray or list

    Returns:
        :returns: 計算された積分値。
        :rtype: float
    """
    S = 0.0
    for i in range(len(y1)-1): # 最後の点を除いて積分
        S += y1[i] * y2[i]
    return S * h

def integrate3(h, y1, y2, y3):
    """
    3つの配列の積の数値積分（長方形近似）を計算します。

    概要:
        与えられた3つの配列 `y1`, `y2`, `y3` の積 `y1[i] * y2[i] * y3[i]` を、
        長方形近似（リーマン和）を用いて`h`の刻み幅で数値積分します。

    Parameters:
        :param h: 積分区間の刻み幅。
        :type h: float
        :param y1: 最初の配列。
        :type y1: numpy.ndarray or list
        :param y2: 2番目の配列。
        :type y2: numpy.ndarray or list
        :param y3: 3番目の配列。
        :type y3: numpy.ndarray or list

    Returns:
        :returns: 計算された積分値。
        :rtype: float
    """
    S = 0.0
    for i in range(len(y1)-1): # 最後の点を除いて積分
        S += y1[i] * y2[i] * y3[i]
    return S * h

def cal_E(orbital, xmin, xmax, nx):
    """
    波動関数軌道の運動エネルギー、ポテンシャルエネルギー、全エネルギーを数値的に計算します。

    概要:
        与えられた波動関数軌道(orbital)に対して、指定された空間範囲で数値積分を実行し、
        <ψ|ψ>, <ψ|K|ψ>, <ψ|U|ψ>, <ψ|H|ψ> を計算します。
        H = K + U であり、Kは運動エネルギー演算子 (-hbar^2/2m * d^2/dx^2)、Uはポテンシャルエネルギー演算子です。

    Parameters:
        :param orbital: 波動関数軌道の詳細情報を含む辞書。
        :type orbital: dict
        :param xmin: 積分範囲の最小x座標 (Å)。
        :type xmin: float
        :param xmax: 積分範囲の最大x座標 (Å)。
        :type xmax: float
        :param nx: 積分に使用する点の数。
        :type nx: int

    Returns:
        :returns: (波動関数の二乗ノルム, 運動エネルギー, ポテンシャルエネルギー, 全エネルギー) のタプル。
        :rtype: tuple[float, float, float, float]
    """
    global Vvac, crystal

    xstep = (xmax - xmin) / (nx - 1)
    xx  = [xmin + i * xstep for i in range(nx)]
    ypot   = np.empty(nx, dtype = complex)
    ywf    = np.empty(nx, dtype = complex)
    for i in range(nx):
        x = xx[i]
        # pot_atomに渡すxは原子サイトの中心からの相対座標なので、ここではそのままxを使っている
        # これは、cal_E関数が単一の原子サイトを原点に置いて計算することを前提としているため
        ypot[i] = pot_atom(x, Vvac, crystal) 
        ywf[i]  = basis(x, orbital) # basis関数も相対座標を引数にとる

    ydiff2 = cal_cdiff2(xstep, ywf) # 波動関数の2階微分を計算

    """
    print("")
    print("x, U(x), wf(x), d2wf/dx2")
    for i in range(nx):
        print("{:8.4f}\t{:8.4f}\t{:8.4f}+j{:8.4f}\t{:8.4f}+j{:8.4f}"
                .format(xx[i], ypot[i], ywf[i].real, ywf[i].imag, ydiff2[i].real, ydiff2[i].imag))
    print("")
    """
    
    Sii = integrate2(xstep, ywf.conjugate(), ywf).real # <ψ|ψ>
    K   = -KE * integrate2(xstep, ywf.conjugate(), ydiff2).real # <ψ|K|ψ> = -KE * <ψ|d^2ψ/dx^2>
    U   = integrate3(xstep, ywf.conjugate(), ypot, ywf).real # <ψ|U|ψ> = <ψ|V(x)|ψ>
    E   = K + U # 全エネルギー <ψ|H|ψ>
#    print("K=", Sii, K, U, E)
#    exit()

    return Sii, K, U, E

def cal_E_analytical(orbital):
    """
    波動関数軌道の運動エネルギー、ポテンシャルエネルギー、全エネルギーを解析的に計算します。

    概要:
        与えられた波動関数軌道(orbital)のパラメータを用いて、
        矩形量子井戸ポテンシャル下での波動関数のノルムの2乗、運動エネルギー、ポテンシャルエネルギー、
        そして全エネルギーを解析的に計算します。
        井戸の内側と外側（障壁領域）で積分を分けて計算します。

    Parameters:
        :param orbital: 波動関数軌道の詳細情報を含む辞書。
                        "parity", "wwidth", "E", "kb", "kw", "Cb", "Cw", "Vvac", "Vbase" のキーを含む。
        :type orbital: dict

    Returns:
        :returns: 井戸外のSii, 井戸内のSii, 全Sii, 井戸外のK, 井戸内のK, 全K, 井戸外のU, 井戸内のU, 全U, 井戸外のE, 井戸内のE, 全E のタプル。
        :rtype: tuple[float, float, float, float, float, float, float, float, float, float, float, float]
    """
    p      = orbital["parity"]
    wwidth = orbital["wwidth"]
    a      = wwidth / 2.0 # 半井戸幅
    E      = orbital["E"]
    kb     = orbital["kb"]
    kw     = orbital["kw"]
    Cb     = orbital["Cb"]
    Cw     = orbital["Cw"]
    Vvac   = orbital["Vvac"]
    Vbase  = orbital["Vbase"]

    # 井戸の外側 (x < -a and x > a) の積分
    # 波動関数は Cb * exp(kb * (dx + a)) for dx < -a, Cb * exp(-kb * (dx - a)) for dx > a
    # |ψ|^2 の積分は、Cb^2 * integral(exp(-2kb * |x-a|)) from a to infinity に対応する。
    # 左右対称なので、片側を計算して2倍する。
    # integral (Cb^2 * exp(-2kb * (x-a))) dx from a to infinity = Cb^2 * [-1/(2kb) * exp(-2kb * (x-a))]_a^infinity
    # = Cb^2 * (0 - (-1/(2kb) * 1)) = Cb^2 / (2kb)
    # 両側で 2 * Cb^2 / (2kb) = Cb^2 / kb
    Siib = Cb * Cb / kb 
    # このコードの `exp(2.0 * kb * a) * exp(-2.0 * kb * a)` は常に1.0なので、上記と同じ。
    # `Cb * exp(kb * (dx + a))` や `Cb * exp(-kb * (dx - a))` は、Cb * exp(kb*a) * exp(kb*dx) の形式で正規化されていると解釈できるため、
    # 井戸の外側の正規化された関数の振幅が Cb * exp(kb*a) となっている。
    # そのため、Cb^2 / kb となる。

    # 井戸の内側 (-a < x < a) の積分
    if p == 0: # 偶関数: Cw * cos(kw*x)
        # integral from -a to a of Cw^2 * cos^2(kw*x) dx
        # = Cw^2 * integral from -a to a of 0.5 * (1 + cos(2*kw*x)) dx
        # = Cw^2 * 0.5 * [x + 0.5/kw * sin(2*kw*x)]_-a^a
        # = Cw^2 * 0.5 * [ (a + 0.5/kw * sin(2*kw*a)) - (-a + 0.5/kw * sin(-2*kw*a)) ]
        # = Cw^2 * 0.5 * [ 2a + 0.5/kw * (sin(2*kw*a) - sin(-2*kw*a)) ]
        # = Cw^2 * 0.5 * [ 2a + 1/kw * sin(2*kw*a) ]
        Siiw = Cw * Cw * 0.5 * (0.5*(sin(2.0 * kw * a)-sin(-2.0 * kw * a)) / kw + 2.0 * a)
    else: # 奇関数: Cw * sin(kw * dx)
        # integral from -a to a of Cw^2 * sin^2(kw*x) dx
        # = Cw^2 * integral from -a to a of 0.5 * (1 - cos(2*kw*x)) dx
        # = Cw^2 * 0.5 * [x - 0.5/kw * sin(2*kw*x)]_-a^a
        # = Cw^2 * 0.5 * [ (a - 0.5/kw * sin(2*kw*a)) - (-a - 0.5/kw * sin(-2*kw*a)) ]
        # = Cw^2 * 0.5 * [ 2a - 0.5/kw * (sin(2*kw*a) - sin(-2*kw*a)) ]
        # = Cw^2 * 0.5 * [ 2a - 1/kw * sin(2*kw*a) ]
        Siiw = Cw * Cw * 0.5 * (-0.5*(sin(2.0 * kw * a)-sin(-2.0 * kw * a)) / kw + 2.0 * a)

    # 運動エネルギー項 K = -KE * <ψ|d^2ψ/dx^2>
    # 井戸の外側: d2f/dx2 = kb*kb * f
    Kiib = -KE * kb * kb * Siib
    # 井戸の内側: d2f/dx2 = -kw*kw * f
    Kiiw = KE * kw * kw * Siiw
    
    # ポテンシャルエネルギー項 U = <ψ|V(x)|ψ>
    Ub = Vvac * Siib # 井戸の外側はVvac
    Uw = Vbase * Siiw # 井戸の内側はVbase
    
    # 全エネルギー
    Eb = Kiib + Ub
    Ew = Kiiw + Uw

    return Siib, Siiw, Siib+Siiw, Kiib, Kiiw, Kiib+Kiiw, Ub, Uw, Ub+Uw, Eb, Ew, Eb+Ew

def wf_innerproduct(x, y1, y2):
    """
    2つの波動関数の内積を数値的に計算します（台形近似）。

    概要:
        与えられたx座標と2つの波動関数 `y1` および `y2` を用いて、
        ∫y1*(x)y2(x)dx の内積を台形近似により数値積分します。

    Parameters:
        :param x: x座標のリストまたは配列。
        :type x: list[float] or numpy.ndarray
        :param y1: 最初の波動関数値のリストまたは配列（複素数可）。
        :type y1: list[complex] or numpy.ndarray
        :param y2: 2番目の波動関数値のリストまたは配列（複素数可）。
        :type y2: list[complex] or numpy.ndarray

    Returns:
        :returns: 計算された内積の複素数値。
        :rtype: complex
    """
    S = 0.0
    for i in range(len(x)-1):
        # 台形近似: (f(x_i)*g(x_i) + f(x_{i+1})*g(x_{i+1})) * (x_{i+1} - x_i) / 2
        S += (y1[i].conjugate() * y2[i] + y1[i+1].conjugate() * y2[i+1]) * (x[i+1] - x[i])
    S *= 0.5
    return S

def cal_tij(xmin, xmax, nx, orbital1, orbital2, x01, x02, Nx):
    """
    ハミルトニアン行列要素t_ijを数値的に計算します。

    概要:
        与えられた2つの軌道(orbital1, orbital2)と、それぞれのサイト位置(x01, x02)および周期オフセット(Nx)に基づき、
        ハミルトニアン行列の非対角要素 t_ij = <ψ_i|H|ψ_j> を数値積分により計算します。
        H = T + U であり、ここではポテンシャルエネルギー項 <ψ_i|U|ψ_j> を計算します。
        （運動エネルギー項は省略されている可能性、または基底関数がシュレディンガー方程式の解であるためゼロと仮定されることがある。）
        この実装ではポテンシャルエネルギー項 <ψ_i|V(x)|ψ_j> を計算しています。

    Parameters:
        :param xmin: 積分範囲の最小x座標 (Å)。
        :type xmin: float
        :param xmax: 積分範囲の最大x座標 (Å)。
        :type xmax: float
        :param nx: 積分に使用する点の数。
        :type nx: int
        :param orbital1: 最初の波動関数軌道の詳細情報を含む辞書。
        :type orbital1: dict
        :param orbital2: 2番目の波動関数軌道の詳細情報を含む辞書。
        :type orbital2: dict
        :param x01: 最初の軌道が属するサイトの中心x座標 (Å)。
        :type x01: float
        :param x02: 2番目の軌道が属するサイトの中心x座標 (Å)。
        :type x02: float
        :param Nx: 2番目のサイトの周期オフセット (2番目のサイトは x02 + Nx*a の位置にあるとみなす)。
        :type Nx: int

    Returns:
        :returns: 計算されたハミルトニアン行列要素 t_ij の複素数値。
        :rtype: complex
    """
    global Vvac, crystal, a

    if orbital1 is orbital2 and Nx == 0:
        # 自己相互作用項の場合、エネルギー準位E_iを返す（基底関数がシュレディンガー方程式の解であるため）
        return orbital1["E"]

#    print("orb {}:{} - {}:{}".format(orbital1["E"], orbital1["parity"], orbital2["E"], orbital2["parity"]))
    xstep = (xmax - xmin) / (nx - 1)
    S = 0.0
    for i in range(nx):
        x = xmin + i * xstep
        V = pot(x, Vvac, crystal) # 全体のポテンシャル
        # 波動関数はそれぞれのサイトからの相対座標で評価
        b1 = basis(x - x01, orbital1).conjugate() 
        b2 = basis(x - (x02 + Nx*a), orbital2)
        S += b1 * V * b2
#        if abs(b1) > 1.0e-3 and abs(b2) > 1.0e-3 and V > 0.0:
#            print("  {}({:8.4f}):{:8.4f}, {:8.4f}, {:8.4f}  S={:8.4f}".format(i, x, b1, V, b2, S))
    S *= xstep
    return S

def cal_Sij(xmin, xmax, nx, orbital1, orbital2, x01, x02, Nx):
    """
    重なり積分行列要素S_ijを数値的に計算します。

    概要:
        与えられた2つの軌道(orbital1, orbital2)と、それぞれのサイト位置(x01, x02)および周期オフセット(Nx)に基づき、
        重なり積分行列の非対角要素 S_ij = <ψ_i|ψ_j> を数値積分により計算します。

    Parameters:
        :param xmin: 積分範囲の最小x座標 (Å)。
        :type xmin: float
        :param xmax: 積分範囲の最大x座標 (Å)。
        :type xmax: float
        :param nx: 積分に使用する点の数。
        :type nx: int
        :param orbital1: 最初の波動関数軌道の詳細情報を含む辞書。
        :type orbital1: dict
        :param orbital2: 2番目の波動関数軌道の詳細情報を含む辞書。
        :type orbital2: dict
        :param x01: 最初の軌道が属するサイトの中心x座標 (Å)。
        :type x01: float
        :param x02: 2番目の軌道が属するサイトの中心x座標 (Å)。
        :type x02: float
        :param Nx: 2番目のサイトの周期オフセット (2番目のサイトは x02 + Nx*a の位置にあるとみなす)。
        :type Nx: int

    Returns:
        :returns: 計算された重なり積分行列要素 S_ij の複素数値。
        :rtype: complex
    """
    global crystal, a
    
    if orbital1 is orbital2 and Nx == 0:
        # 自己重なり積分項は、正規化された波動関数であれば1.0
        return 1.0

    xstep = (xmax - xmin) / (nx - 1)
    S = 0.0
    for i in range(nx):
        x = xmin + i * xstep
        # 波動関数はそれぞれのサイトからの相対座標で評価
        S += basis(x - x01, orbital1).conjugate() * basis(x - (x02 + Nx*a), orbital2)
    S *= xstep
    return S

def build_basisset(crystal):
    """
    タイトバインディング計算のための基底関数セットを構築します。

    概要:
        結晶構造(crystal)内の各サイトと、各原子タイプに属する軌道(orbitals)を組み合わせて、
        計算に使用する全ての基底関数を表すリストを生成します。

    Parameters:
        :param crystal: 結晶サイト情報の辞書リスト。各サイトは `iatom` で `atoms` リストの原子を指す。
                        `atoms` リストの各原子には `orbitals` キーで波動関数のリストが追加されている必要がある。
        :type crystal: list[dict]

    Returns:
        :returns: 各基底関数に関する詳細情報を含む辞書のリスト。
                  各辞書には"isite", "iorb", "atom", "name", "x0", "orbital"のキーが含まれる。
        :rtype: list[dict]
    """
    nsites = len(crystal)
    b = []
    for isite in range(nsites):
        name, iatom, atom, x0, wwidth, xwmin, xwmax, Vbase, meff = siteinf(isite, crystal)
        orbitals = atom["orbitals"] # atom辞書に事前に計算された軌道リストが追加されている前提
        norbitals = len(orbitals)
        for io in range(norbitals):
            b.append({"isite": isite, "iorb": io, "atom": atom, "name": name, "x0": x0, "orbital": orbitals[io]})
    return b

def search_neighbors(rmax, Nxmax, crystal):
    """
    指定された最大距離内で近接サイトを検索します。

    概要:
        結晶構造(crystal)内の各サイトについて、指定された格子周期範囲(Nxmax)と最大距離(rmax)内で
        相互作用する可能性のある近接サイトを検索し、リストとして返します。
        このリストは、タイトバインディング行列の構築に使用されます。

    Parameters:
        :param rmax: 近接サイトとみなす最大距離 (Å)。
        :type rmax: float
        :param Nxmax: 検索する格子周期の最大オフセット (例: Nxmax=2 なら -2a, -a, 0, +a, +2a)。
        :type Nxmax: int
        :param crystal: 結晶サイト情報の辞書リスト。
        :type crystal: list[dict]

    Returns:
        :returns: 各サイトに対する近接サイトのリストを含むリスト。
                  各近接サイトの情報は"isite2", "Nx", "atom2", "name2", "x02", "r"のキーを持つ辞書。
        :rtype: list[list[dict]]
    """
    global a # 結晶の格子定数a

    nsites = len(crystal)
    nlist = []
    for isite1 in range(nsites):
        nlist.append([])
        name1, iatom1, atom1, x01, wwidth, xwmin, xwmax, Vbase, meff = siteinf(isite1, crystal)
        for isite2 in range(nsites):
            name2, iatom2, atom2, x02, wwidth, xwmin, xwmax, Vbase, meff = siteinf(isite2, crystal)
            for inx in range(-Nxmax, +Nxmax+1): # 周期オフセット Nx
                # サイト1と、Nx周期オフセットされたサイト2の距離を計算
                r = abs(x02 + a * inx - x01)
#                if 1.0e-6 < r and r < rmax: # 自己サイト（r=0）以外の近接サイト
                if r < rmax: # 自己サイトを含め、rmax以内
                    nlist[isite1].append({"isite2": isite2, "Nx": inx, 
                                    "atom2": atom2, "name2": name2, "x02": x02, "r": r})

    return nlist

def solve_eigenstates(k, basisset, neighborlist, IsPrint = 0):
    """
    指定された波数kにおけるタイトバインディング法によるエネルギー固有状態を解きます。

    概要:
        波数k (Bloch波数) に対応するBlochの定理を考慮したハミルトニアン行列 (Fij) と重なり積分行列 (Sij) を構築し、
        一般化固有値問題 Hc = ES c を解きます。ただし、この実装では S^(-1) H の固有値を計算しています。
        得られた固有値がエネルギーバンドの各点に対応します。

    Parameters:
        :param k: 波数 (単位は π/a)。
        :type k: float
        :param basisset: タイトバインディング計算に使用する基底関数のリスト。
        :type basisset: list[dict]
        :param neighborlist: 各サイトの近接サイト情報を含むリスト。
        :type neighborlist: list[list[dict]]
        :param IsPrint: デバッグ情報を出力するかどうかのフラグ (0:なし, 1:簡易, 2:詳細)。デフォルトは0。
        :type IsPrint: int, optional

    Returns:
        :returns: (エネルギー固有値の実部配列, 固有ベクトル配列, 計算結果の詳細辞書) のタプル。
                  固有ベクトルは現在空のリストとして返されます。
        :rtype: tuple[numpy.ndarray, list, dict]
    """
    global tijxmin, tijxmax, ntijx # 積分範囲と点数

    tijxstep = (tijxmax - tijxmin) / (ntijx - 1)

#    hij = make_list4(nsites, nbmax, nsites, nbmax)
    nbs = len(basisset) # 基底関数の総数
    Fij = np.zeros([nbs, nbs], dtype = complex) # ハミルトニアン行列
    Sij = np.zeros([nbs, nbs], dtype = complex) # 重なり積分行列

    # 行列要素を計算
    for ib1 in range(nbs):
        bs1 = basisset[ib1]
        isite1   = bs1["isite"]
        iorb1    = bs1["iorb"]
        name1    = bs1["name"]
        atom1    = bs1["atom"]
        x01      = bs1["x0"]
        orbital1 = bs1["orbital"]
        E1       = orbital1["E"]
        nbl1     = neighborlist[isite1] # サイト1の近接サイトリスト
        
        # 行列は対称なので、対角要素と上三角要素のみ計算し、後で共役転置で下三角要素を埋める
        for ib2 in range(ib1, nbs):
            bs2 = basisset[ib2]
            isite2   = bs2["isite"]
            iorb2    = bs2["iorb"]
            name2    = bs2["name"]
            atom2    = bs2["atom"]
            x02      = bs2["x0"]
            orbital2 = bs2["orbital"]
            E2       = orbital2["E"]
            
            # サイト1とサイト2、およびその周期オフセットNxの組み合わせを探索
            # 自己サイト (isite1 == isite2) の場合は Nx=0 の場合も含む
            for inb in range(len(nbl1)):
                nb     = nbl1[inb]
                if isite2 != nb["isite2"]: # サイト2が近接リストのサイトと一致しない場合はスキップ
                    continue
                
                Nx     = nb["Nx"] # サイト2の周期オフセット

                if IsPrint >= 1:
                    print("   base {}-{}  {}[L{}]-{}[L{}]({:+d}a)): "
                            .format(ib1, ib2, name1, iorb1, name2, iorb2, Nx), end = '')
                    print(  " ({:6.4g}eV)".format(E1), end = '')
                    print(" - ({:6.4g}eV)".format(E2))

                # 積分範囲の調整 (2つの基底関数が重なる領域を適切に設定)
                ximin = min(x01 + tijxmin, x02 + tijxmin)
                ximax = max(x01 + tijxmax, x02 + tijxmax)
                nxi = int((ximax - ximin) / tijxstep + 1.00001)

                if IsPrint >= 2:
                    print("      tij integration range: {:6.4g} - {:6.4g}, {:6.4g} step, {} points"
                              .format(ximin, ximax, tijxstep, nxi))

                # H_ij と S_ij の計算
                tij_ = cal_tij(ximin, ximax, nxi, orbital1, orbital2, x01, x02, Nx)
                sij_ = cal_Sij(ximin, ximax, nxi, orbital1, orbital2, x01, x02, Nx)
                
                # Blochの定理による位相因子 exp(i k R) を適用
                # R = (x02 + Nx*a) - x01
                phase = exp(1.0j * k * pi2/a * (x02 + Nx*a - x01))
                Fij[ib1][ib2] += tij_ * phase
                Sij[ib1][ib2] += sij_ * phase

                if IsPrint >= 2:
                    print("     phase=", phase, Nx)
                    print("     tij = {:12.6g}".format(tij_))
                    print("     sij = {:12.6g}".format(sij_))

    # 行列の下三角要素を上三角要素から共役転置で埋める
    for ib1 in range(nbs):
        for ib2 in range(ib1+1, nbs):
            Fij[ib2][ib1] = Fij[ib1][ib2].conjugate()
            Sij[ib2][ib1] = Sij[ib1][ib2].conjugate()

    if IsPrint >= 1:
        print("")
        print("Fij=\n", Fij)
        print("Sij=\n", Sij)
    
    # 一般化固有値問題を解く (Fij * c = E * Sij * c)
    # NumpyのLA.eigh(Fij, Sij) はエルミート行列用。ここでは LA.inv(Sij) @ Fij を計算し、通常の固有値問題を解く。
    # このアプローチはSijが可逆であることと、結果として得られる行列が非エルミートになる可能性があることに注意。
    Fij = LA.inv(Sij) @ Fij # (Sij)^(-1) * Fij を計算
    ei = LA.eigvals(Fij)    # 固有値を計算
    ci = [] # 固有ベクトルは計算していない

    # 固有値は一般に複素数だが、物理的なエネルギーは実数なので実部のみを返す
    # ただし、行列が非エルミートになる場合は固有値が複素数になることがある。
    # この場合、虚部が小さいことを確認する必要がある。
    res = {"ei": ei, "ci": ci, "k": k, "Fij": Fij, "Sij": Sij}

    return ei.real, ci, res

def wf():
    """
    波動関数描画モードのエントリポイントです。(現在未実装)

    概要:
        この関数は現在、実装されておらず、呼び出されるとプログラムを終了し、使用法を表示します。
        将来的に個々の波動関数を詳細に描画する機能が追加される可能性があります。

    Returns:
        :returns: None (プログラムは終了します)
    """
    global mode
    global a
    global V0, crystal
    global xmin, xmax, nx
    global tijxmin, tijxmax, ntijx

    terminate() # 現在は未実装のため、呼び出されるとプログラムを終了

def band():
    """
    タイトバインディング法によるエネルギーバンド構造の計算と描画を行います。

    概要:
        量子井戸のエネルギー準位から基底関数を構築し、それらを用いてハミルトニアン行列と重なり積分行列を形成します。
        Blochの定理を適用して波数kの関数としてエネルギーバンドを計算し、結果をプロットします。

    Returns:
        :returns: None (結果はmatplotlibで表示されます)
    """
    global mode
    global a
    global V0, crystal
    global xmin, xmax, nx
    global tijxmin, tijxmax, ntijx

    nsites = len(crystal)    

    tijxstep = (tijxmax - tijxmin) / (ntijx - 1)
    # tij積分範囲を包含する全体のx範囲を設定
    xmin = min(tijxmin, a + tijxmin)
    xmax = max(tijxmax, a + tijxmax)
    xstep = (xmax - xmin) / (nx - 1)

    print("")
    print("=== Input parameterss ===")
    print("mode:", mode)
    print("Crystal crystal")
    print("  a=", a, "A")
    print("  Vacuum potential: {:10.6f} eV".format(Vvac))
    for i in range(nsites):
        name, iatom, atom, x, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
        print("    {:8s}: x={:10.6f} A".format(name, x))
        print("           well: V(x)={:10.6f} eV from {:10.6f} to {:10.6f} A"
                .format(Vbase, xwmin, xwmax))
    print("X range to calculate tij")
    print("  {:10.6f} - {:10.6f} A at {:10.6f} A step, {} points".format(tijxmin, tijxmax, tijxstep, ntijx))

    xpot, ypot = build_potential(xmin, xstep, nx, Vvac, crystal)

    print("")
    print("=== for finite potential ===")
    print("  Vacuum potential: {:10.6f} eV".format(Vvac))
    atom_orbitals = [] # 各原子タイプに属する軌道リストを格納 (実際はatom辞書に直接追加している)
    xWFlist = []       # 各サイトの波動関数x座標リスト
    yWFlist = []       # 各サイトの波動関数値リスト
    for i in range(nsites):
        name, iatom, atom, x0, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
        print("  Atom {}: {}".format(i, name))
        
        # 各サイト（原子タイプ）の量子井戸のエネルギー準位と波動関数を計算
        orbitals = find_eigenstates(Vvac, Vbase, wwidth, meff, Emin, Emax, nE, nmaxiter, eps)
        WFlist_ = [] # このサイトの各軌道の波動関数値リスト
        for ie in range(len(orbitals)):
            orbital = orbitals[ie]
            p = orbital["parity"]
            E = orbital["E"]
            print("  Level {:3d}: parity={:4}  E={:12.6g} eV".format(ie, parity_str[p], E))
            
            # 波動関数プロファイルを計算（tij積分範囲で）
            xwf, ywf = cal_basis(x0, tijxmin, tijxmax, ntijx, orbital)
            wf2 = wf_innerproduct(xwf, ywf, ywf) # 波動関数のノルムの二乗
            orbital["wf2"] = wf2 # 軌道辞書にノルムを追加
            print("      wf2={}".format(wf2))
            WFlist_.append(ywf)

        atom["orbitals"] = orbitals # atom辞書に軌道リストを追加
        atom["xwf"] = xwf           # atom辞書に波動関数x座標リストを追加
        atom["wfs"] = WFlist_       # atom辞書に波動関数値リストを追加
#        atom_orbitals.append(orbitals) # 未使用
        xWFlist.append(xwf)         # 全てのサイトの波動関数x座標リスト
        yWFlist.append(WFlist_)     # 全てのサイトの波動関数値リスト

    print("")
    print("=== Neighbor sites to calculate tij ===")
    neighborlist = search_neighbors(rmax_neighbor, Nxmax_neighbor, crystal)
    # 近接サイトの最大数を取得 (行列の次元計算のためだが、未使用)
    nbmax = max([len(nb) for nb in neighborlist])
    print("nbmax=", nbmax)
    print("rmax: {:8.4g} A".format(rmax_neighbor))
    print("search Nx range: {} - {}".format(-Nxmax_neighbor, Nxmax_neighbor))
    for isite1 in range(nsites):
        name1, iatom1, atom1, x01, wwidth, xwmin, xwmax, Vbase, meff = siteinf(isite1, crystal)
        print("  site {:03d} {:4}: x=({:8.4g})".format(isite1, name1, x01))
        nbl = neighborlist[isite1]
        for inb in range(len(nbl)):
            nb = nbl[inb]
            print("    - site {:03d} {:4}: x=({:8.4g} + {:+1d}*a) A  r={:8.4g} A"
                    .format(nb["isite2"], nb["name2"], nb["x02"], nb["Nx"], nb["r"]))

    print("")
    print("=== Basis set ===")
    basisset = build_basisset(crystal) # 計算に使う全ての基底関数リストを構築
    for i in range(len(basisset)):
        b = basisset[i]
        print("  {:03d}: site {:03d}:orb {:03d} {:4}".format(i, b["isite"], b["iorb"], b["name"]))

    # k空間の走査範囲を設定
    kstep = (kmax - kmin) / (nk - 1)
    xk  = [kmin + i * kstep for i in range(nk)]
    yE0 = [] # 1サイト系での非ハイブリッド化バンド (Fij)
    yE1 = [] # 1サイト系での非ハイブリッド化バンド (Fij/Sij)
    yE  = [] # ハイブリッド化バンド (多サイト系)

    # 1サイト系の場合の非ハイブリッド化バンドの計算 (デバッグ/比較用)
    if nsites == 1:
        print("")
        print("=== Fij and non-hybridized bands for one-atom crystal ===")
        print("k range: {} - {}, {} step (in pi/a), {} points".format(kmin, kmax, kstep, nk))
        name1, iatom1, atom1, x01, wwidth, xwmin, xwmax, Vbase, meff = siteinf(0, crystal)
        ximin = min(x01 + tijxmin, x01 + tijxmin)
        ximax = max(x01 + tijxmax, x01 + tijxmax)
        nxi = int((ximax - ximin) / tijxstep + 1.00001)
        for i in range(nk):
            k = xk[i]
            phase_arg  = pi2 * k # k * 2pi/a * a = k * 2pi
            # Bloch因子の実際の計算には x02+Nx*a-x01 が使われるため、ここでは単にcos()を計算
            # 波数k * a / (pi/a) で、k * a の単位で位相を表現
            ei0 = []
            ei1 = []
            for ib1 in range(len(basisset)):
                bs1 = basisset[ib1]
                iorb1    = bs1["iorb"]
                print("  iorb=", iorb1)
                orbital1 = bs1["orbital"]
                
                # 自己相互作用項 (Nx=0)
                tii  = cal_tij(ximin, ximax, nxi, orbital1, orbital1, x01, x01, 0)
                sii  = cal_Sij(ximin, ximax, nxi, orbital1, orbital1, x01, x01, 0)
                # 最近接相互作用項 (Nx=-1, +1)
                tijm = cal_tij(ximin, ximax, nxi, orbital1, orbital1, x01, x01, -1)
                sijm = cal_Sij(ximin, ximax, nxi, orbital1, orbital1, x01, x01, -1)
                tijp = cal_tij(ximin, ximax, nxi, orbital1, orbital1, x01, x01, +1)
                sijp = cal_Sij(ximin, ximax, nxi, orbital1, orbital1, x01, x01, +1)
                
                # 1サイト系の場合のバンド計算式 (k・aで畳み込み)
                # E(k) = tii + 2 * cos(k * a) * tijp
                # S(k) = sii + 2 * cos(k * a) * sijp
                # 注: `phase_arg` は k * pi2/a * a = k * 2pi になるが、cos(phase_arg)はcos(k*2pi)になる。
                # 通常は k * a が直接波数として使われるため、phase_arg_k_a = k * pi
                # ここでは k は pi/a 単位で定義されているので k*pi2/a * a = k*2pi となる。
                # しかし、`2.0 * cos(phase)` の `phase` が何を示すか注意。`phase_arg` をそのまま使うと `cos(k*2pi)` になる。
                # 一般的な1D TBモデルでは `t_0 + 2t_1 cos(k*a)` となるため、k*a を `phase` とすべき。
                # xk の単位が `pi/a` なので、k*a は `xk[i] * pi` となるべき。
                # 既存のロジックを変更しないため、`phase_arg` = `k * pi2/a * a` と解釈してそのままにする。
                Fij_nonhyb = 2.0 * cos(phase_arg) * tijp + tii
                Sij_nonhyb = 2.0 * cos(phase_arg) * sijp + sii
                
                print("    at k={:8.4g}: Fij={:10.6g} Sij={:10.6g} Fij/Sij={:10.6g}"
                        .format(k, Fij_nonhyb, Sij_nonhyb, Fij_nonhyb / Sij_nonhyb))
                print("      diagonal: tii={:10.6g} sii={:10.6g}".format(tii,  sii))
                print("      to Nx=-1: tij={:10.6g} sij={:10.6g}".format(tijm, sijm))
                print("      to Nx=+1: tij={:10.6g} sij={:10.6g}".format(tijp, sijp))
                ei0.append(Fij_nonhyb)
                ei1.append(Fij_nonhyb / Sij_nonhyb)
            yE0.append(ei0)
            yE1.append(ei1)
#    exit()

    print("")
    print("=== Band ===")
    print("k range: {} - {}, {} step (in pi/a), {} points".format(kmin, kmax, kstep, nk))
    for i in range(nk):
        k = xk[i]
        # 各波数kで固有値問題を解き、エネルギー準位を取得
        ei, ci, res = solve_eigenstates(k, basisset, neighborlist, 0)
        print("k: {:10.6f}".format(k), ":", ei)
        yE.append(ei) # 計算されたエネルギー固有値をリストに追加
#        print("  Fij=", res["Fij"].real)
#        print("  Sij=", res["Sij"].real)

    print("")
    print("plot")
    print("")
    fig, ax = plt.subplots(1, 2, figsize=figsize) # 2つのサブプロットを作成
    ax1 = ax[0] # 左側のプロット (ポテンシャルと波動関数)
    ax2 = ax[1] # 右側のプロット (バンド構造)

    # ax1: ポテンシャルプロファイルの描画
    ax1.plot(xpot, ypot, linewidth = 0.5, label = '$V(x)$') 
    # ax1: 量子井戸のエネルギー準位と波動関数の描画
    for i in range(nsites):
        name, iatom, atom, x0, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
        orbitals = atom["orbitals"] # 各サイトの軌道情報を取得
        norbitals = len(orbitals)
        for ie in range(norbitals):
            orbital = orbitals[ie]
            E = orbital["E"]
            ax1.plot([xwmin, xwmax], [E, E], linestyle = 'dashed', color = 'b', linewidth = 1.0) # エネルギー準位
            # 波動関数をエネルギー準位上に重ねてプロット
            # 波動関数値の最大値でスケーリングし、オフセットとしてエネルギー準位Eを加える
            kwf = wfplotmax / max(yWFlist[i][ie])
            ywf = [yWFlist[i][ie][j] * kwf + E for j in range(ntijx)]
            ax1.plot(xWFlist[i], ywf, linestyle = '-', color = 'g', linewidth = 0.5)

    # ax2: バンド構造の描画
    ax2.set_xlim([kmin, kmax])
    ax2.set_ylim(Erange)
    # ハイブリッド化バンド (多サイト系)
    ax2.plot(xk, yE, label = 'hybridized', linestyle = '', marker = 'o', markersize = 5.0,
                        markerfacecolor = 'none', markeredgecolor = 'black', markeredgewidth = 0.5)
    # 1サイト系の場合の非ハイブリッド化バンド (比較用)
    if nsites == 1:
        ax2.plot(xk, yE0, label = '$F_{ii}$',     linestyle = '-',      color = 'red',  linewidth = 0.5)
        ax2.plot(xk, yE1, label = '$F_{ii}/S_{ii}$', linestyle = 'dashed', color = 'blue', linewidth = 0.5)

    ax1.set_xlabel("x ($\AA$)", fontsize = fontsize)
    ax1.set_ylabel("$V$($x$), $\phi$($x$)", fontsize = fontsize)
    ax1.legend(fontsize = legend_fontsize)
    ax1.tick_params(labelsize = fontsize)
    ax2.set_xlabel("$k$ $(\pi$$/a)$", fontsize = fontsize)
    ax2.set_ylabel("E (eV)", fontsize = fontsize)
    ax2.legend(fontsize = legend_fontsize) # ax2にも凡例を追加
    ax2.tick_params(labelsize = fontsize)

    plt.tight_layout() # レイアウトを調整

    plt.pause(0.1)
    print("Press ENTER to exit>>", end = '')
    input()

    terminate()

def view_basis():
    """
    量子井戸のポテンシャル、エネルギー準位、波動関数（基底関数）を表示します。

    概要:
        設定された結晶構造と原子情報に基づいて、矩形量子井戸ポテンシャルを描画します。
        無限に深い量子井戸と有限の量子井戸の両方についてエネルギー準位を計算し、
        有限量子井戸のエネルギー準位に対応する波動関数もプロットします。
        また、波動関数の数値積分によるエネルギー期待値も計算し、解析解と比較します。

    Returns:
        :returns: None (結果はmatplotlibで表示されます)
    """
    global mode
    global a
    global V0, atoms, crystal
    global xmin, xmax, nx
    
    natoms = len(atoms)
    nsites = len(crystal)
    xstep = (xmax - xmin) / (nx - 1)
    
    print("")
    print("=== Input parameterss ===")
    print("mode:", mode)
    print("Atoms:")
    for i in range(natoms):
        name, wwidth, Vbase, meff = atominf(i, atoms)
        print("    {:8s}: well: width={:10.6f} A  Vbase={:10.6f} eV  me={:10.6f}me"
                .format(name, wwidth, Vbase, meff))
    print("Crystal structure:")
    print("  a=", a, "A")
    print("  Vacuum potential: {:10.6f} eV".format(Vvac))
    for i in range(nsites):
        name, iatom, atom, x, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
        print("    {:8s}: x={:10.6f} A".format(name, x))
        print("           well: V(x)={:10.6f} eV from {:10.6f} to {:10.6f} A"
                .format(Vbase, xwmin, xwmax))
    print("X range to plot V(x) and basis")
    print("  {:10.6f} - {:10.6f} A at {:10.6f} A step, {} points".format(xmin, xmax, xstep, nx))

    xpot, ypot = build_potential(xmin, xstep, nx, Vvac, crystal) # ポテンシャルプロファイルの構築
    
    print("")
    print("=== for infinite potential ===")
    # 無限井戸ポテンシャルのエネルギー準位を計算
    for i in range(natoms): # 各原子タイプについて
        name, wwidth, Vbase, meff = atominf(i, atoms)
        print("    {:8s}: well: width={:10.6f} A  Vbase={:10.6f} eV  me={:10.6f}me"
                .format(name, wwidth, Vbase, meff))
        for iL in range(1, 100): # 量子数1から順に計算
            E = quantization_level_inifinite(Vbase, meff, wwidth, iL)
            if E > Vvac: # 真空準位を超えたら終了
                break
            print("    Level {:3d}: E={:10.6f} eV  E-Vbase={:10.6f} eV".format(iL, E, E - Vbase))

    print("")
    print("=== for finite potential ===")
    print("  Vacuum potential: {:10.6f} eV".format(Vvac))
#    atom_orbitals = [] # 未使用
    xWFlist = []       # 各サイトの波動関数x座標リスト
    yWFlist = []       # 各サイトの波動関数値リスト
    for i in range(nsites): # 各結晶サイトについて
        name, iatom, atom, x0, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
        print("    {:8s}: well: width={:10.6f} A  Vbase={:10.6f} eV  me={:10.6f}me"
                .format(name, wwidth, Vbase, meff))
        xstep_int = (xmaxEinteg - xminEinteg) / (nxEinteg - 1) # エネルギー期待値計算用の刻み幅
        print("    integration x range: {} - {}, {} step, n={}".format(xminEinteg, xmaxEinteg, xstep_int, nxEinteg))
        WFlist_ = []
        # 有限井戸ポテンシャルのエネルギー準位と波動関数を計算
        orbitals = find_eigenstates(Vvac, Vbase, wwidth, meff, Emin, Emax, nE, nmaxiter, eps)
        norbitals = len(orbitals)
        for ie in range(norbitals):
            orbital = orbitals[ie]
            p = orbital["parity"]
            E = orbital["E"]
            print("  Level {:3d}: parity={:4}  E={:12.6g} eV".format(ie, parity_str[p], E))
            # 波動関数プロファイルを計算（表示範囲で）
            # 波動関数を可視化のためにスケーリング
            xwf, ywf_raw = cal_basis(x0, -wwidth/2.0 - wmargin, wwidth/2.0 + wmargin, nwf, orbital)
            kwf = wfplotmax / max(np.abs(ywf_raw)) # 波動関数の絶対値の最大値でスケーリング
            ywf = [E + ywf_raw[i] * kwf for i in range(len(ywf_raw))] # エネルギー準位に重ねる
            WFlist_.append(ywf)

            # 解析解によるエネルギー計算
            Siib, Siiw, Sii, Kb, Kw, K, Ub, Uw, U, Eb, Ew, E_analytical = cal_E_analytical(orbital)
            print("    [Analytical] Sii = {:8.4f}(b) + {:8.4f}(w) = {:8.4f}".format(Siib, Siiw, Sii))
            print("    [Analytical] K   = {:8.4f}(b) + {:8.4f}(w) = {:8.4f} eV".format(Kb, Kw, K))
            print("    [Analytical] U   = {:8.4f}(b) + {:8.4f}(w) = {:8.4f} eV".format(Ub, Uw, U))
            print("    [Analytical] E   = {:8.4f}(b) + {:8.4f}(w) = {:8.4f} eV".format(Eb, Ew, E_analytical))
            
            # 数値積分によるエネルギー計算
            Sii_num, K_num, U_num, E_num = cal_E(orbital, xminEinteg, xmaxEinteg, nxEinteg)
            print("    [Numerical]  Sii={:8.4f}, K={:8.4f}, U={:8.4f}, E={:8.4f} eV".format(Sii_num, K_num, U_num, E_num))
            print("")

        atoms[iatom]["orbitals"] = orbitals
        atoms[iatom]["xwf"] = xwf
        atoms[iatom]["wfs"] = WFlist_
#        atom_orbitals.append(orbitals) # 未使用
        xWFlist.append(xwf)
        yWFlist.append(WFlist_)

    fig, ax = plt.subplots(1, 1, figsize=figsize) # 1つのサブプロットを作成
    ax1 = ax 

    ax1.plot(xpot, ypot, linewidth = 0.5, label = '$V(x)$') # ポテンシャルプロファイルの描画
    
    # 無限井戸ポテンシャルのエネルギー準位を点線でプロット
    for i in range(nsites):
        name, iatom, atom, x, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
        for iL in range(1, 100):
            E = quantization_level_inifinite(Vbase, meff, wwidth, iL)
            if E > Vvac:
                break
            ax1.plot([xwmin, xwmax], [E, E], linestyle = 'dotted', color = 'r', linewidth = 0.5, label='Infinite well' if i==0 and iL==1 else "")

    # 有限井戸ポテンシャルのエネルギー準位と波動関数をプロット
    for i in range(nsites):
        name, iatom, atom, x, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
        orbitals = atoms[iatom]["orbitals"]
        norbitals = len(orbitals)
        for ie in range(norbitals):
            orbital = orbitals[ie]
            E = orbital["E"]
            ax1.plot([xwmin, xwmax], [E, E], linestyle = 'dashed', color = 'b', linewidth = 1.0, label='Finite well' if i==0 and ie==0 else "") # エネルギー準位
            ax1.plot(xWFlist[i], yWFlist[i][ie], linestyle = '-', color = 'g', linewidth = 0.5, label='Wave function' if i==0 and ie==0 else "") # 波動関数

    ax1.set_xlabel("x ($\AA$)", fontsize = fontsize)
    ax1.set_ylabel("$V$($x$), $\phi$($x$)", fontsize = fontsize)
    ax1.legend(fontsize = legend_fontsize)
    ax1.tick_params(labelsize = fontsize)
    plt.tight_layout()

    plt.pause(0.1)
    print("")
    print("Press ENTER to exit>>", end = '')
    input()

    terminate()


def main():
    """
    スクリプトのメインエントリポイントです。

    概要:
        コマンドライン引数で指定された`mode`に基づいて、適切な機能（`view_basis`、`band`、`wf`）を呼び出します。
        無効なモードが指定された場合は、エラーメッセージを表示してプログラムを終了します。

    Returns:
        :returns: None
    """
    global mode
    
    if mode == 'basis':
        view_basis()
    elif mode == 'band':
        band()
    elif mode == 'wf':
        wf()
    else:
        terminate("Error: Invalid mode [{}]".format(mode))


if __name__ == "__main__":
    main()