Let the ellipse be in standard position, with (fraction-free) equation
$$b^2 x^2 + a^2 y^2 = a^2 b^2$$
Let our equilateral triangle have circumcenter $(p,q)$ and circumradius $r$. Note that maximizing the area of the triangle is equivalent to maximizing $r$.
For some angle $\theta$ ---actually, for three choices of $\theta$--- the vertices of the triangle have coordinates
$$
(p,q) + r \; \mathrm{cis}\theta \qquad (p,q)+r\;\mathrm{cis}\left(\theta+120^{\circ}\right) \qquad (p,q) + r\;\mathrm{cis}\left(\theta-120^{\circ}\right)
$$
where I abuse the notation "$\mathrm{cis}\theta$" to indicate the vector $(\cos\theta,\sin\theta)$.
Substituting these coordinates into the ellipse equation gives a system of three equations in four parameters $p$, $q$, $r$, $\theta$. I used Mathematica's Resultant[] function to help me eliminate $r$ and $\theta$, arriving at a huge polynomial equation in $p$ and $q$. One factor of the polynomial gives rise to this equation:
$$p^2 b^2\left( a^2+3b^2 \right)^2 + q^2 a^2\left(3a^2+b^2\right)^2= a^2b^2\left(a^2-b^2\right)^2$$
This says that the family of circumcenters $(p,q)$ lie on their own ellipse! We can therefore write
$$p = \frac{a\left(a^2-b^2\right)}{a^2+3b^2}\cos\phi \qquad q = \frac{b\left(a^2-b^2\right)}{3a^2+b^2}\sin\phi$$
for some $\phi$. Back-substituting into the system of equations gives this formula for $r$:
$$r = \frac{4 a b \sqrt{a^2\left(a^2+3b^2\right)^2-\left(a-b\right)^3\left(a+b\right)^3 \cos^2\phi}}{\left(a^2 + 3 b^2\right)\left(3a^2+b^2\right)}$$
The maximum value, $R$, is attained when $\cos\phi = 0$, so
$$R := \frac{4a^2b}{3a^2+b^2}$$
The area of the maximal triangle is
$$\frac{3\sqrt{3}}{4}R^2 = \frac{12a^4 b^2\sqrt{3}}{\left(3a^2+b^2\right)^2}$$
Perhaps-unsurprisingly, the corresponding triangles are centered horizontally within the ellipse, with a vertex at either the top or bottom of the minor axis.
Now, I should point out that my big $pq$ polynomial has other factors, namely, $p$ itself, $q$ itself, and a giant I'll call $f$.
One can verify that the cases $p=0$ and $q=0$ lead to the same results as above. (Specifically, they correspond to the respective cases $\cos\phi=0$ and $\sin\phi=0$.) Intuitively, if a circle's center lies on an axis of the ellipse, then the points of intersection with the ellipse have reflective symmetry over that axis. If there are four distinct points (or two, or none), then we cannot choose three to be the vertices of our equilateral triangle; consequently, there must be only three points of intersection, with one of them on the axis, serving as the point of tangency for the circle and ellipse.
As for the case $f=0$ ... I'll just irresponsibly call it extraneous. (The method of resultants tends to spawn such things.)
Well, since the distances form a Pythagorean triple the choice was not that random. You are on the right track and reflection is a great idea, but you need to take it a step further.
Check that in the (imperfect) drawing below $\triangle RBM$, $\triangle AMQ$, $\triangle MPC$ are equilateral, since they each have two equal sides enclosing angles of $\frac{\pi}{3}$. Furthermore, $S_{\triangle ARM}=S_{\triangle QMC}=S_{\triangle MBP}$ each having sides of length 3,4,5 respectively (sometimes known as the Egyptian triangle as the ancient Egyptians are said to have known the method of constructing a right angle by marking 12 equal segments on the rope and tying it on the poles to form a triangle; all this long before the Pythagoras' theorem was conceived)
By construction the area of the entire polygon $ARBPCQ$ is $2S_{\triangle ABC}$
On the other hand
$$ARBPCQ= S_{\triangle AMQ}+S_{\triangle MPC}+S_{\triangle RBM}+3S_{\triangle ARM}\\=\frac{3^2\sqrt{3}}{4}+\frac{4^2\sqrt{3}}{4}+\frac{5^2\sqrt{3}}{4}+3\frac{1}{2}\cdot 3\cdot 4 = 18+\frac{25}{2}\sqrt{3}$$
Hence
$$S_{\triangle ABC}= 9+\frac{25\sqrt{3}}{4}$$
Best Answer
As the figure shows the distances between the base of the tower $P$ and each of the three vertices $A,B,C$ of equilateral $\triangle ABC$ is
$ PA = h \cot \alpha = h p $
$ PB = h \cot \beta = h q $
$ PC = h \cot \gamma = h r $
Now using Barycentric coordinates, and taking $A$ to be the origin, we can express point $P$ in terms of $B $ and $C$ as follows
$ P = c_1 B + c_2 C $
so that
$ AP = P = c_1 B + c_2 C , BP = (c_1 - 1) B + c_2 C , CP = c_1 B + (c_2 - 1) C $
I'll assume that the side length $a = 1$. Taking the magnitude of these vectors, and noting that $|B| = |C| = 1$ and that $B \cdot C = \dfrac{1}{2} $, then we can write the following three equations
$c_1^2 + c_2^2 + c_1 c_2 = h^2 p^2 \hspace{10pt} (*) $
$(c_1 - 1)^2 + c_2^2 + (c_1 - 1) c_2 = h^2 q^2 $
$ c_1^2 + (c_2 - 1)^2 + c_1 (c_2 - 1) = h^2 r^2 $
We want to eliminate $c_1, c_2$, so subtract each pair of equations, i.e. (1)-(2) , and (1) - (3), we get
$ 2 c_1 + c_2 = h^2 (p^2 - q^2) + 1 $
$ c_1 + 2 c_2 = h^2 (p^2 - r^2) +1$
Solving this $2 \times 2$ system, gives us
$ c_1 = \dfrac{1}{3} ( h^2 (p^2 - 2 q^2 + r^2 ) + 1 )$
$ c_2 = \dfrac{1}{3} ( h^2 (p^2 - 2 r^2 + q^2 ) + 1 )$
substituting $c_1, c_2$ into Eq. $(*)$ and expanding, gives us
$ h^4 K + h^2 L + M = 0 $
where
$K = p^4 + q^4 + r^4 - p^2 q^2 - p^2 r^2 - q^2 r^2 $
$ L = - ( p^2 + q^2 + r^2 ) $
$ M = 1 $
Combining the $h^2 $ coefficients, we get the quadratic equation (in $h^2$)
From this, using the quadratic formula, the solution is
$ h^2 = \dfrac{ -L \pm \sqrt{ L^2 - 4 K } }{ 2 K } $
we have
$ - L = p^2 + q^2 + r^2 $
and
$L^2 - 4 K = p^4 + q^4 + r^4 + 2 p^2 q^2 + 2 p^2 r^2 + 2 q^2 r^2 - 4 (p^4 + q^4 + r^4 - p^2 q^2 - p^2 r^2 - q^2 r^2) $
And this reduces to
$ L^2 - 4 K = 3 ( 2 p^2 q^2 + 2 p^2 r^2 + 2 q^2 r^2 - p^4 - q^4 - r^4 ) $
Hence,
$ h^2 = \dfrac{ p^2 + q^2 + r^2 \pm \sqrt{3} \sqrt{2 p^2 q^2 + 2 p^2 r^2 + 2 q^2 r^2 - p^4 - q^4 - r^4} }{ 2 (p^4 + q^4 + r^4 - p^2 q^2 - p^2 r^2 - q^2 r^2)}$