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

kronig_penney_refactored.py をダウンロード

kronig_penney_refactored.py
kronig_penney_refactored.py
   1import argparse
   2import sys
   3import traceback
   4from types import SimpleNamespace
   5
   6import matplotlib.pyplot as plt
   7import numpy as np
   8from numpy import cos, cosh, exp, sin, sinh, sqrt
   9import numpy.linalg as LA
  10
  11
  12"""
  13Kronig-Penneyモデルによる1次元バンド計算スクリプト。
  14
  15概要:
  16    矩形ポテンシャルを用いたKronig-Penneyモデルに基づき、電子のエネルギーバンド構造、
  17    波動関数、およびデルタ関数グラフを計算しプロットします。
  18詳細説明:
  19    このスクリプトは、以下の3つの主要なモードで動作します。
  20    1. delta(E) グラフのプロット (graphモード)
  21    2. バンド構造のプロット (bandモード)
  22    3. 波動関数のプロット (wfモード)
  23    各モードでは、ポテンシャルの幅、高さ、格子定数などのパラメータを調整できます。
  24主な機能:
  25    - Kronig-Penney方程式の解に基づくエネルギーデルタ関数の計算。
  26    - エネルギーバンドの計算とプロット。
  27    - 選択されたエネルギー準位に対応する波動関数の計算とプロット。
  28関連リンク:
  29    kronig_penney_refactored_usage
  30"""
  31
  32
  33# =============================================================================
  34# Constants
  35# =============================================================================
  36PI = 3.14159265358979323846
  37PI2 = 2.0 * PI
  38
  39H = 6.6260755e-34
  40HBAR = 1.05459e-34
  41ELEMENTARY_CHARGE = 1.60218e-19
  42ELECTRON_MASS = 9.1093897e-31
  43
  44DEFAULT = SimpleNamespace(
  45    mode="graph",
  46    a=5.4064,          # A
  47    bwidth=0.5,        # A
  48    bpot=10.0,         # eV
  49    kg=0.0,            # pi/a
  50    emin=0.0,          # eV
  51    emax=9.5,          # eV
  52    nE=51,
  53    nEsearch=51,
  54    eps=1.0e-8,
  55    nmaxiter=100,
  56    dump=0.0,
  57    kmin=-0.5,         # pi/a
  58    kmax=0.5,          # pi/a
  59    nk=21,
  60    erange_min=0.0,    # eV
  61    erange_max=10.0,   # eV
  62    nMaxLevel=15,
  63    xwmin=0.0,         # A
  64    xwmax=None,        # A; None -> 3.0 * a
  65    nxw=101,
  66    kw=0.0,            # pi/a
  67    iLevel=0,
  68    figsize_w=6.0,
  69    figsize_h=8.0,
  70    wf_figsize_w=16.0,
  71    wf_figsize_h=4.0,
  72    fontsize=12,
  73    legend_fontsize=8,
  74    show=1,
  75    save=0,
  76    output="kronig_penney.png",
  77)
  78
  79
  80# =============================================================================
  81# argparse
  82# =============================================================================
  83def parse_args() -> argparse.Namespace:
  84    """
  85    概要:
  86        コマンドライン引数を解析します。
  87    詳細説明:
  88        Kronig-Penneyモデルの計算およびプロットに必要なパラメータをコマンドラインから受け取ります。
  89        従来の引数形式もサポートしています。xwmax が指定されない場合は、格子定数 a の3倍がデフォルト値となります。
  90    引数:
  91        なし。
  92    戻り値:
  93        :returns: 解析されたコマンドライン引数を格納したNamespaceオブジェクト。
  94        :rtype: argparse.Namespace
  95    """
  96    parser = argparse.ArgumentParser(
  97        description="Kronig-Penney model: graph, band, and wavefunction plotter.",
  98        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
  99    )
 100
 101    # 従来の位置引数形式も維持するため、mode と legacy_args を受ける。
 102    parser.add_argument("mode", nargs="?", default=DEFAULT.mode, choices=["graph", "band", "wf"])
 103    parser.add_argument("legacy_args", nargs="*")
 104
 105    parser.add_argument("--a", type=float, default=DEFAULT.a)
 106    parser.add_argument("--bwidth", type=float, default=DEFAULT.bwidth)
 107    parser.add_argument("--bpot", type=float, default=DEFAULT.bpot)
 108
 109    parser.add_argument("--kg", type=float, default=DEFAULT.kg)
 110    parser.add_argument("--Emin", dest="emin", type=float, default=DEFAULT.emin)
 111    parser.add_argument("--Emax", dest="emax", type=float, default=DEFAULT.emax)
 112    parser.add_argument("--nE", type=int, default=DEFAULT.nE)
 113
 114    parser.add_argument("--nEsearch", type=int, default=DEFAULT.nEsearch)
 115    parser.add_argument("--eps", type=float, default=DEFAULT.eps)
 116    parser.add_argument("--nmaxiter", type=int, default=DEFAULT.nmaxiter)
 117    parser.add_argument("--dump", type=float, default=DEFAULT.dump)
 118
 119    parser.add_argument("--kmin", type=float, default=DEFAULT.kmin)
 120    parser.add_argument("--kmax", type=float, default=DEFAULT.kmax)
 121    parser.add_argument("--nk", type=int, default=DEFAULT.nk)
 122
 123    parser.add_argument("--ErangeMin", dest="erange_min", type=float, default=DEFAULT.erange_min)
 124    parser.add_argument("--ErangeMax", dest="erange_max", type=float, default=DEFAULT.erange_max)
 125    parser.add_argument("--nMaxLevel", type=int, default=DEFAULT.nMaxLevel)
 126
 127    parser.add_argument("--kw", type=float, default=DEFAULT.kw)
 128    parser.add_argument("--iLevel", type=int, default=DEFAULT.iLevel)
 129    parser.add_argument("--xwmin", type=float, default=DEFAULT.xwmin)
 130    parser.add_argument("--xwmax", type=float, default=DEFAULT.xwmax)
 131    parser.add_argument("--nxw", type=int, default=DEFAULT.nxw)
 132
 133    parser.add_argument("--figsizeW", dest="figsize_w", type=float, default=DEFAULT.figsize_w)
 134    parser.add_argument("--figsizeH", dest="figsize_h", type=float, default=DEFAULT.figsize_h)
 135    parser.add_argument("--wfFigsizeW", dest="wf_figsize_w", type=float, default=DEFAULT.wf_figsize_w)
 136    parser.add_argument("--wfFigsizeH", dest="wf_figsize_h", type=float, default=DEFAULT.wf_figsize_h)
 137    parser.add_argument("--fontsize", type=int, default=DEFAULT.fontsize)
 138    parser.add_argument("--legendFontsize", dest="legend_fontsize", type=int, default=DEFAULT.legend_fontsize)
 139
 140    parser.add_argument("--show", type=int, choices=[0, 1], default=DEFAULT.show)
 141    parser.add_argument("--save", type=int, choices=[0, 1], default=DEFAULT.save)
 142    parser.add_argument("--output", type=str, default=DEFAULT.output)
 143
 144    args = parser.parse_args()
 145    apply_legacy_args(args)
 146    if args.xwmax is None:
 147        args.xwmax = 3.0 * args.a
 148    args.nEsearch = args.nE if args.nEsearch == DEFAULT.nEsearch else args.nEsearch
 149    return args
 150
 151
 152def apply_legacy_args(args: argparse.Namespace) -> None:
 153    """
 154    概要:
 155        旧形式のコマンドライン引数を適用します。
 156    詳細説明:
 157        modeの後に続く位置引数として指定された値を解析し、argsオブジェクトの対応する属性に設定します。
 158        これにより、スクリプトの古いバージョンとの互換性が保たれます。
 159    引数:
 160        :param args: コマンドライン引数を格納したNamespaceオブジェクト。
 161        :type args: argparse.Namespace
 162    戻り値:
 163        なし。
 164    """
 165    values = args.legacy_args
 166    if not values:
 167        return
 168
 169    if len(values) >= 1:
 170        args.a = float(values[0])
 171    if len(values) >= 2:
 172        args.bwidth = float(values[1])
 173    if len(values) >= 3:
 174        args.bpot = float(values[2])
 175
 176    if args.mode == "graph":
 177        if len(values) >= 4:
 178            args.kg = float(values[3])
 179        if len(values) >= 5:
 180            args.emin = float(values[4])
 181        if len(values) >= 6:
 182            args.emax = float(values[5])
 183        if len(values) >= 7:
 184            args.nE = int(values[6])
 185    elif args.mode == "band":
 186        if len(values) >= 4:
 187            args.kmin = float(values[3])
 188        if len(values) >= 5:
 189            args.kmax = float(values[4])
 190        if len(values) >= 6:
 191            args.nk = int(values[5])
 192    elif args.mode == "wf":
 193        if len(values) >= 4:
 194            args.kw = float(values[3])
 195        if len(values) >= 5:
 196            args.iLevel = int(values[4])
 197        if len(values) >= 6:
 198            args.xwmin = float(values[5])
 199        if len(values) >= 7:
 200            args.xwmax = float(values[6])
 201        if len(values) >= 8:
 202            args.nxw = int(values[7])
 203
 204
 205# =============================================================================
 206# utility functions
 207# =============================================================================
 208def round01(x: float, a: float) -> tuple[float, int]:
 209    """
 210    概要:
 211        値を周期 a で区間 [0, a) に丸め、周期数を返します。
 212    詳細説明:
 213        入力値 x を周期 a で正規化し、周期内の値 x0 と、何周期分シフトしたかを示す整数 n を計算します。
 214        x >= 0 の場合は n = floor(x / a)、x < 0 の場合は n = floor(x / a) - 1 となります。
 215        結果として x0 は常に [0, a) の範囲に収まります。
 216    引数:
 217        :param x: 丸める対象の浮動小数点数。
 218        :type x: float
 219        :param a: 周期を表す浮動小数点数。
 220        :type a: float
 221    戻り値:
 222        :returns: 周期内の値 x0 と周期数 n のタプル。
 223        :rtype: tuple
 224    """
 225    if x >= 0.0:
 226        n = int(x / a)
 227    else:
 228        n = int(x / a) - 1
 229    x0 = x - n * a
 230    return x0, n
 231
 232
 233def validate_args(args: argparse.Namespace) -> None:
 234    """
 235    概要:
 236        コマンドライン引数の値を検証します。
 237    詳細説明:
 238        各引数が物理的に妥当な範囲内にあるか、または計算に必要な最小値を満たしているかを確認します。
 239        例えば、幅やポテンシャル高さが正であること、格子定数 a がポテンシャル幅 bwidth より大きいことなどをチェックします。
 240    引数:
 241        :param args: コマンドライン引数を格納したNamespaceオブジェクト。
 242        :type args: argparse.Namespace
 243    戻り値:
 244        なし。
 245    例外:
 246        :raises ValueError: 引数の値が不正な場合に発生します。
 247    """
 248    if args.bwidth <= 0.0:
 249        raise ValueError("--bwidth must be positive.")
 250    if args.a <= args.bwidth:
 251        raise ValueError("--a must be larger than --bwidth.")
 252    if args.bpot <= 0.0:
 253        raise ValueError("--bpot must be positive.")
 254    if args.nE < 2:
 255        raise ValueError("--nE must be >= 2.")
 256    if args.nEsearch < 2:
 257        raise ValueError("--nEsearch must be >= 2.")
 258    if args.nk < 2:
 259        raise ValueError("--nk must be >= 2.")
 260    if args.nxw < 2:
 261        raise ValueError("--nxw must be >= 2.")
 262    if args.mode == "wf" and args.iLevel < 0:
 263        raise ValueError("--iLevel must be >= 0.")
 264
 265
 266def create_context(args: argparse.Namespace) -> SimpleNamespace:
 267    """
 268    概要:
 269        解析された引数からコンテキストオブジェクトを作成します。
 270    詳細説明:
 271        argparse.Namespaceオブジェクトの引数を元に、計算やプロットに必要な全てのパラメータを
 272        SimpleNamespaceオブジェクトとして集約します。これにより、関数の引数渡しが簡潔になります。
 273        特に、ポテンシャルの幅 b と井戸の幅 w を計算して格納します。
 274    引数:
 275        :param args: コマンドライン引数を格納したNamespaceオブジェクト。
 276        :type args: argparse.Namespace
 277    戻り値:
 278        :returns: 計算コンテキストを格納したSimpleNamespaceオブジェクト。
 279        :rtype: types.SimpleNamespace
 280    """
 281    return SimpleNamespace(
 282        mode=args.mode,
 283        a=args.a,
 284        bwidth=args.bwidth,
 285        bpot=args.bpot,
 286        w=args.a - args.bwidth,
 287        b=args.bwidth,
 288        V0=args.bpot,
 289        kg=args.kg,
 290        emin=args.emin,
 291        emax=args.emax,
 292        nE=args.nE,
 293        nEsearch=args.nEsearch,
 294        eps=args.eps,
 295        nmaxiter=args.nmaxiter,
 296        dump=args.dump,
 297        kmin=args.kmin,
 298        kmax=args.kmax,
 299        nk=args.nk,
 300        erange=(args.erange_min, args.erange_max),
 301        nMaxLevel=args.nMaxLevel,
 302        kw=args.kw,
 303        iLevel=args.iLevel,
 304        xwmin=args.xwmin,
 305        xwmax=args.xwmax,
 306        nxw=args.nxw,
 307        figsize=(args.figsize_w, args.figsize_h),
 308        wf_figsize=(args.wf_figsize_w, args.wf_figsize_h),
 309        fontsize=args.fontsize,
 310        legend_fontsize=args.legend_fontsize,
 311        show=args.show,
 312        save=args.save,
 313        output=args.output,
 314    )
 315
 316
 317def save_show_close(fig: plt.Figure, ctx: SimpleNamespace) -> None:
 318    """
 319    概要:
 320        Matplotlibの図を保存、表示、閉じます。
 321    詳細説明:
 322        コンテキストオブジェクトの設定 (ctx.save, ctx.show, ctx.output) に応じて、
 323        生成されたMatplotlibの図 (Figureオブジェクト) をファイルに保存し、画面に表示し、
 324        最終的に閉じます。
 325    引数:
 326        :param fig: 保存、表示、閉じる対象のMatplotlib Figureオブジェクト。
 327        :type fig: matplotlib.pyplot.Figure
 328        :param ctx: 計算コンテキストを格納したSimpleNamespaceオブジェクト。
 329        :type ctx: types.SimpleNamespace
 330    戻り値:
 331        なし。
 332    """
 333    if ctx.save:
 334        fig.savefig(ctx.output, dpi=300, bbox_inches="tight")
 335        print(f"Saved figure: {ctx.output}")
 336    if ctx.show:
 337        plt.show()
 338    plt.close(fig)
 339
 340
 341# =============================================================================
 342# core functions
 343# =============================================================================
 344def compute_potential_value(x: float, ctx: SimpleNamespace) -> float:
 345    """
 346    概要:
 347        指定された位置 x におけるポテンシャルの値を計算します。
 348    詳細説明:
 349        Kronig-Penneyモデルにおける矩形ポテンシャル関数を定義します。
 350        x が周期 a のセル内でポテンシャル幅 bwidth の範囲内にある場合、ポテンシャルの高さ ctx.bpot を返し、
 351        それ以外の場合は 0.0 を返します。
 352    引数:
 353        :param x: ポテンシャルの値を評価するx座標。
 354        :type x: float
 355        :param ctx: 計算コンテキストを格納したSimpleNamespaceオブジェクト。
 356        :type ctx: types.SimpleNamespace
 357    戻り値:
 358        :returns: xにおけるポテンシャルの値。
 359        :rtype: float
 360    """
 361    xred, _ = round01(x, ctx.a)
 362    if ctx.a - ctx.bwidth <= xred < ctx.a:
 363        return ctx.bpot
 364    return 0.0
 365
 366
 367def compute_potential_profile(xmin: float, xstep: float, n: int, ctx: SimpleNamespace) -> tuple[np.ndarray, np.ndarray]:
 368    """
 369    概要:
 370        指定された範囲とステップでポテンシャルプロファイルを計算します。
 371    詳細説明:
 372        xminから始まり、xstep間隔でn個の点におけるポテンシャル値を計算し、x座標配列とポテンシャル値配列を生成します。
 373        各点のポテンシャル値は compute_potential_value を呼び出して決定されます。
 374    引数:
 375        :param xmin: プロファイルの開始x座標。
 376        :type xmin: float
 377        :param xstep: 各点間のxステップサイズ。
 378        :type xstep: float
 379        :param n: 計算する点の総数。
 380        :type n: int
 381        :param ctx: 計算コンテキストを格納したSimpleNamespaceオブジェクト。
 382        :type ctx: types.SimpleNamespace
 383    戻り値:
 384        :returns: x座標のnumpy.ndarrayと対応するポテンシャル値のnumpy.ndarrayのタプル。
 385        :rtype: tuple
 386    """
 387    xpot = xmin + np.arange(n) * xstep
 388    ypot = np.array([compute_potential_value(x, ctx) for x in xpot], dtype=float)
 389    return xpot, ypot
 390
 391
 392def compute_delta(E: float, k: float, w: float, b: float, V0: float) -> float:
 393    """
 394    概要:
 395        Kronig-Penneyモデルのデルタ関数 (delta(E)) の値を計算します。
 396    詳細説明:
 397        与えられたエネルギー E、波数 k、井戸の幅 w、ポテンシャル障壁の幅 b、
 398        およびポテンシャルの高さ V0 に基づいて、Kronig-Penneyモデルのデルタ関数を評価します。
 399        この関数は、電子の透過・反射特性を表し、delta(E) の値が -1 から 1 の間に収まるエネルギーが
 400        許容エネルギーバンドに対応します。
 401        計算には物理定数 HBAR, ELECTRON_MASS, ELEMENTARY_CHARGE を使用します。
 402    引数:
 403        :param E: 電子のエネルギー (eV)。
 404        :type E: float
 405        :param k: 電子の波数 (pi/a 単位)。
 406        :type k: float
 407        :param w: 井戸の幅 (A)。
 408        :type w: float
 409        :param b: ポテンシャル障壁の幅 (A)。
 410        :type b: float
 411        :param V0: ポテンシャルの高さ (eV)。
 412        :type V0: float
 413    戻り値:
 414        :returns: デルタ関数の値。
 415        :rtype: float
 416    """
 417    alpha = sqrt(2.0 * ELECTRON_MASS * E * ELEMENTARY_CHARGE) / HBAR
 418    beta = sqrt(2.0 * ELECTRON_MASS * (V0 - E) * ELEMENTARY_CHARGE) / HBAR
 419    ka = k * PI2
 420    alphaw = alpha * w * 1.0e-10
 421    betab = beta * b * 1.0e-10
 422
 423    delta = (
 424        (beta * beta - alpha * alpha) / 2.0 / alpha / beta * sin(alphaw) * sinh(betab)
 425        + cos(alphaw) * cosh(betab)
 426        - cos(ka)
 427    )
 428    return float(delta)
 429
 430
 431def compute_boundary_matrix(k: float, E: float, w: float, b: float, V0: float) -> np.ndarray:
 432    """
 433    概要:
 434        波動関数の境界条件を記述する4x4行列を計算します。
 435    詳細説明:
 436        Kronig-Penneyモデルにおける各領域 (井戸と障壁) での波動関数の係数間の関係を記述する行列を構築します。
 437        この行列は、井戸と障壁の界面における波動関数とその導関数の連続性条件に基づいています。
 438        電子のエネルギー E、波数 k、井戸の幅 w、障壁の幅 b、およびポテンシャルの高さ V0 が入力として必要です。
 439    引数:
 440        :param k: 電子の波数 (pi/a 単位)。
 441        :type k: float
 442        :param E: 電子のエネルギー (eV)。
 443        :type E: float
 444        :param w: 井戸の幅 (A)。
 445        :type w: float
 446        :param b: ポテンシャル障壁の幅 (A)。
 447        :type b: float
 448        :param V0: ポテンシャルの高さ (eV)。
 449        :type V0: float
 450    戻り値:
 451        :returns: 境界条件を記述する4x4の複素数行列。
 452        :rtype: numpy.ndarray
 453    """
 454    alpha = sqrt(2.0 * ELECTRON_MASS * E * ELEMENTARY_CHARGE) / HBAR
 455    beta = sqrt(2.0 * ELECTRON_MASS * (V0 - E) * ELEMENTARY_CHARGE) / HBAR
 456    ka = k * PI2
 457    lambda_ = exp(1.0j * ka)
 458    alphaw = alpha * w * 1.0e-10
 459    betab = beta * b * 1.0e-10
 460    alpha *= 1.0e-10
 461    beta *= 1.0e-10
 462
 463    matrix = np.empty((4, 4), dtype=complex)
 464    matrix[0, 0] = matrix[0, 1] = 1.0
 465    matrix[0, 2] = matrix[0, 3] = -1.0
 466    matrix[1, 0] = 1.0j * alpha
 467    matrix[1, 1] = -1.0j * alpha
 468    matrix[1, 2] = -beta
 469    matrix[1, 3] = beta
 470    matrix[2, 0] = exp(1.0j * alphaw)
 471    matrix[2, 1] = exp(-1.0j * alphaw)
 472    matrix[2, 2] = -lambda_ * exp(-betab)
 473    matrix[2, 3] = -lambda_ * exp(betab)
 474    matrix[3, 0] = 1.0j * alpha * exp(1.0j * alphaw)
 475    matrix[3, 1] = -1.0j * alpha * exp(-1.0j * alphaw)
 476    matrix[3, 2] = -lambda_ * beta * exp(-betab)
 477    matrix[3, 3] = lambda_ * beta * exp(betab)
 478    return matrix
 479
 480
 481def compute_wave_coefficients(k: float, E: float, w: float, b: float, V0: float) -> list[complex]:
 482    """
 483    概要:
 484        波動関数の係数を計算します。
 485    詳細説明:
 486        周期ポテンシャル中の電子の波動関数は、Blochの定理により周期部分と平面波部分に分けられます。
 487        この関数は、境界条件行列を解くことで、井戸と障壁領域における波動関数の未知の係数 A, B, C, D を計算します。
 488        特に、井戸領域の左側から右側への入射波の係数を 1.0 と仮定し、残りの3つの係数を連立一次方程式の解として求めます。
 489    引数:
 490        :param k: 電子の波数 (pi/a 単位)。
 491        :type k: float
 492        :param E: 電子のエネルギー (eV)。
 493        :type E: float
 494        :param w: 井戸の幅 (A)。
 495        :type w: float
 496        :param b: ポテンシャル障壁の幅 (A)。
 497        :type b: float
 498        :param V0: ポテンシャルの高さ (eV)。
 499        :type V0: float
 500    戻り値:
 501        :returns: 波動関数の4つの複素数係数 (A, B, C, D) のリスト。
 502        :rtype: list
 503    """
 504    matrix = compute_boundary_matrix(k, E, w, b, V0)
 505
 506    a0 = 1.0
 507    matrix3 = np.empty((3, 3), dtype=complex)
 508    vector3 = np.empty((3, 1), dtype=complex)
 509
 510    matrix3[0, 0] = matrix[1, 1]
 511    matrix3[0, 1] = matrix[1, 2]
 512    matrix3[0, 2] = matrix[1, 3]
 513    matrix3[1, 0] = matrix[2, 1]
 514    matrix3[1, 1] = matrix[2, 2]
 515    matrix3[1, 2] = matrix[2, 3]
 516    matrix3[2, 0] = matrix[3, 1]
 517    matrix3[2, 1] = matrix[3, 2]
 518    matrix3[2, 2] = matrix[3, 3]
 519
 520    vector3[0, 0] = -a0 * matrix[1, 0]
 521    vector3[1, 0] = -a0 * matrix[2, 0]
 522    vector3[2, 0] = -a0 * matrix[3, 0]
 523
 524    ai = LA.solve(matrix3, vector3)
 525    return [a0, ai[0, 0], ai[1, 0], ai[2, 0]]
 526
 527
 528def check_wave_coefficients(
 529    ci: list[complex],
 530    kw: float,
 531    E: float,
 532    w: float,
 533    b: float,
 534    V0: float,
 535    eps: float,
 536    is_print: int = 0,
 537) -> None:
 538    """
 539    概要:
 540        計算された波動関数係数が境界条件を満たしているか検証します。
 541    詳細説明:
 542        計算された波動関数の係数 ci を用いて、境界条件行列 Mij との積 Mij @ ci が
 543        ほぼゼロになっていることを確認します。これは、係数が正しく計算されているかの健全性チェックです。
 544        積の絶対値の最大値が許容誤差 eps を超える場合、RuntimeErrorを発生させます。
 545    引数:
 546        :param ci: 波動関数の複素数係数のリスト。
 547        :type ci: list
 548        :param kw: 電子の波数 (pi/a 単位)。
 549        :type kw: float
 550        :param E: 電子のエネルギー (eV)。
 551        :type E: float
 552        :param w: 井戸の幅 (A)。
 553        :type w: float
 554        :param b: ポテンシャル障壁の幅 (A)。
 555        :type b: float
 556        :param V0: ポテンシャルの高さ (eV)。
 557        :type V0: float
 558        :param eps: 許容誤差。
 559        :type eps: float
 560        :param is_print: 結果を標準出力に表示するかどうかを示すフラグ (0:表示しない, 1:表示する)。
 561        :type is_print: int
 562    戻り値:
 563        なし。
 564    例外:
 565        :raises RuntimeError: 係数が境界条件を満たさない場合に発生します。
 566    """
 567    matrix = compute_boundary_matrix(kw, E, w, b, V0)
 568    values = matrix @ np.asarray(ci, dtype=complex)
 569
 570    if is_print:
 571        for i, coef in enumerate(ci):
 572            print(f"  ci[{i}] = {coef.real:12.4g}+j{coef.imag:12.4g}")
 573        for i, value in enumerate(values):
 574            print(f"  abs(Mij@ci[{i}]) = {abs(value)} {eps}")
 575
 576    vmax = float(np.max(np.abs(values)))
 577    if vmax > eps:
 578        raise RuntimeError(f"Mij @ ci is not zero: abs(Mij@ci)={vmax} > eps={eps}")
 579
 580
 581def compute_refined_energy(
 582    E0: float,
 583    E1: float,
 584    k: float,
 585    w: float,
 586    b: float,
 587    V0: float,
 588    ctx: SimpleNamespace,
 589    is_print: int = 0,
 590) -> tuple[float | None, float | None, float | None]:
 591    """
 592    概要:
 593        ニュートン法を用いてエネルギー準位を精密化します。
 594    詳細説明:
 595        Kronig-Penneyモデルのデルタ関数 compute_delta がゼロに近づくエネルギー E を、
 596        初期推定値 E0 と E1 からニュートン法を適用して見つけます。
 597        許容誤差 ctx.eps または最大反復回数 ctx.nmaxiter に達するまで反復計算を行います。
 598        収束しない場合はNoneを返します。
 599    引数:
 600        :param E0: エネルギーの初期推定値1 (eV)。
 601        :type E0: float
 602        :param E1: エネルギーの初期推定値2 (eV)。
 603        :type E1: float
 604        :param k: 電子の波数 (pi/a 単位)。
 605        :type k: float
 606        :param w: 井戸の幅 (A)。
 607        :type w: float
 608        :param b: ポテンシャル障壁の幅 (A)。
 609        :type b: float
 610        :param V0: ポテンシャルの高さ (eV)。
 611        :type V0: float
 612        :param ctx: 計算コンテキストを格納したSimpleNamespaceオブジェクト。
 613        :type ctx: types.SimpleNamespace
 614        :param is_print: 結果を標準出力に表示するかどうかを示すフラグ (0:表示しない, 1:表示する)。
 615        :type is_print: int
 616    戻り値:
 617        :returns: 精密化されたエネルギー、エネルギー変化 dE、およびデルタ関数の値のタプル。
 618                  収束しない場合は (None, None, None) を返します。
 619        :rtype: tuple
 620    """
 621    delta0 = compute_delta(E0, k, w, b, V0)
 622    delta1 = compute_delta(E1, k, w, b, V0)
 623
 624    E2 = E1
 625    dE = 0.0
 626    delta2 = delta1
 627
 628    for _ in range(ctx.nmaxiter):
 629        diff = (delta1 - delta0) / (E1 - E0)
 630        if diff >= 0.0:
 631            diff += ctx.dump
 632        else:
 633            diff = -(abs(diff) + ctx.dump)
 634
 635        dE = -delta1 / diff
 636        E2 = E1 + dE
 637        delta2 = compute_delta(E2, k, w, b, V0)
 638
 639        if abs(dE) < ctx.eps:
 640            if is_print:
 641                print(f"  converged at E = {E2:12.6g} with dE = {dE:12.6g}  delta = {delta2:12.6g}")
 642            return E2, dE, delta2
 643
 644        E0 = E1
 645        E1 = E2
 646        delta0 = delta1
 647        delta1 = delta2
 648
 649    print(f"  Not converged for {ctx.nmaxiter} iterations.")
 650    print(f"    E = {E2:12.6g} with dE = {dE:12.6g}  delta = {delta2:12.6g}")
 651    return None, None, None
 652
 653
 654def compute_energy_levels(
 655    emin: float,
 656    emax: float,
 657    nEsearch: int,
 658    k: float,
 659    w: float,
 660    b: float,
 661    V0: float,
 662    ctx: SimpleNamespace,
 663) -> tuple[list[float], list[list[complex]]]:
 664    """
 665    概要:
 666        指定されたエネルギー範囲内で許容されるエネルギー準位を計算します。
 667    詳細説明:
 668        エネルギー範囲 emin から emax を nEsearch 個のステップで探索し、
 669        compute_delta 関数の符号が反転する点を特定します。
 670        これらの点を初期値として compute_refined_energy を呼び出し、
 671        各エネルギー準位を精密化します。対応する波動関数の係数も計算して返します。
 672    引数:
 673        :param emin: 探索開始エネルギー (eV)。
 674        :type emin: float
 675        :param emax: 探索終了エネルギー (eV)。
 676        :type emax: float
 677        :param nEsearch: エネルギー探索点の数。
 678        :type nEsearch: int
 679        :param k: 電子の波数 (pi/a 単位)。
 680        :type k: float
 681        :param w: 井戸の幅 (A)。
 682        :type w: float
 683        :param b: ポテンシャル障壁の幅 (A)。
 684        :type b: float
 685        :param V0: ポテンシャルの高さ (eV)。
 686        :type V0: float
 687        :param ctx: 計算コンテキストを格納したSimpleNamespaceオブジェクト。
 688        :type ctx: types.SimpleNamespace
 689    戻り値:
 690        :returns: 許容されるエネルギー準位のリストと、それぞれのエネルギー準位に対応する波動関数係数のリストのタプル。
 691        :rtype: tuple
 692    """
 693    estep = (emax - emin) / (nEsearch - 1)
 694    previous_delta = None
 695    iband = 0
 696    energy_list = []
 697    coefficient_list = []
 698
 699    for iE in range(nEsearch):
 700        E = emin + iE * estep
 701        if E == 0.0:
 702            continue
 703        if V0 <= E:
 704            break
 705
 706        delta = compute_delta(E, k, w, b, V0)
 707
 708        if previous_delta is None:
 709            previous_delta = delta
 710            continue
 711
 712        if previous_delta * delta < 0.0:
 713            previous_delta = delta
 714            refined_E, dE, delta0 = compute_refined_energy(E - estep, E, k, w, b, V0, ctx, is_print=0)
 715            print(f"  E[{iband}]={refined_E:12.6g} eV  dE={dE:12.6g} delta={delta0:12.6g}")
 716
 717            if refined_E is not None:
 718                energy_list.append(refined_E)
 719                coefficient_list.append(compute_wave_coefficients(k, refined_E, w, b, V0))
 720                iband += 1
 721
 722    return energy_list, coefficient_list
 723
 724
 725def compute_wavefunction(ci: list[complex], x: float, kw: float, E: float, w: float, b: float, V0: float) -> complex:
 726    """
 727    概要:
 728        指定された位置 x における波動関数の値を計算します。
 729    詳細説明:
 730        計算された波動関数係数 ci を使用して、Blochの定理に基づく波動関数を評価します。
 731        位置 x を周期 a で正規化し、それが井戸領域か障壁領域かに応じて適切な波動関数形式を適用します。
 732        結果は複素数値の波動関数となります。
 733    引数:
 734        :param ci: 波動関数の複素数係数のリスト (A, B, C, D)。
 735        :type ci: list
 736        :param x: 波動関数の値を評価するx座標 (A)。
 737        :type x: float
 738        :param kw: 電子の波数 (pi/a 単位)。
 739        :type kw: float
 740        :param E: 電子のエネルギー (eV)。
 741        :type E: float
 742        :param w: 井戸の幅 (A)。
 743        :type w: float
 744        :param b: ポテンシャル障壁の幅 (A)。
 745        :type b: float
 746        :param V0: ポテンシャルの高さ (eV)。
 747        :type V0: float
 748    戻り値:
 749        :returns: xにおける波動関数の複素数値。
 750        :rtype: complex
 751    例外:
 752        :raises ValueError: 内部でのx座標の正規化が不正な場合に発生します。
 753    """
 754    a = w + b
 755    xmin = -b
 756    xmax = w
 757    x0, n_period = round01(x, a)
 758
 759    if x0 < -xmin:
 760        x0 += a
 761    if x0 >= xmax:
 762        x0 -= a
 763    if not xmin <= x0 < xmax:
 764        raise ValueError(f"x0 out of range: x={x:8.4g} {n_period} x0={x0:8.4g} w={w:8.4g} b={b:8.4g}")
 765
 766    alpha = sqrt(2.0 * ELECTRON_MASS * E * ELEMENTARY_CHARGE) / HBAR * 1.0e-10
 767    beta = sqrt(2.0 * ELECTRON_MASS * (V0 - E) * ELEMENTARY_CHARGE) / HBAR * 1.0e-10
 768    phase0 = PI2 / a * kw * x0
 769    kph0 = exp(1.0j * phase0)
 770
 771    if xmin <= x0 < 0.0:
 772        phi = ci[2] * exp(beta * x0) + ci[3] * exp(-beta * x0)
 773        periodic_part = phi / kph0
 774    else:
 775        phi = ci[0] * exp(1.0j * alpha * x0) + ci[1] * exp(-1.0j * alpha * x0)
 776        periodic_part = phi / kph0
 777
 778    return exp(1.0j * PI2 / a * kw * x) * periodic_part + 0.0j
 779
 780
 781def normalize_coefficients(
 782    ci: list[complex],
 783    E: float,
 784    kw: float,
 785    xstep: float,
 786    ctx: SimpleNamespace,
 787) -> list[complex]:
 788    """
 789    概要:
 790        波動関数の係数を正規化します。
 791    詳細説明:
 792        波動関数の全空間での確率密度積分 (abs(psi(x))^2 の積分) が1になるように、
 793        波動関数の係数 ci を調整します。積分は周期 a の範囲で数値的に行われます。
 794        これにより、物理的な解釈が可能な波動関数が得られます。
 795    引数:
 796        :param ci: 正規化する前の波動関数係数のリスト。
 797        :type ci: list
 798        :param E: 電子のエネルギー (eV)。
 799        :type E: float
 800        :param kw: 電子の波数 (pi/a 単位)。
 801        :type kw: float
 802        :param xstep: 積分に使用するxステップサイズ。
 803        :type xstep: float
 804        :param ctx: 計算コンテキストを格納したSimpleNamespaceオブジェクト。
 805        :type ctx: types.SimpleNamespace
 806    戻り値:
 807        :returns: 正規化された波動関数係数のリスト。
 808        :rtype: list
 809    """
 810    nxintg = int(ctx.a / xstep + 1.0001)
 811    xintgstep = ctx.a / (nxintg - 1)
 812    charge = 0.0
 813
 814    for i in range(nxintg):
 815        x = i * xintgstep
 816        yval = compute_wavefunction(ci, x, kw, E, ctx.w, ctx.b, ctx.V0)
 817        charge += yval * yval.conjugate()
 818
 819    charge = charge.real * xintgstep
 820    coefficient = 1.0 / sqrt(charge)
 821
 822    print("integ(|psi(x)|^2) = ", charge)
 823    print("Normalization coefficient = ", coefficient)
 824
 825    return [c * coefficient for c in ci]
 826
 827
 828def compute_graph_data(ctx: SimpleNamespace) -> tuple[list[float], list[float]]:
 829    """
 830    概要:
 831        delta(E) グラフプロット用のデータを計算します。
 832    詳細説明:
 833        コンテキストオブジェクト (ctx) で指定されたエネルギー範囲 (emin から emax) と
 834        ステップ数 (nE) に基づいて、各エネルギー E におけるデルタ関数 (compute_delta) の値を計算します。
 835        このデータは、許容エネルギーバンドを視覚化するために使用されます。
 836    引数:
 837        :param ctx: 計算コンテキストを格納したSimpleNamespaceオブジェクト。
 838        :type ctx: types.SimpleNamespace
 839    戻り値:
 840        :returns: エネルギー値のリストと対応するデルタ関数の値のリストのタプル。
 841        :rtype: tuple
 842    """
 843    estep = (ctx.emax - ctx.emin) / (ctx.nE - 1)
 844    xE = []
 845    yD = []
 846
 847    for i in range(1, ctx.nE):
 848        E = ctx.emin + i * estep
 849        if ctx.V0 <= E:
 850            break
 851        xE.append(E)
 852        yD.append(compute_delta(E, ctx.kg, ctx.w, ctx.b, ctx.V0))
 853
 854    return xE, yD
 855
 856
 857def compute_band_data(ctx: SimpleNamespace) -> tuple[list[float], np.ndarray, int]:
 858    """
 859    概要:
 860        バンド構造プロット用のデータを計算します。
 861    詳細説明:
 862        コンテキストオブジェクト (ctx) で指定された波数範囲 (kmin から kmax) と
 863        ステップ数 (nk) に基づいて、各波数 k における許容エネルギー準位を計算します。
 864        これらのエネルギー準位は、compute_energy_levels 関数を用いて探索され、
 865        バンド構造としてプロットされます。
 866    引数:
 867        :param ctx: 計算コンテキストを格納したSimpleNamespaceオブジェクト。
 868        :type ctx: types.SimpleNamespace
 869    戻り値:
 870        :returns: 波数 k のリスト、各波数におけるエネルギー準位のnumpy.ndarray、および見つかった最大バンドレベル数のタプル。
 871        :rtype: tuple
 872    """
 873    kstep = (ctx.kmax - ctx.kmin) / (ctx.nk - 1)
 874    xk = [ctx.kmin + i * kstep for i in range(ctx.nk)]
 875    yE = np.zeros((ctx.nMaxLevel, ctx.nk))
 876    n_max_band = 0
 877
 878    for ik, k in enumerate(xk):
 879        print(f"at k={k:8.4g}")
 880        energy_list, _ = compute_energy_levels(0.0, ctx.V0, ctx.nEsearch, k, ctx.w, ctx.b, ctx.V0, ctx)
 881        n_energy = len(energy_list)
 882        n_max_band = max(n_max_band, n_energy)
 883
 884        for iband in range(min(n_energy, ctx.nMaxLevel)):
 885            yE[iband, ik] = energy_list[iband]
 886
 887    return xk, yE, n_max_band
 888
 889
 890def compute_wavefunction_data(ctx: SimpleNamespace) -> tuple[np.ndarray, np.ndarray, list[float], np.ndarray, np.ndarray]:
 891    """
 892    概要:
 893        波動関数プロット用のデータを計算します。
 894    詳細説明:
 895        コンテキストオブジェクト (ctx) で指定された波数 (kw) とエネルギー準位インデックス (iLevel) に基づいて、
 896        対応する波動関数とその確率密度を計算します。
 897        また、ポテンシャルプロファイルも計算し、波動関数と共にプロットできるように準備します。
 898        波動関数は計算後に正規化されます。
 899    引数:
 900        :param ctx: 計算コンテキストを格納したSimpleNamespaceオブジェクト。
 901        :type ctx: types.SimpleNamespace
 902    戻り値:
 903        :returns: x座標、ポテンシャル値、エネルギー準位のリスト、波動関数の複素数値、確率密度のタプル。
 904        :rtype: tuple
 905    例外:
 906        :raises IndexError: 指定されたiLevelが利用可能なエネルギー準位の範囲外である場合に発生します。
 907    """
 908    xwstep = (ctx.xwmax - ctx.xwmin) / (ctx.nxw - 1)
 909    xplot, ypot = compute_potential_profile(ctx.xwmin, xwstep, ctx.nxw, ctx)
 910
 911    energy_list, coefficient_list = compute_energy_levels(0.0, ctx.V0, ctx.nEsearch, ctx.kw, ctx.w, ctx.b, ctx.V0, ctx)
 912    if ctx.iLevel >= len(energy_list):
 913        raise IndexError(f"iLevel={ctx.iLevel} is out of range. Found {len(energy_list)} energy levels.")
 914
 915    E = energy_list[ctx.iLevel]
 916    ci = coefficient_list[ctx.iLevel]
 917
 918    print("")
 919    print("=== Calculate wave function ===")
 920    print("Energy levels:", energy_list, "eV")
 921    print(f"at k = {ctx.kw}")
 922    print(f"{ctx.iLevel}-th energy level")
 923    print(f"  E = {E:12.6g} eV")
 924    print_coefficients(ci)
 925
 926    sumci = abs(ci[0] + ci[1] - ci[2] - ci[3])
 927    alpha = sqrt(2.0 * ELECTRON_MASS * E * ELEMENTARY_CHARGE) / HBAR * 1.0e-10
 928    beta = sqrt(2.0 * ELECTRON_MASS * (ctx.V0 - E) * ELEMENTARY_CHARGE) / HBAR * 1.0e-10
 929    print(f"  sum(ci) = {sumci:12.4e}")
 930    print(f"  alpha = {alpha:12.6g} A^-1")
 931    print(f"  beta  = {beta:12.6g} A^-1")
 932
 933    print("")
 934    print("Normalization")
 935    ci = normalize_coefficients(ci, E, ctx.kw, xwstep, ctx)
 936    print_coefficients(ci)
 937
 938    ywf = np.array(
 939        [compute_wavefunction(ci, ctx.xwmin + i * xwstep, ctx.kw, E, ctx.w, ctx.b, ctx.V0) for i in range(ctx.nxw)],
 940        dtype=complex,
 941    )
 942    charge = np.array([(value * value.conjugate()).real for value in ywf], dtype=float)
 943    return xplot, ypot, energy_list, ywf, charge
 944
 945
 946def print_coefficients(ci: list[complex]) -> None:
 947    """
 948    概要:
 949        波動関数の係数を整形して標準出力に表示します。
 950    詳細説明:
 951        与えられた波動関数係数 ci (通常は A, B, C, D に対応) を、
 952        実部と虚部に分けて読みやすい形式でコンソールに出力します。
 953    引数:
 954        :param ci: 波動関数の複素数係数のリスト。
 955        :type ci: list
 956    戻り値:
 957        なし。
 958    """
 959    names = ["A", "B", "C", "D"]
 960    for name, value in zip(names, ci):
 961        print(f"  {name} = {value.real:12.4g}+j{value.imag:12.4g}")
 962
 963
 964# =============================================================================
 965# plot functions
 966# =============================================================================
 967def plot_graph(ctx: SimpleNamespace) -> None:
 968    """
 969    概要:
 970        Kronig-Penneyモデルのdelta(E)グラフをプロットします。
 971    詳細説明:
 972        与えられたコンテキストオブジェクト (ctx) のパラメータに基づいて、
 973        エネルギー E に対するデルタ関数 (compute_delta) の値を計算し、プロットします。
 974        デルタ関数の絶対値が1を超える領域は許容エネルギーバンドに対応しません。
 975        生成された図は、ctx.save と ctx.show の設定に応じて保存または表示されます。
 976    引数:
 977        :param ctx: 計算コンテキストを格納したSimpleNamespaceオブジェクト。
 978        :type ctx: types.SimpleNamespace
 979    戻り値:
 980        なし。
 981    """
 982    estep = (ctx.emax - ctx.emin) / (ctx.nE - 1)
 983
 984    print("")
 985    print("=== Input parameters ===")
 986    print("mode:", ctx.mode)
 987    print("a=", ctx.a, "A")
 988    print(f"  barrier: w={ctx.b} A  h={ctx.V0} eV")
 989    print(f"  well   : w={ctx.w} A  h={0.0} eV")
 990    print(f"Energy range: {ctx.emin} - {ctx.emax}, {estep} eV step  {ctx.nE} points")
 991    print(f"at k = {ctx.kg}")
 992    print("")
 993
 994    xE, yD = compute_graph_data(ctx)
 995
 996    fig = plt.figure(figsize=ctx.figsize)
 997    ax = fig.add_subplot(1, 1, 1)
 998    ax.plot(xE, yD)
 999    ax.set_xlim([ctx.emin, ctx.emax])
1000    ax.axhline(0.0, linestyle="dashed", linewidth=0.5)
1001    ax.set_xlabel("E (eV)", fontsize=ctx.fontsize)
1002    ax.set_ylabel("delta", fontsize=ctx.fontsize)
1003    ax.tick_params(labelsize=ctx.fontsize)
1004    plt.tight_layout()
1005
1006    save_show_close(fig, ctx)
1007
1008
1009def plot_band(ctx: SimpleNamespace) -> None:
1010    """
1011    概要:
1012        Kronig-Penneyモデルのバンド構造をプロットします。
1013    詳細説明:
1014        与えられたコンテキストオブジェクト (ctx) のパラメータに基づいて、
1015        波数 k に対する許容エネルギー準位 (バンド) を計算し、プロットします。
1016        これにより、電子が占めることのできるエネルギー領域と、バンドギャップを視覚化します。
1017        生成された図は、ctx.save と ctx.show の設定に応じて保存または表示されます。
1018    引数:
1019        :param ctx: 計算コンテキストを格納したSimpleNamespaceオブジェクト。
1020        :type ctx: types.SimpleNamespace
1021    戻り値:
1022        なし。
1023    """
1024    kstep = (ctx.kmax - ctx.kmin) / (ctx.nk - 1)
1025
1026    print("")
1027    print("=== Input parameters ===")
1028    print("mode:", ctx.mode)
1029    print("a=", ctx.a, "A")
1030    print(f"potential: w={ctx.bwidth} A  h={ctx.bpot} eV")
1031    print(f"k range: {ctx.kmin} - {ctx.kmax} at {kstep} step, {ctx.nk} points")
1032    print("")
1033
1034    xk, yE, n_max_band = compute_band_data(ctx)
1035
1036    fig = plt.figure(figsize=ctx.figsize)
1037    ax = fig.add_subplot(1, 1, 1)
1038    ax.set_xlim([-0.5, 0.5])
1039    ax.set_ylim(ctx.erange)
1040
1041    for i_level in range(n_max_band):
1042        ax.plot(
1043            xk,
1044            yE[i_level],
1045            linestyle="",
1046            marker="o",
1047            markersize=5.0,
1048            markerfacecolor="none",
1049            markeredgewidth=0.5,
1050        )
1051
1052    ax.set_xlabel("$k$ $(\\pi$$/a)$", fontsize=ctx.fontsize)
1053    ax.set_ylabel("E (eV)", fontsize=ctx.fontsize)
1054    ax.tick_params(labelsize=ctx.fontsize)
1055    plt.tight_layout()
1056
1057    save_show_close(fig, ctx)
1058
1059
1060def plot_wavefunction(ctx: SimpleNamespace) -> None:
1061    """
1062    概要:
1063        Kronig-Penneyモデルの波動関数とその確率密度、およびポテンシャルプロファイルをプロットします。
1064    詳細説明:
1065        与えられたコンテキストオブジェクト (ctx) のパラメータに基づいて、
1066        特定の波数 k とエネルギー準位 iLevel における波動関数 (実部、虚部) とその確率密度 (charge) を計算します。
1067        これらのデータは、計算されたポテンシャルプロファイルと共に一つの図にプロットされ、
1068        電子の局在化とポテンシャルの関係を視覚化します。
1069        生成された図は、ctx.save と ctx.show の設定に応じて保存または表示されます。
1070    引数:
1071        :param ctx: 計算コンテキストを格納したSimpleNamespaceオブジェクト。
1072        :type ctx: types.SimpleNamespace
1073    戻り値:
1074        なし。
1075    """
1076    xwstep = (ctx.xwmax - ctx.xwmin) / (ctx.nxw - 1)
1077
1078    print("")
1079    print("=== Input parameters ===")
1080    print("mode:", ctx.mode)
1081    print("a=", ctx.a, "A")
1082    print(f"Wave function to be plotted: k = {ctx.kw}  iLevel = {ctx.iLevel}")
1083    print(f"x range: {ctx.xwmin} - {ctx.xwmax} at {xwstep} step, {ctx.nxw} points")
1084    print(f"potential: w={ctx.bwidth} A  h={ctx.bpot} eV")
1085    print("")
1086    print(f"at k={ctx.kw:8.4g}")
1087
1088    xplot, ypot, _, ywf, charge = compute_wavefunction_data(ctx)
1089
1090    fig = plt.figure(figsize=ctx.wf_figsize)
1091    ax_wave = fig.add_subplot(1, 1, 1)
1092    ax_potential = ax_wave.twinx()
1093
1094    ax_potential.set_xlim([ctx.xwmin, ctx.xwmax])
1095    ax_potential.plot(xplot, ypot, linewidth=0.5, label="U(x)")
1096    ax_potential.axhline(0.0, linestyle="dashed", linewidth=0.5)
1097
1098    ax_wave.set_xlim([ctx.xwmin, ctx.xwmax])
1099    ax_wave.plot(xplot, ywf.real, linewidth=1.5, label="real")
1100    ax_wave.plot(xplot, ywf.imag, linewidth=1.5, label="imaginary")
1101    ax_wave.plot(xplot, charge, linewidth=0.5, label="charge")
1102    ax_wave.axhline(0.0, linestyle="dashed", linewidth=0.5)
1103
1104    ax_potential.set_xlabel("x (A)", fontsize=ctx.fontsize)
1105    ax_potential.set_ylabel("U(x)", fontsize=ctx.fontsize)
1106    ax_wave.set_xlabel("x (A)", fontsize=ctx.fontsize)
1107    ax_wave.set_ylabel("$\\Psi$($x$)", fontsize=ctx.fontsize)
1108
1109    handler1, label1 = ax_potential.get_legend_handles_labels()
1110    handler2, label2 = ax_wave.get_legend_handles_labels()
1111    ax_wave.legend(
1112        handler1 + handler2,
1113        label1 + label2,
1114        loc=2,
1115        borderaxespad=0.0,
1116        fontsize=ctx.legend_fontsize,
1117    )
1118
1119    ax_potential.tick_params(labelsize=ctx.fontsize)
1120    ax_wave.tick_params(labelsize=ctx.fontsize)
1121    plt.tight_layout()
1122
1123    save_show_close(fig, ctx)
1124
1125
1126# =============================================================================
1127# run / main
1128# =============================================================================
1129def run(args: argparse.Namespace) -> None:
1130    """
1131    概要:
1132        Kronig-Penneyモデルの計算とプロットのメイン処理を実行します。
1133    詳細説明:
1134        まずコマンドライン引数を検証し、その後コンテキストオブジェクトを作成します。
1135        args.mode の値に応じて、delta(E)グラフ、バンド構造、または波動関数のいずれかのプロット関数を呼び出します。
1136    引数:
1137        :param args: コマンドライン引数を格納したNamespaceオブジェクト。
1138        :type args: argparse.Namespace
1139    戻り値:
1140        なし。
1141    例外:
1142        :raises ValueError: args.mode が認識されない値である場合に発生します。
1143    """
1144    validate_args(args)
1145    ctx = create_context(args)
1146
1147    if ctx.mode == "graph":
1148        plot_graph(ctx)
1149    elif ctx.mode == "band":
1150        plot_band(ctx)
1151    elif ctx.mode == "wf":
1152        plot_wavefunction(ctx)
1153    else:
1154        raise ValueError(f"Invalid mode: {ctx.mode}")
1155
1156
1157def main() -> None:
1158    """
1159    概要:
1160        スクリプトのエントリポイントです。
1161    詳細説明:
1162        コマンドライン引数を解析し、run関数を呼び出してKronig-Penneyモデルの計算とプロットを実行します。
1163        例外が発生した場合は、エラーメッセージを表示してプログラムを終了します。
1164    引数:
1165        なし。
1166    戻り値:
1167        なし。
1168    """
1169    args = parse_args()
1170    run(args)
1171
1172
1173if __name__ == "__main__":
1174    try:
1175        main()
1176    except Exception as exc:
1177        print("")
1178        print(f"Error: {exc}")
1179        traceback.print_exc()
1180        input("\nPress ENTER to terminate>>\n")
1181        sys.exit(1)