SPECTRAL SIMULATION CONSOLE

import numpy as np
from scipy.linalg import eigvalsh

N       = 128          # matrix dimension
GAMMA_C = 50          # zeta zero cutoff  γ

# Construct GUE ensemble seeded by Riemann zeros
rng   = np.random.default_rng(seed=0xANOR)
H     = rng.standard_normal((N, N)) + 1j*rng.standard_normal((N, N))
H     = (H + H.conj().T) / np.sqrt(2 * N)

# Fold in zeta zeros up to cutoff γ
zeros = load_riemann_zeros(cutoff=GAMMA_C)
for k, gamma in enumerate(zeros):
    H[k % N, k % N] += gamma / (2*np.pi)

spectrum = eigvalsh(H)
density  = neutronic_profile(spectrum, bins=256)