Earth section paths

Earth section paths are plane curves defined by the intersection of an earth ellipsoid and a plane (ellipsoid plane sections). Common examples include the great ellipse (containing the center of the ellipsoid) and normal sections (containing an ellipsoid normal direction). Earth section paths are useful as approximate solutions for geodetic problems, the direct and inverse calculation of geographic distances. The rigorous solution of geodetic problems involves skew curves known as geodesics.

Plane section of an ellipsoid

Inverse problem

The inverse problem for earth sections is: given two points, ${\displaystyle P_{1}}$ and ${\displaystyle P_{2}}$ on the surface of the reference ellipsoid, find the length, ${\displaystyle s_{12}}$, of the short arc of a spheroid section from ${\displaystyle P_{1}}$ to ${\displaystyle P_{2}}$ and also find the departure and arrival azimuths (angle from true north) of that curve, ${\displaystyle \alpha _{1}}$ and ${\displaystyle \alpha _{2}}$. The figure to the right illustrates the notation used here. Let ${\displaystyle P_{k}}$ have geodetic latitude ${\displaystyle \phi _{k}}$ and longitude ${\displaystyle \lambda _{k}}$ (k=1,2). This problem is best solved using analytic geometry in earth-centered, earth-fixed (ECEF) Cartesian coordinates. Let ${\displaystyle R_{1}=\mathrm {ECEF} (P_{1})}$ and ${\displaystyle R_{2}=\mathrm {ECEF} (P_{2})}$ be the ECEF coordinates of the two points, computed using the geodetic to ECEF transformation discussed here.

This illustrates the notation used for the geodetic problems discussed here.

Section plane

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 \mathbf {\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}}$ is ${\displaystyle \mathbf {V_{0}} =\mathbf {R_{1}} -\mathbf {R_{0}} }$, and the vector from ${\displaystyle R_{1}}$ to ${\displaystyle R_{2}}$ is ${\displaystyle \mathbf {V_{1}} =\mathbf {R_{2}} -\mathbf {R_{1}} }$. Therefore, ${\displaystyle \mathbf {\hat {N}} =\mathrm {unit} (\mathbf {V_{0}} \times \mathbf {V_{1}} )}$), where ${\displaystyle \mathrm {unit} (\mathbf {V} )}$ is the unit vector in the direction of ${\displaystyle \mathbf {V} }$. The orientation convention used here is that ${\displaystyle \mathbf {\hat {N}} }$ points to the left of the path. If this is not the case then redefine ${\displaystyle \mathbf {V_{0}} =-\mathbf {V_{0}} }$. Finally, the parameter d for the plane may be computed using the dot product of ${\displaystyle \mathbf {\hat {N}} }$ with a vector from the origin to any point on the plane, such as ${\displaystyle R_{1}}$, i.e. ${\displaystyle d=\mathbf {\hat {N}} \cdot \mathbf {R_{1}} }$. The equation of the plane (in vector form) is thus ${\displaystyle \mathbf {\hat {N}} \cdot \mathbf {R} =d}$, where ${\displaystyle \mathbf {R} }$ is the position vector of ${\displaystyle (x,y,z)}$.

Azimuth

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}} =(-\sin \lambda ,\cos \lambda ,0)}$, a unit vector pointing north is ${\displaystyle \mathbf {\hat {n}} =(-\sin \phi \cos \lambda ,-\sin \phi \sin \lambda ,\cos \phi )}$, and a unit vector pointing up is ${\displaystyle \mathbf {\hat {u}} =(\cos \phi \cos \lambda ,\cos \phi \sin \lambda ,\sin \phi )}$. A vector tangent to the path is: ${\displaystyle \mathbf {t} =\mathbf {\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 =\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}}$.

Section ellipse

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 or numerical integration. 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 ${\displaystyle p=0}$ then the section is a horizontal circle of radius ${\textstyle a{\sqrt {1-{\frac {d^{2}}{b^{2}}}}}}$, which has no solution if ${\displaystyle |d|>b}$.

If ${\displaystyle p>0}$ then Gilbertson^[1] showed that the ECEF coordinates of the center of the ellipse is ${\textstyle {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 ${\textstyle a^{*}=a{\sqrt {1-{\frac {d^{2}}{C}}}}}$, in the direction ${\textstyle \mathbf {{\hat {i}}^{*}} =\left({\frac {m}{p}},{\frac {-l}{p}},0\right)}$, and the semi-minor axis is ${\textstyle b^{*}={\frac {b}{\sqrt {C}}}a^{*}}$, in the direction ${\textstyle \mathbf {{\hat {j}}^{*}} =\left({\frac {ln}{p}},{\frac {mn}{p}},-p\right)}$, which has no solution if ${\displaystyle |d|>{\sqrt {C}}}$.

Arc Length

The above referenced paper provides a derivation for an arc length formula involving the central angle and powers of ${\displaystyle e^{2}}$ to compute the arc length to millimeter accuracy, where ${\textstyle e^{2}=1-\left({\frac {b^{*}}{a^{*}}}\right)^{2}}$. That arc length formula may be rearranged and put into the form: ${\displaystyle s_{12}=s(\theta _{2})-s(\theta _{1})}$, where ${\displaystyle s(\theta )=b^{*}({C_{0}}\theta +{C_{2}}\sin(2\theta )+{C_{4}}\sin(4\theta )+{C_{6}}\sin(6\theta ))}$ and the coefficients are

${\displaystyle C_{0}=1.0+e^{2}(1/4+13e^{2}/64+45e^{4}/256+2577e^{6}/16384)}$

${\displaystyle C_{2}=e^{2}(1/8+3e^{2}/32+95e^{4}/1024+385e^{6}/4096)}$

${\displaystyle C_{4}=-e^{4}(1/256+5e^{2}/1024+19e^{4}/16384)}$

${\displaystyle C_{6}=-e^{6}(15/3072+35e^{2}/4096)}$

To compute the central angle, let ${\displaystyle P}$ be any point on the section ellipse and ${\displaystyle R=\mathrm {ECEF} (P)}$. Then ${\displaystyle \mathbf {V} =\mathbf {R} -\mathbf {R_{c}} }$ is a vector from the center of the ellipse to the point. The central angle ${\displaystyle \theta }$ is the angle from the semi-major axis to ${\displaystyle \mathbf {V} }$. Letting ${\displaystyle \mathbf {\hat {V}} =\mathrm {unit} (\mathbf {V} )}$, we have ${\displaystyle \theta =\operatorname {atan2} (\mathbf {\hat {V}} \cdot \mathbf {{\hat {j}}^{*}} ,\mathbf {\hat {V}} \cdot \mathbf {{\hat {i}}^{*}} )}$. In this way we obtain ${\displaystyle \theta _{1}}$ and ${\displaystyle \theta _{2}}$.

On the other hand it's possible to use Meridian arc formulas in the more general case provided that the section ellipse parameters are used rather than the spheroid parameters. One such rapidly convergent series is given in Series in terms of the parametric latitude. If we use ${\displaystyle \varepsilon }$ to denote the spheroid eccentricity, i.e. ${\textstyle \varepsilon ^{2}=1-\left({\frac {b}{a}}\right)^{2}}$, then ${\displaystyle e^{8}}$ ≤ ${\displaystyle \varepsilon ^{8}}$ ≅ 1.8×10⁻⁹. Similarly the third flattening of the section ellipse is bounded by the corresponding value for the spheroid, and for the spheroid we have ${\displaystyle n^{3}}$ ≅ 4.4×10⁻⁹, and ${\displaystyle n^{4}}$ ≅ 7.3×10⁻¹². Therefore it may suffice to ignore terms beyond ${\displaystyle B_{6}}$ in the parametric latitude series. To apply ${\textstyle s(\beta )={\frac {a^{*}+b^{*}}{2}}(B_{0}\beta +B_{2}\sin 2\beta +B_{4}\sin 4\beta +B_{6}\sin 6\beta )}$ in the current context requires converting the central angle to the parametric angle using ${\displaystyle \beta =\tan ^{-1}\left(\tan \theta /((1-f)\right)}$, and using the section ellipse third flattening. Whichever method is used, care must be taken when using ${\displaystyle \theta _{1}}$ & ${\displaystyle \theta _{2}}$ or ${\displaystyle \beta _{1}}$ & ${\displaystyle \beta _{2}}$ to ensure that the shorter arc connecting the 2 points is used.

Direct problem

The direct problem is given ${\displaystyle {P_{1}}}$, the distance ${\displaystyle s_{12}}$, and departure azimuth ${\displaystyle \alpha _{1}}$, find ${\displaystyle {P_{2}}}$ and the arrival azimuth ${\displaystyle \alpha _{2}}$.

Section plane

The answer to this problem depends the choice of ${\displaystyle \mathbf {V_{0}} }$. i.e. on the type of section. Observe that ${\displaystyle \mathbf {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). Having made such a choice, and considering orientation proceed as follows. 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}}}$. The normal vector ${\displaystyle \mathbf {\hat {N}} =\mathrm {unit} (\mathbf {V_{0}} \times \mathbf {\hat {t}} _{1}}$), together with ${\displaystyle \mathbf {P_{1}} }$ defines the plane. In other words, the tangent takes the place of the chord since the destination is unknown.

Locate arrival point

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. If the arc length, ${\displaystyle s_{12}}$ is given then the problem is to find the corresponding change in the central angle ${\displaystyle \theta _{12}}$, so that ${\displaystyle \theta _{2}=\theta _{1}+\theta _{12}}$ and the position can be calculated. Assuming that we have a series that gives ${\displaystyle s=s(\theta )}$ then what we seek now is ${\displaystyle \theta _{2}=s^{-1}(s_{1}+s_{12})}$. The inverse of the central angle arc length series above may be found on page 8a of Rapp, Vol. 1,^[2] who credits Ganshin.^[3] An alternative to using the inverse series is using Newton's method of successive approximations to ${\displaystyle \theta _{12}}$. The inverse meridian problem for the ellipsoid provides the inverse to Bessel's arc length series in terms of the parametric angle. Before the inverse series can be used, the parametric angle series must be used to compute the arc length from the semi-major axis to ${\displaystyle P_{1}}$, ${\textstyle s_{1}=s(\beta _{1})={\frac {a^{*}+b^{*}}{2}}(B_{0}\beta _{1}+B_{2}\sin 2\beta _{1}+B_{4}\sin 4\beta _{1}+B_{6}\sin 6\beta _{1})}$. Once ${\displaystyle s_{1}}$ is known apply the inverse formula to obtain${\displaystyle \beta _{2}=\beta (s_{1}+s_{12})=\mu _{2}+B'_{2}\sin 2\mu _{2}+B'_{4}\sin 4\mu _{2}+B'_{6}\sin 6\mu _{2}}$, where ${\displaystyle \mu _{2}=2(s_{1}+s_{12})/(B_{0}(a^{*}+b^{*}))}$. Rectangular coordinates in the section plane are ${\displaystyle x_{2}=a^{*}\cos \beta _{2},y_{2}=b^{*}\sin \beta _{2}}$. So an ECEF vector may be computed using ${\displaystyle \mathbf {V_{2}} =\mathbf {R_{c}} +(x_{2}\mathbf {{\hat {i}}^{*}} +y_{2}\mathbf {{\hat {j}}^{*}} )}$. Finally calculate geographic coordinates via ${\displaystyle P_{2}=\mathrm {Geo} (V_{2})}$ using Bowring's 1985 algorithm,^[4] or the algorithm here.

Azimuth

Azimuth may be obtained by the same method as the indirect problem: ${\displaystyle \mathbf {t_{2}} =\mathbf {\hat {N}} \times \mathbf {{\hat {u}}_{2}} }$ and ${\displaystyle {\alpha _{2}}=\operatorname {atan2} (\mathbf {t_{2}} \cdot \mathbf {{\hat {e}}_{2}} ,\mathbf {t_{2}} \cdot \mathbf {{\hat {n}}_{2}} )}$.

Examples

Shows the geodesic deviation for various sections connecting New York to Paris

The great ellipse

Main article: Great ellipse

The great ellipse is the curve formed by intersecting the ellipsoid with a plane through its center. Therefore, to use the method above, just let ${\displaystyle R_{0}}$ be the origin, so that ${\displaystyle \mathbf {V_{0}} =\mathbf {R_{1}} }$ (the position vector of ${\displaystyle R_{1}}$). This method avoids the esoteric and sometimes ambiguous formulas of spherical trigonometry, and provides an alternative to the formulas of Bowring.^[5] The shortest path between two points on a spheroid is known as a geodesic. Such paths are developed using differential geometry. The equator and meridians are great ellipses that are also geodesics^[a]. The maximum difference in length between a great ellipse and the corresponding geodesic of length 5,000 nautical miles is about 10.5 meters. The lateral deviation between them may be as large as 3.7 nautical miles. A normal section connecting the two points will be closer to the geodesic than the great ellipse, unless the path touches the equator.

On the WGS84 ellipsoid, the results for the great elliptic arc from New York, ${\displaystyle \phi _{1}}$ = 40.64130°, ${\displaystyle \lambda _{1}}$ = -73.77810° to Paris, ${\displaystyle \phi _{2}}$ = 49.00970°, ${\displaystyle \lambda _{2}}$= 2.54800° are:

${\displaystyle \alpha _{1}}$ = 53.596810°, ${\displaystyle \alpha _{2}}$ = 111.537138° and ${\displaystyle s_{12}}$ = 5849159.753 (m) = 3158.293603 (nm). The corresponding numbers for the geodesic are:

${\displaystyle \alpha _{1}}$ = 53.511007°, ${\displaystyle \alpha _{2}}$ = 111.626714° and ${\displaystyle s_{12}}$ = 5849157.543 (m) = 3158.292410 (nm).

To illustrate the dependence on section type for the direct problem, let the departure azimuth and trip distance be those of the geodesic above, and use the great ellipse to define the direct problem. In this case the arrival point is ${\displaystyle \phi _{2}}$ = 49.073057°, ${\displaystyle \lambda _{2}}$= 2.586154°, which is about 4.1 nm from the arrival point in Paris defined above. Of course using the departure azimuth and distance from the great ellipse indirect problem will properly locate the destination, ${\displaystyle \phi _{2}}$ = 49.00970°, ${\displaystyle \lambda _{2}}$= 2.54800°, and the arrival azimuth ${\displaystyle \alpha _{2}}$ = 111.537138°.

Shows the geodesic deviation for various sections connecting Sydney to Bangkok

Normal sections

A normal section at ${\displaystyle P_{1}}$ is determined by letting ${\displaystyle \mathbf {V_{0}} =\mathbf {\hat {u}} _{1}}$ (the surface normal at ${\displaystyle P_{1}}$). Another normal section, known as the reciprocal normal section, results from using the surface normal at ${\displaystyle P_{2}}$. Unless the two points are both on the same parallel or the same meridian, the reciprocal normal section will be a different path than the normal section. The above approach provides an alternative to that of others, such as Bowring.^[7] The importance of normal sections in surveying as well as a discussion of the meaning of the term line in such a context is given in the paper by Deakin, Sheppard and Ross.^[8]

On the WGS84 ellipsoid, the results for the normal section from New York, ${\displaystyle \phi _{1}}$ = 40.64130°, ${\displaystyle \lambda _{1}}$ = -73.77810° to Paris, ${\displaystyle \phi _{2}}$ = 49.00970°, ${\displaystyle \lambda _{2}}$= 2.54800° are:

${\displaystyle \alpha _{1}}$ = 53.521396°, ${\displaystyle \alpha _{2}}$ = 111.612516° and ${\displaystyle s_{12}}$ = 5849157.595 (m) = 3158.292438 (nm). The results for the reciprocal normal section from New York to Paris are:

${\displaystyle \alpha _{1}}$ = 53.509422°, ${\displaystyle \alpha _{2}}$ = 111.624483° and ${\displaystyle s_{12}}$ = 5849157.545 (m) = 3158.292411 (nm).

The maximum difference in length between a normal section and the corresponding geodesic of length 5,000 nautical miles is about 6.0 meters. The lateral deviation between them may be as large as 2.8 nautical miles.

To illustrate the dependence on section type for the direct problem, let the departure azimuth and trip distance be those of the geodesic above, and use the surface normal at NY to define the direct problem. In this case the arrival point is ${\displaystyle \phi _{2}}$ = 49.017378°, ${\displaystyle \lambda _{2}}$= 2.552626°, which is about 1/2 nm from the arrival point defined above. Of course, using the departure azimuth and distance from the normal section indirect problem will properly locate the destination in Paris. Presumably the direct problem is used when the arrival point is unknown, yet it is possible to use whatever vector ${\displaystyle \mathbf {V_{0}} }$ one pleases. For example, using the surface normal at Paris, ${\displaystyle \mathbf {\hat {u}} _{2}}$, results in an arrival point of ${\displaystyle \phi _{2}}$ = 49.007778°, ${\displaystyle \lambda _{2}}$= 2.546842°, which is about 1/8 nm from the arrival point defined above. Using the surface normal at Reykjavik (while still using the departure azimuth and trip distance of the geodesic to Paris) will have you arriving about 347 nm from Paris, while the normal at Zürich brings you to within 5.5 nm.

The search for a section that's closer to the geodesic led to the next two examples.

Shows how the geodesic deviation varies with azimuth for sections originating at 20° latitude.

The mean normal section

The mean normal section from ${\displaystyle P_{1}}$ to ${\displaystyle P_{2}}$ is determined by letting ${\displaystyle \mathbf {V_{0}} =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.The maximum difference in length between the mean normal section and the corresponding geodesic of length 5,000 nautical miles is about 0.5 meters. The lateral deviation between them is no more than about 0.8 nautical miles. For paths of length 1000 nautical miles the length error is less than a millimeter, and the worst case lateral deviation is about 4.4 meters. Continuing the example from New York to Paris on WGS84 gives the following results for the mean normal section:

${\displaystyle \alpha _{1}}$ = 53.515409°, ${\displaystyle \alpha _{2}}$ = 111.618500° and ${\displaystyle s_{12}}$ = 5849157.560 (m) = 3158.292419 (nm).

Shows the geodesic deviation for various 5000nm normal sections from the equator.

The midpoint normal section

The midpoint normal section from ${\displaystyle P_{1}}$ to ${\displaystyle P_{2}}$ is determined by letting ${\displaystyle \mathbf {V_{0}} }$ = the surface normal at the midpoint of the geodesic from ${\displaystyle P_{1}}$ to ${\displaystyle P_{2}}$. This path is only slightly closer to the geodesic that the mean normal section. The maximum difference in length between a midpoint normal section and the corresponding geodesic of length 5,000 nautical miles is about 0.3 meters. The worst case lateral deviation between them is about 0.3 nautical miles.

Finishing the example from New York to Paris on WGS84 gives the following results for the geodesic midpoint normal section: ${\displaystyle \alpha _{1}}$ = 53.506207°, ${\displaystyle \alpha _{2}}$ = 111.627697° and ${\displaystyle s_{12}}$ = 5849157.545 (m) = 3158.292411 (nm).

Discussion

All of the section paths used in the charts to the right were defined using the indirect method above. In the third and fourth charts the terminal point was defined using the direct algorithm for the geodesic with the given distance and initial azimuth. On each of the geodesics some points were selected, the nearest point on the section plane was located by vector projection, and the distance between the two points computed. This distance is described as the lateral deviation from the geodesic, or briefly geodesic deviation, and is displayed in the charts on the right. The alternative of finding the corresponding point on the section path and computing geodesic distances would produce slightly different results.

The first chart is typical of mid-latitude cases where the great ellipse is the outlier. The normal section associated with the point farthest from the equator is a good choice for these cases.

The second example is longer and is typical of equator crossing cases, where the great ellipse beats the normal sections. However, the two normal sections deviate on opposite sides of the geodesic, making the mean normal section a good choice here.

The third chart shows how the geodesic deviations vary with initial geodesic azimuth originating from 20 degrees north latitude. The worst case deviation for normal sections of 5000 nautical miles length is about 2.8 nm and occurs at initial geodesic azimuth of 132° from 18° north latitude (48° azimuth for south latitude).

The fourth chart is what the third chart looks like when departing from the equator. On the equator there are more symmetries since sections at 90° and 270° azimuths are also geodesics. Consequently the fourth chart shows only 7 distinct lines out of the 24 with 15 degree spacing. Specifically, the lines at azimuths 15, 75, 195 and 255 coincide, as do the lines at 105, 165, 285, and 345 on the other side as the inner most (other than the geodesics). Next farthest coincident lines from the four geodesic lines are at azimuths 30, 60, 210, and 240 on one side and 120, 150, 300, and 330 on the other side. The outer most lines are at azimuths 45, and 225 on one side and 135 and 315 on the other. As the departure point moves north the lines at azimuths 90 and 270 are no longer geodesics, and other coincident lines separate and fan out until 18° latitude where the maximum deviation is attained. Beyond this point the deviations contract like a Japanese fan as the initial point proceeds north. So that by 84° latitude the maximum deviation for normal sections is about 0.25 nm.

The midpoint normal section is (almost) always a good choice.

Intersections

Let two section planes be given: ${\displaystyle \mathbf {\hat {N}} _{1}\cdot \mathbf {R} =d_{1}}$, and ${\displaystyle \mathbf {\hat {N}} _{2}\cdot \mathbf {R} =d_{2}}$. Assuming that the two planes are not parallel, the line of intersection is on both planes. Hence orthogonal to both normals, i.e. in the direction of ${\displaystyle \mathbf {N_{3}} =\mathbf {\hat {N}} _{1}\times \mathbf {\hat {N}} _{2}}$ (there is no reason to normalize ${\displaystyle \mathbf {N_{3}} }$).

Since ${\displaystyle \mathbf {\hat {N}} _{1}}$ and ${\displaystyle \mathbf {\hat {N}} _{2}}$ are not collinear ${\displaystyle \mathbf {\hat {N}} _{1}}$, ${\displaystyle \mathbf {\hat {N}} _{2}}$, ${\displaystyle \mathbf {N_{3}} }$ is a basis for ${\displaystyle \mathbb {R} ^{3}}$. Therefore, there exist constants ${\displaystyle C_{1}}$ and ${\displaystyle C_{2}}$ such that the line of intersection of the 2 planes is given by ${\displaystyle R=C_{1}\mathbf {\hat {N}} _{1}+C_{2}\mathbf {\hat {N}} _{2}+t\mathbf {N_{3}} }$, where t is an independent parameter.

Since this line is on both section planes, it satisfies both: ${\displaystyle C_{1}+C_{2}(\mathbf {\hat {N}} _{1}\cdot \mathbf {\hat {N}} _{2})=d_{1}}$, and ${\displaystyle C_{1}(\mathbf {\hat {N}} _{1}\cdot \mathbf {\hat {N}} _{2})+C_{2}=d_{2}}$.

Solving these equations for ${\displaystyle {C_{1}}}$ and ${\displaystyle {C_{2}}}$ gives ${\displaystyle C_{1}[1-(\mathbf {\hat {N}} _{1}\cdot \mathbf {\hat {N}} _{2})^{2}]=d_{1}-d_{2}(\mathbf {\hat {N}} _{1}\cdot \mathbf {\hat {N}} _{2})}$, and ${\displaystyle C_{2}[1-(\mathbf {\hat {N}} _{1}\cdot \mathbf {\hat {N}} _{2})^{2}]=d_{2}-d_{1}(\mathbf {\hat {N}} _{1}\cdot \mathbf {\hat {N}} _{2})}$.

Define the "dihedral angle", ${\displaystyle \nu }$, by ${\displaystyle \cos \nu ={\mathbf {\hat {N}} _{1}}\cdot {\mathbf {\hat {N}} _{2}}}$. Then ${\displaystyle C_{1}={\frac {(d_{1}-d_{2}\cos \nu )}{\sin ^{2}\nu }}}$ , and ${\displaystyle C_{2}={\frac {(d_{2}-d_{1}\cos \nu )}{\sin ^{2}\nu }}}$.

On the intersection line we have ${\displaystyle \mathbf {R} =\mathbf {R_{0}} +t\mathbf {N_{3}} }$, where ${\displaystyle \mathbf {R_{0}} =C_{1}\mathbf {\hat {N}} _{1}+C_{2}\mathbf {\hat {N}} _{2}}$. Hence: ${\displaystyle x=x_{0}+tl_{3}}$, ${\displaystyle y=y_{0}+tm_{3}}$, and ${\displaystyle z=z_{0}+tn_{3}}$, where ${\displaystyle x_{0}=C_{1}l_{1}+C_{2}l_{2}}$, ${\displaystyle y_{0}=C_{1}m_{1}+C_{2}m_{2}}$, and ${\displaystyle z_{0}=C_{1}n_{1}+C_{2}n_{2}}$, ${\displaystyle \mathbf {\hat {N}} _{i}=(l_{i},m_{i},n_{i})}$, for i=1,2, and ${\displaystyle \mathbf {N_{3}} =(l_{3},m_{3},n_{3})}$.

To find the intersection of this line with the earth, plug the line equations into ${\displaystyle {\frac {x^{2}}{a^{2}}}+{\frac {y^{2}}{a^{2}}}+{\frac {z^{2}}{b^{2}}}=1}$, to get ${\displaystyle At^{2}+2Bt+C=0}$, where ${\displaystyle A=l_{3}^{2}+m_{3}^{2}+{\frac {a^{2}}{b^{2}}}n_{3}^{2}}$, ${\displaystyle B=x_{0}l_{3}+y_{0}m_{3}+{\frac {a^{2}}{b^{2}}}z_{0}n_{3}}$, ${\displaystyle C=x_{0}^{2}+y_{0}^{2}+{\frac {a^{2}}{b^{2}}}z_{0}^{2}-a^{2}}$.

Therefore, the line intersects the earth at ${\displaystyle t={\frac {-B\pm {\sqrt {{B}^{2}-AC}}}{A}}}$. If ${\displaystyle B^{2}<AC}$, then there is no intersection. If ${\displaystyle B^{2}=AC}$, then the line is tangent to the earth at ${\displaystyle t=-B/A}$ (i.e. the sections intersect at that single point).

Observe that ${\displaystyle A\neq 0}$ since ${\displaystyle \mathbf {\hat {N}} _{1}}$ and ${\displaystyle \mathbf {\hat {N}} _{2}}$ are not collinear. Plugging t into ${\displaystyle \mathbf {R} =\mathbf {R_{0}} +t\mathbf {N_{3}} }$, gives the points of intersection of the earth sections.

Example

Find where a section from New York to Paris, intersects the Greenwich meridian. The plane of the prime meridian may be described by ${\displaystyle \mathbf {\hat {N}} =(0,1,0)}$ and ${\displaystyle d=0}$. The results are as follows:

Intersections

Section	Latitude
Great Ellipse	49.634970°
Normal	49.637377°
Mean Normal	49.637568°
Reciprocal	49.637759°
Midpoint	49.637862°

Extreme latitudes and longitudes

The maximum (or minimum) latitude is where the section ellipse intersections a parallel at a single point. To set up the problem, let ${\displaystyle \mathbf {{\hat {N}}_{1}} =(l,m,n)}$, ${\displaystyle d_{1}=d}$ be the given section plane. The parallel is ${\displaystyle \mathbf {{\hat {N}}_{2}} =(0,0,1)}$, ${\displaystyle d_{2}=z_{0}}$, where ${\displaystyle z_{0}}$ is to be determined so that there is only one intersection point. Applying the intersection method above results in ${\displaystyle \mathbf {N_{3}} =\mathbf {\hat {N}} _{1}\times \mathbf {\hat {N}} _{2}=(m,-l,0)}$, ${\displaystyle {\mathbf {\hat {N}} _{1}}\cdot {\mathbf {\hat {N}} _{2}}=n}$, ${\textstyle C_{1}={\frac {1}{p^{2}}}(d-nz_{0})}$, and ${\textstyle C_{2}={\frac {1}{p^{2}}}(z_{0}-nd)}$, since ${\displaystyle 1-n^{2}=l^{2}+m^{2}=p^{2}}$. The resulting linear equations become ${\displaystyle x=x_{0}+tm}$, ${\displaystyle y=y_{0}-tl}$, and ${\displaystyle z=z_{0}}$, where ${\displaystyle x_{0}=C_{1}l}$, ${\displaystyle y_{0}=C_{1}m}$, and ${\displaystyle z_{0}}$ is to be determined. The resulting quadratic coefficients are ${\displaystyle A=m^{2}+l^{2}=p^{2}}$, ${\displaystyle B=mx_{0}-ly_{0}=lmC_{1}-lmC_{1}=0}$, ${\displaystyle C=x_{0}^{2}+y_{0}^{2}+{\frac {a^{2}}{b^{2}}}z_{0}^{2}-a^{2}=p^{2}C_{1}^{2}+{\frac {a^{2}}{b^{2}}}z_{0}^{2}-a^{2}={\frac {1}{p^{2}}}(d-nz_{0})^{2}+{\frac {a^{2}}{b^{2}}}z_{0}^{2}-a^{2}}$. Therefore the intersection will result in only one solution if ${\displaystyle B^{2}=AC}$, but since ${\displaystyle B=0}$ and ${\displaystyle A>0}$^[b], the critical equation becomes ${\displaystyle C=0}$. This equation may be rearranged and put into the form ${\displaystyle Ez_{0}^{2}-2Fz_{0}+G=0}$, where ${\displaystyle E={\frac {a^{2}}{b^{2}}}p^{2}+n^{2}}$, ${\displaystyle F=nd}$, and ${\displaystyle G=d^{2}-a^{2}p^{2}}$. Therefore, ${\textstyle z_{0}={\frac {F\pm {\sqrt {{F}^{2}-EG}}}{E}}}$ provides the distance from the origin of the desired parallel planes. Plugging ${\displaystyle z_{0}}$ into ${\displaystyle C_{1}}$ gives the values for ${\displaystyle x_{0}}$ and ${\displaystyle y_{0}}$. Recall that ${\displaystyle t=-B/A=0}$ so ${\displaystyle x=x_{0}}$, ${\displaystyle y=y_{0}}$ are the remaining coordinates of the intersections. The geographic coordinates may then be computed using the ECEF_to_Geo conversion.

The same method may be applied to meridians to find extreme longitudes, but the results are not easy to interpret due to the modular nature of longitude. However, the results can always be verified using the following approach.

The simpler approach is to compute the end points of the minor and major axes of the section ellipse using ${\displaystyle \mathbf {R} =\mathbf {R_{c}} \pm b^{*}\mathbf {\hat {j}} ^{*}}$, and ${\displaystyle \mathbf {R} =\mathbf {R_{c}} \pm a^{*}\mathbf {\hat {i}} ^{*}}$, and then converting to geographic coordinates. It may be worth mentioning here that the line of intersection of two planes consists of the set of fixed points, hence the rotation axis, of a coordinate rotation that maps one plane onto the other.

For the New York to Paris example the results are:

Section	Minor Axis Point 1	Minor Axis Point 2	Major Axis Point 1	Major Axis Point 2
Great Ellipse	${\displaystyle \phi _{1}}$ = 52.418061°, ${\displaystyle \lambda _{1}}$ = -25.123079°	${\displaystyle \phi _{2}}$ = -52.418061°, ${\displaystyle \lambda _{2}}$ = 154.876921°	${\displaystyle \phi _{1}}$ = 0.000000°, ${\displaystyle \lambda _{1}}$ = 64.876921°	${\displaystyle \phi _{2}}$ = 0.000000°, ${\displaystyle \lambda _{2}}$ = -115.123079°
Normal	${\displaystyle \phi _{1}}$ = 52.433790°, ${\displaystyle \lambda _{1}}$ = -25.154863°	${\displaystyle \phi _{2}}$ = -52.739188°, ${\displaystyle \lambda _{2}}$ = 154.845137°	${\displaystyle \phi _{1}}$ = -0.093365°, ${\displaystyle \lambda _{1}}$ = 64.723898°	${\displaystyle \phi _{2}}$ = -0.093365°, ${\displaystyle \lambda _{2}}$ = -115.033623°
Mean Normal	${\displaystyle \phi _{1}}$ = 52.435039°, ${\displaystyle \lambda _{1}}$ = -25.157380°	${\displaystyle \phi _{2}}$ = -52.764681°, ${\displaystyle \lambda _{2}}$ = 154.842620°	${\displaystyle \phi _{1}}$ = -0.100746°, ${\displaystyle \lambda _{1}}$ = 64.711732°	${\displaystyle \phi _{2}}$ = -0.100746°, ${\displaystyle \lambda _{2}}$ = -115.026491°
Reciprocal	${\displaystyle \phi _{1}}$ = 52.436288°, ${\displaystyle \lambda _{1}}$ = -25.159896°	${\displaystyle \phi _{2}}$ = -52.790172°, ${\displaystyle \lambda _{2}}$ = 154.840104°	${\displaystyle \phi _{1}}$ = -0.108122°, ${\displaystyle \lambda _{1}}$ = 64.699565°	${\displaystyle \phi _{2}}$ = -0.108122°, ${\displaystyle \lambda _{2}}$ = -115.019357°
Midpoint	${\displaystyle \phi _{1}}$ = 52.436959°, ${\displaystyle \lambda _{1}}$ = -25.161247°	${\displaystyle \phi _{2}}$ = -52.803863°, ${\displaystyle \lambda _{2}}$ = 154.838753°	${\displaystyle \phi _{1}}$ = -0.112082°, ${\displaystyle \lambda _{1}}$ = 64.693029°	${\displaystyle \phi _{2}}$ = -0.112082°, ${\displaystyle \lambda _{2}}$ = -115.015522°

See also

• Geographical distance#Altitude correction
• Normal-section azimuth

Notes

1. Equatorial paths are geodesics up to a point. For example the geodesic connecting two points that are 180° apart on the equator is a meridian path over a pole, whereas the equator is still a great ellipse. In fact there are infinitely many great ellipses in this case, only two of which are geodesics. For short arcs the geodesic and great ellipse coincide. So at what point does it change? Rapp calculates the answer to be 179° 23' 38.18182".^[6] At that point the geodesic starts to move off of the equator, and by 180° it is all the way to a pole.
2. Otherwise the section is a parallel, so there is nothing to solve, since all latitudes are the same.

References

1. Gilbertson, Charles (Spring 2012). "Earth Section Paths". Navigation. 59 (1): 1–7. doi:10.1002/navi.2.
2. Rapp, R. H. (1991), Geometric geodesy, part I, Ohio State Univ., hdl:1811/24333
3. Gan'shin, V. V. (1969) [1967].Geometry of the Earth Ellipsoid. Translated by Willis, J. M. St. Louis: Aeronautical Chart and Information Center. doi:10.5281/zenodo.32854. OCLC 493553. Translation from Russian of Геометрия земного эллипсоида (Moscow, 1967)
4. Bowring, B.R. (1985). "The accuracy of geodetic latitude and height equations". Survey Review. 28 (218): 202–206. Bibcode:1985SurRv..28..202B. doi:10.1179/sre.1985.28.218.202.
5. Bowring, B.R. (1984). "The Direct and Inverse Solutions for the Great Elliptic Line on the Reference Ellipsoid". Bulletin Géodésique. 58 (1): 101–108. Bibcode:1984BGeod..58..101B. doi:10.1007/BF02521760. S2CID 123161737.
6. Rapp, R. H. (1993), Geometric geodesy, part II, Ohio State Univ., hdl:1811/24409
7. Bowring, B.R. (1971). "The normal section -- forward and inverse formulae at any distance". Survey Review. XXI (161): 131–136. Bibcode:1971SurRv..21..131B. doi:10.1179/sre.1971.21.161.131.
8. Deakin, R. E.; Sheppard, S. W.; Ross, R. (2011). "The Black-Allan Line Revisited" (PDF). 24th Victorian Regional Survey Conference, Shepparton, 1–3 April 2011. Archived from the original (PDF) on 5 January 2012. Retrieved 3 February 2012.

Further reading

• Helmert, Friedrich Robert (1964-01-01). "Mathematical and Physical Theories of Higher Geodesy, Part 1, Preface and the Mathematical Theories". Zenodo. doi:10.5281/zenodo.32050. Retrieved 2022-04-17.
• Jordan, Wilhelm; Eggert, Otto (1962-01-01). "Jordan's Handbook of Geodesy, Vol. 3, 2nd half". Zenodo. doi:10.5281/zenodo.35316. Retrieved 2022-04-17.