以下は、bose_condensation.py のソースコードを解析し、Sphinx（MyST）で問題なくビルド可能なように記述されたMarkdownドキュメントです。
bose_condensation.py のドキュメント
概要
このドキュメントは、Pythonプログラム bose_condensation.py のソースコードを解析し、その機能、構造、および使用方法について説明します。
モジュールの目的と機能
bose_condensation.py は、理想ボーズ気体におけるボーズ・アインシュタイン凝縮現象をシミュレーションし、関連する物理量を計算およびプロットするためのモジュールです。
概要
このモジュールは、理想ボーズ気体におけるボーズ・アインシュタイン凝縮現象をシミュレーションし、関連する物理量（化学ポテンシャル、粒子数、凝縮粒子数など）を計算およびプロットします。
詳細説明
コマンドライン引数により、動作モード（Fs-α 関数のプロットまたは化学ポテンシャルと粒子数の計算）と計算範囲を設定できます。主要な機能として、ボーズ・アインシュタイン積分の計算、与えられた粒子数密度に対する化学ポテンシャルの決定（二分法を使用）、そしてそれらの結果のグラフ表示が含まれています。
関連情報
コード内のdocstringに 関連リンク: :doc:bose_condensation_usage`` と記述されています。これは、bose_condensation_usage という名前の別のドキュメントファイルへの参照を示唆していますが、その詳細な内容やリンク先のパスは本コードからは確認できません。
非標準ライブラリ
bose_condensation.py は以下の非標準Pythonライブラリを使用しています。これらはプログラムの実行前にインストールが必要です。
numpy: 数値計算を効率的に行うためのライブラリです。scipy: 科学技術計算のためのライブラリで、ここでは数値積分 (scipy.integrate) と数値最適化 (scipy.optimize) の機能が使用されています。matplotlib: グラフ描画のためのライブラリで、計算結果の可視化に使用されています。
物理定数
このモジュールでは、物理学の基本的な定数が定義されています。
PI: 円周率 (3.14159265358979323846)H: プランク定数 (6.6260755e-34J s)HBAR: 換算プランク定数 (1.05459e-34J s)C: 真空中の光速 (2.99792458e8m/s)E: 素電荷 (1.60218e-19C)KB: ボルツマン定数 (1.380658e-23J K^-1)ME: 電子の質量 (9.1093897e-31kg)MP: プロトンの質量 (1.6726231e-27kg)MN: 中性子の質量 (1.67495e-27kg)
デフォルトパラメータ
プログラムの実行時に使用されるデフォルトの物理パラメータと計算範囲です。
m_def: 粒子の質量 (6.64e-27kg, ヘリウム4の質量に相当)N_def: 粒子数密度 (2.18e22cm^-3, 液体ヘリウムの密度に相当)Tmin_def: 温度範囲の最小値 (3.00K)Tmax_def: 温度範囲の最大値 (4.00K)Tstep_def: 温度の刻み幅 (0.01K)alphamin_def:alpha(-mu / kB / T) 範囲の最小値 (0.0)alphamax_def:alpha範囲の最大値 (2.0)alphastep_def:alphaの刻み幅 (0.02)EPS_DEF: 収束判定に使用される微小値 (1.0e-12)NMAXITER: 二分法の最大反復回数 (100)
関数
Gamma(sigma: float) -> float
ガンマ関数を計算します。
概要
階乗関数の実数・複素数への一般化であるガンマ関数を再帰的に計算します。
詳細説明
この実装は再帰を使用し、sigma が 1 や 0.5 に近い特殊なケースを処理します。sigma が 0.5 未満の場合、関数はエラーメッセージを表示してプログラムを終了します。
引数
sigma:float: ガンマ関数を計算する入力値。
戻り値
float: 与えられたsigmaに対するガンマ関数の計算結果。
IntegFunc(y: float, sigma: float, alpha: float) -> float
ボーズ・アインシュタイン積分 Fs_sigma(alpha) の被積分関数を定義します。
概要
ボーズ・アインシュタイン分布関数に関連する積分の一部を形成する被積分関数を計算します。
引数
y:float: 積分変数。sigma:float: 関数の次数。alpha:float: 無次元パラメータ (-mu/(kBT))。
戻り値
float: 被積分関数の計算結果。注:
y + alpha > 100.0の場合、指数関数のオーバーフローを避けるために0.0を返します。
Fsalpha(sigma: float, alpha: float, Emax: float = 10.0, eps: float = 1.0e-8) -> float
ボーズ・アインシュタイン積分 Fs_sigma(alpha) を数値的に計算します。
概要
ボーズ・アインシュタイン関数 Fs_sigma(alpha) を scipy.integrate.quad を用いた数値積分で計算します。
引数
sigma:float: 関数の次数。alpha:float: 無次元パラメータ (-mu/(kBT))。Emax:float, オプション: 積分の最大値。デフォルトは10.0。eps:float, オプション: 相対誤差の許容範囲。デフォルトは1.0e-8。
戻り値
float: 計算されたボーズ・アインシュタイン積分の値。
Ne_func(EF: float, T: float, mass: float, Emax: float = 10.0, eps: float = 1.0e-8) -> float
特定の化学ポテンシャル (EF) と温度 (T) におけるボーズ粒子の密度を計算します。
引数
EF:float: 化学ポテンシャル。T:float: 温度。mass:float: 粒子の質量。Emax:float, オプション:Fsalpha()に渡される積分の最大値。デフォルトは10.0。eps:float, オプション:Fsalpha()に渡される相対誤差の許容範囲。デフォルトは1.0e-8。
戻り値
float: 粒子数密度を[cm^-3]単位で返します。
CalEF(T: float, target_N: float, mass: float, EFmin: float, EFmax: float, eps_val: float) -> float
二分法を用いて、指定された温度と粒子数密度に対する化学ポテンシャルを決定します。
引数
T:float: 温度。target_N:float: 目標とする粒子数密度。mass:float: 粒子の質量。EFmin:float: 化学ポテンシャルの探索範囲の最小値。EFmax:float: 化学ポテンシャルの探索範囲の最大値。eps_val:float: 収束判定に使用される許容誤差。
戻り値
float: 計算された化学ポテンシャルの値。初期範囲が根を囲まない場合、0を返しエラーメッセージを出力します。また、最大反復回数NMAXITERを超過した場合は警告メッセージを出力します。
ExecFs(params)
Fs-α 曲線と粒子数密度 N のプロットを実行します。
引数
params:dict: 以下のキーを含む辞書:'Tmin','Tmax','Tstep': 温度範囲。'alphamin','alphamax','alphastep':alpha(-mu / kB / T) の範囲。'mass': 粒子の質量。'target_N': 目標粒子数密度。'zeta32':Fsalpha(1.5, 0.0)の値。
動作
複数の温度に対する粒子数密度
Neとalphaの関係をプロットします。目標粒子数密度を破線で表示します。
プロット後、ユーザーの入力を待って終了します。
ExecMu(params)
化学ポテンシャル (mu) および関連するパラメータの計算とプロットを実行します。
引数
params:dict: 以下のキーを含む辞書:'Tmin','Tmax','Tstep': 温度範囲。'mass': 粒子の質量。'target_N': 目標粒子数密度。'zeta32':Fsalpha(1.5, 0.0)の値。'eps':CalEF()に渡される収束許容誤差。
動作
指定された温度範囲で、ボーズ・アインシュタイン凝縮の臨界温度
Tcを計算します。T < Tcの場合、化学ポテンシャルEFは0.0となり、凝縮粒子数N0を計算します。T >= Tcの場合、CalEF()を使用して化学ポテンシャルEFを決定し、近似値と比較します。計算結果を標準出力に表形式で出力し、化学ポテンシャル (
alpha) と温度の関係をプロットします。プロット後、ユーザーの入力を待って終了します。
main()
スクリプトの主要な実行フローを制御します。
概要
コマンドライン引数を解析し、それに基づいて ExecFs() または ExecMu() 関数を呼び出します。
動作
デフォルトパラメータを
params辞書に設定します。sys.argvからコマンドライン引数を解析し、実行モード ('Fs'または'mu') および関連する計算範囲パラメータを更新します。引数が不足している場合、使用方法を表示して終了します。Gamma(1.5)とFsalpha(1.5, 0.0)(リーマン・ゼータ関数zeta(3/2)に相当) を計算し、臨界温度Tcを導出します。これらの定数を標準出力に出力します。modeの値に応じてExecFs(params)またはExecMu(params)を呼び出します。
スクリプトの実行方法とコマンドライン引数
このスクリプトは、コマンドライン引数によって動作モードと計算範囲を制御します。
実行形式
python bose_condensation.py <mode> [parameters...]
コマンドライン引数
<mode>: 必須。実行モードを指定します。mu: 化学ポテンシャル (mu) と関連パラメータの計算およびプロットを実行します。Fs:Fs-α(ボーズ・アインシュタイン積分) 曲線と粒子数密度Nのプロットを実行します。
mu モードの場合
構文:
python bose_condensation.py mu Tmin Tmax Tstep
パラメータ:
Tmin:float, オプション: 温度の最小値 (K)。デフォルトはTmin_def(3.00)。Tmax:float, オプション: 温度の最大値 (K)。デフォルトはTmax_def(4.00)。Tstep:float, オプション: 温度の刻み幅 (K)。デフォルトはTstep_def(0.01)。
Fs モードの場合
構文:
python bose_condensation.py Fs Tmin Tmax Tstep alphamin alphamax alphastep
パラメータ:
Tmin:float, オプション: 温度の最小値 (K)。デフォルトはTmin_def(3.00)。Tmax:float, オプション: 温度の最大値 (K)。デフォルトはTmax_def(4.00)。Tstep:float, オプション: 温度の刻み幅 (K)。デフォルトはTstep_def(0.01)。alphamin:float, オプション:alphaの最小値。デフォルトはalphamin_def(0.0)。alphamax:float, オプション:alphaの最大値。デフォルトはalphamax_def(2.0)。alphastep:float, オプション:alphaの刻み幅。デフォルトはalphastep_def(0.02)。
引数不足の場合の出力例
引数が不足している場合、以下のような使用方法のメッセージが標準出力に表示され、プログラムは終了します。
Usage: python bose_condensation.py mu Tmin Tmax Tstep
Usage: python bose_condensation.py Fs Tmin Tmax Tstep alphamin alphamax alphastep
出力
スクリプトは、実行モードに応じて標準出力とグラフィカルなプロットを生成します。
標準出力
初期パラメータと臨界温度 Tc の表示
スクリプトの開始時に、設定された質量、粒子数密度、ガンマ関数値、ボーズ・アインシュタイン積分値、および計算された臨界温度 Tc が標準出力に表示されます。
例:
mass = 6.64e-27 kg, N = 2.18e+22 cm^-3
G(1.5) = 0.886227, F3/2(0) = 2.612375
Calculated Tc = 3.125799 K
ExecFs() 実行時の出力
各温度における T、lambda_T (熱的ドブロイ波長)、および Nmax (最大粒子数密度) が標準出力に表示されます。
例:
T = 4 K lambda_T = 0.443048 nm Nmax = 1.0575e+23 cm-3
T = 3.99 K lambda_T = 0.44358 nm Nmax = 1.0567e+23 cm-3
ExecMu() 実行時の出力
温度 T、熱的ドブロイ波長 lambda_T、化学ポテンシャル EF、EF の近似値 Fapprox、計算された粒子数密度 N'、最大粒子数密度 Nmax、および凝縮粒子数 N0 (T < Tc の場合のみ) が表形式で標準出力に出力されます。
例:
T(K) lambda_T(nm) EF(eV) Fapprox(eV) N'(cm-3) Nmax(cm-3) N0(cm-3)
4.00000 0.4430 -1.077e-07 -1.045e-07 2.180e+22 1.058e+23
3.99000 0.4436 -1.052e-07 -1.020e-07 2.180e+22 1.057e+23
3.12580 0.5042 -1.096e-11 -2.362e-14 2.180e+22 2.180e+22
3.11580 0.5050 0.0 0.0 2.174e+22 2.180e+22 6.004e+19
プロット出力
matplotlib.pyplot を使用してグラフが生成されます。
Fsモード: 「Particle Density vs Alpha」というタイトルで、異なる温度における粒子数密度Neとalpha(-mu / kB / T) の関係がプロットされます。目標粒子数密度は赤色の破線で示されます。muモード: 「Chemical Potential vs Temperature」というタイトルで、化学ポテンシャル (alpha=-mu / kB / T) と温度Tの関係がプロットされます。計算値と近似値が比較表示されます。
いずれのプロットも、表示後にユーザーが ENTER キーを押すまでプログラムは終了しません (input("Press ENTER to exit>>") )。
ソースコード
"""
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()