Suppose we have two random variables
$$R\sim U[1-\varepsilon,1]\;\;\;\;\; \Theta\sim U[0,2\pi],$$
and a third random variable
$$X=g(R,\Theta)=R\cos\Theta.$$
What is the density $p_X(x)$?
The section Transformations of Random Variables on Random contains a few results, the most general being the Distribution Function Method, since the Change of Variables Method does not seem to apply as $g$ is neither one-to-one (although can be made by restricting $\Theta$ and using symmetry) or strictly increasing/decreasing.
One idea is to introduce an auxiliary random variable $Y=R\sin\Theta$, define the bivariate transform $[X,Y] = [g_1(R,\Theta), g_2(R,\Theta)] = [R\cos\Theta, R\sin\Theta]$ and then solve
$$p_{X,Y}(x,y) = p_{R,\Theta}\left(\sqrt{x^2+y^2},\arctan\left(\frac{x}{y} \right) \right)|J_h(x,y)|$$
where $h(X,Y)$ is the inverse transform of $g$, and then marginalize over $Y$. Although I remain uncertain whether this formula is applicable for the transform $h$. Is there a more direct or alternative method to find the marginal density $p_X$?
I created a few plots of the marginal of $X$ by simulation to visualize what its distribution looks like for various values of $\varepsilon$.
Assuming $\delta=1-\epsilon\ge0$, for $x\ge0$ we have \begin{align} \def\pro#1{\textsf{Pr}\left(#1\right)} \pro{X\ge x}&=\pro{R\cos\Theta\ge x}\\ &=\pro{\cos\Theta\ge\frac xR}\\ &=\int_\delta^1\pro{\cos\Theta\ge\frac xr}\mathrm dr\\ &=\begin{cases} \frac1{\epsilon\pi}\int_\delta^1\arccos\frac xr\mathrm dr&0\le x\le\delta\\ \frac1{\epsilon\pi}\int_x^1\arccos\frac xr\mathrm dr&\delta\le x\le1\\ 0&x\ge1 \end{cases}\\ &=\begin{cases} \frac1{\epsilon\pi}\left(\arccos x-x\log\left(\sqrt{1-x^2}+1\right)-\delta\arccos\frac x\delta+x\log\left(\sqrt{\delta^2-x^2}+\delta\right)\right)&0\le x\le\delta\\ \frac1{\epsilon\pi}\left(\arccos x-x\log\left(\sqrt{1-x^2}+1\right)+x\log x\right)&\delta\le x\le1\\ 0&x\ge1\;. \end{cases} \end{align}
The density is the negative derivative of this probability with respect to $x$:
$$ p_X(x)=\begin{cases} \frac1{\epsilon\pi}\left(\log\left(\sqrt{1-x^2}+1\right)-\log\left(\sqrt{\delta^2-x^2}+\delta\right)\right)&0\le x\le\delta\\ \frac1{\epsilon\pi}\left(\log\left(\sqrt{1-x^2}+1\right)-\log x\right)&\delta\le x\le1\\ 0&x\ge1\;. \end{cases} $$
If you don't need the cumulative distribution function, you can avoid the detour of integrating the arccosine by applying the derivative with respect to $x$ directly to the integrals.
Here's a plot of the density for the four intermediate cases you simulated, $\epsilon=0.1,0.2,0.5,0.7$: