(Se supone que esto es un comentario a la respuesta de ja72, pero se hizo demasiado largo).
La respuesta de ja72 menciona la fórmula
$$\tan\,2\varphi = \frac{c}{a-b}$$
Considero que es un poco derrochador evaluar una arctangente y posteriormente pasarla como argumento de una unción trigonométrica mucho más tarde; para evitarlo, podemos utilizar la fórmula del ángulo doble para la tangente y la fórmula "Citardauq" en tándem para obtener la relación
$$\tan\,\varphi=\frac{c}{a-b+(\mathrm{sign}\,c)\sqrt{c^2+(a-b)^2}}$$
Denotando esta expresión como $t$ podemos sustituirlo por la matriz de rotación
$$\frac1{\sqrt{1+t^2}}\begin{pmatrix}1&-t\\t&1\end{pmatrix}$$
La respuesta de ja72 ya daba las fórmulas de los ejes; recuerda que $(p\sin\,u,q\cos\,u)$ y $(p\cos\,u,q\sin\,u)$ son la misma elipse atravesada de forma diferente, por lo que puedes rotar posteriormente cualquiera de las dos expresiones que elijas con la matriz de rotación que te di antes.