2D wave front set: simple planar shapes (in progress)

Let \(h>0\) be the semiclassical parameter.
The FBI (Fourier-Bros-Iagolnitzer) transform of an \(n\)-dimensional tempered distribution \(u \in \mathscr{S}^\prime(\mathbb{R}^n)\) is defined by
\[
T_hu(x;\xi)
:=
2^{-n/2}
(\pi h)^{-3n/4}
\int_{\mathbb{R}^n}
e^{-i(x-y)\cdot\xi/h-\lvert{x-y}\rvert^2/2h}
u(y)
dy,
\]
where \((x,\xi) \in T^\ast\mathbb{R}^n\simeq\mathbb{R}^n\times\mathbb{R}^n\).
It is well-known that if \(u\) is independent of \(h>0\),
then the wave front set of \(u\) is characterized by its FBI transform.
In fact it follows that \((x_0,\xi_0) \in T^\ast\mathbb{R}^n\setminus0\)
does not belong to \(\operatorname{WF}(u)\) if and only if
\[
T_hu(x;\xi)=\mathcal{O}(h^\infty)
\quad
(h \downarrow0)
\]
uniformly near \((x_0,\xi_0)\). Note that a wave front set is conic in \(\xi\) and independent of the size of \(\xi\). We are trying to detect the wave front sets of 2D distributions using MATLAB and the FBI transform. For \(u(x,y) \in \mathscr{S}^\prime(\mathbb{R}^2)\) we compute
\[
T_hu(x,y;\lambda\cos\theta,\lambda\sin\theta)
=
2^{-1}(\pi h)^{-3/2}
\iint_{\mathbb{R}^2}
e^{i\lambda(p\cos\theta+q\sin\theta)/h-(p^2+q^2)/2h}
u(x+p,y+q)
dpdq
\]
where \(\lambda\) is the size of the frequency.