(* Special Schwartz function.  The Gaussian is its own Fourier transform and its integral is 1. *)
φ[x_]=Exp[-Pi x^2]; 

(* Real part of bilinear character. *)
ψ[x_,ξ_]=Cos[2Pi x ξ];

(* Dampened Gaussian; the t here is (Sqrt[t]) from the handout. *)
Clear[t];
g[t_,ξ_]=Exp[-Pi t^2 ξ^2];

(* Its Fourier transform. *)
ϕ[t_,ξ_]=g[1/t,ξ]/t;

FirstIntegrand[x_,ξ_]=φ[x]ψ[x,ξ];

SecondIntegrand[x_,ξ_,t_]=φ[x]g[t,ξ]ψ[x,ξ];

xmax=0.8; ximax=5;

(* Real part of naive integrand for the double integral.  The x axis runs along the front of the figure, the ξ axis along the side.

Note that the front cross section φ(x)cos(2(Subscript[πξ, aax])x)is dampened while the side cross section cos(2(Subscript[πx, aax])ξ) is not. 

For fixed ξ the cross sectional x integral should be the Gaussian of ξ and is a classical integral.  For fixed x the cross sectional ξ integral should be ∞ if x is 0 and 0 otherwise but this is a distribution.
*)

Plot3D[FirstIntegrand[x,ξ],{x,-xmax,xmax},{ξ,-ximax,ximax},Mesh->None,PlotPoints->50,PlotRange->All]

(* Absolute value of naive integrand for the double integral.  This figure shows that the naive double integral does not converge absolutely.  It needs to be dampened in the ξ direction. 
*)

Plot3D[φ[x],{x,-xmax,xmax},{ξ,-ximax,ximax},Mesh->None,PlotPoints->50,PlotRange->All]

(* Successively less ξ dampened integrands to make the double integral absolutely convergent.  Also the two outer cross sectional integrands in each case.  The first cross sectional integrand is the function of ξ obtained by integrating in the x direction. We expect it to converge monotonically up to the original Schwartz function φ (also shown in its graphs) as t shrinks to 0.  The second cross sectional integrand is the function of x obtained by integrating in the ξ drection.  We expect it to converge to a spike as t shrinks to 0. 
*)

t=1/2;
Plot3D[SecondIntegrand[x,ξ,t],{x,-xmax,xmax},{ξ,-ximax,ximax},Mesh->None,PlotPoints->50,PlotRange->All]
Plot[{φ[ξ],φ[ξ]g[t,ξ]},{ξ,-ximax,ximax},PlotRange->All]
Plot[φ[x]ϕ[t,x],{x,-xmax,xmax},PlotRange->All]

t=1/5;
Plot3D[SecondIntegrand[x,ξ,t],{x,-xmax,xmax},{ξ,-ximax,ximax},Mesh->None,PlotPoints->50,PlotRange->All]
Plot[{φ[ξ],φ[ξ]g[t,ξ]},{ξ,-ximax,ximax},PlotRange->All]
Plot[φ[x]ϕ[t,x],{x,-xmax,xmax},PlotRange->All]

t=1/10;
Plot3D[SecondIntegrand[x,ξ,t],{x,-xmax,xmax},{ξ,-ximax,ximax},Mesh->None,PlotPoints->50,PlotRange->All]
Plot[{φ[ξ],φ[ξ]g[t,ξ]},{ξ,-ximax,ximax},PlotRange->All]
Plot[φ[x]ϕ[t,x],{x,-xmax,xmax},PlotRange->All]

t=1/25;
Plot3D[SecondIntegrand[x,ξ,t],{x,-xmax,xmax},{ξ,-ximax,ximax},Mesh->None,PlotPoints->50,PlotRange->All]
Plot[{φ[ξ],φ[ξ]g[t,ξ]},{ξ,-ximax,ximax},PlotRange->All]
Plot[φ[x]ϕ[t,x],{x,-xmax,xmax},PlotRange->All]