Haar measure on SO(n)

import numpy as np
import numpy.linalg as la
import scipy.stats as ss
import matplotlib.pyplot as plt


d = 20 # dimesion of matrices
n = 1000 # number of matrices
h = 500 # number of Histogram bins for plotting only


evals = np.zeros((d,n))
# seps = np.zeros((d-1,n))
for i in range(n):
    # a = ss.special_ortho_group.rvs(d)
    tmp = np.random.randn(d,d)
    q, r = la.qr(tmp)
    a = q @ np.diag(np.sign(np.diag(r)))
    # a[:,0] = a[:,0] * np.sign(la.det(a)) # (un)comment this line to toggle between SO(n) and O(n)
    # Haar measure:
    b,_ = la.eig(a)
    evals[:,i] = np.sort(np.angle(b))


# %%
seps = np.diff(evals,axis=0)*d/(2*np.pi)

plt.hist(evals.flatten(), bins=h)
plt.title(f'Dimension: {d}, Number: {n}, HistogramBins: {h}')
plt.show()


plt.hist(seps.flatten(), bins=h)
plt.title(f'Dimension: {d}, Number: {n}, HistogramBins: {h}')
plt.show()