理想気体の状態方程式

import numpy as np
import matplotlib.pyplot as plt

R = 8.314  # 気体定数 [J/(mol·K)]

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# --- (1) ボイル・シャルル・アボガドロの法則 ---
# ボイルの法則: pV = const (T固定)
V = np.linspace(0.5, 5, 200)
temperatures = [200, 300, 400, 500]
n = 1.0

for T in temperatures:
    p = n * R * T / V  # 圧力 [Pa]（Vはリットル相当にスケーリング）
    axes[0].plot(V, p, linewidth=2, label=f'$T = {T}$ K')

axes[0].set_xlabel('Volume $V$ [arb. units]')
axes[0].set_ylabel('Pressure $p$ [arb. units]')
axes[0].set_title("Boyle's Law: $pV = nRT$")
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(0.5, 5)
axes[0].set_ylim(0, 8000)

# --- (2) 理想気体 vs ファンデルワールス気体 (CO2) ---
a_CO2 = 0.3640    # Pa·m^6/mol^2
b_CO2 = 4.267e-5  # m^3/mol
T_list = [250, 304.2, 400]  # 304.2 Kは臨界温度
colors_vdw = ['blue', 'red', 'green']

Vm = np.linspace(8e-5, 5e-4, 500)  # モル体積 [m^3/mol]

for T, color in zip(T_list, colors_vdw):
    # 理想気体
    p_ideal = R * T / Vm
    # ファンデルワールス
    p_vdw = R * T / (Vm - b_CO2) - a_CO2 / Vm**2

    axes[1].plot(Vm * 1e6, p_ideal / 1e6, color=color, linestyle='--', alpha=0.5)
    axes[1].plot(Vm * 1e6, p_vdw / 1e6, color=color, linewidth=2,
                label=f'$T = {T}$ K')

axes[1].set_xlabel('Molar Volume $V_m$ [cm$^3$/mol]')
axes[1].set_ylabel('Pressure $p$ [MPa]')
axes[1].set_title('Ideal Gas (dashed) vs van der Waals ($CO_2$)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
axes[1].set_ylim(0, 20)

# 臨界点の計算
Vc = 3 * b_CO2
pc = a_CO2 / (27 * b_CO2**2)
Tc = 8 * a_CO2 / (27 * R * b_CO2)
axes[1].plot(Vc * 1e6, pc / 1e6, 'ko', markersize=8, zorder=5)
axes[1].annotate(f'臨界点\n$T_c={Tc:.1f}$ K', (Vc * 1e6, pc / 1e6),
                textcoords="offset points", xytext=(30, -10), fontsize=9,
                arrowprops=dict(arrowstyle='->', color='black'))

# --- (3) 分子速度分布（マクスウェル-ボルツマン分布） ---
kB = 1.381e-23  # ボルツマン定数

def maxwell_boltzmann(v, m, T):
    """マクスウェル-ボルツマン速度分布"""
    return 4 * np.pi * (m / (2 * np.pi * kB * T))**1.5 * v**2 * np.exp(-m * v**2 / (2 * kB * T))

gases = {
    'He ($M=4$)': {'M': 4e-3, 'color': 'blue'},
    '$N_2$ ($M=28$)': {'M': 28e-3, 'color': 'red'},
    '$CO_2$ ($M=44$)': {'M': 44e-3, 'color': 'green'},
}

v = np.linspace(0, 2500, 500)
T = 300  # K

for name, props in gases.items():
    m = props['M'] / 6.022e23  # 1分子の質量
    f_v = maxwell_boltzmann(v, m, T)
    v_rms = np.sqrt(3 * kB * T / m)

    axes[2].plot(v, f_v * 1e3, color=props['color'], linewidth=2, label=name)
    axes[2].axvline(x=v_rms, color=props['color'], linestyle='--', alpha=0.5)

axes[2].set_xlabel('Speed $v$ [m/s]')
axes[2].set_ylabel('Probability density $f(v)$ [×10$^{-3}$]')
axes[2].set_title(f'Maxwell-Boltzmann Distribution ($T = {T}$ K)')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('ideal_gas_equation.png', dpi=150, bbox_inches='tight')
plt.show()

# 数値結果
print("=== 臨界点（CO2、ファンデルワールス方程式から） ===")
print(f"臨界温度 Tc = {Tc:.1f} K（実測値: 304.2 K）")
print(f"臨界圧力 pc = {pc/1e6:.2f} MPa（実測値: 7.38 MPa）")
print(f"臨界モル体積 Vc = {Vc*1e6:.1f} cm³/mol")

print("\n=== 二乗平均速度 (T = 300 K) ===")
for name, props in gases.items():
    v_rms = np.sqrt(3 * R * T / props['M'])
    print(f"{name}: v_rms = {v_rms:.1f} m/s")