Given a quadric of Cartesian equation:
$$
\small
a_{11}\,x^2+a_{22}\,y^2+a_{33}\,z^2+2\,a_{12}\,x\,y + 2\,a_{13}\,x\,z + 2\,a_{23}\,y\,z + 2\,a_{14}\,x + 2\,a_{24}\,y + 2\,a_{34}\,z + a_{44} = 0
$$
the first thing is to build the respective real and symmetric matrices:
$$
A_1 \equiv
\begin{bmatrix}
a_{11} & a_{12} & a_{13} & a_{14} \\
a_{12} & a_{22} & a_{23} & a_{24} \\
a_{13} & a_{23} & a_{33} & a_{34} \\
a_{14} & a_{24} & a_{34} & a_{44} \\
\end{bmatrix},
\quad \quad \quad
A_2 \equiv
\begin{bmatrix}
a_{11} & a_{12} & a_{13} \\
a_{12} & a_{22} & a_{23} \\
a_{13} & a_{23} & a_{33} \\
\end{bmatrix}.
$$
So, whenever:
$$
\det(A_1) \ne 0,
\quad \quad \quad
\det(A_2) \ne 0
$$
it's a non-degenerate quadric whose center is the solution of the system of linear equations:
$$
\begin{cases}
a_{11}\,x + a_{12}\,y + a_{13}\,z + a_{14} = 0 \\
a_{12}\,x + a_{22}\,y + a_{23}\,z + a_{24} = 0 \\
a_{13}\,x + a_{23}\,y + a_{33}\,z + a_{34} = 0 \\
\end{cases}
\quad \quad \Rightarrow \quad \quad
C=(x_c,y_c,z_c).
$$
In particular, if $\det(A_1)<0$ and the eigenvalues of $A_2$ have the same sign, it's a real ellipsoid.
So, if the eigenvalues of $A_2$ are $\lambda_1,\,\lambda_2,\,\lambda_3$ and the orthonormal eigenvectors are $\mathbf{u}_1,\,\mathbf{u}_2,\,\mathbf{u}_3$:
$$
\boxed{\mathcal{E}(u,v) = C +
\frac{\mathbf{u}_1}{\sqrt{|\lambda_1|}}\,\cos u\,\cos v +
\frac{\mathbf{u}_2}{\sqrt{|\lambda_2|}}\,\cos u\,\sin v +
\frac{\mathbf{u}_3}{\sqrt{|\lambda_3|}}\,\sin u\;}
$$
with $-\frac{\pi}{2} \le u \le \frac{\pi}{2}$ and $0 \le v < 2\pi$ is a parametrization of the real ellipsoid considered.
Therefore, the principal ellipses are trivially obtained by setting:
$$
\boxed{\mathcal{E}_1(t) = \mathcal{E}\left(0,t\right),
\quad \quad \quad
\mathcal{E}_2(t) = \mathcal{E}\left(t,0\right),
\quad \quad \quad
\mathcal{E}_3(t) = \mathcal{E}\left(t,\frac{\pi}{2}\right)\;}
$$
where, in all three parametrizations, we have $0 \le t < 2\pi$.
ADDENDUM:
Given the quadric of Cartesian equation:
$$
\left(15+4\sqrt{2}\right)x^2+\left(15-4\sqrt{2}\right)y^2+13\,z^2 + 4\,x\,y + \left(8+4\sqrt{2}\right)x\,z + \left(8-4\sqrt{2}\right)y\,z - 1 = 0
$$
it's easy to verify that it's a real ellipsoid with center $O$ and the eigen/values/vectors of $A_2$ are:
$$
\begin{aligned}
& \lambda_1 = 25,
\quad \quad \;\;
\mathbf{u}_1 = \left(\sqrt{2}+1,\,\sqrt{2}-1,\,\sqrt{2}\right); \\
& \lambda_2 = 9,
\quad \quad \quad
\mathbf{u}_2 = \left(\sqrt{2}-2,\,0,\,1\right); \\
& \lambda_3 = 9,
\quad \quad \quad
\mathbf{u}_3 = \left(2\sqrt{2}-3,\,1,\,0\right); \\
\end{aligned}
$$
but since $\lambda_2=\lambda_3$ not always $\mathbf{u}_2\perp\mathbf{u}_3$; in fact, it's $\mathbf{u}_2\cdot\mathbf{u}_3 = 10-7\sqrt{2}\ne 0$.
So, in these cases it's necessary to apply the Gram–Schmidt process, by which we obtain:
$$
\begin{aligned}
& \mathbf{u}_1 = \left(\frac{2+\sqrt{2}}{4},\,\frac{2-\sqrt{2}}{4},\,\frac{1}{2}\right); \\
& \mathbf{u}_2 = \left(\frac{\sqrt{2}-2}{\sqrt{7-4\sqrt{2}}},\,0,\,\frac{1}{\sqrt{7-4\sqrt{2}}}\right); \\
& \mathbf{u}_3 = \left(\frac{2\sqrt{2}-3}{2\sqrt{74-52\sqrt{2}}},\,\frac{7-4\sqrt{2}}{2\sqrt{74-52\sqrt{2}}},\,\frac{7\sqrt{2}-10}{2\sqrt{74-52\sqrt{2}}}\right); \\
\end{aligned}
$$
which are eigenvectors of $A_2$ of unit length and orthogonal to each other, as desired.
In all those cases where the eigenvectors are already orthogonal, it's sufficient to normalize them!