(***** This is an extremely simple script to show how any state can be approximated by smearing boundary operators over a finite time-band. Here the extent of the time-band is taken to be [0,1]
in units where the light-crossing time of AdS is π. This can be changed by changing the "window" function below. The window function used below has the advantage that its Fourier transform dies off very rapidly for large positive and negative frequencies. On the other hand, it has the disadvantage that the Fourier transform is not known analytically.  For this reason, with the window function below, it is difficult to increase N_{cut} or nostates to much larger values than 100. This also means that the accuracy of the reconstruction is limited to about 1%.  We invite users of this script to modify the window function, and inform us if they obtain interesting results and can increase the numerical accuracy significanty. *****)

window[x_] := If[(x > 0 && x < 1), Exp[1/(x (x - 1))], 0]
windowf[n_?NumericQ] := windowf[n] = NIntegrate[window[x] Exp[I n x], {x, 0, 1}, MaxRecursion -> 20]
state[Δ_,  l_, d_, Ncut_,m_] := Table[windowf[2 n + Δ + l + m]*Sqrt[G[n,l,d,Δ]], {n, 0, Ncut}](*computes the state created by \int O(t) e^{i m t} window[t] d t, where O has an expansion as O(t) = \sum_n O_n e^{i (2 n + l + \Delta) t} + h.c *)
G[n_, l_,d_,  Δ_] := G[n,l,d,Δ] = 4 π^d Gamma[Δ + n + l] Gamma[Δ + n + 1 - d/2]/(Gamma[d/2] Gamma[n + 1] Gamma[Δ] Gamma[Δ + 1 - d/2] Gamma[d/2 + n + l])
eneigstate[n_, Ncut_] := Module[{temp,i}, temp=Table[0, {i,0, Ncut}]; temp[[n+1]] = 1; temp];(*This simply puts a "1" in the (n+1)th place. The idea is to compare states with a fixed angular momentum, so that they are then labelled by n.  *)

Plot[window[x], {x, 0, π}, PlotRange -> All]

Plot[Abs[windowf[k]], {k, -150, 150}, PlotRange -> All]
(***** This step may take some time, but is worth running since it initializes the windowf function apart from producing the graph ********)

(******** defines an error function ************)
errfunc[dim_, ang_, Δ_,en_,  Ncutoff_, nostates_] := Module[{rsbasistates,i, var, trialstate}, rsbasistates = Table[state[Δ, ang, dim, Ncutoff, i]/Sqrt[Norm[state[Δ, ang, dim, Ncutoff, i]]], {i, -nostates, nostates}]; var = LeastSquares[Transpose[rsbasistates], eneigstate[en, Ncutoff]];trialstate = Sum[var[[i]] rsbasistates[[i]], {i, 1, Length[rsbasistates]}];
Norm[trialstate - eneigstate[en, Ncutoff]]]

d = 4;
Δ = 4;
En=3;
ang = 2;
(******** The state we are trying to approximate is |E=2*En+ang+\Delta, l=ang>. You can change these values. *********)

(****** This produces an estimate of the error ******)
errtable =ParallelTable[errfunc[d, ang, Δ, En, 100, nostates],{nostates, 1, 30}]

ListPlot[errtable]