The equation $$ax^2+2hxy+by^2=c$$ represents a conic rotated about the axis (depending on $a,b,h,c$) with centre $(0,0)$.
You can always rotate the axis to remove the $xy$ term. This makes it easy to find the focus, length of major axis, etc.
To remove $xy$ term without shifting origin, lets rotate the axis by an angle $\theta$.
Let $X,Y$ represent the new coordinate system. Then,
$$x=X\cos\theta-Y\sin\theta$$
$$y=X\sin\theta+Y\cos\theta$$
Substituting and letting coefficient of $XY=0$, you get $$\tan2\theta=\frac{2h}{a-b}$$
For example, lets find $e$ of the following equation (which represents an ellipse): $$4x^2+2y^2+2\sqrt{3}xy=1$$
Raotate the axis by $$\tan2\theta=\frac{2h}{a-b}=\sqrt{3}$$
$$\theta=\frac{\pi}{6}$$
Now, substitute $$x=X\frac{\sqrt{3}}{2}-Y\frac{1}{2}$$
$$y=Y\frac{\sqrt{3}}{2}+X\frac{1}{2}$$
You will get $$5X^2+Y^2=5$$
Now you can find eccentricity using the formula you specified. (Rotating the ellipse or hyperbola will not change the eccentricity and the length of major\minor axis.
This is the graph for the two equations:
(Black represents the new ellipse)
Parameterization of an ellipse by angle from the center is
$$
\gamma(\phi)=(\cos(\phi),\sin(\phi))\frac{ab}{\sqrt{a^2\sin^2(\phi)+b^2\cos^2(\phi)}}
$$
and
$$
\left|\gamma'(\phi)\right|=ab\sqrt{\frac{a^4\sin^2(\phi)+b^4\cos^2(\phi)}{\left(a^2\sin^2(\phi)+b^2\cos^2(\phi)\right)^3}}
$$
and integrating $\left|\gamma'(\phi)\right|$ gets extremely messy.
So instead, we use $\theta$, where
$$
b\tan(\theta)=a\tan(\phi)
$$
Then
$$
\gamma(\theta)=(a\cos(\theta),b\sin(\theta))
$$
and
$$
\left|\gamma'(\theta)\right|=\sqrt{a^2\sin^2(\theta)+b^2\cos^2(\theta)}
$$
Now, integrating $\left|\gamma'(\theta)\right|$ is a lot simpler.
$$
\int\left|\gamma'(\theta)\right|\,\mathrm{d}\theta
=b\,\mathrm{EllipticE}\left(\theta,\frac{b^2-a^2}{b^2}\right)
$$
However, to go from from arc length to angle, we still need to invert the Elliptic integral.
Solution of the problem given using Mathematica
Here we get the solution to 30 places.
In[]= Phi[a_,b_,s_,opts:OptionsPattern[]] := Block[ {t}, t+ArcTan[(b-a)Tan[t]/(a+b Tan[t]^2)] /. FindRoot[b EllipticE[t,(b^2-a^2)/b^2] == s, {t, 0}, opts]]
In[]= Phi[4, 2.9`32, 3.31`32, WorkingPrecision->30]
Out[]= 0.87052028743193111752524449959
Thus, $\phi\doteq0.87052028743193111752524449959\text{ radians}$.
Let me describe the algorithm a bit.
FindRoot[b EllipticE[t,(b^2-a^2)/b^2] == s, {t, 0}, opts]
inverts the elliptic integral to get $\theta$ from $a$, $b$, and $s$.
Now, we want to find $\phi$ so that $b\tan(\theta)=a\tan(\phi)$. However, simply using $\tan^{-1}\left(\frac ba\tan(\theta)\right)$ will only return a result in $\left(-\frac\pi2,\frac\pi2\right)$. To get the correct value, we use the relation
$$
\begin{align}
\tan(\phi-\theta)
&=\frac{\tan(\phi)-\tan(\theta)}{1+\tan(\phi)\tan(\theta)}\\
&=\frac{\frac ba\tan(\theta)-\tan(\theta)}{1+\frac ba\tan(\theta)\tan(\theta)}\\
&=\frac{b\tan(\theta)-a\tan(\theta)}{a+b\tan(\theta)\tan(\theta)}\\
\phi
&=\theta+\tan^{-1}\left(\frac{(b-a)\tan(\theta)}{a+b\tan^2(\theta)}\right)
\end{align}
$$
This is why we have
t+ArcTan[(b-a)Tan[t]/(a+b Tan[t]^2)]
Best Answer
Parametrize the ellipse as $$(x(t), y(t)) = (w \cos t, h \sin t), \quad 0 \le t < 2\pi, \quad w > h > 0. \tag{1}$$ Then the point of tangency corresponds to a value of $t$, say $t = t_0$, such that $x(t_0) = h$, or $$t_0 = \arccos \frac{h}{w}. \tag{2}$$ The slope $m(t)$ of the tangent line at $t = t_0$ is given by $$m(t_0) = \left[\frac{dy/dt}{dx/dt}\right]_{t=t_0} = \left[-\frac{h}{w} \cot t \right]_{t=t_0} = -\frac{h^2}{w \sqrt{w^2 - h^2}}. \tag{3}$$ The angle of the tangent line is the arctangent of the slope; i.e. $$\tan \theta = m(t_0). \tag{4}$$ So to obtain the semimajor axis $w$ in terms of a given fixed semiminor axis $h$ and desired tangent angle $\theta$, we solve the equation $(4)$ after substituting in $(3)$, giving $$w = \frac{h}{\sqrt{2}} \sqrt{1 + \sqrt{1 + 4 \cot^2 \theta}}.$$