%% Wigner surmises

ws1 = @(s) pi*s/2 .* exp(-s.^2 * pi/4);
ws2 = @(s) 32/pi^2*s.^2 .* exp(-s.^2 * 4/pi);
ws4 = @(s) 2^18/3^6/pi^3 * s.^4 .* exp(-s.^2 * 64/9/pi);

%% Plot the empirical eigenvalue spacing distribution

RUNS = 10000;
N = 100;
a = round(0.3*N);
b = round(0.7*N);

spacings = zeros(N-1, RUNS);

scLaw = @(x) sqrt(4-x.^2)/2/pi;

for j = 1:RUNS
    X = GUE(N);
    sp = eig(X);
    sps = sort(sp);
    %sps = sps(2:2:end);  % for GSE only
    spacings(:,j) = N * diff( sps .* scLaw(sps) );
end;

bulkSpacings = spacings(a:b, :);

%%
figure;
ttt = bulkSpacings * diag(1./mean(bulkSpacings));
bS = ttt(1:end);
%histogram(bSO, 100, 'Normalization', 'pdf');
[n, xout] = hist(bS, 100);
bar(xout, n / length(bS) / (xout(2)-xout(1)));
xlim([0,3]);
hold on;

ss = linspace(0,3);
plot(ss, ws2(ss), 'r-');

%histogram(bulkSpacingsU(1:end), 100, 'Normalization', 'pdf');

%% Circular ensembles

RUNS = 1000;
N = 50;

sp = zeros(N, RUNS);

for j = 1:RUNS
    spec = eig(CUE(N));
    angles = sort(angle(spec));
    sp(:,j) = [diff(angles); angles(1) + 2*pi - angles(end)] * N / (2*pi);
end;

figure;
sp1 = sp(1:end);
%histogram(sp(1:end), 25, 'Normalization', 'pdf');
[n,xout] = hist(sp1, 100);
bar(xout, n / length(sp1) / (xout(2)-xout(1)));
hold on;
xlim([0,3]);

ss = linspace(0,3);
plot(ss, ws2(ss), 'k-');