熱力学第1法則（エネルギー保存）

import numpy as np
import matplotlib.pyplot as plt

R = 8.314  # 気体定数 [J/(mol·K)]
n = 1.0    # 物質量 [mol]
cv = 3/2 * R  # 単原子理想気体の定容モル比熱
cp = 5/2 * R  # 単原子理想気体の定圧モル比熱

T1, T2 = 300, 500  # 温度 [K]

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

# --- (1) 定容過程と定圧過程のp-V図 ---
# 初期状態: p1, V1
p1 = 101325  # 1 atm [Pa]
V1 = n * R * T1 / p1  # 理想気体

# 定容過程（V = const）
V_cv = V1
p_cv_1 = p1
p_cv_2 = n * R * T2 / V_cv

# 定圧過程（p = const）
p_cp = p1
V_cp_1 = V1
V_cp_2 = n * R * T2 / p_cp

# プロット
V_range = np.linspace(V1 * 0.8, V_cp_2 * 1.2, 200)

# 等温線
p_T1 = n * R * T1 / V_range
p_T2 = n * R * T2 / V_range

axes[0].plot(V_range * 1000, p_T1 / 1000, 'b--', alpha=0.3, label=f'$T = {T1}$ K')
axes[0].plot(V_range * 1000, p_T2 / 1000, 'r--', alpha=0.3, label=f'$T = {T2}$ K')

# 定容過程
axes[0].plot([V_cv * 1000, V_cv * 1000], [p_cv_1 / 1000, p_cv_2 / 1000],
            'g-', linewidth=2.5, label='定容過程')
axes[0].annotate('', xy=(V_cv * 1000, p_cv_2 / 1000),
                xytext=(V_cv * 1000, p_cv_1 / 1000),
                arrowprops=dict(arrowstyle='->', color='green', lw=2))

# 定圧過程
axes[0].plot([V_cp_1 * 1000, V_cp_2 * 1000], [p_cp / 1000, p_cp / 1000],
            'm-', linewidth=2.5, label='定圧過程')
axes[0].annotate('', xy=(V_cp_2 * 1000, p_cp / 1000),
                xytext=(V_cp_1 * 1000, p_cp / 1000),
                arrowprops=dict(arrowstyle='->', color='purple', lw=2))

axes[0].set_xlabel('Volume $V$ [L]')
axes[0].set_ylabel('Pressure $p$ [kPa]')
axes[0].set_title('$p$-$V$ Diagram')
axes[0].legend(fontsize=8)
axes[0].grid(True, alpha=0.3)

# --- (2) エネルギー収支の棒グラフ ---
# 定容過程
Q_cv = n * cv * (T2 - T1)
W_cv = 0
dU_cv = Q_cv - W_cv

# 定圧過程
Q_cp = n * cp * (T2 - T1)
W_cp = p_cp * (V_cp_2 - V_cp_1)
dU_cp = Q_cp - W_cp

x = np.arange(2)
width = 0.25

bars1 = axes[1].bar(x - width, [Q_cv, Q_cp], width, label='$Q$（熱）', color='red', alpha=0.7)
bars2 = axes[1].bar(x, [dU_cv, dU_cp], width, label='$\\Delta U$（内部エネルギー変化）', color='blue', alpha=0.7)
bars3 = axes[1].bar(x + width, [W_cv, W_cp], width, label='$W$（仕事）', color='green', alpha=0.7)

axes[1].set_xticks(x)
axes[1].set_xticklabels(['定容過程', '定圧過程'])
axes[1].set_ylabel('Energy [J]')
axes[1].set_title('Energy Balance: $\\Delta U = Q - W$')
axes[1].legend(fontsize=8)
axes[1].grid(True, alpha=0.3, axis='y')

# 数値を棒の上に表示
for bar_group in [bars1, bars2, bars3]:
    for bar in bar_group:
        height = bar.get_height()
        if height > 0:
            axes[1].text(bar.get_x() + bar.get_width()/2., height,
                        f'{height:.0f}', ha='center', va='bottom', fontsize=8)

# --- (3) 温度に対する内部エネルギーとエンタルピー ---
T_range = np.linspace(200, 600, 200)
U = n * cv * T_range     # U(T) = ncvT（定数を除く）
H = n * cp * T_range     # H(T) = ncpT（定数を除く）

axes[2].plot(T_range, U, 'b-', linewidth=2, label=f'$U = nc_vT$ ($c_v = {cv:.1f}$ J/(mol·K))')
axes[2].plot(T_range, H, 'r-', linewidth=2, label=f'$H = nc_pT$ ($c_p = {cp:.1f}$ J/(mol·K))')
axes[2].fill_between(T_range, U, H, alpha=0.15, color='green', label='$pV = nRT$')
axes[2].axvline(x=T1, color='gray', linestyle='--', alpha=0.5)
axes[2].axvline(x=T2, color='gray', linestyle='--', alpha=0.5)
axes[2].set_xlabel('Temperature $T$ [K]')
axes[2].set_ylabel('Energy [J]')
axes[2].set_title('$U(T)$ and $H(T)$ for Ideal Gas')
axes[2].legend(fontsize=8)
axes[2].grid(True, alpha=0.3)

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

# 数値結果の表示
print("=== 定容過程 ===")
print(f"Q_V = {Q_cv:.1f} J")
print(f"W   = {W_cv:.1f} J")
print(f"ΔU  = {dU_cv:.1f} J")
print(f"検証: Q - W = {Q_cv - W_cv:.1f} = ΔU ✓")

print("\n=== 定圧過程 ===")
print(f"Q_p = {Q_cp:.1f} J")
print(f"W   = {W_cp:.1f} J")
print(f"ΔU  = {dU_cp:.1f} J")
print(f"検証: Q - W = {Q_cp - W_cp:.1f} = ΔU ✓")
print(f"ΔH  = Q_p = {Q_cp:.1f} J")