Consider the hyperbola:
$$
\frac{x^2}{a^2}-\frac{y^2}{b^2}=1\tag{1}
$$
Given $f_1,f_2$, the two foci of a hyperbola, one property of a hyperbola is that there is a constant, $\Delta\text{ distance}$, so that any point on the hyperbola, p, satisfies
$$
\Big||p-f_1|-|p-f_2|\Big|=\Delta\text{ distance}\tag{2}
$$
Consider a point at the intersection of the hyperbola and the line between the foci. The $\Delta\text{ distance}$ given in $(2)$ is the distance between the two branches of the hyperbola.
$\hspace{3.3cm}$
Using $(1)$, we get that the distance between the two branches of the hyperbola is $2a$. Therefore,
$$
\Delta\text{ distance}=2a\tag{3}
$$
Consider a point at an infinite distance on the upper right branch of the hyperbola. Because $\triangle gpf_2$ is essentially isosceles, the $\Delta\text{ distance}$ given in $(2)$ is
$$
\Delta\text{ distance}=|f_1-f_2|\cos(\theta)\tag{4}
$$
Using $(1)$, we get that
$$
\begin{align}
\tan(\theta)&=\lim_{x,y\to\infty}\frac yx=\frac ba\\
\cos(\theta)&=\frac{a}{\sqrt{a^2+b^2}}\tag{5}
\end{align}
$$
Combining $(3)$, $(4)$, and $(5)$, we get
$$
|f_1-f_2|=2\sqrt{a^2+b^2}\tag{6}
$$
Thus, if $c$ is the distance from the center of the hyperbola to each of the foci, then $(6)$ gives
$$
c^2=a^2+b^2\tag{7}
$$
[Transferring the discussion in comments to an answer.]
Under a general affine transformation, the image of the center of an ellipse is the center of the ellipse’s image, so that’s easily computed. Unfortunately, the principal axes aren’t usually mapped to the axes of the image, although their images are conjugate diameters of the resulting ellipse. There are a couple of approaches that you might take to determine the semiaxis lengths and rotation angle of the image.
Writing the general conic equation $Ax^2+Bxy+Cy^2+Dx+Ey+F=0$ in matrix form $$\mathbf x^TQ\mathbf x = \begin{bmatrix}x&y&1\end{bmatrix} \begin{bmatrix}A&\frac B2&\frac D2\\\frac B2&C&\frac E2\\\frac D2&\frac E2&F\end{bmatrix} \begin{bmatrix}x\\y\\1\end{bmatrix} =0, \tag1$$ under the point transformation $\mathbf x' = M\mathbf x$ the matrix $Q$ transforms into $M^{-T}QM^{-1}$. (You can easily derive this by substituting for $\mathbf x$ in (1) and rearranging.) So, you could compute this transformation, extract the new coefficients of the conic equation from the result, and then use standard methods and formulas to extract the parameters that you’re interested in.
You can save yourself some work by looking only at the quadratic part of the equation: the shape of the ellipse is only affected by the linear part of the transformation, so you need only focus on the ellipse with equation $${(x\cos\theta+y\sin\theta)^2\over a^2}+{(x\sin\theta-y\cos\theta)^2\over b^2}=1, \tag 2$$ i.e., decompose the affine transformation into the composition of a translation by $(-h,-k)$, a linear transformation, and then another translation. (Note that I’ve corrected the sign error in your original equation.) In matrix form, we have $(x,y)Q(x,y)^T=1$ and the ellipse’s matrix again transforms as $M^{-T}QM^{-1}$, except that now we’re working with $2\times2$ matrices. Apply the transformation and extract the rotation and semiaxis lengths from the resulting coefficients.
One could, with a lot of tedious work or assistance of a symbolic algebra program, work out a general formula for the resulting parameters, but it’s going to be very ugly. If you’re planning to code this, it’s probably best to take it step-by-step instead of trying to transcribe some complex formula.
There are a couple of alternative approaches that might be interesting, though they don’t lead to simple general formulas, either, and might not be any less work than the above. First notice that you can parameterize your original ellipse as $\mathbf c + \mathbf u\cos t + \mathbf v\sin t$, where $\mathbf c = (h,k)$, $\mathbf u = a(\cos\theta,\sin\theta)$ and $\mathbf v=b(-\sin\theta,\cos\theta)$. Focusing again on the linear part of the transformation, the new ellipse is a translate of $\mathbf r(t) = M\mathbf u\cos t+M\mathbf v\sin t$. The tangent to an ellipse at a vertex is orthogonal to the corresponding principal axis, so a way to find these vertices is to solve $\mathbf r(t)\cdot\mathbf r'(t)=0$ for $t$. Unless you ended up with a circle, there will be two pairs of diametrically opposite solutions. The semiaxis lengths are then just the two lengths $\|\mathbf r(t)\|$ and the rotation angle is also easily computed.
Another possible approach takes advantage of the fact that the original ellipse (again after translation back to the origin) is itself obtained from the unit circle by scaling and rotation, so that this ellipse when further transformed is obtained from the unit circle via some composite transformation given by the matrix $M'=MRS$, with $$R = \begin{bmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{bmatrix} \text{ and } S = \begin{bmatrix}a&0\\0&b\end{bmatrix}.$$ The SVD of $M'$ refactors it into the product of a rotation, a diagonal matrix (i.e., scaling along orthogonal directions) and another rotation. The half-axis lengths of the ellipse are therefore given by the central diagonal matrix, and, since the first rotation leaves the unit circle unchanged, the ellipse’s principal axis directions are the columns of the left-hand rotation matrix in the SVD. The rotation angle is easily recovered from them.
Best Answer
As long as $d$ is sufficiently small (where "sufficiently small" depends on the eccentricity of the ellipse), you can proceed as follows:
Start at point $P_0$, and set off around the ellipse in steps of (Euclidean) length $d$, leaving point $P_i$ at the $i$th step. When you reach or pass the original point $P_0$, leave a point $P_n$ there, and stop.
If $P_n$ and $P_0$ coincide, we are done, Otherwise, decrease $d$ continuously until they do. Now we have $n$ equally spaced points on the ellipse. And we can repreat this procedure to find $n+1$ equally spaced points, and so on.
There are two things that can go wrong: