Flux integral of Gauss law Consider a point charge enclosed by some surface, using spherical coordinates, and taking $\hat a$ to be the unit vector in the direction of the surface element, flux is
$$\oint\vec E\cdot d\vec A   =  \frac q{4\pi}\oint\frac{\hat r}{r^2}\cdot  r^2d\phi \sin\theta d\theta\hat a \\ ‎\\=\frac q{4\pi}\oint\sin\theta d\theta d\phi (\hat r\cdot\hat a)$$
How to show that this integral is equal to $4\pi$ when the charge is inside and $0$ otherwise? I know that this can be written in terms of solid angle, but my acquaintance with that topic is very low, thus it would be appreciated if one refrains from using that to answer the question.
PS: I am working in HL system which for our purposes means that $\epsilon_0=1$ (though this doesn't affect the question at hand since it's about solving the integral, but is said just for the sake of completeness).
 A: PRIMER:
Note that in regions that are absent of charge, $\nabla \cdot E=0$.  Let $V$ be a region, bounded by the surface $S$, that contains a point charge.  Furthermore, let $V'$ be a region bounded by the surface $S'$, contained in $V$.
If the point charge is in $V'$, then it is outside of the region $V\setminus V'$ and from Gauss's law and the Divergence Theorem we have
$$\begin{align}
\int_{V\setminus V'}0\,dV&=\int_{V\setminus V'}\nabla\cdot\vec E\,dV\\\\
&=\oint_{S} \vec E\cdot \hat n\,dS-\oint_{S'} \vec E\cdot \hat n\,dS
\end{align}$$
from which we conclude
$$\oint_{S} \vec E\cdot \hat n\,dS=\oint_{S'} \vec E\cdot \hat n\,dS$$
That is, the value of the surface integral is independent of the surface over which the integration is taken.
If the point charge is in $V\setminus V'$, then $\nabla \cdot E=0$ in $V'$ and $\oint_{S'}\vec E\cdot \hat n\,dS=0$.
We now will show explicitly that if $S'$ is any sphere that contains a point charge, that $\oint_{S'}\vec E\cdot \hat n\,dS=q/\varepsilon_0$ (We could have even chosen a sphere for which the point charge is at the origin.).  To that end we proceed.

The electric field $\vec E(\vec r')$ due to a point charge at $\vec r$ at points on the surface of a sphere, centered at $0$ and with radius $a$ is
$$\vec E(\vec r')=\frac{q(a\hat r'-\vec r)}{4\pi\varepsilon_0|\vec r-a\hat r'|^3}$$
So, we wish to compute the integral
$$\oint_{|\vec r'|=a} \vec E(\vec r')\cdot \hat n'\,dS'=q\int_0^{2\pi}\int_0^\pi \left(\frac{a\hat r'-\vec r}{4\pi \epsilon_0|\vec r-a\hat r'|^3}\right)\cdot\hat r'\,a^2\,\sin(\theta')\,d\theta'\,d\phi'$$
Due to the spherical symmetry, without loss of generality, we align the polar axis so that the point charge lies at $\vec r=r\hat z$.  Then,
$$\begin{align}
\oint_{|\vec r'|=a} \vec E(\vec r')\cdot \hat n'\,dS'&=\frac{qa^2}{2\varepsilon_0}\int_0^\pi \frac{a-r\cos(\theta')}{(r^2+a^2-2ar\cos(\theta'))^{3/2}}\sin(\theta')\,d\theta'\\\\
&=\frac{qa^2}{2\varepsilon_0}\int_{-1}^1 \frac{a-rx}{(r^2+a^2-2arx)^{3/2}}\,dx\tag1
\end{align}$$

The integral on the right-hand side of $(1)$ is
$$\begin{align}\int_{-1}^1 \frac{a-rx}{(r^2+a^2-2arx)^{3/2}}\,dx&=-\left.\left(\frac{r-ax}{a^2(r^2+a^2-2arx)^{1/2}}\right)\right|_{-1}^1\\\\
&=\frac{r+a}{a^2(r+a)}-\frac{r-a}{a^2|r-a|}\\\\
&=\begin{cases}
0&, r>a\\\\
\frac 2{a^2}&,r<a\tag2
\end{cases}
\end{align}$$
Using $(2)$ in $(1)$ we find the coveted result
$$\oint_{|\vec r'|=a} \vec E(\vec r')\cdot \hat n'\,dS'=\begin{cases}\frac{q}{\varepsilon_0}&,r<a\\\\
0&,r>a\end{cases}$$
as was to be shown!
A: This answer is only helpful if you are familiar (and comfortable) with the divergence theorem and Dirac delta distribution.
By the divergence theorem we have
$$
\oint_{\partial V} \vec E\cdot d\vec A = \int_V \nabla\cdot\vec E\, dV.
$$
It is very useful to know that $\nabla\cdot\left(\frac{\hat r}{r^2}\right)= 4\pi \delta(\vec r)$. Here is a link. This gives the divergence of $\vec E$ of a charge at position $\vec r_0$ to be $\nabla \cdot \vec E= \frac{q}{4\pi}\times 4\pi\,\times \delta(\vec r- \vec r_0)$. So we have:
$$
\oint_{\partial V} \vec E\cdot d\vec A = q\int_V \delta(\vec r- \vec r_0)\, dV=\begin{cases}q&\vec r_0 \in V\\0&\vec r_0 \notin V\end{cases}.
$$
