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

transfer_matrix.py をダウンロード

transfer_matrix.py

transfer_matrix.py

"""
伝達行列法を用いて1次元の波動関数と透過確率を計算するモジュール。

このスクリプトは、多層構造における電子の1次元波動関数と透過確率を、
伝達行列法（Transfer Matrix Method）を用いて計算します。
ポテンシャルプロファイルは、Excelファイルから読み込むか、
または多重量子井戸（Multiple Quantum Well, MQW）モデルに基づいて生成することができます。
計算結果はグラフとして表示され、透過確率データはExcelファイルに保存されます。

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

import os
import sys
import numpy as np
from numpy import sqrt, exp, log, sin, cos, tan, cosh, sinh
import numpy.linalg as LA 
import csv
import pandas as pd
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";


#========================
# global configuration
#========================
mode = 'wf'   # wf|tr

#========================
# potential
#========================
pottype = ''
#pottype = 'potential_defect.xlsx'
#pottype = 'mqw'

#barrierwidth  =  0.5  # A
barrierwidth  =  1.0  # Si
#wellwidth     = 10.0  # A 
wellwidth     = 5.4064 - barrierwidth  # Si
barrierheight =  10.0  # eV
nbarriers     = 10

# for wave functin plot
Ez0 = 0.1 # eV

# for transmission probability
Emin = 0.01 # eV
#Emax =  1.0 # eV
Emax =  9.5 # Si
nE   =  1001

# z range
zmin = -20.0 # zL, in angstrom
#zmax = 200.0 # zR, in angstrom
zmax = 70.0 # Si
nz   =  201

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


#==============================================
# fundamental functions
#==============================================
def pfloat(str):
    """
    文字列を浮動小数点数に安全に変換します。

    `float()` 関数と同様に文字列を浮動小数点数に変換しますが、
    変換に失敗した場合はエラーを発生させずにNoneを返します。

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

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

    `int()` 関数と同様に文字列を整数に変換しますが、
    変換に失敗した場合はエラーを発生させずにNoneを返します。

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

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

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

    :param position: 取得する引数の位置 (0-indexed)。
    :type position: int
    :param defval: 引数が存在しない場合に返すデフォルト値。
    :type defval: any, optional
    :returns: 指定された位置の引数、またはデフォルト値。
    :rtype: str or any
    """
    try:
        return sys.argv[position]
    except:
        return defval

def getfloatarg(position, defval = None):
    """
    コマンドライン引数を浮動小数点数として安全に取得します。

    `sys.argv` から指定された位置の引数を取得し、`pfloat` 関数を使って
    浮動小数点数に変換します。引数が存在しない、または変換できない場合は
    デフォルト値を返します。

    :param position: 取得する引数の位置 (0-indexed)。
    :type position: int
    :param defval: 引数が存在しない、または変換できない場合に返すデフォルト値。
    :type defval: float or None, optional
    :returns: 変換された浮動小数点数、またはデフォルト値。
    :rtype: float or None
    """
    return pfloat(getarg(position, defval))

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

    `sys.argv` から指定された位置の引数を取得し、`pint` 関数を使って
    整数に変換します。引数が存在しない、または変換できない場合は
    デフォルト値を返します。

    :param position: 取得する引数の位置 (0-indexed)。
    :type position: int
    :param defval: 引数が存在しない、または変換できない場合に返すデフォルト値。
    :type defval: int or None, optional
    :returns: 変換された整数、またはデフォルト値。
    :rtype: int or None
    """
    return pint(getarg(position, defval))

def usage():
    """
    プログラムの正しい使用方法を標準出力に表示します。

    プログラムの実行モード (`wf` または `tr`) に応じたコマンドライン引数の例を示します。
    """
    global mode, Ez, nz

    print("")
    print("Usage: Variables in () are optional")
    print("  python {} (pottype wf nz Ez0)".format(sys.argv[0]))
    print("    ex: python {} (wf {} {} {})".format(sys.argv[0], pottype, nz, Ez0))
    print("  python {} (tr nz Ez0 Emin Emax nE)".format(sys.argv[0]))
    print("    ex: python {} (tr {} {} {} {} {} {})".format(sys.argv[0], pottype, nz, Ez0, Emin, Emax, nE))

def terminate(message = None):
    """
    指定されたメッセージを表示し、プログラムを終了します。

    エラーメッセージなどを表示した後、`usage()` 関数を呼び出して使用方法を表示し、
    システムを終了します。

    :param message: 終了時に表示するメッセージ。
    :type message: str, optional
    """
    print("")
    if message is not None:
        print("")
        print(message)
        print("")

    usage()
    print("")
    exit()


def IsInBarrier(z):
    """
    指定されたz座標がポテンシャル障壁内にあるかを判定します。

    グローバル変数 `wellwidth`, `barrierwidth`, `nbarriers` に基づき、
    多重量子井戸構造における障壁領域にzが属するかをチェックします。

    :param z: 判定するz座標 (Å)。
    :type z: float
    :returns: zが障壁内にある場合は1、それ以外の場合は0。
    :rtype: int
    """
    global wellwidth, barrierwidth, barrierheight, nbarriers

    w = wellwidth + barrierwidth
    n = int(z / w)
    if n < nbarriers and 0.0 <= z - n * w < barrierwidth:
        return 1
    else:
        return 0

def U(z):   # eV
    """
    指定されたz座標におけるポテンシャルエネルギーを返します。

    `IsInBarrier` 関数を使用してz座標が障壁内にあるかを判断し、
    障壁内であれば `barrierheight` (eV) を、そうでなければ0.0 eVを返します。

    :param z: ポテンシャルエネルギーを評価するz座標 (Å)。
    :type z: float
    :returns: zにおけるポテンシャルエネルギー (eV)。
    :rtype: float
    """
    global wellwidth, barrierwidth, barrierheight, nbarriers

    if IsInBarrier(z):
        return barrierheight
    else:
        return 0.0

def build_U(pottype, xz = None):
    """
    ポテンシャルエネルギーと有効質量プロファイルを作成または読み込みます。

    `pottype` が 'mqw' の場合、多重量子井戸モデルに基づいてポテンシャル `yU` と
    有効質量 `ym` を生成します。それ以外の場合、指定されたExcelファイルから
    ポテンシャルと有効質量のデータを読み込みます。

    :param pottype: ポテンシャルの種類 ('mqw' または Excelファイルパス)。
    :type pottype: str
    :param xz: 'mqw'タイプの場合に使用されるz座標の配列。Noneの場合は生成される。
    :type xz: numpy.ndarray, optional
    :returns: (データ点の数, z座標配列, ポテンシャルエネルギー配列, 有効質量配列)。
    :rtype: tuple[int, numpy.ndarray, numpy.ndarray, numpy.ndarray]
    """
    if pottype == 'mqw':
        print()
        print(f"Build multiple quantum wells potential")
        nz = len(xz)
        yU = np.array([U(z) for z in xz])
        ym = np.array([meff(z) for z in xz])
        return nz, xz, yU, ym
    else:
        print()
        if not os.path.exists(pottype):
            print(f"Potential file [{pottype}] does not exists.")
            exit()

        print(f"Read potential from [{pottype}]")
        df = pd.read_excel(pottype)
        nz = len(df)
        xz = df.iloc[:, 0].to_numpy()
        yU = df.iloc[:, 1].to_numpy()
        ym = df.iloc[:, 2].to_numpy()
        return nz, xz, yU, ym
    
def meff(z):
    """
    指定されたz座標における有効質量を返します。

    z座標に基づいて、井戸領域または障壁領域の有効質量（電子質量meの倍数）を返します。
    この関数はハードコードされた領域に依存します。

    :param z: 有効質量を評価するz座標 (Å)。
    :type z: float
    :returns: zにおける有効質量（電子質量 me の倍数）。
    :rtype: float
    """
    mwell    = 1.0
    mbarrier = 1.0
    if 0.0 <= z < 5.0:
        return mbarrier
    elif 25.0 <= z < 30.0:
        return mbarrier
    else:
        return mwell


def cal_wf(xz, yU, ym, Ez):
    """
    伝達行列法を用いて1次元波動関数と透過確率を計算します。

    与えられたz座標、ポテンシャル、有効質量、およびエネルギーに対して、
    各セクションの波数 `kz`、振幅 `A` と `B`、および波動関数 `Psi` を計算します。
    最終的に透過確率 `T` も算出します。

    :param xz: z座標の配列 (Å)。
    :type xz: numpy.ndarray
    :param yU: 各z座標におけるポテンシャルエネルギーの配列 (eV)。
    :type yU: numpy.ndarray
    :param ym: 各z座標における有効質量の配列 (電子質量 me の倍数)。
    :type ym: numpy.ndarray
    :param Ez: 電子の入射エネルギー (eV)。
    :type Ez: float
    :returns: (波数配列, A振幅配列, B振幅配列, 波動関数配列, 透過確率)。
    :rtype: tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, float]
    """
    nz = len(xz)
    ykz  = np.empty(nz, dtype = complex)
    for i in range(nz):
        z = xz[i]
        kz = sqrt(2.0 * ym[i] * me / hbar / hbar * (Ez - yU[i])*e + 0.0j) * 1.0e-10
        ykz[i] = kz

    Ai  = np.empty(nz, dtype = complex)
    Bi  = np.empty(nz, dtype = complex)
    Psi = np.empty(nz, dtype = complex)
    Ai[0] = 1.0
    Bi[0] = 0.0
    for i in range(1, nz):
#        print(i, xz[i], Ez, ym[i], ykz[i])
        if ykz[i] == 0.0:
            # 波数kzが0の場合、数値的な安定性を保つため微小な値に置き換える
            # これによりゼロ除算を防ぎ、伝達行列の計算を続行できるようにする
            ykz[i] = 1.0e-10
        mk = ym[i] / ym[i-1] * ykz[i-1] / ykz[i]
        ap = 0.5 * (1.0 + mk)
        am = 0.5 * (1.0 - mk)
        P = exp(1.0j * (ykz[i-1] - ykz[i]) * xz[i])
        Q = exp(1.0j * (ykz[i-1] + ykz[i]) * xz[i])
        Ai[i] = ap * P * Ai[i-1] + am / Q * Bi[i-1]
        Bi[i] = am * Q * Ai[i-1] + ap / P * Bi[i-1]
        Psi[i] = Ai[i] * exp(1.0j * ykz[i] * xz[i]) + Bi[i] * exp(-1.0j * ykz[i] * xz[i])

    # 最後のステップでの透過確率の計算
    # 透過確率はAiとBiの比、および有効質量と波数に依存する
    T = ym[nz-1] / ym[0] * ykz[0] / ykz[nz-1] * Ai[0] / Ai[nz-1]
    T = pow(abs(T), 2)

    # 波動関数を規格化（最終的なPsiの振幅をkで割る）
    # ただし、この部分はPsi[i]の計算ループ外にあるため、Psi配列全体を規格化するには修正が必要
    # 現在の実装ではPsi[i]の最後の要素のみが更新されているが、Ai, Biも同様に規格化されている。
    k = sqrt(pow(abs(Ai[nz-1]), 2) + pow(abs(Bi[nz-1]), 2))
    # 以下は最後の要素のみを規格化するコードであることに注意
    Ai[i] /= k
    Bi[i] /= k
    Psi[i] /= k

    return ykz, Ai, Bi, Psi, T

def analytical_check(Ez, a):
    """
    無限に深い井戸型ポテンシャルの解析解との比較情報を提供します。

    指定されたエネルギー `Ez` と井戸幅 `a` を用いて、自由粒子の波数と、
    無限に深い井戸型ポテンシャルの基底状態エネルギーを計算し、表示します。
    これは、計算結果の妥当性を確認するための参考情報です。

    :param Ez: 比較に使用するエネルギー (eV)。
    :type Ez: float
    :param a: 無限井戸の幅 (Å)。
    :type a: float
    """
    print("")
    print("=== Analyatical results to check ===")
    kz = sqrt(Ez * e / (hbar * hbar / 2.0 / me)) * 1.0e-10 # A^-1
    print("Ez = {} eV".format(Ez))
    print("  kz = {} A^-1".format(kz*1.0e-10))

    print("=== Infinit potential well ===")
    print("Well width: {} A".format(a))
    k = 2.0 * pi / (2.0 * a * 1.0e-10) 
    E = hbar * hbar / 2.0 / me * k * k / e # eV
    print("  E0={:12.6g} eV".format(E))

def tr():
    """
    エネルギーに対する透過確率を計算し、結果をプロットします。

    コマンドライン引数またはグローバル設定に基づき、指定されたエネルギー範囲で
    透過確率を計算します。計算されたポテンシャルプロファイルと透過確率は
    Excelファイルに保存され、複数のグラフとして表示されます。
    """
    global mode
    global zmin, zmax, nz
    global Emin, Emax, nE
    global Ez0

    zstep = (zmax - zmin) / (nz - 1)
    Estep = (Emax - Emin) / (nE - 1)

    analytical_check(Ez0, 20.0)

    print("")
    print("=== Input parameterss ===")
    print("zmin(zL)=", zmin, "A")
    print("zmax(zL)=", zmax, "A")
    print("nz=", nz)
    print("zstep=", zstep)
    print("Ez0={} eV".format(Ez0))
    print("Erange: {} - {} eV, {} step, nE={}".format(Emin, Emax, Estep, nE))
    
    # z座標をzmaxからzminへと逆順に設定
    xz = np.array([zmax - i * zstep for i in range(nz)])
    nz, xz, yU, ym = build_U(pottype, xz)

    # エネルギー範囲を対数スケールで設定
    Elogstep = (log(Emax) - log(Emin)) / (nE - 1)
    xE = np.array([exp(log(Emin) + i * Elogstep) for i in range(nE)])
    
    potential_xlsx = 'potential.xlsx'
    print()
    print("Save potential to [{potential_xlsx}]")
    df = pd.DataFrame(np.array([xz, yU, ym]).T, columns = ["z (A)", "U (eV)", "m* (me)"])
    df.to_excel(potential_xlsx, index = False)

    print("")
    print("=== Wave function at Ez0 = {} eV ===".format(Ez0))
    # 特定のエネルギーEz0での波動関数と透過確率を計算
    ykz0, Ai0, Bi0, Psi0, T0 = cal_wf(xz, yU, ym, Ez0)
    print("Transmission probability at {} eV: {}".format(Ez0, T0))

    transmission_xlsx = 'transmission.xlsx'
    print("")
    print("=== Transmission probability vs Energy")
    yT = []
    print("{:10}\t{:10}".format("E (eV)", "T"))
    nprint_data = 100
    nskip = int(len(xE) / nprint_data)
    # 各エネルギーに対する透過確率を計算
    for i, E in enumerate(xE):
        ykz, Ai, Bi, Psi, T = cal_wf(xz, yU, ym, E)
        yT.append(T)
        if i % nskip == 0:
            print("{:10.6g}\t{:14.6g}".format(E, T))
    df = pd.DataFrame(np.array([xE, yT]).T, columns = ["E (eV)", "T"])
    df.to_excel(transmission_xlsx, index = False)


    print("")
    print("plot")
    fig = plt.figure(figsize = figsize)
    ax1 = fig.add_subplot(2, 2, 1)
# ax2をax1に関連させる
    ax2 = ax1.twinx()
    ax3 = fig.add_subplot(2, 2, 2)
    ax4 = fig.add_subplot(2, 2, 3)
    ax5 = fig.add_subplot(2, 2, 4)

    # プロット: ポテンシャルと有効質量
    ax1.plot(xz, yU, label = 'U(z)',         linewidth = 0.5, color = 'red')
    ax2.plot(xz, ym, label = 'm$_{eff}$(z)', linewidth = 0.5, color ='blue')
    
    # プロット: 波数kzの虚部と実部
    ykzr = [ykz0[i].real for i in range(nz)]
    ykzi = [ykz0[i].imag for i in range(nz)]
    ax3.plot(xz, ykzr, label = 'kz(real) (A$^{-1}$)', linewidth = 0.5, marker = 'o', markersize = 2)
    ax3.plot(xz, ykzi, label = 'kz(imag) (A$^{-1}$)', linewidth = 0.5)

    # プロット: 透過確率
    ax4.plot(xE, yT, label = 'T',    linewidth = 0.5, color = 'red', marker = 'o', markersize = 0.5)
    
    # プロット: 波動関数
    yPsir = [Psi0[i].real for i in range(nz)]
    yPsii = [Psi0[i].imag for i in range(nz)]
    yPsia = [pow(abs(Psi0[i]), 2) for i in range(nz)]
    ax5.plot(xz, yPsir, label = '$\Psi$(real)', linewidth = 0.5)
    ax5.plot(xz, yPsii, label = '$\Psi$(imag)', linewidth = 0.5)
    ax5.plot(xz, yPsia, label = '$|\Psi|^2$', linewidth = 0.5)

    # ラベル設定
    ax1.set_xlabel("z (A)", fontsize = fontsize)
    ax1.set_ylabel("U(z)", fontsize = fontsize)
    ax2.set_ylabel("m*(z)", fontsize = fontsize)
    ax3.set_xlabel("z (A)", fontsize = fontsize)
    ax3.set_ylabel("kz (A$^{-1}$)", fontsize = fontsize)
    ax4.set_xlabel("E (eV)", fontsize = fontsize)
    ax4.set_ylabel("Transmission probability", fontsize = fontsize)
    ax5.set_xlabel("z (A)", fontsize = fontsize)
    ax5.set_ylabel("$\Psi$(z)", fontsize = fontsize)

    # 軸範囲設定
    ax1.set_xlim([zmin, zmax])
    ax2.set_xlim([zmin, zmax])
    ax1.set_ylim([min(yU), max(yU) * 1.1])
    ax2.set_ylim([min(ym), max(ym) * 1.3])
    ax3.set_xlim([zmin, zmax])
#    ax4.set_xlim([min(0.0, Emin), Emax])
    ax4.set_ylim([0.0, 1.1])
    ax5.set_xlim([zmin, zmax])

# 凡例をまとめて出力する
    handler1, label1 = ax1.get_legend_handles_labels()
    handler2, label2 = ax2.get_legend_handles_labels()
    ax1.legend(handler1 + handler2, label1 + label2, loc = 2, borderaxespad = 0.0, fontsize = legend_fontsize)
    ax3.legend(fontsize = legend_fontsize)
    ax4.legend(fontsize = legend_fontsize)
    ax5.legend(fontsize = legend_fontsize)

    # 目盛り設定
    ax1.tick_params(labelsize = fontsize)
    ax2.tick_params(labelsize = fontsize)
    ax3.tick_params(labelsize = fontsize)
    ax4.tick_params(labelsize = fontsize)
    ax5.tick_params(labelsize = fontsize)
    plt.tight_layout()

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

    terminate()

def wf():
    """
    特定のエネルギーにおける波動関数を計算し、結果をプロットします。

    コマンドライン引数またはグローバル設定に基づき、指定された単一のエネルギー `Ez0` に対して
    波動関数 `Psi` を計算します。計算されたポテンシャルプロファイルと波動関数は
    複数のグラフとして表示されます。
    """
    global mode
    global zmin, zmax, nz
    global Ez0

    analytical_check(Ez0, 20.0)

    zstep = (zmax - zmin) / (nz - 1)

    print("")
    print("=== Input parameterss ===")
    print("zmin(zL)=", zmin, "A")
    print("zmax(zL)=", zmax, "A")
    print("nz=", nz)
    print("zstep=", zstep)
    print("Ez0=", Ez0, "eV")

    # z座標をzmaxからzminへと逆順に設定
    xz = np.array([zmax - i * zstep for i in range(nz)])
    nz, xz, yU, ym = build_U(pottype, xz)

    print("")
    print("=== Calculate wave function")
    ykz, Ai, Bi, Psi, T = cal_wf(xz, yU, ym, Ez0)

    print("Transmission probability from z = {} to {} A: {:14.6g}".format(xz[nz-1], xz[0], T))

    print("")
    print("plot")
    fig = plt.figure(figsize = figsize)
    ax1 = fig.add_subplot(2, 2, 1)
# ax2をax1に関連させる
    ax2 = ax1.twinx()
    ax3 = fig.add_subplot(2, 2, 2)
    ax4 = fig.add_subplot(2, 2, 3)
    ax5 = fig.add_subplot(2, 2, 4)

    # プロット: ポテンシャルと有効質量
    ax1.plot(xz, yU, label = 'U(z)',    linewidth = 0.5, color = 'red')
    ax2.plot(xz, ym, label = 'm$_{eff}$(z)', linewidth = 0.5, color ='blue')
    
    # プロット: 波数kzの虚部と実部
    ykzr = [ykz[i].real for i in range(nz)]
    ykzi = [ykz[i].imag for i in range(nz)]
    ax3.plot(xz, ykzr, label = 'kz(real) (A$^{-1}$)', linewidth = 0.5, marker = 'o', markersize = 2)
    ax3.plot(xz, ykzi, label = 'kz(imag) (A$^{-1}$)', linewidth = 0.5)
    
    # プロット: Ai, Biの絶対値
    yAir = [Ai[i].real for i in range(nz)]
    yAii = [Ai[i].imag for i in range(nz)]
    yAia = [abs(Ai[i]) for i in range(nz)]
    yBir = [Bi[i].real for i in range(nz)]
    yBii = [Bi[i].imag for i in range(nz)]
    yBia = [abs(Bi[i]) for i in range(nz)]
#    ax4.plot(xz, yAir, label = 'Ai(real)', linewidth = 0.5)
#    ax4.plot(xz, yAii, label = 'Ai(imag)', linewidth = 0.5)
    ax4.plot(xz, yAia, label = 'Ai(abs)', linewidth = 0.5)
#    ax4.plot(xz, yBii, label = 'Bi(imag)', linewidth = 0.5)
#    ax4.plot(xz, yBir, label = 'Bi(real)', linewidth = 0.5)
    ax4.plot(xz, yBia, label = 'Bi(abs)', linewidth = 0.5) # yBiiではなくyBiaをプロットする
    
    # プロット: 波動関数
    yPsir = [Psi[i].real for i in range(nz)]
    yPsii = [Psi[i].imag for i in range(nz)]
    yPsia = [pow(abs(Psi[i]), 2) for i in range(nz)]
    ax5.plot(xz, yPsir, label = '$\Psi$(real)', linewidth = 0.5)
    ax5.plot(xz, yPsii, label = '$\Psi$(imag)', linewidth = 0.5)
    ax5.plot(xz, yPsia, label = '$|\Psi|^2$', linewidth = 0.5)

    # ラベル設定
    ax1.set_xlabel("z (A)", fontsize = fontsize)
    ax1.set_ylabel("U(z)", fontsize = fontsize)
    ax2.set_ylabel("m*(z)", fontsize = fontsize)
    ax3.set_xlabel("z (A)", fontsize = fontsize)
    ax3.set_ylabel("kz (A$^{-1}$)", fontsize = fontsize)
    ax4.set_xlabel("z (A)", fontsize = fontsize)
    ax4.set_ylabel("Ai, Bi", fontsize = fontsize)
    ax5.set_xlabel("z (A)", fontsize = fontsize)
    ax5.set_ylabel("$\Psi$(z)", fontsize = fontsize)

    # 軸範囲設定
    ax1.set_xlim([zmin, zmax])
    ax2.set_xlim([zmin, zmax])
    ax3.set_xlim([zmin, zmax])
    ax4.set_xlim([zmin, zmax])
    ax5.set_xlim([zmin, zmax])
    ax1.set_ylim([min(yU), max(yU) * 1.1])
    ax2.set_ylim([min(ym), max(ym) * 1.3])
    # kzの軸範囲は計算された実部と虚部の最小値と最大値に基づいて調整
    ax3.set_ylim([min(min(ykzr), min(ykzi) - 0.1), max(max(ykzr) * 1.1, max(ykzi) * 1.3)]) # 実部と虚部の両方を考慮

# 凡例をまとめて出力する
    handler1, label1 = ax1.get_legend_handles_labels()
    handler2, label2 = ax2.get_legend_handles_labels()
    ax1.legend(handler1 + handler2, label1 + label2, loc = 2, borderaxespad = 0.0, fontsize = legend_fontsize)
    ax3.legend(fontsize = legend_fontsize)
    ax4.legend(fontsize = legend_fontsize)
    ax5.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()


if __name__ == "__main__":
#==============================================
# update default values by startup arguments
#==============================================
    argv = sys.argv

    mode    = getarg     (1, mode)
    pottype = getarg     (2, pottype)
    nz      = getintarg  (3, nz)
    Ez0     = getfloatarg(4, Ez0)

    if pottype == "":
        exit()

    if mode == 'wf':
        pass
    elif mode == 'tr':
        Emin  = getfloatarg(5, Emin)
        Emax  = getfloatarg(6, Emax)
        nE    = getintarg(7, nE)
    else:
        terminate("Error: Invalide mode [{}]".format(mode))

    if mode == 'wf':
        wf()
    elif mode == 'tr':
        tr()
    else:
        terminate("Error: Invalid mode [{}]".format(mode))