"""
Bose-Einstein凝縮の計算と分析を行うモジュール
概要:
このモジュールは、理想ボーズ気体におけるBose-Einstein凝縮現象をシミュレーションし、
関連する物理量（化学ポテンシャル、粒子数、凝縮粒子数など）を計算およびプロットします。
詳細説明:
コマンドライン引数により、動作モード（Fs-α関数のプロットまたは化学ポテンシャルと粒子数の計算）と
計算範囲を設定できます。主要な機能として、Bose-Einstein積分の計算、
与えられた粒子数密度に対する化学ポテンシャルの決定（二分法を使用）、
そしてそれらの結果のグラフ表示が含まれています。
関連リンク: :doc:`bose_condensation_usage`
"""
import sys
import time
from math import exp, sqrt
import numpy as np
from scipy import integrate # 数値積分関数 integrateを読み込む
from scipy import optimize # newton関数はscipy.optimizeモジュールに入っている
from matplotlib import pyplot as plt
#===================
# 物理定数
#===================
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-1 (J/K)
ME = 9.1093897e-31 # kg
MP = 1.6726231e-27 # kg
MN = 1.67495e-27 # kg
#===================
# デフォルトパラメータ
#===================
m_def = 6.64e-27 # kg, 4He
N_def = 2.18e22 # cm^-3, Liq He
# Temperature range
Tmin_def = 3.00 # K
Tmax_def = 4.00
Tstep_def = 0.01
# alpha = -mu / kB / T range
alphamin_def = 0.0
alphamax_def = 2.0
alphastep_def = 0.02
# 収束条件
EPS_DEF = 1.0e-12
NMAXITER = 100
[ドキュメント]
def Gamma(sigma: float) -> float:
"""
ガンマ関数を計算します。
概要: 階乗関数の実数・複素数への一般化であるガンマ関数を再帰的に計算します。
詳細説明:
この実装は再帰を使用し、`sigma`が1や0.5に近い特殊なケースを処理します。
`sigma`が0.5未満の場合、関数はエラーメッセージを表示してプログラムを終了します。
:param sigma: float: ガンマ関数を計算する入力値。
:returns: float: 与えられた`sigma`に対するガンマ関数の計算結果。
"""
if abs(sigma - 1.0) < 1.0e-6:
return 1.0
if abs(sigma - 0.5) < 1.0e-6:
return sqrt(PI)
if sigma < 0.5 - 1.0e-6:
print("Gamma: Abnormal argument sigma = ", sigma)
sys.exit()
return (sigma - 1.0) * Gamma(sigma - 1.0)
[ドキュメント]
def IntegFunc(y: float, sigma: float, alpha: float) -> float:
"""
Bose-Einstein積分 F_sigma(alpha) の被積分関数を定義します。
概要: Bose-Einstein分布関数に関連する積分の一部を形成する被積分関数を計算します。
:param y: float: 積分変数。
:param sigma: float: 関数の次数。
:param alpha: float: 無次元パラメータ（-mu/(kBT)）。
:returns: float: 被積分関数の計算結果。
"""
if y + alpha > 100.0:
return 0.0
# y=0付近での発散を避けるため微小値を加算する等の処理はquadに任せる
return pow(y, sigma - 1.0) / (exp(y + alpha) - 1.0)
[ドキュメント]
def Fsalpha(sigma: float, alpha: float, Emax: float = 10.0, eps: float = 1.0e-8) -> float:
"""
Bose-Einstein積分 F_sigma(alpha) を数値的に計算します。
概要: Bose-Einstein関数 F_sigma(alpha) を数値積分を用いて計算します。
"""
# 0付近の特異性に対応するため points=[0] を指定する場合もあるが、quadのデフォルトで概ね動作する
ret = integrate.quad(lambda y: IntegFunc(y, sigma, alpha), 0.0, Emax, epsrel = eps)
return 1.0 / Gamma(sigma) * ret[0]
[ドキュメント]
def Ne_func(EF: float, T: float, mass: float, Emax: float = 10.0, eps: float = 1.0e-8) -> float:
"""
特定の化学ポテンシャル (EF) と温度 (T) におけるボーズ粒子の密度を計算します。
"""
lambdaT = H / sqrt(2.0 * PI * mass * KB * T)
alpha = -EF / (KB * T / E)
# 粒子数密度を [cm^-3] 単位に変換 (* 1.0e-6)
density = 1.0 / pow(lambdaT, 3.0) * Fsalpha(1.5, alpha, Emax = Emax, eps = eps) * 1.0e-6
return density
[ドキュメント]
def CalEF(T: float, target_N: float, mass: float, EFmin: float, EFmax: float, eps_val: float) -> float:
"""
二分法を用いて、指定された温度と粒子数密度に対する化学ポテンシャルを決定します。
"""
dQmin = Ne_func(EFmin, T, mass) - target_N
dQmax = Ne_func(EFmax, T, mass) - target_N
if dQmin * dQmax > 0.0:
print(f"Error: Initial range [{EFmin}, {EFmax}] does not bracket the root (dQmin*dQmax > 0)")
return 0
EFhalf = (EFmin + EFmax) / 2.0
for i in range(NMAXITER):
EFhalf = (EFmin + EFmax) / 2.0
dQhalf = Ne_func(EFhalf, T, mass) - target_N
if abs(EFmin - EFhalf) < eps_val and abs(EFmax - EFhalf) < eps_val:
break
if dQmin * dQhalf < 0.0:
EFmax = EFhalf
dQmax = dQhalf
else:
EFmin = EFhalf
dQmin = dQhalf
else:
print("Warning: Convergence did not reach max iterations")
return EFhalf
[ドキュメント]
def ExecFs(params):
"""Fs-α曲線と粒子数密度Nのプロットを実行します。"""
Tmin, Tmax, Tstep = params['Tmin'], params['Tmax'], params['Tstep']
alphamin, alphamax, alphastep = params['alphamin'], params['alphamax'], params['alphastep']
mass, target_N, zeta32 = params['mass'], params['target_N'], params['zeta32']
nT = int((Tmax - Tmin) / Tstep + 1.00001)
nalpha = int((alphamax - alphamin) / alphastep + 1.00001)
fig, ax1 = plt.subplots()
ax1.set_xlabel("alpha = -mu / kB / T")
ax1.set_ylabel("Ne (cm-3)")
for i in range(nT):
T = Tmax - i * Tstep
lambdaT = H / sqrt(2.0 * PI * mass * KB * T)
Nmax = 1.0 / pow(lambdaT, 3.0) * zeta32 * 1.0e-6
print(f"T = {T:g} K lambda_T = {lambdaT * 1e9:g} nm Nmax = {Nmax:g} cm-3")
xalpha, yNe = [], []
for ia in range(nalpha):
alpha = alphamin + ia * alphastep
EF = -alpha * KB * T / E
ne = Ne_func(EF, T, mass)
xalpha.append(alpha)
yNe.append(ne)
ax1.plot(xalpha, yNe, label=f'Ne({T:g} K)', linewidth=0.8)
ax1.plot([alphamin, alphamax], [target_N, target_N], color='red', linewidth=1.0, linestyle='dashed', label='Target N')
ax1.set_xlim([alphamin, alphamax])
ax1.legend()
plt.title("Particle Density vs Alpha")
plt.show(block=False)
input("Press ENTER to exit>>")
[ドキュメント]
def ExecMu(params):
"""化学ポテンシャル (mu) および関連するパラメータの計算とプロットを実行します。"""
Tmin, Tmax, Tstep = params['Tmin'], params['Tmax'], params['Tstep']
mass, target_N, zeta32, eps_val = params['mass'], params['target_N'], params['zeta32'], params['eps']
Tc = H * H / 2.0 / PI / mass / KB * pow(target_N * 1.0e6 / zeta32, 2.0/3.0)
nT = int((Tmax - Tmin) / Tstep + 1.00001)
xT, yEF, yEFapprox = [], [], []
print("\n T(K) \tlambda_T(nm)\t EF(eV)\t Fapprox(eV)\t N'(cm-3)\t Nmax(cm-3)\t N0(cm-3)")
for i in range(nT):
T = Tmax - i * Tstep
lambdaT = H / sqrt(2.0 * PI * mass * KB * T)
Nmax = 1.0 / pow(lambdaT, 3.0) * zeta32 * 1.0e-6
if T < Tc:
EF = 0.0
EFapprox = 0.0
Ncal = target_N * pow(T / Tc, 1.5)
N0 = target_N - Ncal
print(f"{T:8.5f}\t{lambdaT*1e9:8.4g}\t{EF:12.4g}\t{EFapprox:12.4g}\t{Ncal:12.4e}\t{Nmax:12.4e}\t{N0:12.4e}")
else:
EFmin, EFmax = -2.0, -1.0e-100
EF = CalEF(T, target_N, mass, EFmin, EFmax, eps_val)
t_rel = (T - Tc) / Tc
Aapprox = 9.0 / 16.0 / PI * zeta32 * zeta32 * t_rel * t_rel
EFapprox = -Aapprox * KB * T / E
Ncal = Ne_func(EF, T, mass)
print(f"{T:8.5f}\t{lambdaT*1e9:8.4g}\t{EF:12.4g}\t{EFapprox:12.4g}\t{Ncal:12.4e}\t{Nmax:12.4e}")
xT.append(T)
yEF.append(-EF / (KB * T / E))
yEFapprox.append(-EFapprox / (KB * T / E))
fig, ax1 = plt.subplots()
ax1.plot(xT, yEF, label='alpha (calculated)')
ax1.plot(xT, yEFapprox, label='alpha (approx)', linestyle='--', color='red')
ax1.set_xlabel("T (K)")
ax1.set_ylabel("alpha = -mu/(kBT)")
ax1.set_xlim([Tmin, Tmax])
ax1.legend()
plt.title("Chemical Potential vs Temperature")
plt.show(block=False)
input("Press ENTER to exit>>")
[ドキュメント]
def main():
"""スクリプトの主要な実行フローを制御します。"""
# 初期値の設定
params = {
'mass': m_def, 'target_N': N_def,
'Tmin': Tmin_def, 'Tmax': Tmax_def, 'Tstep': Tstep_def,
'alphamin': alphamin_def, 'alphamax': alphamax_def, 'alphastep': alphastep_def,
'eps': EPS_DEF
}
mode = 'Fs'
# コマンドライン引数の処理
argv = sys.argv
if len(argv) > 1:
mode = argv[1]
if mode == 'mu':
if len(argv) >= 3: params['Tmin'] = float(argv[2])
if len(argv) >= 4: params['Tmax'] = float(argv[3])
if len(argv) >= 5: params['Tstep'] = float(argv[4])
elif mode == 'Fs':
if len(argv) >= 3: params['Tmin'] = float(argv[2])
if len(argv) >= 4: params['Tmax'] = float(argv[3])
if len(argv) >= 5: params['Tstep'] = float(argv[4])
if len(argv) >= 6: params['alphamin'] = float(argv[5])
if len(argv) >= 7: params['alphamax'] = float(argv[6])
if len(argv) >= 8: params['alphastep'] = float(argv[7])
else:
print("Usage: python bose_condensation.py mu Tmin Tmax Tstep")
print("Usage: python bose_condensation.py Fs Tmin Tmax Tstep alphamin alphamax alphastep")
return
# 基本定数の計算
params['zeta32'] = Fsalpha(1.5, 0.0)
print(f"mass = {params['mass']} kg, N = {params['target_N']} cm^-3")
print(f"G(1.5) = {Gamma(1.5):.6f}, F3/2(0) = {params['zeta32']:.6f}")
Tc = H * H / 2.0 / PI / params['mass'] / KB * pow(params['target_N']*1.0e6 / params['zeta32'], 2.0/3.0)
print(f"Calculated Tc = {Tc:.6f} K\n")
if mode == 'Fs':
ExecFs(params)
elif mode == 'mu':
ExecMu(params)
if __name__ == '__main__':
main()