Earth section paths

To define the section plane select any third point ${\displaystyle R_{0}}$ not on the line from ${\displaystyle R_{1}}$ to ${\displaystyle R_{2}}$. Choosing ${\displaystyle R_{0}}$ to be on the surface normal at ${\displaystyle P_{1}}$ will define the normal section at ${\displaystyle P_{1}}$. If ${\displaystyle R_{0}}$ is the origin then the earth section is the great ellipse. (The origin would be co-linear with 2 antipodal points so a different point must be used in that case). Since there are infinitely many choices for ${\displaystyle R_{0}}$, the above problem is really a class of problems (one for each plane). Let ${\displaystyle R_{0}}$ be given. To put the equation of the plane into the standard form, ${\displaystyle lx+my+nz=d}$, where ${\displaystyle \ l^{2}+m^{2}+n^{2}=1}$, requires the components of a unit vector, ${\displaystyle {\hat {N}}=(l,m,n)}$, normal to the section plane. These components may be computed as follows: The vector from ${\displaystyle R_{0}}$ to ${\displaystyle R_{1}}$ has components ${\displaystyle V_{0}=R_{1}-R_{0}}$, and the vector from ${\displaystyle R_{1}}$ to ${\displaystyle R_{2}}$ has components ${\displaystyle V_{1}=R_{2}-R_{1}}$. Therefore, ${\displaystyle {\hat {N}}}$ = ${\displaystyle unit(V_{0}}$×${\displaystyle V_{1}}$), where ${\displaystyle unit(V)}$ is the unit vector in the direction of ${\displaystyle V}$. The orientation convention used here is that ${\displaystyle {\hat {N}}}$ points to the left of the path. If this is not the case then redefine ${\displaystyle V_{0}}$ = -${\displaystyle V_{0}}$. Finally, the parameter d for the plane may be computed using the dot product of ${\displaystyle {\hat {N}}}$ with a vector from the origin to any point on the plane, such as ${\displaystyle R_{1}}$, i.e. d = ${\displaystyle {\hat {N}}\cdot {R_{1}}}$. The equation of the plane (in vector form) is thus ${\displaystyle {\hat {N}}}$ ⋅ ${\displaystyle R}$ = d, where ${\displaystyle R}$ is the position vector of (x, y, z).

Examination of the ENU to ECEF transformation reveals that the ECEF coordinates of a unit vector pointing east at any point on the ellipsoid is: ${\displaystyle \mathbf {\hat {e}} }$=${\displaystyle (-\sin \lambda ,\cos \lambda ,0)}$, a unit vector pointing north is ${\displaystyle \mathbf {\hat {n}} }$=${\displaystyle (-\sin \phi \cos \lambda ,-\sin \phi \sin \lambda ,\cos \phi )}$, and a unit vector pointing up is ${\displaystyle \mathbf {\hat {u}} }$=${\displaystyle (\cos \phi \cos \lambda ,\cos \phi \sin \lambda ,\sin \phi )}$. A vector tangent to the path is: ${\displaystyle \mathbf {t} ={\hat {N}}\times \mathbf {\hat {u}} }$ so the east component of ${\displaystyle \mathbf {t} }$ is ${\displaystyle \mathbf {t} \cdot \mathbf {\hat {e}} }$, and the north component is ${\displaystyle \mathbf {t} \cdot \mathbf {\hat {n}} }$. Therefore, the azimuth may be obtained from a two-argument arctangent function, ${\displaystyle \alpha }$=${\displaystyle \operatorname {atan2} (\mathbf {t} \cdot \mathbf {\hat {e}} ,\mathbf {t} \cdot \mathbf {\hat {n}} )}$. Use this method at both ${\displaystyle P_{1}}$ and ${\displaystyle P_{2}}$ to get ${\displaystyle \alpha _{1}}$ and ${\displaystyle \alpha _{2}}$.

The (non-trivial) intersection of a plane and ellipsoid is an ellipse. Therefore, the arc length, ${\displaystyle s_{12}}$, on the section path from ${\displaystyle P_{1}}$ to ${\displaystyle P_{2}}$ is an elliptic integral that may be computed to any desired accuracy using a truncated series. Before this can be done the ellipse must be defined and the limits of integration computed. Let the ellipsoid given by ${\displaystyle {\frac {x^{2}}{a^{2}}}+{\frac {y^{2}}{a^{2}}}+{\frac {z^{2}}{b^{2}}}=1}$, and let ${\displaystyle p={\sqrt {l^{2}+m^{2}}}}$. If p=0 then the section is a horizontal circle of radius ${\displaystyle a{\sqrt {1-{\frac {d^{2}}{b^{2}}}}}}$, which has no solution if ${\displaystyle |d|>b}$.

If p>0 then Gilbertson[1] showed that the ECEF coordinates of the center of the ellipse is ${\displaystyle {R_{c}}={\frac {d}{C}}(la^{2},ma^{2},nb^{2})}$, where ${\displaystyle C=a^{2}p^{2}+b^{2}n^{2}}$,

the semi-major axis is ${\displaystyle a^{*}=a{\sqrt {1-{\frac {d^{2}}{C}}}}}$, in the direction ${\displaystyle \mathbf {\hat {i^{*}}} =({\frac {m}{p}},{\frac {-l}{p}},0)}$, and the semi-minor axis is ${\displaystyle b^{*}={\frac {b}{\sqrt {C}}}a^{*}}$, in the direction ${\displaystyle \mathbf {\hat {j^{*}}} =({\frac {ln}{p}},{\frac {mn}{p}},-p)}$, which has no solution if ${\displaystyle |d|>{\sqrt {C}}}$.

The polar form relative to center for the equation of an ellipse is ${\displaystyle R(\theta )={\frac {b^{*}}{\sqrt {1-e^{2}\cos ^{2}\theta }}}}$, where ${\displaystyle e^{2}=1-{\frac {(b^{*})^{2}}{(a^{*})^{2}}}}$, relates to the ellipse eccentricity, not the spheroid eccentricity (see ellipse). Let P be a point on the ellipse and ${\displaystyle R=ECEF(P)}$, then the vector from ${\displaystyle {R_{c}}}$ to ${\displaystyle R}$ has components ${\displaystyle R-R_{c}}$. Using an argument similar to the one for azimuth above, let ${\displaystyle \mathbf {\hat {w}} =unit(R-R_{c})}$, then ${\displaystyle \cos \theta =\mathbf {\hat {w}} \cdot \mathbf {\hat {i^{*}}} }$, and ${\displaystyle \sin \theta =\mathbf {\hat {w}} \cdot \mathbf {\hat {j^{*}}} }$, and ${\displaystyle \theta =\operatorname {atan2} (\sin \theta ,\cos \theta )}$. In this way we obtain the central angles ${\displaystyle \theta _{1}}$ and ${\displaystyle \theta _{2}}$ corresponding to ${\displaystyle P_{1}}$ and ${\displaystyle P_{2}}$ respectively. Care must be taken to ensure that ${\displaystyle \theta _{1}}$ ≤ ${\displaystyle \theta _{2}}$ ≤ ${\displaystyle \theta _{1}+\pi }$. Then the arc length along the ellipse is given by ${\displaystyle s_{12}}$ =${\displaystyle \int _{\theta _{1}}^{\theta _{2}}{\sqrt {R^{2}+(R^{'})^{2}}}d\theta .}$ Substituting ${\displaystyle R(\theta )}$ above into this formula, performing the indicated operations, using one more term than Gilbertson's expression and regrouping, results in ${\displaystyle s_{12}=ArcLength(\theta _{1},\theta _{2})=b^{*}({c_{0}}\Delta \theta +{c_{1}}\Delta s2+{c_{2}}\Delta s4+{c_{3}}\Delta s6)}$, where

${\displaystyle {\begin{aligned}\Delta \theta &=\theta _{2}-\theta _{1},\\[6pt]\Delta s2&=\sin(2\theta _{2})-\sin(2\theta _{1}),\\[6pt]\Delta s4&=\sin(4\theta _{2})-\sin(4\theta _{1}),\\[6pt]\Delta s6&=\sin(6\theta _{2})-\sin(6\theta _{1}),\\[6pt]{c_{0}}&=1+e^{2}(4096+3328e^{2}+2880e^{4})/16384,\\[6pt]{c_{1}}&=e^{2}(512+384e^{2}+380e^{4})/4096,\\[6pt]{c_{2}}&=-e^{4}(64+80e^{2})/16384,\\[6pt]{c_{3}}&=-60e^{6}/12288.\\[6pt]\end{aligned}}}$

Alternatively, expansions for the Meridian arc may be used here by replacing the spheroid eccentricity with the section ellipse eccentricity.

Construct the tangent vector at ${\displaystyle {P_{1}}}$, ${\displaystyle \mathbf {\hat {t_{1}}} =\mathbf {\hat {n_{1}}} \cos {\alpha _{1}}+\mathbf {\hat {e_{1}}} \sin {\alpha _{1}}}$, where ${\displaystyle \mathbf {\hat {n_{1}}} }$ and ${\displaystyle \mathbf {\hat {e_{1}}} }$ are unit vectors pointing north and east (respectively) at ${\displaystyle {P_{1}}}$. Pick a vector, ${\displaystyle V_{0}}$, to define the section plane, paying attention to orientation. Observe that ${\displaystyle V_{0}}$ must not be in span{${\displaystyle \mathbf {\hat {n_{1}}} ,\mathbf {\hat {e_{1}}} }$} (otherwise the plane would be tangent to the earth at ${\displaystyle {P_{1}}}$, so no path would result). The normal vector ${\displaystyle {\hat {N}}}$ = ${\displaystyle unit(V_{0}}$×${\displaystyle \mathbf {\hat {t_{1}}} }$), together with ${\displaystyle {P_{1}}}$ defines the plane.

This is a 2-d problem in span{${\displaystyle \mathbf {\hat {i^{*}}} ,\mathbf {\hat {j^{*}}} }$}, which will be solved with the help of the arc length formula above. The basic approach is to use Newton-Raphson iteration to arrive at ${\displaystyle {P_{2}}}$. The basis of the estimate is that the position vector of any point on the section ellipse may be expressed in terms of the position vector of the center and the central angle as ${\displaystyle V={V_{c}}+{r(\theta )}(\mathbf {\hat {i^{*}}} \cos \theta +\mathbf {\hat {j^{*}}} \sin \theta )}$. To get an initial estimate of ${\displaystyle {\theta _{2}}}$, let ${\displaystyle {V_{1}}={R_{1}}-R{_{c}}}$, ${\displaystyle {\theta _{1}}}$=Central_Angle${\displaystyle ({V_{1}})}$ (see the arc length section above), ${\displaystyle {r_{1}}=r({\theta _{1}})}$, ${\displaystyle \Delta \theta ={\frac {\delta }{r_{1}}}}$.

Now initialize ${\displaystyle \theta _{2}}$ = ${\displaystyle \theta _{1}+\Delta \theta }$, and iterate the following steps:

${\displaystyle {\begin{aligned}s&=ArcLength({\theta _{1}},{\theta _{2}}),\\[6pt]Err&=\delta -s,\\[6pt]s'({\theta })&={\frac {b^{*}}{(1-e^{2}\cos ^{2}\theta )}}{\sqrt {\frac {(1-(2-e^{2})e^{2}\cos ^{2}\theta }{1-e^{2}\cos ^{2}\theta }}},\\[6pt]\Delta \theta &={\frac {Err}{s'({\theta _{2}})}},\\[6pt]\theta _{2}&=\theta _{2}+\Delta \theta ,\\[6pt]\end{aligned}}}$

exit when ${\displaystyle abs(\Delta \theta )<10^{-12}}$

No more than three iterations are usually necessary, although nearly antipodal cases can be problematic. Finally, let ${\displaystyle {V_{2}}={V_{c}}+{r(\theta _{2})}(\mathbf {\hat {i^{*}}} \cos \theta _{2}+\mathbf {\hat {j^{*}}} \sin \theta _{2})}$, and ${\displaystyle {P_{2}}}$ = ECEF_to_Geo${\displaystyle ({V_{2}})}$ using Bowring's 1985 algorithm,[2] or the algorithm here.

Alternatively, inversion of the arc length series may be used to avoid iterations.

Azimuth may be obtained by the same method as the indirect problem: ${\displaystyle {\alpha _{2}}}$=${\displaystyle \operatorname {atan2} (\mathbf {t_{2}} \cdot \mathbf {\hat {e_{2}}} ,\mathbf {t_{2}} \cdot \mathbf {\hat {n_{2}}} )}$, where the subscript 2 indicates evaluation of the associated quantity at ${\displaystyle P_{2}}$.

Let ${\displaystyle R_{0}}$ be the origin, so that ${\displaystyle {V_{0}}}$ = the position vector of ${\displaystyle R_{1}}$. The above approach provides an alternative to that of others, such as Bowring.[3]

The normal section at ${\displaystyle P_{1}}$ is determined by letting ${\displaystyle {V_{0}}}$ = ${\displaystyle \mathbf {\hat {u_{1}}} }$ (the surface normal at ${\displaystyle P_{1}}$). The above approach provides an alternative to that of others, such as Bowring.[4]

The mean normal section from ${\displaystyle P_{1}}$ to ${\displaystyle P_{2}}$ is determined by letting ${\displaystyle {V_{0}}}$ = ${\displaystyle 0.5(\mathbf {\hat {u_{1}}} +\mathbf {\hat {u_{2}}} )}$. This is a good approximation to the geodesic from ${\displaystyle P_{1}}$ to ${\displaystyle P_{2}}$ for aviation or sailing.

A class of sections may be imagined by rotating ${\displaystyle \mathbf {\hat {u_{1}}} }$ about the chord connecting ${\displaystyle P_{1}}$ and ${\displaystyle P_{2}.}$ All of these may be solved with the single approach above.

on an earth section path may be found by dropping the subscripts on the given section; ${\displaystyle {\hat {N_{1}}}=(l,m,n)}$, ${\displaystyle d_{1}=d}$, and setting ${\displaystyle {\hat {N_{2}}}=(0,0,1)}$, so that ${\displaystyle {\hat {N_{3}}}=(m,-l,0)}$. Then solve for ${\displaystyle d_{2}=z_{2}}$ such that ${\displaystyle B^{2}=AC}$.

Since ${\displaystyle B=0}$, and ${\displaystyle A\neq 0}$, we must have ${\displaystyle C=0}$. Plugging t into ${\displaystyle R}$ = ${\displaystyle R_{0}+t{\hat {N_{3}}}}$, gives the points of intersection of the earth sections. Alternatively, just set ${\displaystyle R={R_{c}}\pm b^{*}\mathbf {\hat {j^{*}}} }$.

on an earth section path may be found by dropping the subscripts on the given section; ${\displaystyle {\hat {N_{1}}}=(l,m,n)}$, ${\displaystyle d_{1}=d}$, and setting ${\displaystyle {\hat {N_{2}}}=(-\sin \theta ,\cos \theta ,0)}$, where ${\displaystyle \theta }$ is the longitude to be solved for such that ${\displaystyle B^{2}=AC}$.

Alternatively, just set ${\displaystyle R={R_{c}}\pm a^{*}\mathbf {\hat {i^{*}}} }$.