Let $\Omega$ be the solid angle we seek. It is actually not that hard to compute $\Omega$ by evaluating your integral in polar coordinates (see method 1 below).
Method 1 - direct integration in polar coordinate.
Rotate the figure $180^\circ$ with respect to the center. By symmetry, $\Omega$ is equal to
$$\begin{align}
\int_0^{\frac{\pi}{3}}\int_0^{\frac{a}{2*\sqrt{3}\cos\theta}} \frac{6a rdrd\theta}{\sqrt{a^2+r^2}^3}
= \int_0^{\frac{\pi}{3}} \left[\frac{-6a}{\sqrt{a^2+r^2}}\right]_0^{\frac{a}{2*\sqrt{3}\cos\theta}} d\theta
= 2 \pi - 6\int_0^{\frac{\pi}{3}} \frac{d\theta}{\sqrt{1+\frac{1}{12\cos^2\theta}}}
\end{align}$$
We can evaluate the integral on RHS as
$$
\int_0^{\frac{\pi}{3}}\frac{\cos\theta}{\sqrt{\frac{1}{12} + \cos^2\theta}}
= \int_0^{\frac{\pi}{3}} \frac{d\sin\theta}{\sqrt{\frac{13}{12}-\sin^2\theta}}
= \sin^{-1}\left(\sqrt{\frac{12}{13}}\sin\frac{\pi}{3}\right)$$
As a result,
$$\Omega = 2\pi - 6\sin^{-1}\frac{3}{\sqrt{13}} \approx 0.3864229676956112$$
Method 2 - Oosterom and Strackee's formula for solid angle.
Alternatively, one can construct three vectors
$$
\vec{p} = \begin{bmatrix}\frac{a}{\sqrt{3}} \\ 0 \\ a \end{bmatrix},\quad
\vec{q} = \begin{bmatrix}-\frac{a}{2\sqrt{3}} \\ \frac{a}{2} \\ a \end{bmatrix}
\quad\text{ and }\quad
\vec{r} = \begin{bmatrix}-\frac{a}{2\sqrt{3}} \\ -\frac{a}{2} \\ a \end{bmatrix},
$$
normalize them and throw them
$\displaystyle\;\hat{p} = \frac{\vec{p}}{|\vec{p}|},
\hat{q} = \frac{\vec{q}}{|\vec{r}|},
\hat{r} = \frac{\vec{q}}{|\vec{r}|}$ to following formula of solid angle by
Oosterom and Strackee:
$$\tan\frac{\Omega}{2} = \frac{|\hat{p} \cdot ( \hat{q} \times \hat{r} )|}{1 + \hat{p}\cdot\hat{q} + \hat{q}\cdot\hat{r} + \hat{r}\cdot\hat{p}}$$
and get a different form of same answer.
$$\Omega = 2\tan^{-1}\left(\frac{\frac{9}{16}}{1 + 3\times \frac{5}{8}}\right)
= 2\tan^{-1}\frac{9}{46} \approx 0.3864229676956112$$
Method 3 - external angle formula for geodesic polygon.
Note - I added this because I may reuse the Lemma below somewhere else.
Feel free to ignore this.
Change coordinates so that the triangle $T$ at hand is living on the plane
$\mathcal{P} = \{ (x,y,z) : z = a \}$ and the reference point for computing
the solid angle is located at origin.
Consider following projection $\eta : \mathcal{P} \to S^2$:
$$\mathcal{P} \ni (x,y,a) \quad\mapsto\quad \frac{1}{\sqrt{x^2+y^2+a^2}}(x,y,a) \in \mathcal{S}^2$$
The $\Omega$ we want is simply the area of $\eta(T)$ on $S^2$. Under such a projection, straight lines on $\mathcal{P}$ get mapped to geodesics, i.e. arcs of great circles, on $S^2$. This implies $\eta(T)$ is a geodesic polygon.
Gauss-Bonnet theorem tells us
$$\Omega = \verb/Area/(\eta(T)) = 2\pi - \sum_i \theta_i
\quad\text{ where } \theta_i \text{ is the external angle at } i^{th} \text{ vertex}.$$
What we need to do is figure out the external angle for $\eta(T)$.
Parametrize $\mathcal{P}$ by polar coordinates and $S^2$ by spherical polar coordinates
$$\begin{array}{rcl}
[ 0, \infty) \times [ 0, 2\pi) \ni (r,\phi)
&\mapsto&
(x,y,z) = (r\cos\phi,r\sin\phi, a) \in \mathcal{P}\\
[ 0, \pi ) \times [ 0, 2\pi ) \ni (\theta, \phi )
&\mapsto&
(x,y,z) = (\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\phi) \in S^2
\end{array}
$$
Under projection $\eta$, the point with spherical polar coordinate $(\theta,\phi)$ on $S^2$ corresponds to the point with polar coordinate $(r,\phi)$ with $r = a\tan\theta$ on $\mathcal{P}$. Notice the metric on $S^2$ is given by
$$\begin{align}
d\theta^2 + \sin\theta^2 d\phi^2
&= \frac{(d\tan\theta)^2}{(1+\tan^2\theta)^2} + \frac{\tan^2\theta}{1+\tan^2\theta}d\phi^2
= \frac{a^2 dr^2}{(a^2+r^2)^2} + \frac{r^2 d\phi^2}{a^2+r^2}\\
&= \frac{1}{a^2+r^2}\left( \frac{dr^2}{1+(r/a)^2} + r^2d\phi^2\right)
\end{align}
$$
This means line segments parallel to and perpendicular to the radial directions
on $\mathcal{P}$ are scaled differently. Using this info, it is easy to derive
following
Lemma Let $\gamma$ be a curve in $\mathcal{P}$ passing through $p = (r\cos\phi,r\sin\phi,a )$ and making an angle $\alpha$ with the radial direction. If $\beta$ is the angle between $\eta(\gamma)$ and the meridian passing through $\eta(p)$, then
$$\tan\beta = \sqrt{1+\frac{r^2}{a^2}} \tan\alpha$$
Apply this lemma to our problem,
- At any vertex of $T$, the side of $T$ is making an angle $\frac{\pi}{6}$ with the radial direction.
- This implies its image under $\eta$ is making an angle $\tan^{-1}\left(\sqrt{1 + \frac13}\tan\frac{\pi}{6}\right) = \tan^{-1}\frac{2}{3}$ with the corresponding meridian.
- As a result, the external angle for $\eta(T)$ is equal to
$\pi - 2\tan^{-1}\frac23$.
From this, we find
$$\Omega = 2\pi - 3(\pi - 2\tan^{-1}\frac23) = 6\tan^{-1}\frac23 - \pi
\approx 0.3864229676956121$$
Yet another form of the same answer.
