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

bose_condensation.py をダウンロード

bose_condensation.py
bose_condensation.py
  1"""
  2Bose-Einstein凝縮の計算と分析を行うモジュール
  3
  4概要:
  5    このモジュールは、理想ボーズ気体におけるBose-Einstein凝縮現象をシミュレーションし、
  6    関連する物理量（化学ポテンシャル、粒子数、凝縮粒子数など）を計算およびプロットします。
  7
  8詳細説明:
  9    コマンドライン引数により、動作モード（Fs-α関数のプロットまたは化学ポテンシャルと粒子数の計算）と
 10    計算範囲を設定できます。主要な機能として、Bose-Einstein積分の計算、
 11    与えられた粒子数密度に対する化学ポテンシャルの決定（二分法を使用）、
 12    そしてそれらの結果のグラフ表示が含まれています。
 13
 14関連リンク: :doc:`bose_condensation_usage`
 15"""
 16
 17import sys
 18import time
 19from math import exp, sqrt
 20import numpy as np
 21from scipy import integrate         # 数値積分関数 integrateを読み込む
 22from scipy import optimize          # newton関数はscipy.optimizeモジュールに入っている
 23from matplotlib import pyplot as plt
 24
 25#===================
 26# 物理定数
 27#===================
 28PI   = 3.14159265358979323846
 29H    = 6.6260755e-34    # Js
 30HBAR = 1.05459e-34      # Js
 31C    = 2.99792458e8     # m/s
 32E    = 1.60218e-19      # C
 33KB   = 1.380658e-23     # JK-1 (J/K)
 34ME   = 9.1093897e-31    # kg
 35MP   = 1.6726231e-27    # kg
 36MN   = 1.67495e-27      # kg
 37
 38#===================
 39# デフォルトパラメータ
 40#===================
 41m_def = 6.64e-27      # kg, 4He  
 42N_def = 2.18e22       # cm^-3, Liq He
 43
 44# Temperature range
 45Tmin_def  =  3.00     # K
 46Tmax_def  =  4.00
 47Tstep_def =  0.01
 48
 49# alpha = -mu / kB / T range
 50alphamin_def  = 0.0
 51alphamax_def  = 2.0
 52alphastep_def = 0.02
 53
 54# 収束条件
 55EPS_DEF   = 1.0e-12
 56NMAXITER = 100
 57
 58def Gamma(sigma: float) -> float:
 59    """
 60    ガンマ関数を計算します。
 61
 62    概要: 階乗関数の実数・複素数への一般化であるガンマ関数を再帰的に計算します。
 63    詳細説明:
 64        この実装は再帰を使用し、`sigma`が1や0.5に近い特殊なケースを処理します。
 65        `sigma`が0.5未満の場合、関数はエラーメッセージを表示してプログラムを終了します。
 66    :param sigma: float: ガンマ関数を計算する入力値。
 67    :returns: float: 与えられた`sigma`に対するガンマ関数の計算結果。
 68    """
 69    if abs(sigma - 1.0) < 1.0e-6:
 70        return 1.0
 71    if abs(sigma - 0.5) < 1.0e-6:
 72        return sqrt(PI)
 73    if sigma < 0.5 - 1.0e-6:
 74        print("Gamma: Abnormal argument sigma = ", sigma)
 75        sys.exit()
 76    return (sigma - 1.0) * Gamma(sigma - 1.0)
 77
 78def IntegFunc(y: float, sigma: float, alpha: float) -> float:
 79    """
 80    Bose-Einstein積分 F_sigma(alpha) の被積分関数を定義します。
 81
 82    概要: Bose-Einstein分布関数に関連する積分の一部を形成する被積分関数を計算します。
 83    :param y: float: 積分変数。
 84    :param sigma: float: 関数の次数。
 85    :param alpha: float: 無次元パラメータ（-mu/(kBT)）。
 86    :returns: float: 被積分関数の計算結果。
 87    """
 88    if y + alpha > 100.0:
 89        return 0.0
 90    # y=0付近での発散を避けるため微小値を加算する等の処理はquadに任せる
 91    return pow(y, sigma - 1.0) / (exp(y + alpha) - 1.0)
 92
 93def Fsalpha(sigma: float, alpha: float, Emax: float = 10.0, eps: float = 1.0e-8) -> float:
 94    """
 95    Bose-Einstein積分 F_sigma(alpha) を数値的に計算します。
 96
 97    概要: Bose-Einstein関数 F_sigma(alpha) を数値積分を用いて計算します。
 98    """
 99    # 0付近の特異性に対応するため points=[0] を指定する場合もあるが、quadのデフォルトで概ね動作する
100    ret = integrate.quad(lambda y: IntegFunc(y, sigma, alpha), 0.0, Emax, epsrel = eps)
101    return 1.0 / Gamma(sigma) * ret[0]
102
103def Ne_func(EF: float, T: float, mass: float, Emax: float = 10.0, eps: float = 1.0e-8) -> float:
104    """
105    特定の化学ポテンシャル (EF) と温度 (T) におけるボーズ粒子の密度を計算します。
106    """
107    lambdaT = H / sqrt(2.0 * PI * mass * KB * T)
108    alpha = -EF / (KB * T / E)
109    # 粒子数密度を [cm^-3] 単位に変換 (* 1.0e-6)
110    density = 1.0 / pow(lambdaT, 3.0) * Fsalpha(1.5, alpha, Emax = Emax, eps = eps) * 1.0e-6
111    return density
112
113def CalEF(T: float, target_N: float, mass: float, EFmin: float, EFmax: float, eps_val: float) -> float:
114    """
115    二分法を用いて、指定された温度と粒子数密度に対する化学ポテンシャルを決定します。
116    """
117    dQmin = Ne_func(EFmin, T, mass) - target_N
118    dQmax = Ne_func(EFmax, T, mass) - target_N
119
120    if dQmin * dQmax > 0.0:
121        print(f"Error: Initial range [{EFmin}, {EFmax}] does not bracket the root (dQmin*dQmax > 0)")
122        return 0
123
124    EFhalf = (EFmin + EFmax) / 2.0
125    for i in range(NMAXITER):
126        EFhalf = (EFmin + EFmax) / 2.0
127        dQhalf = Ne_func(EFhalf, T, mass) - target_N
128        if abs(EFmin - EFhalf) < eps_val and abs(EFmax - EFhalf) < eps_val:
129            break
130        if dQmin * dQhalf < 0.0:
131            EFmax = EFhalf
132            dQmax = dQhalf
133        else:
134            EFmin = EFhalf
135            dQmin = dQhalf
136    else:
137        print("Warning: Convergence did not reach max iterations")
138    return EFhalf
139
140def ExecFs(params):
141    """Fs-α曲線と粒子数密度Nのプロットを実行します。"""
142    Tmin, Tmax, Tstep = params['Tmin'], params['Tmax'], params['Tstep']
143    alphamin, alphamax, alphastep = params['alphamin'], params['alphamax'], params['alphastep']
144    mass, target_N, zeta32 = params['mass'], params['target_N'], params['zeta32']
145
146    nT = int((Tmax - Tmin) / Tstep + 1.00001)
147    nalpha = int((alphamax - alphamin) / alphastep + 1.00001)
148
149    fig, ax1 = plt.subplots()
150    ax1.set_xlabel("alpha = -mu / kB / T")
151    ax1.set_ylabel("Ne (cm-3)")
152
153    for i in range(nT):
154        T = Tmax - i * Tstep
155        lambdaT = H / sqrt(2.0 * PI * mass * KB * T)
156        Nmax = 1.0 / pow(lambdaT, 3.0) * zeta32 * 1.0e-6
157        print(f"T = {T:g} K  lambda_T = {lambdaT * 1e9:g} nm  Nmax = {Nmax:g} cm-3")
158
159        xalpha, yNe = [], []
160        for ia in range(nalpha):
161            alpha = alphamin + ia * alphastep
162            EF = -alpha * KB * T / E
163            ne = Ne_func(EF, T, mass) 
164            xalpha.append(alpha)
165            yNe.append(ne)
166        
167        ax1.plot(xalpha, yNe, label=f'Ne({T:g} K)', linewidth=0.8)
168
169    ax1.plot([alphamin, alphamax], [target_N, target_N], color='red', linewidth=1.0, linestyle='dashed', label='Target N')
170    ax1.set_xlim([alphamin, alphamax])
171    ax1.legend()
172    plt.title("Particle Density vs Alpha")
173    plt.show(block=False)
174    input("Press ENTER to exit>>")
175
176def ExecMu(params):
177    """化学ポテンシャル (mu) および関連するパラメータの計算とプロットを実行します。"""
178    Tmin, Tmax, Tstep = params['Tmin'], params['Tmax'], params['Tstep']
179    mass, target_N, zeta32, eps_val = params['mass'], params['target_N'], params['zeta32'], params['eps']
180
181    Tc = H * H / 2.0 / PI / mass / KB * pow(target_N * 1.0e6 / zeta32, 2.0/3.0)
182    nT = int((Tmax - Tmin) / Tstep + 1.00001)
183
184    xT, yEF, yEFapprox = [], [], []
185    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)")
186    
187    for i in range(nT):
188        T = Tmax - i * Tstep
189        lambdaT = H / sqrt(2.0 * PI * mass * KB * T)
190        Nmax = 1.0 / pow(lambdaT, 3.0) * zeta32 * 1.0e-6
191        
192        if T < Tc:
193            EF = 0.0
194            EFapprox = 0.0
195            Ncal = target_N * pow(T / Tc, 1.5)
196            N0 = target_N - Ncal
197            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}")
198        else:
199            EFmin, EFmax = -2.0, -1.0e-100
200            EF = CalEF(T, target_N, mass, EFmin, EFmax, eps_val)
201            t_rel = (T - Tc) / Tc
202            Aapprox = 9.0 / 16.0 / PI * zeta32 * zeta32 * t_rel * t_rel
203            EFapprox = -Aapprox * KB * T / E
204            Ncal = Ne_func(EF, T, mass)
205            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}")
206            
207        xT.append(T)
208        yEF.append(-EF / (KB * T / E))
209        yEFapprox.append(-EFapprox / (KB * T / E))
210
211    fig, ax1 = plt.subplots()
212    ax1.plot(xT, yEF, label='alpha (calculated)')
213    ax1.plot(xT, yEFapprox, label='alpha (approx)', linestyle='--', color='red')
214    ax1.set_xlabel("T (K)")
215    ax1.set_ylabel("alpha = -mu/(kBT)")
216    ax1.set_xlim([Tmin, Tmax])
217    ax1.legend()
218    plt.title("Chemical Potential vs Temperature")
219    plt.show(block=False)
220    input("Press ENTER to exit>>")
221
222def main():
223    """スクリプトの主要な実行フローを制御します。"""
224    # 初期値の設定
225    params = {
226        'mass': m_def, 'target_N': N_def,
227        'Tmin': Tmin_def, 'Tmax': Tmax_def, 'Tstep': Tstep_def,
228        'alphamin': alphamin_def, 'alphamax': alphamax_def, 'alphastep': alphastep_def,
229        'eps': EPS_DEF
230    }
231    mode = 'Fs'
232
233    # コマンドライン引数の処理
234    argv = sys.argv
235    if len(argv) > 1:
236        mode = argv[1]
237        if mode == 'mu':
238            if len(argv) >= 3: params['Tmin'] = float(argv[2])
239            if len(argv) >= 4: params['Tmax'] = float(argv[3])
240            if len(argv) >= 5: params['Tstep'] = float(argv[4])
241        elif mode == 'Fs':
242            if len(argv) >= 3: params['Tmin'] = float(argv[2])
243            if len(argv) >= 4: params['Tmax'] = float(argv[3])
244            if len(argv) >= 5: params['Tstep'] = float(argv[4])
245            if len(argv) >= 6: params['alphamin'] = float(argv[5])
246            if len(argv) >= 7: params['alphamax'] = float(argv[6])
247            if len(argv) >= 8: params['alphastep'] = float(argv[7])
248    else:
249        print("Usage: python bose_condensation.py mu Tmin Tmax Tstep")
250        print("Usage: python bose_condensation.py Fs Tmin Tmax Tstep alphamin alphamax alphastep")
251        return
252
253    # 基本定数の計算
254    params['zeta32'] = Fsalpha(1.5, 0.0)
255    print(f"mass = {params['mass']} kg, N = {params['target_N']} cm^-3")
256    print(f"G(1.5) = {Gamma(1.5):.6f}, F3/2(0) = {params['zeta32']:.6f}")
257    
258    Tc = H * H / 2.0 / PI / params['mass'] / KB * pow(params['target_N']*1.0e6 / params['zeta32'], 2.0/3.0)
259    print(f"Calculated Tc = {Tc:.6f} K\n")
260
261    if mode == 'Fs':
262        ExecFs(params)
263    elif mode == 'mu':
264        ExecMu(params)
265
266if __name__ == '__main__':
267    main()