Bethe-Weizsäcker 結合エネルギー曲線, Bethe-Weizsäcker 結合エネルギー曲線の各項の寄与 のplot (Grokによるコード)

plot.py

# %%
import numpy as np
import matplotlib.pyplot as plt

import matplotlib.font_manager as fm
import matplotlib
prop = fm.FontProperties(fname='/System/Library/Fonts/ヒラギノ角ゴシック W3.ttc')
prop.set_weight = 'normal'
matplotlib.rcParams['font.family'] = prop.get_name()
matplotlib.rcParams['font.weight'] = 'normal'
# %%
# 標準定数（観測値にフィットした値）
a_v = 15.5   # 体積項
a_s = 16.8   # 表面項
a_c = 0.72   # クーロン項
a_a = 23     # 非対称項
a_p = 0      # ペアリング項（滑らか曲線のため0）


def binding_energy_per_nucleon(A):
    if A < 2:
        return 0.0

    # 最安定Z（連続値）
    Z = A / (2 * (1 + (a_c / (4 * a_a)) * A**(2/3)))

    volume = a_v
    surface = a_s * A**(-1/3)
    coulomb = a_c * Z * (Z - 1) * A**(-4/3)
    asymmetry = a_a * (A - 2*Z)**2 / A**2

    pairing = 0  # 滑らか曲線

    B_per_A = volume - surface - coulomb - asymmetry + pairing
    return B_per_A


# 質量数Aの範囲
A = np.arange(2, 271)
B_per_A = np.vectorize(binding_energy_per_nucleon)(A)

# 先頭20個の確認（負値なし、滑らか増加）
print("B_per_A（最終修正後）:", B_per_A[:20])

# プロット
plt.figure(figsize=(12, 7))
plt.plot(A, B_per_A, 'b-', linewidth=2, label='Bethe-Weizsäcker 曲線（滑らか版）')

# 観測値の例
obs = {4: 7.07,   # He-4
       12: 7.68,  # C-12
       16: 7.98,  # O-16
       56: 8.79,  # Fe-56
       62: 8.80,  # Ni-62
       238: 7.57}  # U-238
plt.scatter(list(obs.keys()), list(obs.values()), color='red', zorder=5, label='観測値')

# Fe付近強調
plt.axvline(x=56, color='green', linestyle='--', alpha=0.7, label='Fe-56（ピーク）')

plt.title('Bethe-Weizsäcker 結合エネルギー曲線\n（核子1個あたりの結合エネルギー vs 質量数）', fontsize=14)
plt.xlabel('質量数 A', fontsize=12)
plt.ylabel('結合エネルギー/核子 B/A [MeV]', fontsize=12)
plt.grid(True, alpha=0.3)
plt.legend(fontsize=11)
plt.xlim(0, 270)
plt.ylim(6, 9.5)

# 注釈：なぜFeでピークか？
plt.text(60, 8.5, '← ここがピーク！\nFe-56付近で最も安定',
         fontsize=11, color='darkgreen',
         bbox=dict(boxstyle="round,pad=0.3", facecolor="lightgreen", alpha=0.7))

plt.tight_layout()
plt.show()

# %%
# 定数（標準値）
a_v = 15.5
a_s = 16.8
a_c = 0.72
a_a = 23.0
a_p = 0.0  # 滑らか曲線


def terms(A):
    Z = A / (2 * (1 + (a_c / (4 * a_a)) * A**(2/3)))  # 最安定Z
    volume = a_v
    surface = -a_s * A**(-1/3)
    coulomb = -a_c * Z*(Z-1) * A**(-4/3)
    asymmetry = -a_a * (A - 2*Z)**2 / A**2
    pairing = 0.0
    total = volume + surface + coulomb + asymmetry + pairing
    return volume, surface, coulomb, asymmetry, pairing, total


# 計算
A = np.arange(2, 271)
V, S, C, Asym, P, Total = np.vectorize(terms, otypes=[float]*6)(A)

# === プロット ===
plt.figure(figsize=(14, 8))

# 各項
plt.plot(A, V, label='体積項 (+a_v)', color='red', linestyle='-', linewidth=2)
plt.plot(A, S, label='表面項 (-a_s A^{-1/3})', color='orange', linestyle='--')
plt.plot(A, C, label='クーロン項 (-a_c Z(Z-1)/A^{4/3})', color='purple', linestyle='-.')
plt.plot(A, Asym, label='非対称項 (-a_a (A-2Z)^2/A^2)', color='brown', linestyle=':')
plt.plot(A, Total, label='合計 (B/A)', color='blue', linewidth=3)

# 観測値
obs = {4: 7.07, 12: 7.68, 16: 7.98, 56: 8.79, 62: 8.80, 238: 7.57}
plt.scatter(obs.keys(), obs.values(), color='black', zorder=5, s=60, label='観測値')

plt.axvline(56, color='green', linestyle='--', alpha=0.7, label='Fe-56 ピーク')

plt.title('Bethe-Weizsäcker 結合エネルギー曲線の各項の寄与\n（物理モデルに基づく分解）', fontsize=14)
plt.xlabel('質量数 A', fontsize=12)
plt.ylabel('エネルギー寄与 [MeV/核子]', fontsize=12)
plt.grid(True, alpha=0.3)
plt.legend(fontsize=11, loc='lower right')
plt.xlim(0, 270)
plt.ylim(-20, 16)  # 負の項も見えるように
plt.tight_layout()
plt.show()