"""
概要:
H-like 1s軌道エネルギーレベル計算モジュール。
詳細説明:
このモジュールは、H-like 1s軌道のエネルギーレベルを計算するために、1s動径関数とスレーターのX-αポテンシャルを使用します。
ハートリー・フォック・スレーター（HFS）法の基本的な概念に基づき、電子間の相互作用を交換ポテンシャルで近似して扱います。
動径関数の指数係数kaや電子数Neなどを変分法で最適化し、全エネルギーを最小化する計算も可能です。
"""
import os
import sys
from math import exp, sqrt
import numpy as np
from math import log, exp
from numpy import arange
from scipy import integrate         # 数値積分関数 integrateを読み込む
from scipy.interpolate import interp1d
from scipy import optimize
from scipy.optimize import minimize
from matplotlib import pyplot as plt


# constants
pi   = 3.14159265358979323846
h    = 6.6260755e-34    # Js";
hbar = 1.05459e-34      # "Js";
c    = 2.99792458e8     # m/s";
e    = 1.60218e-19      # C";
kB   = 1.380658e-23     # JK<sup>-1</sup>";
me   = 9.1093897e-31    # kg";
e0          = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
e2_4pie0    = 2.30711e-28      # Nm<sup>2</sup>";
a0          = 5.29177e-11 * 1.0e10      # A
HartreeToeV = 27.2116          # eV";
RyToeV      = HartreeToeV / 2.0
pi2 = pi + pi
pi4 = pi2 + pi2


# mode: 'd': debug mode, plot fundamental graphs
# 'g': plot graph
# 'k': sweep ka, 'n': sweep Ne
# 'v': add variational calculation, 'e': energy-based output
mode = 'k'
ELabel = 'E 1s'

# Nuclear and orbital parameters
Z  = 1.0
n  = 1
l  = 0
m  = 0
ka = 1.0   # coefficient of the exponent in R1s
Ne = 1.0  

# Number of electrons, Slater's alpha
alpha = 2.0 / 3.0

# Radius range and integration (quad()) parameters
Rmin  =  0.0
Rmax  = 20.0
nR    = 2001
Rstep = None
Rmaxdata = None
nmaxdiv = 40
epsR = 1.0e-4
eps = 1.0e-8

# Ne mesh to calculate 1st and 2nd derivative of Etot
hparab = 0.01

# iteration parameters
#Nearray = [1.0, 0.8, 0.6]
#Nearray = [1.0, 0.8, 0.6, 0.5, 0.4, 0.3, 0.1, 0.01, 0.001]
Nearray = [1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.07, 0.05, 0.04, 0.03, 0.02, 0.01, 1.0e-3, 1.0e-4, 1.0e-5, 1.0e-6]
kaarray = [0.5, 0.6, 0.65, 0.7, 0.8, 0.9, 0.95, 1.0, 1.05, 1.1, 1.2, 1.4]

#=============================================
# scipy.optimize.minimizeで使うアルゴリズム
#=============================================
#nelder-mead Downhill simplex
#powell Modified Powell
#cg conjugate gradient (Polak-Ribiere method)
#bfgs BFGS法
#newton-cg Newton-CG
#trust-ncg 信頼領域 Newton-CG 法
#dogleg 信頼領域 dog-leg 法
#L-BFGS-B’ (see here)
#TNC’ (see here)
#COBYLA’ (see here)
#SLSQP’ (see here)
#trust-constr’(see here)
#dogleg’ (see here)
#trust-exact’ (see here)
#trust-krylov’ (see here)
method = "nelder-mead"
#method = 'cg'
#method = 'powell'
#method = 'bfgs'

maxiter = 1000
tol    = 1.0e-3
h_diff = 1.0e-3

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: 取得する引数のsys.argv内でのインデックス。
    :type position: int
    :param defval: 引数が存在しない場合に返すデフォルト値。デフォルトはNone。
    :type defval: any, optional
    :returns: 指定された位置の引数、または引数が存在しない場合はdefval。
    :rtype: str or any
    """
    try:
        return sys.argv[position]
    except:
        return defval

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

    詳細説明:
    指定された位置のコマンドライン引数を取得し、pfloat()関数を使用して浮動小数点数に変換します。
    引数が存在しない場合や変換に失敗した場合、エラーを発生させずにデフォルト値を返します。

    :param position: 取得する引数のsys.argv内でのインデックス。
    :type position: int
    :param defval: 引数が存在しない場合や変換に失敗した場合に返すデフォルト値。デフォルトはNone。
    :type defval: float or None, optional
    :returns: 変換された浮動小数点数、または失敗した場合はdefval。
    :rtype: float or None
    """
    return pfloat(getarg(position, defval))

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

    詳細説明:
    指定された位置のコマンドライン引数を取得し、pint()関数を使用して整数に変換します。
    引数が存在しない場合や変換に失敗した場合、エラーを発生させずにデフォルト値を返します。

    :param position: 取得する引数のsys.argv内でのインデックス。
    :type position: int
    :param defval: 引数が存在しない場合や変換に失敗した場合に返すデフォルト値。デフォルトはNone。
    :type defval: int or None, optional
    :returns: 変換された整数、または失敗した場合はdefval。
    :rtype: int or None
    """
    return pint(getarg(position, defval))

def usage(ka = ka, Z = Z, n = n, l = l, m = m):
    """
    概要: プログラムの使用方法を表示します。

    詳細説明:
    コマンドライン引数の構文と利用可能なモードオプションをユーザーに示します。
    デフォルト値のka, Z, n, l, mは現在のグローバル変数から取得されますが、
    直接引数として渡すことも可能です。

    :param ka: 1s動径関数の指数部の係数。デフォルトはグローバル変数ka。
    :type ka: float, optional
    :param Z: 原子番号。デフォルトはグローバル変数Z。
    :type Z: float, optional
    :param n: 主量子数。デフォルトはグローバル変数n。
    :type n: int, optional
    :param l: 方位量子数。デフォルトはグローバル変数l。
    :type l: int, optional
    :param m: 磁気量子数。デフォルトはグローバル変数m。
    :type m: int, optional
    :returns: なし
    :rtype: None
    """
    print("")
    print("Usage: Variables in () are optional")
    print("  (i) python {} mode Z ka Ne".format(sys.argv[0]))
    print("      mode: combination of the following key characters")
    print("               d: debug mode, plot fundamental graphs")
    print("               v: add variational calculations")
    print("               e: output based on energy (default: based on E 1s eigen value)")
    print("               k: sweep ka")
    print("               n: sweep Ne")
    print("               g : plot graph")
    print("     ex: python {} nvg {} {} {}".format(sys.argv[0], Z, ka, Ne))
    print("     ex: python {} k {} {} {}".format(sys.argv[0], Z, ka, Ne))

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

    詳細説明:
    オプションで終了メッセージを表示し、usage()関数を呼び出して使用方法を提示した後、
    プログラムを終了します。

    :param message: 終了時に表示するメッセージ。デフォルトはNone。
    :type message: str, optional
    :returns: なし (プログラムが終了するため)
    :rtype: None
    """
    if message is not None:
        print("")
        print(message)

    usage()
    print("")
    exit()

def updatevars():
    """
    概要: コマンドライン引数に基づいてグローバル変数を更新します。

    詳細説明:
    sys.argvからプログラム起動時の引数を取得し、mode, ELabel, ka, Z, Ne
    などのグローバル変数を更新します。これにより、ユーザーはコマンドラインから計算パラメータを
    指定できます。ELabelはmodeに'e'が含まれるかどうかに応じて設定されます。

    :returns: なし
    :rtype: None
    """
    global mode, ELabel
    global ka, Z, Ne, n, l, m, Ne, Rmax, Rstep

    argv = sys.argv
    if len(argv) == 1:
        terminate()

    mode = getarg( 1, mode)
    Z  = getfloatarg( 2, Z)
    ka = getfloatarg( 3, ka)
    Ne = getfloatarg( 4, Ne)
    
    if 'e' in mode:
        ELabel = 'Etot'
    else:
        ELabel = 'E 1s'
    
# Total charge inside r
yRr = []  # Raidal distributuion function
yQr = []  # Integrated charge inside r: integ_rm(4pi * rm * rm * rho(rm))_[rm = 0, r]

def Rr(ka, Z, n, l, m, r):
    """
    概要: H-like 1s軌道の動径関数 R_nl(r) を計算します。

    詳細説明:
    この関数は、水素様原子の1s軌道に対する動径関数を計算します。
    指数係数 ka と原子番号 Z に依存します。
    正規化条件は integ(4pi*r*r*Rr(r)*Rr(r))[r=0, inf] = 1 です。
    内部でグローバル変数 R1s0 を使用しますが、この変数はこのスクリプト内で定義されていません。

    :param ka: 動径関数の指数部分の係数。
    :type ka: float
    :param Z: 原子番号。
    :type Z: float
    :param n: 主量子数 (1s軌道なので常に1)。
    :type n: int
    :param l: 方位量子数 (1s軌道なので常に0)。
    :type l: int
    :param m: 磁気量子数 (1s軌道なので常に0)。
    :type m: int
    :param r: 核からの距離。
    :type r: float
    :returns: 指定された距離におけるH-like 1s軌道の動径関数の値。
    :rtype: float
    """
    global R1s0

    return R1s0 * pow(ka * Z, 1.5) * exp(-ka * Z * r)

def rho(ka, Z, n, l, m, r):
    """
    概要: H-like 1s軌道の電子密度関数 ρ(r) を計算します。

    詳細説明:
    動径関数 Rr(r) の二乗として電子密度を計算します。
    ρ(r) = Rr(r)^2 です。

    :param ka: 動径関数の指数部分の係数。
    :type ka: float
    :param Z: 原子番号。
    :type Z: float
    :param n: 主量子数 (1s軌道なので常に1)。
    :type n: int
    :param l: 方位量子数 (1s軌道なので常に0)。
    :type l: int
    :param m: 磁気量子数 (1s軌道なので常に0)。
    :type m: int
    :param r: 核からの距離。
    :type r: float
    :returns: 指定された距離における電子密度 ρ(r) の値。
    :rtype: float
    """
    psi = Rr(ka, Z, n, l, m, r)
    return psi * psi

def calQ(R):
    """
    概要: 核からの距離Rの内側に存在する総電荷 (Q(R)) を計算します。

    詳細説明:
    0からRまでの範囲で電子密度 ρ(r) を積分することにより、
    距離Rの内側の電子が形成する電荷を求めます。
    この関数は build_Qr で構築された補間関数 qfunc を使用します。
    qfuncが未構築の場合や、指定されたRが有効範囲外の場合はエラー処理または0を返します。

    :param R: 核からの距離。
    :type R: float
    :returns: 距離Rの内側に存在する総電荷。
    :rtype: float
    """
    global Rmin, Rmaxdata, r, Qr

    if R < Rmin:
        print("Error in Q(r): Given r={} exceeds the R range [{}, {}]".format(R, Rmin, Rmax))
        exit()
    if R > Rmaxdata:
        return 0.0

    return qfunc(R)

def build_Qr(ka, Z, n, l, m):
    """
    概要: 距離r内側の総電荷分布 Q(r) の補間関数を構築します。

    詳細説明:
    rの各点に対して Rr(r) と rho(r) を計算し、
    4pi * r*r * rho(r) を0から r まで積分することで Q(r) を求めます。
    求めた Q(r) のデータポイント yQr を使用して、
    scipy.interpolate.interp1d により三次スプライン補間関数 qfunc を構築し、
    calQ 関数で利用できるようにします。

    :param ka: 動径関数の指数部分の係数。
    :type ka: float
    :param Z: 原子番号。
    :type Z: float
    :param n: 主量子数 (1s軌道なので常に1)。
    :type n: int
    :param l: 方位量子数 (1s軌道なので常に0)。
    :type l: int
    :param m: 磁気量子数 (1s軌道なので常に0)。
    :type m: int
    :returns: なし。グローバル変数yRr, yQr, qfuncを更新します。
    :rtype: None
    """
    global yRr, yQr, qfunc

    yRr = []
    yQr = []
    for i in range(len(r)):
        yRr.append(Rr(ka, Z, n, l, m, r[i]))
        if r[i] == 0.0:
            Q   = 0.0
        else:
            Q, errQ = integrate.quad(lambda r: pi4 * r * r * rho(ka, Z, n, l, m, r), 0, r[i], limit = nmaxdiv, epsrel = eps)
        yQr.append(Ne * Q)
    qfunc = interp1d(r, yQr, kind = 'cubic')

def calTanal(Z = 1.0, ka = 1.0):
    """
    概要: H様原子の動径関数における運動エネルギーの解析解を計算します。

    詳細説明:
    この関数は、原子番号Zと指数係数kaを持つH様原子の1s軌道に対する運動エネルギーの
    解析的な近似値をハートリー単位で計算します。
    T = (ka*Z)^2 / 2 という関係に基づいています。

    :param Z: 原子番号。デフォルトは1.0。
    :type Z: float, optional
    :param ka: 動径関数の指数部分の係数。デフォルトは1.0。
    :type ka: float, optional
    :returns: 運動エネルギーの解析解 (ハートリー単位)。
    :rtype: float
    """
# T = e^2 / 8pi / e0 / a0
    Tanal = Z * Z * ka * ka * e / 8.0 / pi / e0 / (a0*1.0e-10)
    return Tanal

def calUanal(Z = 1.0, ka = 1.0):
    """
    概要: H様原子の動径関数における核電子相互作用ポテンシャルエネルギーの解析解を計算します。

    詳細説明:
    この関数は、原子番号Zと指数係数kaを持つH様原子の1s軌道に対する
    核と電子の間のポテンシャルエネルギーの解析的な近似値をハートリー単位で計算します。
    U = -(ka*Z)^2 という関係に基づいています。

    :param Z: 原子番号。デフォルトは1.0。
    :type Z: float, optional
    :param ka: 動径関数の指数部分の係数。デフォルトは1.0。
    :type ka: float, optional
    :returns: 核電子相互作用ポテンシャルエネルギーの解析解 (ハートリー単位)。
    :rtype: float
    """
# U = -e^2 / 4pi / e0 / a0
    Uanal = -Z * Z * ka * ka * e / 4.0 / pi / e0 / (a0*1.0e-10)
    return Uanal


def integrate_trapezoidal(func, E0, E1, h):
    """
    概要: 台形公式を用いて数値積分を行います。

    詳細説明:
    指定された関数 func を E0 から E1 までの範囲で、
    ステップサイズ h を用いて台形公式で数値積分します。
    この関数は現在使用されていませんが、将来的な高速化のために残されています。

    :param func: 積分する関数。引数を一つ取る関数を想定。
    :type func: callable
    :param E0: 積分範囲の下限。
    :type E0: float
    :param E1: 積分範囲の上限。
    :type E1: float
    :param h: 積分ステップサイズ。
    :type h: float
    :returns: [積分の結果, エラー値 (ここでは常に-1.0)]
    :rtype: list of float
    """
    n = int((E1 - E0) / h + 1.000001)
    h = (E1 - E0) / n
    y = [func(E0 + i * h) for i in range(n)]

    S = 0.5 * (y[0] + y[n-1]) + sum(y[1:n-1])
    return [h * S, -1.0]

def integrate3DR(func, rmin, rmax, limit = 15, epsrel = 1.0e-8):
    """
    概要: 3次元空間における動径方向の積分 (4πr^2 func(r) dr) を行います。

    詳細説明:
    scipy.integrate.quad を使用して、球座標における体積要素 4pi * r*r を考慮した
    動径方向の積分を実行します。

    :param func: 積分する動径関数。引数を一つ取る関数を想定。
    :type func: callable
    :param rmin: 積分範囲の下限。
    :type rmin: float
    :param rmax: 積分範囲の上限。
    :type rmax: float
    :param limit: quad関数の最大サブインターバル数。デフォルトは15。
    :type limit: int, optional
    :param epsrel: quad関数の相対許容誤差。デフォルトは1.0e-8。
    :type epsrel: float, optional
    :returns: [積分の結果, 積分の絶対誤差]
    :rtype: tuple of (float, float)
    """
    I, err = integrate.quad(lambda r: pi4 * r * r * func(r), rmin, rmax, limit, epsrel)
    return I, err

def calUZ(r, Z):
    """
    概要: 核電荷Zによって生成される静電ポテンシャル U_Z(r) を計算します。

    詳細説明:
    点電荷としての核が距離rに生成するクーロンポテンシャル -Z/r を計算します。
    rが非常に小さい場合 (1.0e-3より小さい場合) は、数値的な安定性のために
    -Z/1.0e-3を返します。

    :param r: 核からの距離。
    :type r: float
    :param Z: 原子番号 (核電荷)。
    :type Z: float
    :returns: 核電荷Zによって生成される静電ポテンシャル。
    :rtype: float
    """
    if r < 1.0e-3:
        return -Z / 1.0e-3
    else:
        return -Z / r

def calUrho(r, ka, Z, n, l, m):
    """
    概要: 電子密度 ρ(r) によって距離rに形成される静電ポテンシャル U_ee(r) を計算します。

    詳細説明:
    電子密度 ρ(r) が自身によって生成するポテンシャルを計算します。
    これは以下の2つの部分から構成されます。
    1. 距離rの内側にある電荷 Q(r) によるポテンシャル (Q(r)/r)。
    2. 距離rの外側にある電荷が距離rに生成するポテンシャル (4pi * integ(rm * ρ(rm) drm, rから∞))。

    :param r: 核からの距離。
    :type r: float
    :param ka: 動径関数の指数部分の係数。
    :type ka: float
    :param Z: 原子番号。
    :type Z: float
    :param n: 主量子数。
    :type n: int
    :param l: 方位量子数。
    :type l: int
    :param m: 磁気量子数。
    :type m: int
    :returns: 電子密度 ρ(r) によって距離rに形成される静電ポテンシャル。
    :rtype: float
    """
# Uerho(r) = integ(rho(m) / |r-rm|)_[rm=0, inf]                             [r, rm are vectors]
#          = Q(r) / r + integ[4pi * rm * rm * rho(rm) / rm]_[rm = r, inf])  [r is vector]
    if r == 0.0:
        Uee1 = 0.0
    else:
        Uee1 = calQ(r) / r

    Rmaxint = Rmax
#    Rmaxint = min(-log(eps) / Z / ka, Rmax)
    Uee2, errUee2 = integrate.quad(lambda rm: pi4 * rm * rho(ka, Z, n, l, m, rm), r, Rmaxint, limit = nmaxdiv, epsrel = eps)
    return Uee1 + Uee2
#    return Uee1 + Ne * Uee2

def calT(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps):
    """
    概要: 電子系の運動エネルギーを計算します。

    詳細説明:
    シュレーディンガー方程式の運動エネルギー演算子 (-1/2)∇^2 に基づいて、
    H-like 1s軌道の電子の運動エネルギーを数値積分によって計算します。
    積分の範囲は Rmin から Rmaxint (実質的には∞) です。
    内部でグローバル変数 R1s0 を使用しますが、この変数はこのスクリプト内で定義されていません。

    :param ka: 動径関数の指数部分の係数。
    :type ka: float
    :param Z: 原子番号。
    :type Z: float
    :param n: 主量子数。
    :type n: int
    :param l: 方位量子数。
    :type l: int
    :param m: 磁気量子数。
    :type m: int
    :param Ne: 電子数。
    :type Ne: float
    :param Rmin: 積分範囲の最小半径。
    :type Rmin: float
    :param Rmax: 積分範囲の最大半径（上限）。
    :type Rmax: float
    :param nmaxdiv: scipy.integrate.quad の最大サブインターバル数。
    :type nmaxdiv: int
    :param eps: scipy.integrate.quad の相対許容誤差。
    :type eps: float
    :returns: [運動エネルギー, 積分の絶対誤差]
    :rtype: tuple of (float, float)
    """
    Rmaxint = min(-log(eps) / Z / ka, Rmax)
    T, errT = integrate.quad(lambda r: r * (2.0 - Z * ka * r) * exp(-2.0 * Z * ka * r), 
                    Rmin, Rmaxint, limit = nmaxdiv, epsrel = eps)
    T *= 1.0 / 2.0 * pi4 * Z * ka * R1s0 * R1s0 * pow(ka * Z, 3.0)
    return T, errT

def calUeZ(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps):
    """
    概要: 電子と原子核間の引力ポテンシャルエネルギー U_eZ を計算します。

    詳細説明:
    電子密度 ρ(r) と核電荷 Z の間のクーロン引力によるポテンシャルエネルギーを
    数値積分によって計算します。積分は Rmin から Rmaxint (実質的には∞) です。

    :param ka: 動径関数の指数部分の係数。
    :type ka: float
    :param Z: 原子番号。
    :type Z: float
    :param n: 主量子数。
    :type n: int
    :param l: 方位量子数。
    :type l: int
    :param m: 磁気量子数。
    :type m: int
    :param Ne: 電子数。
    :type Ne: float
    :param Rmin: 積分範囲の最小半径。
    :type Rmin: float
    :param Rmax: 積分範囲の最大半径（上限）。
    :type Rmax: float
    :param nmaxdiv: scipy.integrate.quad の最大サブインターバル数。
    :type nmaxdiv: int
    :param eps: scipy.integrate.quad の相対許容誤差。
    :type eps: float
    :returns: [電子と原子核間の引力ポテンシャルエネルギー, 積分の絶対誤差]
    :rtype: tuple of (float, float)
    """
    Rmaxint = min(-log(eps) / Z / ka, Rmax)
    UeZ, errUeZ = integrate.quad(lambda r: pi4 * r * rho(ka, Z, n, l, m, r), 
                        Rmin, Rmaxint, limit = nmaxdiv, epsrel = eps)
    UeZ *= -Z
    return UeZ, errUeZ

def calUee(mode, ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps):
    """
    概要: 電子間のクーロン反発ポテンシャルエネルギー U_ee を計算します。

    詳細説明:
    電子密度 ρ(r) によって生成される静電ポテンシャル U_rho(r) と ρ(r) の積を
    全空間で積分することにより、電子間のクーロン反発エネルギーを計算します。
    mode に 'e' が含まれる場合、全エネルギー計算モードとして Ne*Ne 倍、
    それ以外の場合は Ne 倍します。

    :param mode: 現在の実行モードを示す文字列（'e'が含まれるかでスケールが変わる）。
    :type mode: str
    :param ka: 動径関数の指数部分の係数。
    :type ka: float
    :param Z: 原子番号。
    :type Z: float
    :param n: 主量子数。
    :type n: int
    :param l: 方位量子数。
    :type l: int
    :param m: 磁気量子数。
    :type m: int
    :param Ne: 電子数。
    :type Ne: float
    :param Rmin: 積分範囲の最小半径。
    :type Rmin: float
    :param Rmax: 積分範囲の最大半径（上限）。
    :type Rmax: float
    :param nmaxdiv: scipy.integrate.quad の最大サブインターバル数。
    :type nmaxdiv: int
    :param eps: scipy.integrate.quad の相対許容誤差。
    :type eps: float
    :returns: [電子間のクーロン反発ポテンシャルエネルギー, 積分の絶対誤差]
    :rtype: tuple of (float, float)
    """
    Rmaxint = min(-log(eps) / Z / ka, Rmax)
    Uee, errUee = integrate.quad(lambda r: pi4 * r * r * calUrho(r, ka, Z, n, l, m) * rho(ka, Z, n, l, m, r),
                         Rmin, Rmaxint, limit = nmaxdiv, epsrel = eps)
    if 'e' in mode:
        Uee *= Ne * Ne
    else:
        Uee *= Ne

    return Uee, errUee

def calLocalUXa(r, ka, Z, n, l, m, Ne):
    """
    概要: スレーターのX-α交換ポテンシャル (局所形式) を計算します。

    詳細説明:
    距離rにおけるスレーターのX-α交換ポテンシャルの局所形式 -3.0 * alpha * (3/(4pi) * ρ(r))^(1/3)
    を計算します。

    :param r: 核からの距離。
    :type r: float
    :param ka: 動径関数の指数部分の係数。
    :type ka: float
    :param Z: 原子番号。
    :type Z: float
    :param n: 主量子数。
    :type n: int
    :param l: 方位量子数。
    :type l: int
    :param m: 磁気量子数。
    :type m: int
    :param Ne: 電子数 (この関数では直接使用されないが引数として存在)。
    :type Ne: float
    :returns: 距離rにおけるX-α交換ポテンシャルの値。
    :rtype: float
    """
    return -3.0 * alpha *  pow(3.0/pi4 * rho(ka, Z, n, l, m, r), 1.0/3.0)

def calUXa(mode, ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps):
    """
    概要: X-α交換エネルギーを計算します。

    詳細説明:
    スレーターのX-α交換ポテンシャル calLocalUXa(r) と電子密度 rho(r) の積を
    全空間で積分することにより、X-α交換エネルギーを計算します。
    mode に 'e' が含まれる場合、全エネルギー計算モードとして Ne^(4/3) 倍、
    それ以外の場合は Ne^(1/3) 倍します。

    :param mode: 現在の実行モードを示す文字列（'e'が含まれるかでスケールが変わる）。
    :type mode: str
    :param ka: 動径関数の指数部分の係数。
    :type ka: float
    :param Z: 原子番号。
    :type Z: float
    :param n: 主量子数。
    :type n: int
    :param l: 方位量子数。
    :type l: int
    :param m: 磁気量子数。
    :type m: int
    :param Ne: 電子数。
    :type Ne: float
    :param Rmin: 積分範囲の最小半径。
    :type Rmin: float
    :param Rmax: 積分範囲の最大半径（上限）。
    :type Rmax: float
    :param nmaxdiv: scipy.integrate.quad の最大サブインターバル数。
    :type nmaxdiv: int
    :param eps: scipy.integrate.quad の相対許容誤差。
    :type eps: float
    :returns: [X-α交換エネルギー, 積分の絶対誤差]
    :rtype: tuple of (float, float)
    """
    Rmaxint = min(-log(eps) / Z / ka, Rmax)
    UXa, errUXa = integrate.quad(lambda r: pi4 * r * r * calLocalUXa(r, ka, Z, n, l, m, Ne) * rho(ka, Z, n, l, m, r), 
                        Rmin, Rmaxint, limit = nmaxdiv, epsrel = eps)
    if 'e' in mode:
        UXa *= pow(Ne, 4.0/3.0)
    else:
        UXa *= pow(Ne, 1.0/3.0)
    
    return UXa, errUXa

def calTotalEnergy(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps):
    """
    概要: H-like 1s軌道の全エネルギーを計算します。

    詳細説明:
    運動エネルギー (T)、核電子引力エネルギー (UeZ)、電子間反発エネルギー (Uee)、
    およびX-α交換エネルギー (UXa) の各項を計算し、それらを合計して全エネルギーを求めます。
    build_Qr を呼び出して Q(r) 分布を事前に構築します。

    :param ka: 動径関数の指数部分の係数。
    :type ka: float
    :param Z: 原子番号。
    :type Z: float
    :param n: 主量子数。
    :type n: int
    :param l: 方位量子数。
    :type l: int
    :param m: 磁気量子数。
    :type m: int
    :param Ne: 電子数。
    :type Ne: float
    :param Rmin: 積分範囲の最小半径。
    :type Rmin: float
    :param Rmax: 積分範囲の最大半径（上限）。
    :type Rmax: float
    :param nmaxdiv: scipy.integrate.quad の最大サブインターバル数。
    :type nmaxdiv: int
    :param eps: scipy.integrate.quad の相対許容誤差。
    :type eps: float
    :returns: [運動エネルギーT, 核電子引力エネルギーUeZ, 電子間反発エネルギーUee, X-α交換エネルギーUXa, 全エネルギーEtot]
    :rtype: tuple of (float, float, float, float, float)
    """
    build_Qr(ka, Z, n, l, m)

#def calT(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps):
    T, errT = calT(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
    UeZ, errUeZ = calUeZ(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
    Uee, errUee = calUee(mode, ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
    UXa, errUXa = calUXa(mode, ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
    Etot = T + UeZ + Uee + UXa

    return T, UeZ, Uee, UXa, Etot

def calTotalEnergyOnly(kav, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps):
    """
    概要: 変分計算のために、特定のka値におけるH-like 1s軌道の全エネルギーを計算します。

    詳細説明:
    この関数は calTotalEnergy と同様に全エネルギーを計算しますが、
    最適化ルーチン (minimize) から呼び出されることを想定し、
    引数として与えられた kav を使用し、最終的な全エネルギー値のみを返します。
    また、電子間反発エネルギー (Uee) とX-α交換エネルギー (UXa) の計算モードは
    常に全エネルギー計算 ('e'モード) として扱われます。

    :param kav: 変分パラメータとして使用される動径関数の指数部分の係数。
    :type kav: float
    :param Z: 原子番号。
    :type Z: float
    :param n: 主量子数。
    :type n: int
    :param l: 方位量子数。
    :type l: int
    :param m: 磁気量子数。
    :type m: int
    :param Ne: 電子数。
    :type Ne: float
    :param Rmin: 積分範囲の最小半径。
    :type Rmin: float
    :param Rmax: 積分範囲の最大半径（上限）。
    :type Rmax: float
    :param nmaxdiv: scipy.integrate.quad の最大サブインターバル数。
    :type nmaxdiv: int
    :param eps: scipy.integrate.quad の相対許容誤差。
    :type eps: float
    :returns: 計算された全エネルギー。
    :rtype: float
    """
    build_Qr(kav, Z, n, l, m)

    T, errT = calT(kav, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
    UeZ, errUeZ = calUeZ(kav, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
    Uee, errUee = calUee('e', kav, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
    UXa, errUXa = calUXa('e', kav, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
    Etot = T + UeZ + Uee + UXa

    return Etot

def calOptimizedTotalEnergy(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps):
    """
    概要: 変分法を用いて全エネルギーを最小化するka値を探索し、そのときの全エネルギーと各エネルギー成分を計算します。

    詳細説明:
    scipy.optimize.minimize 関数を使用して、動径関数の指数係数 ka を変分パラメータとして、
    calTotalEnergyOnly で計算される全エネルギーを最小化します。
    指定された最適化手法 (method) とパラメータ (tol, maxiter) に基づいて探索を行い、
    最小化された ka とそのときの各エネルギー成分、および全エネルギーを返します。

    :param ka: 最適化の初期値となる動径関数の指数部分の係数。
    :type ka: float
    :param Z: 原子番号。
    :type Z: float
    :param n: 主量子数。
    :type n: int
    :param l: 方位量子数。
    :type l: int
    :param m: 磁気量子数。
    :type m: int
    :param Ne: 電子数。
    :type Ne: float
    :param Rmin: 積分範囲の最小半径。
    :type Rmin: float
    :param Rmax: 積分範囲の最大半径（上限）。
    :type Rmax: float
    :param nmaxdiv: scipy.integrate.quad の最大サブインターバル数。
    :type nmaxdiv: int
    :param eps: scipy.integrate.quad の相対許容誤差。
    :type eps: float
    :returns: [運動エネルギーT, 核電子引力エネルギーUeZ, 電子間反発エネルギーUee, X-α交換エネルギーUXa, 最小化された全エネルギーEtot, 最適化されたka値]
    :rtype: tuple of (float, float, float, float, float, float)
    """
    xa = [ka]
    ret = minimize(lambda xa: calTotalEnergyOnly(xa[0], Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps),
                xa, method = method, jac = diff1, tol = tol, 
                callback = lambda xa: callback(xa),
                options = {'maxiter':maxiter, "disp":True})
#       print("ret=", ret)
    if method == 'nelder-mead':
        simplex = ret['final_simplex']
        ai = simplex[0][0]    
        Emin = ret['fun']
    elif method == 'cg':
        ai = ret['x']
        Emin = ret['fun']
    elif method == 'powell':
        ai = ret['x']
        Emin = ret['fun']
    elif method == 'bfgs':
        ai = ret['x']
        ka = xa[0]
        Emin = calTotalEnergyOnly(xa[0], Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)

    ka = ai[0]
    T, UeZ, Uee, UXa, Etot = calTotalEnergy(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)

    return T, UeZ, Uee, UXa, Etot, ka


def diff1(ai):
    """
    概要: 1変数関数の勾配 (1次導関数) を中心差分法で近似します。

    詳細説明:
    scipy.optimize.minimize の勾配情報が必要なメソッド (cgなど) のために、
    パラメータ ai (ここでは ka) の微小変化 h_diff を用いて、
    calTotalEnergyOnly 関数に対する1次導関数を数値的に計算します。

    :param ai: 最適化される変数の現在の値 (ここでは ka の配列)。
    :type ai: numpy.ndarray
    :returns: 各変数に対する1次導関数の配列。
    :rtype: numpy.ndarray
    """
    global Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps
    
    n = len(ai)
    f0 = calTotalEnergyOnly(ai[0], Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
    df = np.empty(n)
    for i in range(n):
        aii = ai
        aii[i] = ai[i] + h_diff
        f1 = calTotalEnergyOnly(aii[0], Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
        df[i] = (f1 - f0) / h_diff
    return df

iter = 0
def callback(xk):
    """
    概要: 最適化の各イテレーションで呼び出され、進行状況を表示します。

    詳細説明:
    scipy.optimize.minimize 関数の callback オプションとして使用されます。
    各イテレーションで現在のパラメータ xk (ここでは ka) と
    それに対応する全エネルギー Etot を表示します。

    :param xk: 現在のイテレーションでの最適化変数の値 (ここでは ka の配列)。
    :type xk: numpy.ndarray
    :returns: なし
    :rtype: None
    """
    global iter
    global Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps

    ka = xk[0]
    Etot = calTotalEnergyOnly(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)

# パラメータと残差二乗和を表示
    print("callback {}: ka={:14.4g} Emin={} eV".format(iter, ka, Etot * HartreeToeV))
    iter += 1


def sweep_Ne(Eanal):
    """
    概要: 電子数Neを掃引し、非最適化および最適化された全エネルギーの変化を評価・プロットします。

    詳細説明:
    Nearray に定義された様々な電子数Neに対して、固定された ka での全エネルギー、
    および変分法で ka を最適化した上での全エネルギーを計算します。
    また、Ne=0.5付近での放物線近似も行い、結果を表形式で出力し、グラフとしてプロットします。
    特に 'v' モードが有効な場合、ka の最適化も行います。

    :param Eanal: 比較のための解析解エネルギー値。
    :type Eanal: float
    :returns: なし
    :rtype: None
    """
    global Rmin, Rstep, nR, Rmaxdata, r
    global ka
    global qfunc

    print("")
    print("Not optimized")

    print("  Parabolic expantion around Ne = 0.5")
    Tp, UeZp, Ueep, UXap, Etotp = calTotalEnergy(ka, Z, n, l, m, 0.5 + hparab, Rmin, Rmax, nmaxdiv, eps)
    T0, UeZ0, Uee0, UXa0, Etot0 = calTotalEnergy(ka, Z, n, l, m, 0.5         , Rmin, Rmax, nmaxdiv, eps)
    Tm, UeZm, Ueem, UXam, Etotm = calTotalEnergy(ka, Z, n, l, m, 0.5 - hparab, Rmin, Rmax, nmaxdiv, eps)
    Etotp *= HartreeToeV
    Etot0 *= HartreeToeV
    Etotm *= HartreeToeV
    print("  Ne={}: E={} eV".format(0.5 + hparab, Etotp))
    print("  Ne={}: E={} eV".format(0.5         , Etot0))
    print("  Ne={}: E={} eV".format(0.5 - hparab, Etotm))
    a2 = (Etotp - 2.0 * Etot0 + Etotm) / hparab / hparab
    a1 = (Etotp - Etotm) / 2.0 / hparab
    print("  E = {} + {} * Ne + {} * Ne^2".format(Etot0, a1, a2))

    xNe = []
    yE1s = []
    yE1sparab = []
    print("")
    print("{:6s}\t{:6s}\t{:6s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}"
            .format('ka', 'Z', 'Ne', 'T in eV', 'U(e-Z)', 'U(e-e)', 'U(Xa)', 
                    ELabel, ELabel + ' (parabolic)'))
    for Ne in Nearray:
        T, UeZ, Uee, UXa, Etot = calTotalEnergy(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
        Eparab = Etot0 + a1 * (Ne - 0.5) + a2 * pow(Ne - 0.5, 2)
        print("{:6.4f}\t{:6.4f}\t{:6.4f}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}"
            .format(ka, Z, Ne, 
                    T * HartreeToeV, UeZ * HartreeToeV, Uee * HartreeToeV, UXa * HartreeToeV, 
                    Etot * HartreeToeV, Eparab))
        xNe.append(Ne)
        yE1s.append(Etot * HartreeToeV)
        yE1sparab.append(Eparab)

    if 'v' in mode:
        yE1sOpt = []
        yE1sOptparab = []
        ykaOpt  = []
        print("")
        print("Optimized for ka by variational principle")

        print("  Parabolic expantion around Ne = 0.5")
        Tp, UeZp, Ueep, UXap, Etotp, kap = calOptimizedTotalEnergy(ka, Z, n, l, m, 0.5 + hparab, Rmin, Rmax, nmaxdiv, eps)
        T0, UeZ0, Uee0, UXa0, Etot0, ka0 = calOptimizedTotalEnergy(ka, Z, n, l, m, 0.5         , Rmin, Rmax, nmaxdiv, eps)
        Tm, UeZm, Ueem, UXam, Etotm, kam = calOptimizedTotalEnergy(ka, Z, n, l, m, 0.5 - hparab, Rmin, Rmax, nmaxdiv, eps)
        Etotp *= HartreeToeV    
        Etot0 *= HartreeToeV
        Etotm *= HartreeToeV
        print("  Ne={}: E={} eV".format(0.5 + hparab, Etotp))
        print("  Ne={}: E={} eV".format(0.5         , Etot0))
        print("  Ne={}: E={} eV".format(0.5 - hparab, Etotm))
        a2 = (Etotp - 2.0 * Etot0 + Etotm) / hparab / hparab
        a1 = (Etotp - Etotm) / 2.0 / hparab
        print("  E = {} + {} * Ne + {} * Ne^2".format(Etot0, a1, a2))

        print("{:6s}\t{:6s}\t{:6s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}"
                .format('ka', 'Z', 'Ne', 'T in eV', 'U(e-Z)', 'U(e-e)', 'U(Xa)', ELabel, ELabel + ' (parabolic)'))
        for Ne in Nearray:
            T, UeZ, Uee, UXa, Etot, ka = calOptimizedTotalEnergy(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
            Eparab = Etot0 + a1 * (Ne - 0.5) + a2 * pow(Ne - 0.5, 2)
            print("{:6.4f}\t{:6.4f}\t{:6.4g}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}t{:10.6g}"
                .format(ka, Z, Ne, 
                       T * HartreeToeV, UeZ * HartreeToeV, Uee * HartreeToeV, UXa * HartreeToeV, 
                       Etot * HartreeToeV, Eparab))

            yE1sOpt.append(Etot * HartreeToeV)
            yE1sOptparab.append(Eparab)
            ykaOpt.append(ka)

#=============================
# Plot graphs
#=============================
    if 'g' not in mode:
        terminate()

    fig = plt.figure(figsize = (12, 8))

    ax1 = fig.add_subplot(1, 2, 1)
    ax2 = fig.add_subplot(1, 2, 2)
#    ax3 = fig.add_subplot(2, 3, 3)
#    ax4 = fig.add_subplot(2, 3, 4)
#    ax5 = fig.add_subplot(2, 3, 5)
    ax1.plot(xNe, yE1s,      label = ELabel + ' (non-optimized)', color = 'black', linewidth = 1.5, marker = 'o', markersize = 1.0)
    ax1.plot(xNe, yE1sparab, label = ELabel + ' (non-opt,parabolic)', color = 'black', linestyle = 'dashed', linewidth = 0.5)
    if 'v' in mode:
        ax1.plot(xNe, yE1sOpt,      label = ELabel + ' (optimized)',     color = 'red',   linewidth = 1.5, marker = 'o', markersize = 1.0)
        ax1.plot(xNe, yE1sOptparab, label = ELabel + ' (opt,parabolic)', color = 'red',   linewidth = 0.5)
    ax1.plot([min(xNe), max(xNe)], [Eanal, Eanal], color = 'red', linestyle = 'dashed')
    ax1.plot([0.5, 0.5], [min(yE1s), max(yE1s)], color = 'red', linestyle = 'dashed')
    ax1.set_xlabel("Ne")
    ax1.set_ylabel(ELabel + " (eV)")
    ax1.legend()
    if 'v' in mode:
        ax2.plot(xNe, ykaOpt, label = 'ka (optimized)', color = 'black', linewidth = 0.5, marker = 'o', markersize = 1.0)
    ax2.plot([min(xNe), max(xNe)], [1.0, 1.0], color = 'red', linestyle = 'dashed')
    if 'v' in mode:
        ax2.plot([0.5, 0.5], [min(ykaOpt), max(ykaOpt)], color = 'red', linestyle = 'dashed')
    ax2.set_xlabel("Ne")
    ax2.set_ylabel("ka (optimized)")
    ax2.legend()

    plt.tight_layout()
    plt.pause(0.1)

    print("Press ENTER to exit>>", end = '')
    input()


def sweep_ka(Eanal):
    """
    概要: 動径関数の指数係数kaを掃引し、全エネルギーの変化を評価・プロットします。

    詳細説明:
    kaarray に定義された様々な ka の値に対して、
    H-like 1s軌道の全エネルギーとその各成分を計算します。
    結果を表形式で出力し、'g' モードが有効な場合はグラフとしてプロットします。

    :param Eanal: 比較のための解析解エネルギー値。
    :type Eanal: float
    :returns: なし
    :rtype: None
    """
    global Rmin, Rstep, nR, Rmaxdata, r
    global qfunc
    global nmaxdiv, eps

    xka = []
    yE1s = []
    print("")
    print("{:6s}\t{:6s}\t{:6s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}"
            .format('ka', 'Z', 'Ne', 'T in eV', 'U(e-Z)', 'U(e-e)', 'U(Xa)', ELabel))
    for kav in kaarray:
        T, UeZ, Uee, UXa, Etot = calTotalEnergy(kav, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
        print("{:6.4f}\t{:6.4f}\t{:6.4f}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}"
            .format(kav, Z, Ne, 
                    T * HartreeToeV, UeZ * HartreeToeV, Uee * HartreeToeV, UXa * HartreeToeV, Etot * HartreeToeV))
        xka.append(kav)
        yE1s.append(Etot * HartreeToeV)

#=============================
# Plot graphs
#=============================
    if 'g' not in mode:
        terminate()

    fig = plt.figure(figsize = (12, 8))

    ax1 = fig.add_subplot(2, 3, 1)
    ax2 = fig.add_subplot(2, 3, 2)
    ax3 = fig.add_subplot(2, 3, 3)
    ax4 = fig.add_subplot(2, 3, 4)
    ax5 = fig.add_subplot(2, 3, 5)
    ax1.plot(xka, yE1s, label = ELabel, color = 'black', marker = 'o')
    ax1.plot([min(xka), max(xka)], [Eanal, Eanal], color = 'red', linestyle = 'dashed')
    ax1.set_xlabel("ka")
    ax1.set_ylabel(ELabel + " (eV)")
    ax1.legend()

    plt.tight_layout()
    plt.pause(0.1)

    print("Press ENTER to exit>>", end = '')
    input()

def debug(Eanal):
    """
    概要: デバッグモードで、H-like 1s軌道の計算結果と関連するグラフを表示します。

    詳細説明:
    現在のパラメータ (ka, Z, Ne など) における運動エネルギー、核電子引力エネルギー、
    電子間反発エネルギー、X-α交換エネルギー、および全エネルギーを計算し、表形式で出力します。
    'v' モードが有効な場合は、ka を最適化した結果も表示します。
    また、'g' モードが有効な場合は、動径関数 Rr(r)、電荷分布 Q(r)、および各種ポテンシャル
    (U(Z), U(Z-rho), U(rho), U(Xa)) のグラフをプロットします。

    :param Eanal: 比較のための解析解エネルギー値。
    :type Eanal: float
    :returns: なし
    :rtype: None
    """
    global ka, Z, n, l, m, Ne
    global Rmin, Rstep, nR, Rmaxdata, r

    ka0, Z0 = ka, Z
    
    print("")
    print("{:6s}\t{:6s}\t{:6s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}"
            .format('ka', 'Z', 'Ne', 'T in eV', 'U(e-Z)', 'U(e-e)', 'U(Xa)', ELabel))
    T, UeZ, Uee, UXa, Etot = calTotalEnergy(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
    print("{:6.4f}\t{:6.4f}\t{:6.4f}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}"
            .format(ka, Z, Ne, 
                    T * HartreeToeV, UeZ * HartreeToeV, Uee * HartreeToeV, UXa * HartreeToeV, Etot * HartreeToeV))

    if 'v' in mode:
        print("")
        print("Optimized for ka by variational principle")

        print("{:6s}\t{:6s}\t{:6s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}"
                .format('ka', 'Z', 'Ne', 'T in eV', 'U(e-Z)', 'U(e-e)', 'U(Xa)', ELabel))
        T, UeZ, Uee, UXa, Etot, ka = calOptimizedTotalEnergy(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
        print("{:6.4f}\t{:6.4f}\t{:6.4f}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}"
            .format(ka, Z, Ne, 
                   T * HartreeToeV, UeZ * HartreeToeV, Uee * HartreeToeV, UXa * HartreeToeV, Etot * HartreeToeV))

#=============================
# Plot graphs
#=============================
    if 'g' not in mode:
        terminate()

    build_Qr(ka0, Z, n, l, m)

    ypotZ    = []  # U(Z) = -Z/r    
    ypotZrho = []  # U(e-Z) for Z = 1
    ypotXa   = []  # Xa potential
    ypotrho  = []  # Electrostatic potential by rho(r)
    for i in range(len(r)):
        if r[i] == 0.0:
            phi = 0.0
        else:
            phi, errpot = integrate.quad(lambda r: pi4 * r * rho(ka, Z, n, l, m, r), 
                                    Rmin, r[i], limit = nmaxdiv, epsrel = eps)
            potZ = -Z / r[i]
        ypotZ.append(calUZ(r[i], Z))
        ypotZrho.append(-Z * Ne * phi)
        ypotrho.append(Ne * calUrho(r[i], ka, Z, n, l, m))
        ypotXa.append(calLocalUXa(r[i], ka, Z, n, l, m, 1.0))

    fig = plt.figure(figsize = (12, 8))

    ax1 = fig.add_subplot(2, 3, 1)
    ax2 = fig.add_subplot(2, 3, 2)
    ax3 = fig.add_subplot(2, 3, 3)
    ax1.plot(r, yRr, label = 'Rr(r)', color = 'black')
    ax1.set_xlim([Rmin, Rmax])
    ax1.set_ylim([0, max(yRr)*1.1])
    ax1.set_xlabel("r (bohr)")
    ax1.set_ylabel("Rr(r)")
    ax1.legend()
    ax2.plot(r, yQr, label = 'Q(r)', color = 'black')
    ax2.set_xlabel("r (bohr)")
    ax2.set_ylabel("Q(r)")
    ax2.legend()
    ax3.plot(r, ypotZ,    label = 'U(Z)', color = 'black')
    ax3.plot(r, ypotZrho, label = 'U(Z-rho)', color = 'red')
    ax3.plot(r, -np.array(ypotrho),  label = '-U(rho)',   color = 'blue')
    ax3.plot(r, ypotXa,   label = 'U(Xa)',    color = 'green')
    ax3.set_ylim([min(ypotXa) * 5.0, 0.0])
    ax3.set_xlabel("r (bohr)")
    ax3.set_ylabel("Potential / Energy (Hartree)")
    ax3.legend()

    plt.tight_layout()
    plt.pause(0.1)

    print("Press ENTER to exit>>", end = '')
    input()


def main():
    """
    概要: プログラムのメインエントリポイントです。

    詳細説明:
    コマンドライン引数を解析し、グローバル変数を更新します。
    水素様原子の解析解エネルギーを計算した後、
    指定されたモードに基づいて異なる計算（デバッグ、ka掃引、Ne掃引）を実行します。
    最後に、動径関数の正規化チェックも行います。

    :returns: なし
    :rtype: None
    """
    global mode
    global ka, Z, n, l, m, Ne
    global Rmin, Rstep, nR, Rmaxdata, r
    global qfunc

    updatevars()

    Rstep = (Rmax - Rmin) / (nR - 1)
    r = [Rmin + i * Rstep for i in range(nR+100)]
    Rmaxdata = max(r)

    print("")
    print("mode: ", mode)
    print("")
    print("Orbital: ka={} Z={} n={} l={} m={}".format(ka, Z, n, l, m))
    print("Ne: ", Ne)
    print("Integration: Rmax=", Rmax)
    print("   Rmax: epsR={}  Rmaxinteg={}".format(epsR, -log(epsR) / Z / ka))

    print("")
    print("Analytical solution")

    Tanal = calTanal(Z, ka)
    Uanal = calUanal(Z, ka)
    Eanal = Tanal + Uanal
    print("T(analytical) = {} eV".format(Tanal))
    print("U(analytical) = {} eV".format(Uanal))
    print("  Etotl(analytical) = {} eV".format(Eanal))

    print("")
    print("Numerical integration")
    Rr2tot, err = integrate.quad(lambda r: r * r * rho(ka, Z, n, l, m, r), Rmin, Rmax, limit = nmaxdiv, epsrel = eps)
    Rr2tot *= pi4
    print("R(r) normalization check: 2pi * integ(r*r * Rr*2)dr = ", Rr2tot)

    if 'd' in mode:
        debug(Eanal)
    if 'k' in mode:
        sweep_ka(Eanal)
    if 'n' in mode:
        sweep_Ne(Eanal)

    terminate()


if __name__ == '__main__':
    usage()
    main()