Comparing analytical and numerical results of 2D Fourier transform

Question

I would like to be able to transform back and forth with Fourier and InverseFourier, however, I can't quite match up the results. A similar question exists for 1D, but I don't see how to the necessary normalising coefficients there are obtained. A related question exists for 2D FFTs, from which I have adapted RotateRight.

I'm trying to compare the FFT of a 2D disc $x^2+y^2\leq a^2$ with the well known analytical result. Defining the circle function as proposed here:

$Assumptions = {Element[{x, y}, Reals], a > 0}
PWE = If[x^2 + y^2 <= a^2, 1, 0] // PiecewiseExpand 
FT = FourierTransform[PWE, {x, y}, {u, v}, FourierParameters -> {1, -1}]

This gives

and the inverse transform

    iFT = InverseFourierTransform[FT, {u, v}, {x, y}, 
   FourierParameters -> {1, -1}] // PiecewiseExpand

Works fine.

Next, I define the disc numerically in realspace in the list Da, and then apply Fourier. I am however certainly doing something wrong in the rotation and scaling of the result, which is encapsulated in the line FFTscaled = RotateRight[FFT, Dimensions[FFT]/2] (L/n)*(L/n) . As we will see in a second, this does not match the analytical result. Here is the full FFT code:

    AbsoluteTiming[
 (* https://mathematica.stackexchange.com/questions/280018/is-there-a-\
built-in-equivalent-to-numpys-fftfreq *)
 fftfreq[n_, d_ : 1] := Mod[Range@n - 1, n, -n/2]/(d n);
 
 (*Define parameters*)
 L = 30;
 a = 1.5;
 n = 200;
 
 
 xReal = Subdivide[-L/2, L/2, n - 1] // N;
 xFreq = fftfreq[n, L/n] 2 Pi // N;
 
 Da = Table[
   If[Sqrt[xReal[[k]]^2 + xReal[[j]]^2] < a, 1.0, 0.0], {k, 1, 
    n}, {j, 1, n}];
 ]

AbsoluteTiming[FFT = Fourier[Da, FourierParameters -> {1, -1}];]
AbsoluteTiming[
 FFTscaled = RotateRight[FFT, Dimensions[FFT]/2] (L/n)*(L/n);]

{ArrayPlot[Da, PlotLabel -> "original"], 
 ArrayPlot[FFTscaled, PlotLabel -> "FFT (rotated, scaled)"]}

The function PrepareForPlot simply transforms the data in a way that is accepted by ListPlot3D and Interpolation:

PrepareForPlot[data_, domain_] := 
     Flatten[Table[{domain[[k]], domain[[j]], data[[k, j]]}, {k, 1, 
        n}, {j, 1, n}], 1]

Now, if we compare a slice of the numerical data with the analytical result, we see some discrepancies. Note that I scaled the frequency data earlier with xFreq = fftfreq[n, L/n] 2 Pi // N and here additionally I rotate it before the plot with RotateRight[xFreq, n/2]:

range = 10;
interp = Interpolation[
   PrepareForPlot[Abs@FFTscaled, RotateRight[xFreq, n/2]]];
Plot[{interp[u, v], (2 a \[Pi] BesselJ[1, a Sqrt[u^2 + v^2]])/Sqrt[
    u^2 + v^2]} /. u -> 0 // Evaluate, {v, -range, range}, 
 PlotRange -> All]

Clearly, there are some discrepancies.

EDIT: Re@FFTscaled produces:

EDIT 2: I add here another code which you can just copy and run. This is minimally altered, should work for general functions. However, as the plot shows at the end, there is a problem. There must be something wrong with my rescaling of the frequency domain...

func[x_, y_] := E^(-Pi (a^2 x^2 + b^2 y^2))
(* result consistent with \
https://adriftjustoffthecoast.wordpress.com/2013/06/06/2d-fourier-\
transform-of-the-unit-disk/ *)
Clear[a, b]
FT[u_, v_] := 
 Evaluate[FourierTransform[func[x, y] // Evaluate, {x, y}, {u, v}, 
    FourierParameters -> {1, -1}] // FullSimplify]
func[x, y] == 
  InverseFourierTransform[FT[u, v], {u, v}, {x, y}, 
   FourierParameters -> {1, -1}] // FullSimplify

AbsoluteTiming[
 (* https://mathematica.stackexchange.com/questions/280018/is-there-a-\
built-in-equivalent-to-numpys-fftfreq *)
 fftfreq[n_, d_ : 1] := Mod[Range@n - 1, n, -n/2]/(d n);
 
 (*Define parameters*)
 L = 30;
 a = 0.5;
 b = 1;
 n = 200;
 
 
 xReal = Subdivide[-L/2, L/2, n - 1] // N;
 xFreq = fftfreq[n, L/n] 2 Pi // N;
 
 Da = Table[FT[xReal[[k]], xReal[[j]]], {k, 1, n}, {j, 1, n}];
 ]

AbsoluteTiming[FFT = Fourier[Da, FourierParameters -> {1, -1}];]
AbsoluteTiming[
 iFFT = InverseFourier[FFT, FourierParameters -> {1, -1}];]
AbsoluteTiming[
 FFTscaled = RotateRight[FFT, Dimensions[FFT]/2] (L/n)*(L/n);]

{ArrayPlot[Da, PlotLabel -> "original"], 
 ArrayPlot[FFTscaled, PlotLabel -> "FFT (rotated, scaled)"], 
 ArrayPlot[iFFT, PlotLabel -> "iFFT"] }

PrepareForPlot[data_, domain_] := 
 Flatten[Table[{domain[[k]], domain[[j]], data[[k, j]]}, {k, 1, 
    n}, {j, 1, n}], 1]

range = 10;
interp = Interpolation[
   PrepareForPlot[Re@FFTscaled, RotateRight[xFreq, n/2]]];
Plot[{interp[u, v], FT[u, v]} /. u -> 0 // Evaluate, {v, -range, 
  range}, PlotRange -> All, PlotLegends -> {"FFT", "analytical"}]