where, $x=r\cos^2\theta , y=r\sin\theta \cos\theta$
I drew the region on the $xy$-plane over which the integral is performed
After the transformation to the $R\Theta$-plane the region changes.
To me it seems the transformed region in the $R\Theta$-plane is unbounded! But I am not sure weather or not I am correct.
Please help me find the transformed region on the $R\Theta$-plane and the Transformed integral itself! So, I could calculate.
The integral can be split into two, each of which being solved using coordinates change differently.
Split the integral as $$ \int_{0}^{2a}\int_{\sqrt{2ax-x^2}}^{\sqrt{4ax-x^2}}(1+\frac{y^2}{x^2})\ dy\ dx {= \underbrace{\int_{0}^{2a}\int_{0}^{\sqrt{4ax-x^2}}(1+\frac{y^2}{x^2})\ dy\ dx}_{I_1} \\-\underbrace{\int_{0}^{2a}\int_0^{\sqrt{2ax-x^2}}(1+\frac{y^2}{x^2})\ dy\ dx}_{I_2}. } $$ For $I_1$, through a change of variables $x\to 2a+x$ we can write $$ I_1 {=\int_{-2a}^{0}\int_{0}^{\sqrt{4a^2-x^2}}(1+\frac{y^2}{(x+2a)^2})\ dy\ dx \\= \int_{\frac{\pi}{2}}^{\pi}\int_{0}^{2a}\left[1+\frac{r^2\sin^2\theta}{(r\cos \theta+2a)^2}\right]rdrd\theta. } $$ Similarly, for $I_2$, through a change of variables $x\to a+x$ we have $$ I_2{= \int_{-a}^{a}\int_0^{\sqrt{a^2-x^2}}(1+\frac{y^2}{(a+x)^2})\ dy\ dx \\= \int_{0}^{\pi}\int_0^{a}\left[1+\frac{r^2\sin^2\theta}{(a+r\cos \theta)^2}\right]rdrd\theta. } $$ The rest of calculations should follow simply, however, I did not go through the details from here on. Good luck!