Double Integration over finite plane. $$\phi(z)=\frac{\sigma}{4\pi\varepsilon_0}\int_{\frac{-a}{2}}^{\frac{a}{2}}\int_{\frac{-a}{2}}^{\frac{a}{2}}\frac{1}{\sqrt{x^2+y^2+z^2}}~dx~dy$$
I'm not sure how to do this integral. For the first integral w.r.t $x$ i tried to substitute $x=\sqrt{y^2+z^2}\sin{\theta}\implies dx=\sqrt{y^2+z^2}\cos\theta~d\theta$.
The integral then becomes:
$$\phi(z)=\frac{\sigma}{4\pi\varepsilon_0}\int_{\frac{-a}{2}}^{\frac{a}{2}}\int_{-\text{?}}^\text{?}1~d\theta~dy$$
But the bounds are $\text{?}=\arcsin{\frac{a}{2\sqrt{y^2+z^2}}}$ since arcsin is an odd function. However this just makes it even harder to solve. So what is the best way to do this integral?
If it helps I will give the context of the question. I am asked to find the strength of the electric field at a height z above the centre of a square sheet with constant charge density $\sigma$ and side lengths $a$.
 A: By symmetry, the electric field along the $z$ axis will have only a $z$ component with 
$$E_z(0,0,z)=\left.-\frac{\partial \phi(x,y,z)}{\partial z}\right|_{(0,0,z)}$$
Therefore, we have
$$E_z(0,0,z)=\frac{\sigma z}{\pi \epsilon_0}\int_0^{a/2}\int_0^{a/2}\frac{1}{(x^2+y^2+z^2)^{3/2}}\,dx\,dy$$
Rather than attempt a transformation to cylindrical coordinates, we proceed here with integrating directly in Cartesian coordinates.  
We evaluate the inner integral by making the substitution $x=\sqrt{y^2+z^2}\tan \theta$.  This yields
$$\begin{align}
\int_0^{a/2}\frac{1}{(x^2+y^2+z^2)^{3/2}}\,dx&=\left.\left(\frac{x}{(y^2+z^2)\sqrt{x^2+y^2+z^2}}\right)\right|_{0}^{a/2}\\\\
&=\frac{a/2}{(y^2+z^2)\sqrt{(a/2)^2+y^2+z^2}}
\end{align}$$
Therefore, we have reduced the expression for the electric field along the $z$ axis to 
$$E_z(0,0,z)=\frac{\sigma z}{\pi \epsilon_0}\int_0^{a/2}\frac{a/2}{(y^2+z^2)\sqrt{(a/2)^2+y^2+z^2}}\,dy$$
To evaluate the remaining integral, we make the standard trigonometric substitution $y=\sqrt{(a/2)^2+z^2}\tan(u)$.  Then, we have
$$\begin{align}
E_z(0,0,z)&=\frac{\sigma z(a/2)}{\pi \epsilon_0}\int_0^{\arctan\left(\frac{a/2}{\sqrt{(a/2)^2+z^2}}\right)}\,\,\frac{\cos(u)}{z^2+(a/2)^2\sin^2(u)}\,du\\\\
&=\frac{\sigma z(a/2)}{\pi \epsilon_0}\int_0^{(a/2)/\sqrt{(a/2)^2+(a/2)^2+z^2}} \frac{1}{z^2+(a/2)^2v^2}\,dv\\\\
&=\frac{\sigma }{\pi \epsilon_0}\arctan\left(\frac{(a/2)^2}{z\sqrt{2(a/2)^2+z^2}}\right)
\end{align}$$

NOTE 1:
We can recover, of course, the potential along the $z$ axis by integrating the electric field.  In this problem, integration by parts facilitates.  This is left as an exercise for the reader. 

NOTE 2:
As $a\to \infty$, the arctangent goes to $\pi/2\,\text{sgn}(z)$, and we recover the familiar result of the electric field from a uniform surface charge on an infinite surface, namely $\vec E(0,0,0^{\pm})=\pm \hat z\,\frac{\sigma}{2\epsilon_0}$.

NOTE 3:
As $z\to 0^{\pm}$, the arctangent goes to $\pi/2 \,\text{sgn}(z)$ and the electric field is $\vec E(0,0,0^{\pm})=\pm \hat z\,\frac{\sigma}{2\epsilon_0}$

NOTE 4:
As $z\to \pm \infty$, the arctangent goes to $\frac{(a/2)^2}{z^2}\,\text{sgn}(z)$ and the electric field is $\vec E(0,0,z\to \pm \infty)=\pm \hat z\,\frac{\sigma a^2}{4\pi \epsilon_0\,z^2}$, which appears as the field from a point charge $q=\sigma a^2$.
