核分裂連鎖反応をPythonで動かす：中性子・臨界質量・元素別比較

import numpy as npimport matplotlib.pyplot as plt
np.random.seed(2)L = 10.0steps = 200
pos = np.array([5.0, 5.0])           # 中央スタートangle = np.random.uniform(0, 2*np.pi)v = np.array([np.cos(angle), np.sin(angle)])trajectory = [pos.copy()]
for _ in range(steps):    pos = pos + v * 0.1    # 壁で反射（簡略化）    if pos[0] < 0 or pos[0] > L:        v[0] *= -1    if pos[1] < 0 or pos[1] > L:        v[1] *= -1    # ときどきランダム散乱    if np.random.random() < 0.05:        angle = np.random.uniform(0, 2*np.pi)        v = np.array([np.cos(angle), np.sin(angle)])    trajectory.append(pos.copy())
trajectory = np.array(trajectory)plt.figure(figsize=(6, 6))plt.plot(trajectory[:, 0], trajectory[:, 1], lw=0.8, color="goldenrod")plt.scatter(*trajectory[0], c="green", s=100, label="start", zorder=5)plt.scatter(*trajectory[-1], c="red", s=100, label="end", zorder=5)plt.xlim(0, L); plt.ylim(0, L)plt.title("Step 1: A single neutron's path")plt.legend(); plt.grid(alpha=0.3)plt.show()

import numpy as np
# 原子核（fissile nuclei）N_NUCLEI = 300L = 10.0nuclei = np.random.rand(N_NUCLEI, 2) * L
# 元素パラメータ（U-235）P_FISSION = 0.95   # 衝突距離内に来たときの分裂確率（規格化）NU = 2.42          # 1分裂で出る中性子の平均数R_HIT = 0.25       # 「衝突」と判定する距離
# 中性子の状態（pos, vel）def make_neutron(pos):    angle = np.random.uniform(0, 2*np.pi)    return [pos.copy(), np.array([np.cos(angle), np.sin(angle)])]
neutrons = [make_neutron(np.array([5.0, 5.0]))]
def step_world(neutrons, nuclei):    new_neutrons = []    survived_nuclei_mask = np.ones(len(nuclei), dtype=bool)    for pos, vel in neutrons:        pos += vel * 0.1        # 箱を出たら消える（漏洩）        if not (0 <= pos[0] <= L and 0 <= pos[1] <= L):            continue        # 原子核との衝突判定        if len(nuclei) > 0:            d2 = np.sum((nuclei - pos)**2, axis=1)            hit = np.argmin(d2)            if d2[hit] < R_HIT**2 and survived_nuclei_mask[hit]:                if np.random.random() < P_FISSION:                    # 核分裂発生 ─ 中性子を ν 個生成（ポアソン）                    n_out = np.random.poisson(NU)                    survived_nuclei_mask[hit] = False                    for _ in range(n_out):                        new_neutrons.append(make_neutron(pos))                    continue   # この中性子は消費された        new_neutrons.append([pos, vel])    return new_neutrons, nuclei[survived_nuclei_mask]
# 200 ステップ走らせるcounts = []for t in range(200):    neutrons, nuclei = step_world(neutrons, nuclei)    counts.append(len(neutrons))    if len(neutrons) == 0:        break

# 規格化したシミュレーション用パラメータELEMENTS = {    "U-235":  {"P": 0.78, "nu": 2.42, "color": "tab:blue"},    "Pu-239": {"P": 1.00, "nu": 2.87, "color": "tab:red"},    "U-233":  {"P": 0.71, "nu": 2.49, "color": "tab:green"},    "U-238":  {"P": 0.00003, "nu": 2.5, "color": "tab:gray"},    "Th-232": {"P": 0.0, "nu": 0.0, "color": "tab:olive"},}
def run(element, T_steps=200):    p = ELEMENTS[element]    nuclei = np.random.rand(N_NUCLEI, 2) * L    neutrons = [make_neutron(np.array([5.0, 5.0]))]    counts = []    for _ in range(T_steps):        neutrons, nuclei = step_world_with(neutrons, nuclei, p["P"], p["nu"])        counts.append(len(neutrons))    return counts
plt.figure(figsize=(10, 5))for el, p in ELEMENTS.items():    plt.semilogy(run(el), label=f"{el} (P={p['P']}, ν={p['nu']})", color=p["color"], linewidth=2)plt.xlabel("Time steps"); plt.ylabel("Neutron count (log)")plt.title("Step 4: Chain reaction by element")plt.legend(loc="upper right"); plt.grid(alpha=0.3, which="both")plt.show()

# 原子核密度を変えて 3 パターンdensities = {"sub-critical (k<1)": 100, "critical (k≈1)": 250, "super-critical (k>1)": 500}
plt.figure(figsize=(10, 5))for label, n in densities.items():    nuclei = np.random.rand(n, 2) * L    neutrons = [make_neutron(np.array([5.0, 5.0]))]    counts = []    for _ in range(200):        neutrons, nuclei = step_world(neutrons, nuclei)        counts.append(len(neutrons))        if len(neutrons) > 5000:            break    plt.semilogy(counts, label=label, linewidth=2)plt.xlabel("Time steps"); plt.ylabel("Neutron count (log)")plt.title("Step 5: Three regimes ─ subcritical, critical, supercritical")plt.legend(); plt.grid(alpha=0.3, which="both")plt.show()

# 制御棒：箱の x=4〜5 の領域に入ったら中性子は吸収def step_with_control_rod(neutrons, nuclei):    new_neutrons = []    survived = np.ones(len(nuclei), dtype=bool)    for pos, vel in neutrons:        pos += vel * 0.1        if not (0 <= pos[0] <= L and 0 <= pos[1] <= L):            continue        if 4 < pos[0] < 5:    # ← 制御棒で吸収            continue        # ... 通常の核分裂判定 ...        new_neutrons.append([pos, vel])    return new_neutrons, nuclei[survived]