Compute the area of a surface, encountering a strange integral

compute the area of this yellow surface, which is actually a paraboloid:
$$x^2+y^2=2az$$(yellow one), cutted by $$(x^2+y^2)^2=2a^2xy$$(blue one)
To compute the part be surrounded in the blue surface, use polor coordinates: $$\begin{cases} x=r\cos\theta\\ y=r\sin\theta\\ z=\frac{r^2}{2a}\\ \end{cases}$$ use Gauss efficient to compute the area: $$\begin{cases} E=(\frac{\partial x}{\partial r})^2+(\frac{\partial y}{\partial r})^2+(\frac{\partial z}{\partial r})^2&=1+\frac{r}{a}\\ F=\frac{\partial x}{\partial r}\frac{\partial x}{\partial\theta}+\frac{\partial y}{\partial r}\frac{\partial y}{\partial\theta}+\frac{\partial z}{\partial r}\frac{\partial z}{\partial\theta}&=0\\ G=(\frac{\partial x}{\partial\theta})^2+(\frac{\partial y}{\partial\theta})^2+(\frac{\partial z}{\partial\theta})^2&=r^2\\ \end{cases}$$ $$S=\iint\limits_D\sqrt{EG-F^2}drd\theta$$ where D becomes $$\left\{ (r,\theta )|\theta \in \left[ 0,\frac{\pi}{2} \right] \cup \left[ \pi ,\frac{3}{2}\pi \right] ,\mathrm{r}\in \left[ 0,\mathrm{a}\sqrt{\sin 2\theta} \right] \right\}$$ since the blue surface can be written as $$r^2=a^2\sin2\theta$$.
I can't find anything wrong until now, but this intergal is extreme complex, and Mathematica give me a non elementary solution. However, the standard solution given by the text book to this question is$$\frac{20-3 \pi}{9} a^{2}$$
How can i get that answer? Or what's wrong with my method?

By homogeneity we may assume $$a=1$$. Let $$z=f(x,y)=\frac{x^2+y^2}{2}$$, then (see surface integral) \begin{align} A_1&=\iint_D\sqrt{1+f_x^2+f_y^2}\,dxdy=\iint_D\sqrt{1+x^2+y^2}\,dxdy\\ &=2\int_{\theta=0}^{\pi/2} \int_{r=0}^{\sqrt{\sin(2\theta)}}\sqrt{1+r^2}\cdot r\,drd\theta\\ &=\frac{2}{3}\int_{\theta=0}^{\pi/2} ((1+\sin(2\theta))^{3/2}-1)\,d\theta\\ &=\frac{2}{3}\int_{\theta=0}^{\pi/2} (\cos(\theta)+\sin(\theta))^{3}\,d\theta-\frac{\pi}{3}\\ &=\frac{20}{9}-\frac{\pi}{3}. \end{align} For the last integral you may expand the cube and use the standard reduction formula. Finally $$A_a=A_1\cdot a^2=\left(\frac{20}{9}-\frac{\pi}{3}\right)a^2.$$
The same follows from $$\begin{cases} E=(\frac{\partial x}{\partial r})^2+(\frac{\partial y}{\partial r})^2+(\frac{\partial z}{\partial r})^2&=1+\color{blue}{\frac{r^2}{a^2}}\qquad \text{(a square is missing)}\\ F=\frac{\partial x}{\partial r}\frac{\partial x}{\partial\theta}+\frac{\partial y}{\partial r}\frac{\partial y}{\partial\theta}+\frac{\partial z}{\partial r}\frac{\partial z}{\partial\theta}&=0\\ G=(\frac{\partial x}{\partial\theta})^2+(\frac{\partial y}{\partial\theta})^2+(\frac{\partial z}{\partial\theta})^2&=r^2 \end{cases}$$ and therefore $$\sqrt{EG-F^2}=\sqrt{1+\frac{r^2}{a^2}}\cdot r$$.
Notice that you missed a square also in the definition of the blue surface: it should be $$(x^2+y^2)^2=2a^2xy$$ so that in cylindrical coordinates we find $$r^2=a^2\sin(2\theta)$$.
I am coming from assumption, that to obtain desired answer you should have cylinder $$(x^2+y^2)^2=2a^2xy$$ in place of brought in title $$x^2+y^2=2a^2xy$$ (blue one).
Then we will have $$\sqrt{1+z'_x{^2}+z'_y{^2}} = \frac{1}{a}\sqrt{a^2+x^2+y^2}$$. After writing double integral it is enough to use usual polar coordinates, so we will have $$S=\frac{4}{a}\int\limits_{0}^{\frac{\pi}{4}}\,d\phi \int\limits_{0}^{a\sqrt{\sin 2\phi}}r\sqrt{a^2+r^2}\,dr = \frac{4}{3} a^2 \left( \int\limits_{0}^{\frac{\pi}{4}}(\sin \phi +\cos \phi)^3\,d \phi - \frac{\pi}{4} \right ) =\\ =\frac{4}{3} a^2 \left( 2\sqrt{2} \int\limits_{0}^{\frac{\pi}{4}} \sin^3 \left( \phi+\frac{\pi}{4} \right )\,d\phi - \frac{\pi}{4}\right ) =\\ =\frac{4}{3} a^2 \left(2\sqrt{2} \int\limits_{0}^{\frac{\pi}{4}} \left( \cos^2 \left( \phi+\frac{\pi}{4} \right )-1 \right )\,d \left( \cos \left( \phi+\frac{\pi}{4} \right )-1 \right ) - \frac{\pi}{4}\right ) = \frac{a^2}{9}(20-3\pi)$$