Computation of the probability density function for $(X,Y) = \sqrt{2 R} ( \cos(\theta), \sin(\theta))$ Let $R$ be a almost surely non-negative continuous random variable with absolutely continuous measure, and $\Theta$ be an independent random variable, uniformly distributed on the interval $[0, 2 \pi)$. I am interested in random vector $(X, Y)$ such that
$$
  X = \sqrt{2 R} \cos(\Theta) \qquad Y = \sqrt{2 R} \sin(\Theta)
$$
In the special case when $R$ is exponential random variable with unit mean, $(X,Y)$ is a pair of independent standard normal random variables. 
I would like to compute the joint probability density function $f_{X,Y}(x,y)$ expressed as an expectation over $R$.
I am able to derive the laws for marginals:
$$
    f_X(x) = \mathbb{E}\left( \frac{[ R > \frac{x^2}{2} ]}{ \pi \sqrt{ 2R - x^2}} \right) \qquad
    f_Y(y) = \mathbb{E}\left( \frac{[ R > \frac{y^2}{2} ]}{ \pi \sqrt{ 2R - y^2}} \right)
$$
but I do not yet see through to finding the joint pdf $f_{X,Y}(x,y)$. 
Any help, hints at derivation, or references is appreciated. Thank you.
 A: As usual, fix a bounded measurable function $\varphi$ and consider
$$
(*)=\mathrm E(\varphi(X,Y))=\mathrm E(\varphi(\sqrt{2R}\cos\Theta,\sqrt{2R}\sin\Theta)),
$$
hence
$$
(*)=\iint[r\gt0,0\lt\theta\lt2\pi]\,\varphi(\sqrt{2r}\cos\theta,\sqrt{2r}\sin\theta)f_R(r)\mathrm dr\frac{\mathrm d\theta}{2\pi}.
$$
The change of variables $x=\sqrt{2r}\cos\theta$ and $y=\sqrt{2r}\sin\theta$, yields $2r=x^2+y^2$ and the Jacobian $\mathrm dr\mathrm d\theta=\mathrm dx\mathrm dy$, hence
$$
(*)=\frac1{2\pi}\iint \varphi(x,y)f_R\left(\frac{x^2+y^2}2\right)\mathrm dx\mathrm dy.
$$
This proves that
$$
\color{red}{f_{X,Y}(x,y)=\frac1{2\pi}f_R\left(\frac{x^2+y^2}2\right)}.
$$
One can recover the marginal densities $f_X=f_Y$ from here by the formula
$$
f_X(x)=f_Y(x)=\int f_{X,Y}(x,y)\mathrm dy=\frac1{2\pi}\int f_R\left(\frac{x^2+y^2}2\right)\mathrm dy,
$$
that is,
$$
f_X(x)=f_Y(x)=\frac1{\pi}\int_{y\gt0} f_R\left(\frac{x^2+y^2}2\right)\mathrm dy.
$$
The change of variable $2r=x^2+y^2$ yields $r\gt\frac12x^2$ and $\mathrm dr=\sqrt{2r-x^2}\mathrm dy$, hence
$$
\color{red}{f_X(x)=f_Y(x)=\frac1{\pi}\int_{x^2/2}^{+\infty}\frac{f_R(r)}{\sqrt{2r-x^2}}\mathrm dr}.
$$
A: You basically defined a spherically symmetric distribution. A spherically symmetric RV can be represented as the product of two independent RV $R$ and $\mathbf U$ where $R$ is almost surely non-negative and $\mathbf U$ is uniform on the unit sphere. After transforming your radius $R$ into $R':=\sqrt{2R}$ you can use the formula for the $n$-dimensional case
$$p(\mathbf x) = \frac{\varrho{(\|\mathbf x\|)}}{\|\mathbf x\|^{n-1} S},$$
where $S=\frac{2\Gamma(1/2)^n}{\Gamma(n/2)}$ is the surface area of the unit sphere in $\mathbb R^n$ and $\varrho$ is the univariate density for $R'$.
There are probably many references for it, but the one I used is
Gupta, A. K., & Song, D. (1997). Lp-norm spherical distribution. Journal of Statistical Planning and Inference, 60(2), 241-260. Elsevier.
which handle a slightly more general case. 
