# -*- coding: utf-8 -*-
"""Chapitre 20 — Simulation stochastique, version corrigée."""
import math
import random

import matplotlib.pyplot as plt


lam = 0.01   # constante radioactive (s^-1)
N0 = 2000    # noyaux présents à t = 0

proba = lam  # pour un pas de 1 s et lam petit : p ≈ lam

ts = [0]
Ns = [N0]
while Ns[len(Ns) - 1] > N0 / 100:
    n_desintegrations = 0
    for i in range(Ns[len(Ns) - 1]):
        if random.random() <= proba:
            n_desintegrations = n_desintegrations + 1
    Ns.append(Ns[len(Ns) - 1] - n_desintegrations)
    ts.append(ts[len(ts) - 1] + 1)

N_exact = []
for t in ts:
    N_exact.append(N0 * math.exp(-lam * t))

plt.plot(ts, Ns, ".", markersize=3, label="simulation (hasard)")
plt.plot(ts, N_exact, "-", label="loi du cours : $N_0\\,e^{-\\lambda t}$")
plt.xlabel("t (s)")
plt.ylabel("N (noyaux restants)")
plt.legend()
plt.grid(True, alpha=0.3)
plt.title("Décroissance radioactive : le hasard redonne la loi")
plt.show()

i_moitie = 0
while Ns[i_moitie] > N0 / 2:
    i_moitie = i_moitie + 1
print("demi-vie simulée  ~", ts[i_moitie], "s")
print("demi-vie théorique =", math.log(2) / lam, "s")