I need help finding the Fourier transform of the function
$$ \rho(\vec{r}) = \alpha \delta_{\vec{r},0} \left(\lambda\lambda' J_1 (\beta |\vec{r}|)Y_1(\beta |\vec{r}|) - \pi^2 J_0 (\beta |\vec{r}|)Y_0(\beta |\vec{r}|) \right), $$
where $\alpha,\beta \in \mathcal{R}$ and $\lambda,\lambda'=\pm1$. The $J_i(x)$ are the $i$-th Bessel functions of the first, the $Y_i(x)$ are the $i$-th Bessel functions of the second kind. $\delta_{\vec{r},0}$ denotes the Kronecker delta.
The vector $\vec{r} = (x,y)^T$ is discrete because I'm working on a discrete set of points. I'm a physicist, so please don't hesitate to ask for more specific information.
I tried something like
$$ \rho(\vec{k}) = \sum_\vec{r} e^{-i\vec{k}\cdot\vec{r}}\rho(\vec{r}), $$
which should in fact not be too difficult because the only vector that contributes to the sum is $\vec{r} = \vec{0}$ due to the Kronecker delta in $\rho(\vec{r})$. The problem is that $\rho(\vec{r})$ diverges at this point. I need some mathematical advice here. Is there anybody who can show me how to compute the Fourier transform of $\rho(\vec{r})$?
I'm almost certain the Fourier Transform doesn't exist.
I'll rewrite your radially symmetric (density?) function, $\rho(\vec{r})$, to be continuous, as opposed to discrete, in 2D polar coordinates:
$$d(r) = {}^2\delta(r) \alpha \left[\lambda \lambda' J_1(\beta r) Y_1(\beta r) - \pi^2 J_0(\beta r) Y_0(\beta r)\right]$$
Where ${}^2\delta(r)$ is the 2D Dirac Delta function at the origin of the polar coordinate system. (For the impulse at the origin, $\theta$ doesn't matter, so it doesn't appear as an argument.) In this analysis, ${}^2\delta(r)$ represents the following limiting sequence of functions (that are asymmetrical about 0):
$${}^2\delta(r) = \lim_{\epsilon \rightarrow 0^+} \dfrac{1}{\pi\epsilon^2} \quad 0 < r < \epsilon$$
so that
$$\int_0^{2\pi} \int_0^{\infty} {}^2\delta(r) \space r \space \mathrm{d}r \space \mathrm{d}\theta = 1$$
To separate the $r$ and $\theta$ portions of ${}^2\delta(r)$, so we can use a 1D Dirac Delta function, I'll use the following limiting sequence of functions (that are asymmetric about 0) for the 1D Dirac Delta function:
$$\delta_a(r) = \lim_{\epsilon \rightarrow 0^+} \dfrac{1}{\epsilon} \quad 0 < r < \epsilon$$
and then note $$\begin{align}\\ \int_0^{\infty} \delta_a(r) \space \mathrm{d}r &= 1\\ &= \int_0^{2\pi} \int_0^{\infty} {}^2\delta(r) \space r \space \mathrm{d}r \space \mathrm{d}\theta\\ &= \int_0^{\infty} {}^2\delta(r) \space 2\pi r \space \mathrm{d}r\\ \end{align}$$
so we have
$${}^2\delta(r) = \dfrac{\delta_a(r)}{2\pi r}$$
$d(r)$ can now be written without any dependence on $\theta$ as
$$d(r) = \dfrac{\delta_a(r)}{2\pi r} \alpha \left[\lambda \lambda' J_1(\beta r) Y_1(\beta r) - \pi^2 J_0(\beta r) Y_0(\beta r)\right]$$
The 2D Fourier Transform of a radially symmetric function can be expressed in terms of the Hankel Transform of order 0:
$$F(\rho) = \mathcal{F}\left\{f(r)\right\} = 2\pi \mathcal{H_0}\left\{f(r)\right\} = 2\pi \int_0^{\infty} f(r) J_0(\rho r) \space r \space \mathrm{d}r$$
When attempting to compute the Fourier Transform
$$\begin{align}\mathcal{F}\left\{d(r)\right\} = D(\rho) &= 2\pi \int_0^{\infty} \dfrac{\delta_a(r)}{2\pi r} \alpha \left[\lambda \lambda' J_1(\beta r) Y_1(\beta r) - \pi^2 J_0(\beta r) Y_0(\beta r)\right] J_0(\rho r) \space r \space \mathrm{d}r\\ \\ &= \alpha \int_0^{\infty} \delta_a(r) \left[\lambda \lambda' J_1(\beta r) Y_1(\beta r) - \pi^2 J_0(\beta r) Y_0(\beta r)\right] J_0(\rho r) \mathrm{d}r\\ \\ &= \alpha\lambda \lambda' \int_0^{\infty} \delta_a(r) J_1(\beta r) Y_1(\beta r) J_0(\rho r) \mathrm{d}r - \alpha \pi^2 \int_0^{\infty} \delta_a(r)J_0(\beta r) Y_0(\beta r) J_0(\rho r) \mathrm{d}r\\ \\ &= \alpha\lambda \lambda' \lim_{r \rightarrow 0+} J_1(\beta r) Y_1(\beta r) J_0(\rho r) - \alpha \pi^2 \lim_{r \rightarrow 0+} J_0(\beta r) Y_0(\beta r) J_0(\rho r) \\ \\ &= \alpha\lambda \lambda'\left(-\dfrac1{\pi}\right) - \alpha \pi^2 (-\infty) \\ \end{align}$$
it diverges.
Note that if the problematic term had converged, the Fourier Transform would simply be an infinite disk of constant height above the $\rho\phi$ plane. That should be no surprise for the transform of a Dirac Delta function.
The problematic term would converge, if there were a factor of $1/\ln(r)$:
$$\lim_{r \rightarrow 0+} \dfrac{J_0(\beta r) Y_0(\beta r)}{\ln(r)} J_0(\rho r) = \dfrac{2}{\pi}$$
I don't know the physical situation you are analyzing, and whether inclusion of such a term could ever be justified.