Geodetic Line at Constant Altitude above the Ellipsoid

Geodetic Line at Constant Altitude above the Ellipsoid

Richard J. Mathar Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands
July 29, 2019

The two-dimensional surface of a bi-axial ellipsoid is characterized by the lengths of its major and minor axes. Longitude and latitude span an angular coordinate system across. We consider the egg-shaped surface of constant altitude above (or below) the ellipsoid surface, and compute the geodetic lines—lines of minimum Euclidean length—within this surface which connect two points of fixed coordinates. This addresses the common “inverse” problem of geodesics generalized to non-zero elevations. The system of differential equations which couples the two angular coordinates along the trajectory is reduced to a single integral, which is handled by Taylor expansion up to fourth power in the eccentricity.

geodesy; ellipsoid; eccentricity; geodetic line; inverse problem
91.10.By, 91.10.Ws, 02.40.Hw
thanks: This work is supported by the NWO VICI grant 639.043.201 “Optical Interferometry: A new Method for Studies of Extrasolar Planets” to A. Quirrenbach.

I Contents

The common three parameters employed to relate Cartesian coordinates to an ellipsoidal surface are the angles of latitude and longitude in a grid on the surface, plus an altitude which is a shortest (perpendicular) distance to the surface. The well-known functional relations (coordinate transformations) are summarized in Section II.

The inverse problem of geodesy is to find the line embedded in the ellipsoid surface which connects two fixed points subject to minimization of its length. We pose the equivalent problem for lines at constant altitude, as if one would ask for the shortest track of the center of a sphere of given radius which rolls on the ellipsoid surface and meets two points of the same, known altitude. In Section III, we formulate this in terms of the generic differential equations of geodesy parametrized by the Christoffel symbols.

In the main Section IV, the coupled system of differential equations of latitude and longitude as a function of path length is reduced to one degree of freedom, here chosen to be the direction along the path at one of the fixed terminal points, measured in the topocentric horizontal system, and dubbed the launching angle. Closed-form expressions of this parameter in terms of the coordinates of the terminal points have not been found; instead, the results are presented as series expansions up to fourth power in the eccentricity.

The standard treatment of this analysis is the projection on an auxiliary sphere; this technique is (almost) completely ignored but for the pragmatic aspect that the case of zero eccentricity is a suitable zeroth-order reference of series expansions around small eccentricities.

Ii Spheroidal Coordinates

ii.1 Surface

The cross section of an ellipsis of equatorial radius with eccentricity in a Cartesian () system is:


The ellipsis defines a geocentric latitude and a geodetic latitude , the latter measured by intersection of the normal to the tangential plane with the equatorial plane (Figure 1).

Figure 1: Major semi-axis , minor semi-axis , straight distance to the center of coordinates , geocentric latitude , geodetic latitude , and their pointing difference .

On the surface of the ellipsoid (Smart, 1949, §IX)Dufour (1958)(Bouasse, 1919, §140):


Supposed denotes the distance from the center of coordinates to the point on the geoid surface, the transformation from to is


defines the difference between the geodetic (astronomical) and the geocentric latitudes,


where the expansion of the denominator has been given in terms of the geometric series of the parameter


The Taylor series in powers of and in powers of are


A flattening factor is also commonly defined Kaplan (2006),


The reference values of the Earth ellipsoid adopted in the WGS84 National Imagery and Mapping Agency (2000) are


The numerical evaluation of (6) in terms of a Fourier series with this parametrization is, in units of radians and arcseconds:


Slightly different values would emerge according to IERS conventions of 2003 McCarthy and Petit (2003).

ii.2 General Altitude

If one moves along the direction of a distance away from the surface of the geoid, the new coordinates relative to (3) are Fukushima (2006); Vermeille (2004); Jones (2004); Pollard (2002); Zhang et al. (2005); Hradilek (1976)


which can be written in terms of a distance ,



Figure 2: A point with Cartesian coordinates has a distance to the equatorial plane, a distance to the polar axis, a distance to the surface of the ellipsoid, and a distance to the polar axis, measured along the local normal to the surface Long (1975). The vector points North at this point.

Rotation of (14) around the polar axis with geographic longitude defines the full 3D transformation between and ,


The only theme of this paper is to generalize the geodetic lines of the literature Dufour (1958); Sodano (1958); Saito (1970); Thomas (1965); Ölander (1952); Bowring (1983) to the case of finite altitude . The physics of gravimetric or potential theory is not involved, only the mathematics of the geometry. It should be noted that points with constant, non-zero do not define a surface of an ellipsoid with effective semi-axes —otherwise the geodetic line could be deduced by mapping the problem onto an equivalent ellipsoidal surface Dufour (1958).

Iii Inverse Problem of Geodesy

iii.1 Topocentric Coordinate System

A line of shortest distance at constant height const between two points 1 and 2 is defined by minimizing the Euclidean distance, the line integral


for some parametrization . The integrand is equivalent to a Lagrange function or .

The gradient of of (15) with respect to and defines the vectors that span the topocentric tangent plane


where a meridional radius of curvature Tienstra (1951)


is defined to simplify the notation. Building squares and dot products computes the three Gauss Fundamental parameters of the surface Reichardt (1957)


Specializing to we get the You formulae You (1999). and provide the principal curvatures along the meridian and azimuth Dorrer (1999), and the coefficients of the metric tensor in the quadratic form of ,


iii.2 Christoffel Symbols

Christoffel symbols are the connection coefficients between differentials of the topocentric axis and of the positions , in a generic definition


This format is matched by first computing the derivative of with respect to and (at ):


and of with respect to and :


The next step splits these two equations to the expanded version of (23),


The eight are extracted by evaluating dot products of the four vector coefficients in (24)–(25) by and ,


The Euler-Lagrange Differential Equations for a stationary Lagrange density (at ) become the differential equations of the geodesic Reichardt (1957), in the generic format


The explicit write-up


simplifies with (28) to


Iv Reduction of the Differential Equations

iv.1 Separation of Angular Variables

Decoupling of the two differential equations (35)–(36) starts with the separation of variables in (35),


Change of the integration variable on the right hand side from to allows to use the underivative of (30)


to generate a first integral


Exponentiation yields




is a constant for each geodesic. Is has been defined with the azimuth and at the start of the line, but could as well be associated with any other point or the end point . is positive for trajectories starting into eastwards direction, negative for the westwards heading, zero for routes to the poles. Moving the square of (40) into (36) yields


To solve this differential equation, we substitute the variable by its projection onto the polar axis,


which implies transformations in the derivatives:


We multiply (42) by , then replace and as noted above,


This is a differential equation with no explicit appearance of the independent variable ,

and the standard way of progressing is the substitution


This is transformed to a linear differential equation by the further substitution , ,


The standard approach is to solve the homogeneous differential equation first,


After division through , the left hand side is easily integrated, and the right hand side (incompletely) decomposed into partial fractions,


Solution of the inhomogeneous differential equation (50) proceeds with the variation of the constant, the ansatz


Back insertion into (50) leads to a first order differential equation for ,

which is decomposed into partial fractions

The ensuing integral over is solved by aid of the substitution ,

Back into (53)—using to indicate placement of any member of an anonymous bag of constants of integration,


The subscript denotes values at the starting point of the curve,


Solving for yields


Compared with the differential version of (22),


we must have . Note that (56) is essentially a write-up for and could be derived quickly by inserting (40) directly into (60).

The square root of (56) is


The sign in front of the square root is to be chosen positive for pieces of the trajectory with (northern direction), negative where (southern direction). The value in the square root may run through zero within one curve, at one point , such that the square root switches sign there (Figure 3). This happens whenever


yields a vanishing discriminant of the square root for some , that is, whenever the difference is sufficiently large to create a point of minimum polar distance along the trajectory that is not one of the terminal points.

Figure 3: Two examples of trajectories of [left, equation (61)] or [right] as a function of . They may or may not pass through a zero of (61) while connecting the starting abscissa with the final abscissa . There are four basic topologies, depending on whether is positive or negative, and depending on whether the sign change in is from to or from to at this point.

iv.2 Launching Direction

So far we have written the bundle of all geodesics through in the format (61), which specifies the change of latitude as a function of distance traveled. The direction at point 1 in the topocentric system of coordinates is represented by . To address the inverse problem of geodesy, that is to pick the particular geodesics which also runs through the terminal point in fulfillment of the Dirichlet boundary conditions, the associated change in longitude, some form of (40),


must obviously get involved. The strategy is to write down or with as a parameter, then to adjust to ensure with what would be called a shooting method that starting at point 1 at that angle eventually passes through point 2. Coupling of to is done by division of (61) and (63),


Separation of both variables yields


For short paths, small , the integral is simply . If one of the cases occurs, where the difference in is too large for a solution, the additional path through the singularity is to be used [with some local extremum in the graph ], and the integral is to be interpreted as . Taking the sign change and the symmetry of the integrand into account, this is , twice the underivative at minus the sum of the underivative at and . In the right plot of Figure 3, the difference is in including or not including the loudspeaker shaped area between and .

is the solution of (62), the positive value of the solution taken if at , the negative value if at . The squared zero of , , is a root of the fourth-order polynomial which emerges by rewriting (62),