Determining the bounds of integration after switching to polar coordinates the integral I am solving is:
$$I=\iint_Sxy\sqrt{1+x^2+y^2}\,dx\,dy$$
where $S$ is a rectangle $0\le x\le1, 0\le y\le1$
after $$x=\rho\cos{\varphi},y=\rho\sin{\varphi}$$
$$dy\,dx=\rho \,d\rho \,d\varphi$$
Where i struggle is finding the bounds for integration after the switch. Drawing the area, i see, for $\varphi\in (0,\pi/4)$ $r$ goes from $0$ to $1/\cos{\varphi}$ and for $\varphi\in(\pi/4,\pi/2)$ $r$ goes from $0$ to $1/\sin{\varphi}$. If i divided the integral into two (four), i obtain:
$$I=\int_0^{\pi/4}\frac{1}{2} \sin{2\varphi} \int_0^{1/\cos{\varphi}} \rho^3\sqrt{1+\rho^2} \, d\rho +\int_{\pi/4}^{\pi/2}\frac{1}{2} \sin{2\varphi} \int_0^{1/\sin{\varphi}} \rho^3\sqrt{1+\rho^2} \, d\rho$$
Is this approach correct or I am doing something terribly wrong?
 A: I can't speak to the computation since I did not check, however integrating over a rectangle using polar coordinates is often a bad idea (even if you see a $x^2+y^2$ in the integrand!). 
This integral in particular is not too bad in cartesian,
$$
\int_0^1\int_0^1 xy\sqrt{1+x^2+y^2}\;\mathrm dy \mathrm dx\\
\stackrel{u=1+x^2+y^2\implies \mathrm du=2y\mathrm dy}{=}
\frac12\int_0^1x\int_1^{2+x^2} \sqrt{u}\mathrm \;du \mathrm dx\\
=\frac13\int_0^1x\left((2+x^2)^{3/2}-1 \right)\mathrm dx\\
=\frac13\int_0^1x(2+x^2)^{3/2}\;\mathrm dx-\frac{1}{3} \int_0^1x\;\mathrm dx 
$$
both of which are easy to evaluate.
A: The integral can be computed by exploiting symmetry, rather than putting a square peg into a round hole (I loved Matthew Leingang's comments above).
$$\begin{eqnarray*} \iint_{(0,1)^2}xy\sqrt{1+x^2+y^2}\,dx \,dy &=& 2 \iint_{0\leq y\leq x\leq 1}xy\sqrt{1+x^2+y^2}\,dx\,dy\\&=&2\iint_{(0,1)^2} x^3 z\sqrt{1+x^2+x^2 z^2}\,dz\,dx\\&=&\frac{2}{3}\int_{0}^{1}x\left((2x^2+1)^{3/2}-(x^2+1)^{3/2}\right)\,dx\\&=&\frac{1}{15}\left[(2x^2+1)^{5/2}-2(x^2+1)^{5/2}\right]_{0}^{1}\\[0.2cm]&=&\color{blue}{\frac{9\sqrt{3}-8\sqrt{2}+1}{15}}.\end{eqnarray*}$$
