Problem
I want to find the geodetic coordinates (longitude, latitude) of the ellipse defined by the intersection between a plane and ellipsoid. Point $P$ is an observer away from Earth, and the ellipse defines the boundary between the visible earth at $P$ and the not-visible portion. That is, a ray of light coming from $P$ will be tangent to all points on the ellipse.
I am assuming that the plane encompassing the elliptical boundary will be normal to the vector $\vec{OP}$. For distances very large (greater than 8 million km), it's safe to assume the ellipse is a great circle (encompasses the origin). Assuming this to be a great circle for the moon will result in a little bit of error, but I'm not sure if I can know $\vec{OA}$ upfront.
Questions
- Did I perform the derivation correctly?
- What if I don't know $d$ upfront, is there any way to solve this?
- In my case, $d$ would be the distance $OA$ in my diagram
- I know everything about points $p$ and $s$
The Solution
Theory
We want to perform an affine map of the ellipsoid and plane to the unit sphere, find a set of orthogonal vectors on the plane of the intersection circle which allow us to get a parametric equation defining the intersection circle, then perform the reverse affine map to get back to the original ellipsoid and plane. From there, we can convert the cartesian coordinates to spherical, and then convert spherical coordinates to geodetic.
My attempt
Affine Map to Unit Sphere
We have the ellipsoid and plane defined by the following equations, respectively.
$$ \frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1, \\ n_xx + n_yy + n_zz = d $$
To map the ellipsoid to the unit sphere we choose $u = \frac{x}{a}, v = \frac{y}{b}, z = \frac{z}{c}$, which yields
$$ u^2 + v^2 + w^2 = 1, \\ n_xax + n_yby + n_zcz = d $$
Choose $m_u = n_xa, m_v = n_yb, m_w = n_zc$ to get the Hesse normal form of the plane after the affine map
$$ m_uu + m_vv + m_ww = \delta $$
The unit normal vector of the new plane is $ \vec{m} = \begin{bmatrix} m_u \\ m_v \\ m_w \end{bmatrix}$, and the vector to the center of the intersection circle is $\vec{e_o} = \delta \vec{m}$. The radius of the circle is $\rho = \sqrt{1 - \delta^2}$. The orthogonal vectors $\vec{e_1}, \vec{e_2}$ that form a basis for the intersection circle are then given by the following.
If $m_w = \pm 1$ (the plane is parallel to the xy-plane),
$$ \vec{e_1} = \begin{bmatrix} \rho \\ 0 \\ 0 \end{bmatrix},\ \vec{e_2} = \begin{bmatrix} 0 \\ \rho \\ 0 \end{bmatrix} $$.
If $m_w \neq \pm 1$ (the plane is not parallel to the xy-plane),
$$ \vec{e_1} = \frac{\rho}{\sqrt{m_u^2 + m_v^2}} \begin{bmatrix} m_v \\ -m_u \\ 0 \end{bmatrix},\ \vec{e_2} = \vec{m} \times \vec{e_1} $$.
The parametric equation of the intersection circle is given by
$$ \vec{u} = \vec{e_o} + \vec{e_1}\cos(t) + \vec{e_2}\sin(t),\ t \in [0, 2\pi) $$
Reverse Affine Map
Now we must reverse the affine transformation to obtain the vector defining the center of the ellipse intersection $\vec{f_0}$ and the orthogonal vectors $\vec{f_1}, \vec{f_2}$ that form a basis on the ellipse. Since
$$ \vec{m} = \begin{bmatrix} m_u \\ m_v \\ m_w \end{bmatrix} = \begin{bmatrix} n_xa \\ n_yb \\ n_zc \end{bmatrix} \\ $$
then,
$$ \begin{align} \vec{e_0} &\mapsto \vec{f_0} = d\vec{m},\\ \rho &\mapsto R = \sqrt{1 - d^2} = \sqrt{1 - (n_xx + n_yy + n_zz)^2)}, \\ \vec{e_1} \mapsto \vec{f_1} &= \begin{cases} \begin{bmatrix}R \\ 0 \\ 0 \end{bmatrix}, & m_w = n_zc = \pm 1 \\ \frac{R}{\sqrt{(n_xa)^2 + (n_yb)^2}}\begin{bmatrix}n_yb \\ -n_xa \\ 0 \end{bmatrix} & m_w = n_zc \neq \pm 1 \\ \end{cases}\\ \vec{e_2} \mapsto \vec{f_2} &= \begin{cases} \begin{bmatrix}0 \\ R \\ 0 \end{bmatrix}, & m_w = n_zc = \pm 1 \\ \vec{m} \times \vec{f_1} & m_w = n_zc \neq \pm 1 \\ \end{cases} \end{align} $$
Therefore, the parametric equation for the intersection ellipse is given by
$$ \vec{x} = \vec{f_0} + \vec{f_1}\cos(t) + \vec{f_2}\sin(t), \ t \in [0, 2\pi) $$
Retrieving Geodetic Coordinates
We are looking for the geodetic coordinates (longitude $\lambda$, latitude $\phi$) of $\vec{x}$. To do this we must convert the cartesian coordinates to spherical coordinates, and then use relationships between spherical and geodetic coordinates. Let $(r, \theta, \psi)$ be the spherical coordinates of $\vec{x}$ where $r$ is the distance from the origin of the ellipsoid, $\theta$ is the polar angle, and $\psi$ is the azimuth. Then,
$$ \begin{align} r &= \sqrt{x^2 + y^2 + z^2}\\ \theta &= \arccos(z/r)\\ \psi &= \begin{cases} \arctan(y/x), & x, y \neq 0\\ \text{undefined}, & x = y = 0 \\ \end{cases} \end{align} $$
The relationship between the spherical coordinates and geodetic coordinates are given by
$$ \begin{align} \lambda &= \psi \\ \phi &= \frac{\pi}{2} - \theta \end{align} $$
And we are done.