(* Load Noah Graham's Mathematica package of spheroidal functions *)
Import["http://community.middlebury.edu/~ngraham/Spheroidal6.m"]
(* In Slepian's notation, Ω=πt/2, T=2. Hence, c=ΩT/2=πt/2 *)
(* Linear Prolate Functions are Angular Spheroidal Functions taking m=0 *)
T=2;
LinearProlate[n_,t_,z_]:=SpheroidalPS[n,0,z,π t/2]
(* Normalization factor for LPF.See Eq.(6.6) of this thesis *)
u2[n_,t_]:=Integrate[LinearProlate[n,t,z]^2,{z,-T/2,T/2}]
(* Normalized LPF assumed in Odlyzko's paper.See Eq.(6.7) of this thesis. *)
f[n_,t_,z_]:=LinearProlate[n,t,z]/Sqrt[u2[n,t]]
(* 1. By (24) (27) in Slepian's paper, the eigenvalue associated to f (defined above) is given by ((((2c )/π)[(Subscript[R, 0n]^((1)))(c,1)])^2), where (Subscript[R, 0n]^((1)))(c,1) is the Radial spheroidal function. Since c=ΩT/2=πt/2, we have ((2c )/π)=t. *)
(* SpheroidalS1 outputs an imaginary number with 0 imaginary part. Use Re to take out the real part. *)
γ[n_,t_]:=t Re[SpheroidalS1[n,0,1,π t/2]]^2

(* N0 is the truncation point when evaluate the infinite product and infinite sum in Eq.(6.8) of this thesis. The true value is when N0->∞. *)
N0=4; 
(* Define the spacing distribution function. *)
k=0; (* (k+1)st smallest eigenvalue of those that are larger than it *)
J=Subsets[Range[0,N0],{k+2}];
Ne=Length[J];
(*GUE predicted distribution.See Eq.(6.9) of this thesis*)
pk[t_]:=pk[t]=((16 )/(t^2)) Product[1-γ[m,t],{m,0,N0}]Sum[Product[γ[J[[v,s]],t]/(1-γ[J[[v,s]],t]) ,{s,1,k+2}]
Sum[f[J[[v,s1]],t,1]^2  f[J[[v,s2]],t,-1]^2 Boole[Mod[J[[v,s1]]+J[[v,s2]],2]==1] Boole[s1<s2], {s1,1,k+2},{s2,1,k+2}],
{v,1,Ne}];

(*Take the nearest neighbor distribution for example*)

(* Calculate the GUE predicted distribution  *)
pdf={};
For[t=0.05,t≤3,t=t+0.05,AppendTo[pdf,N[pk[t]]]];

(*Plot the GUE predicted distribution *)
xl=Table[1/20*i,{i,60}];
ltemp={};
For[i=1,i≤60,i++,AppendTo[ltemp,{xl[[i]],pdf[[i]]}]] 
p1=ListPlot[ltemp,PlotRange->{{0,3},{0,1.2}},Frame->True, Joined->True]

(* Empirical distribution *)
(* NZ: the number of data read from data set*)
NZ=10000;
(* This is the directory where I put the data file *)
SetDirectory["D:\\Data"];
data=ToExpression@ReadList["1e28.zeros.1000_10001001.txt",Record,NZ,RecordSeparators->" "];
(* Hiary's data set separates the integer part and the decimal part *)
(* Using Reap and Sow is faster than AppendTo *)
Zeros=Flatten[Reap[Do[Sow[N[(data[[i+1]]+data[[i]])]],{i,1,NZ-1,2}]][[2]]];

(*Normalize the spacing.See Eq.(6.12) of this thesis*)
height=10^28;
(Subscript[δ, n])={};
(Subscript[δ, n])=Flatten[Reap[Do[Sow[N[(Zeros[[i+1]]-Zeros[[i]]) Log[(height+Zeros[[i]])/(2 Pi)]/(2 Pi)]],{i,1,NZ/2-1,1}]][[2]]];

(* Draw the histogram of distribution *)
𝒟=HistogramDistribution[(Subscript[δ, n]),{Range[0.025,3.025,0.05]}];
p2=DiscretePlot[#[𝒟,x],{x,0.05,3,0.05},PlotLabel->#, PlotStyle->Red]&/@{PDF}

(* Compare in the same plot and add legends *)
Show[p1,p2,Epilog->Inset[Framed[Column[{PointLegend[{Red},{"    Empirical data"}],LineLegend[{Blue},{"GUE prediction"}]}],RoundingRadius->5],Scaled[{0.8,0.85}]]]