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
###### Abstract

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
###### pacs:
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:

 ρ2p=ρ2e(1−e2);x2ρ2e+z2ρ2p=1. (1)

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).

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

 tanϕ=zxρ2eρ2p;zx=tanϕ′=ρ2pρ2etanϕ. (2)

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

 x=ρcosϕ′=ρecosϕ√1−e2sin2ϕ,z=ρsinϕ′=ρe(1−e2)sinϕ√1−e2sin2ϕ. (3)

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

 tanv=e2sin(2ϕ)2(1−e2sin2ϕ)=msin(2ϕ)1+mcos(2ϕ)=msin(2ϕ)−m2sin(2ϕ)cos(2ϕ)+m3sin(2ϕ)cos2(2ϕ)+⋯; (4)

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

 m≡e22−e2. (5)

The Taylor series in powers of and in powers of are

 v = m1+msin(2ϕ)+m2(2+m)6(1+m)3sin3(2ϕ)+m2(5+25m+15m2+3m3)40(1+m)5sin5(2ϕ)+⋯ (6) = sin(2ϕ)m−12sin(4ϕ)m2+13sin(6ϕ)m3−14sin(8ϕ)m4+15sin(10ϕ)m5+⋯. (7)

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

 f=ρe−ρpρe=1−√1−e2;e2=f(2−f). (8)

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

 f=1/298.257223563,ρe=6378137.0 m. (9)

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

 v = 0.0033584338sin(2ϕ)−0.56395388×10−5sin(4ϕ)+0.12626678×10−7sin(6ϕ)+⋯ (10) = 692.′′72669sin(2ϕ)−1.′′16324sin(4ϕ)+0.′′00260sin(6ϕ)+⋯. (11)

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)

 x=ρcosϕ′+hcosϕ;z=ρsinϕ′+hsinϕ, (12)

which can be written in terms of a distance ,

 N(ϕ)≡ρe√1−e2sin2ϕ (13)

as

 x=(N+h)cosϕ;z=[N(1−e2)+h]sinϕ. (14)

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

 ⎛⎜⎝xyz⎞⎟⎠=⎛⎜ ⎜⎝[N(ϕ)+h]cosϕcosλ[N(ϕ)+h]cosϕsinλ[N(ϕ)(1−e2)+h]sinϕ⎞⎟ ⎟⎠. (15)

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

 S=∫21√dx2+dy2+dz2=∫21√ds2≡∫21Lds (16)

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

 eλ=⎛⎜⎝−[N(ϕ)+h]cosϕsinλ[N(ϕ)+h]cosϕcosλ0⎞⎟⎠;eϕ=⎛⎜⎝−sinϕcosλ[M(ϕ)+h]−sinϕsinλ[M(ϕ)+h]cosϕ[M(ϕ)+h]⎞⎟⎠, (17)

where a meridional radius of curvature Tienstra (1951)

 M(ϕ)≡ρe(1−e2)(1−e2sin2ϕ)3/2=N(ϕ)1−e21−e2sin2ϕ (18)

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

 e2λ=E = (N+h)2cos2ϕ; (19) eλ⋅eϕ=F = 0; (20) e2ϕ=G = [N1−e21−e2sin2ϕ+h]2=(M+h)2. (21)

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 ,

 S=∫√Edλ2+2Fdλdϕ+Gdϕ2. (22)

### iii.2 Christoffel Symbols

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

 deϵ=∑α,βeαΓαβϵdxβ. (23)

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

 deλ=⎛⎜⎝−(N+h)cosϕcosλ−(N+h)cosϕsinλ0⎞⎟⎠dλ+⎛⎜⎝(M+h)sinϕsinλ−(M+h)sinϕcosλ0⎞⎟⎠dϕ (24)

and of with respect to and :

 deϕ=⎛⎜⎝(M+h)sinϕsinλ−(M+h)sinϕcosλ0⎞⎟⎠dλ+⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝[Ne2cos2ϕ1−e2sin2ϕ−3Ne2(1−e2)sin2ϕ(1−e2sin2ϕ)2−(N+h)]cosϕcosλ[Ne2cos2ϕ1−e2sin2ϕ−3Ne2(1−e2)sin2ϕ(1−e2sin2ϕ)2−(N+h)]cosϕsinλ[Ne2(1−e2)(3cos2ϕ(1−e2sin2ϕ)2−sin2ϕ1−e2sin2ϕ)−[N(1−e2)+h]]sinϕ⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠dϕ. (25)

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

 deλ = (eλΓλλλ+eϕΓϕλλ)dλ+(eλΓλϕλ+eϕΓϕϕλ)dϕ; (26) deϕ = (eλΓλλϕ+eϕΓϕλϕ)dλ+(eλΓλϕϕ+eϕΓϕϕϕ)dϕ. (27)

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

 Γλλλ=Γϕϕλ=Γϕλϕ=Γλϕϕ=0; (28)
 Γϕλλ=(N+h)cosϕsinϕ1−e2sin2ϕh(1−e2sin2ϕ)+N(1−e2); (29)
 Γλϕλ=Γλλϕ=−sinϕh(1−e2sin2ϕ)+N(1−e2)(N+h)(1−e2sin2ϕ)cosϕ; (30)
 Γϕϕϕ=3Ne2(1−e2)sinϕcosϕ1[h(1−e2sin2ϕ)+N(1−e2)](1−e2sin2ϕ). (31)

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

 d2xϵds2+∑μνΓϵμνdxμdsdxνds=0. (32)

The explicit write-up

 d2λds2+Γλλλdλdsdλds+Γλλϕdλdsdϕds+Γλϕλdϕdsdλds+Γλϕϕdϕdsdϕds=0; (33)
 d2ϕds2+Γϕλλdλdsdλds+Γϕλϕdλdsdϕds+Γϕϕλdϕdsdλds+Γϕϕϕdϕdsdϕds=0, (34)

simplifies with (28) to

 d2λds2+2Γλλϕdλdsdϕds = 0; (35) d2ϕds2+Γϕλλ(dλds)2+Γϕϕϕ(dϕds)2 = 0. (36)

## 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),

 d2λds2dλds=−2Γλϕλdϕds. (37)

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

 ∫Γϕλϕ(ϕ)dϕ=log{[N(ϕ)+h]cosϕ}+const (38)

to generate a first integral

 logdλds=−2log[(N+h)cosϕ]+const. (39)

Exponentiation yields

 dλds=c31(N+h)2cos2ϕ, (40)

where

 c3≡dλds∣1(N1+h)2cos2ϕ1 (41)

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

 d2ϕds2+sinϕ(N+h)3cos3ϕ1−e2sin2ϕh(1−e2sin2ϕ)+N(1−e2)c23+3Ne2(1−e2)sinϕcosϕ[h(1−e2sin2ϕ)+N(1−e2)](1−e2sin2ϕ)(dϕds)2=0. (42)

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

 τ≡sinϕ, (43)

which implies transformations in the derivatives:

 dτds=cosϕdϕds; (44)
 d2τds2=−sinϕdϕdsdϕds+cosϕd2ϕds2; (45)
 cosϕd2ϕds2=d2τds2+τ(dϕds)2=d2τds2+τcos2ϕ(dτds)2=d2τds2+τ1−τ2(dτds)2. (46)

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

 d2τds2+τ1−τ2(dτds)2+τ(N+h)3(1−τ2)1−e2τ2h(1−e2τ2)+N(1−e2)c23+3Ne2(1−e2)τ[h(1−e2τ2)+N(1−e2)](1−e2τ2)(dτds)2=0;
 d2τds2+τ1−τ2h(1−e2τ2)2+N(1−4e2τ2+2e2+4e4τ2−3e4)[h(1−e2τ2)+N(1−e2)][1−e2τ2](dτds)2+τ(N+h)3(1−τ2)1−e2τ2h(1−e2τ2)+N(1−e2)c23=0. (47)

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

 (1−τ2)d2τds2+h(1−e2τ2)2+N(1−e2)(1+3e2−4e2τ2)[h(1−e2τ2)+N(1−e2)][1−e2τ2]τ(dτds)2+τc23(N+h)31−e2τ2h(1−e2τ2)+N(1−e2)=0,

and the standard way of progressing is the substitution

 dτds≡p;d2τds2=pdpdτ; (48)
 (1−τ2)pdpdτ+h(1−e2τ2)2+N(1−e2)(1+3e2−4e2τ2)[h(1−e2τ2)+N(1−e2)][1−e2τ2]τp2+τc23(N+h)31−e2τ2h(1−e2τ2)+N(1−e2)=0. (49)

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

 12(1−τ2)dPdτ+h(1−e2τ2)2+N(1−e2)(1+3e2−4e2τ2)[h(1−e2τ2)+N(1−e2)][1−e2τ2]τP+τc23(N+h)31−e2τ2h(1−e2τ2)+N(1−e2)=0. (50)

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

 dPdτ=−2τh(1−e2τ2)2+N(1−e2)(1+3e2−4e2τ2)[1−τ2][h(1−e2τ2)+N(1−e2)][1−e2τ2]P. (51)

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

 logP = −2∫τh(1−e2τ2)2+N(1−e2)(1+3e2−4e2τ2)[1−τ2][h(1−e2τ2)+N(1−e2)][1−e2τ2]dτ (52) = ∫−2τ1−τ2dτ+3∫−2e2τ1−e2τ2dτ+6he2∫τh(1−e2τ2)+N(1−e2)dτ = log(1−τ2)+3log(1−e2τ2)−2log[h(1−e2τ2)3/2+ρe(1−e2)]+const.
 P=const⋅(1−τ2)(1−e2τ2)3[h(1−e2τ2)3/2+ρe(1−e2)]2=const⋅(1−τ2)(1−e2τ2)2[h(1−e2τ2)+N(1−e2)]2=const⋅(1−τ2)G(τ).

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

 P=c(τ)⋅(1−τ2)(1−e2τ2)2[h(1−e2τ2)+N(1−e2)]2. (53)

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

 dc(τ)dτ1−τ2G(τ)=−2τc23(1−τ2)(N+h)31−e2τ2h(1−e2τ2)+N(1−e2),

which is decomposed into partial fractions

 dc(τ)dτ=c23[eN2eτ1−e2τ21(N+h)3+−2τ1−τ21(N+h)2]11−τ2.

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

 c(τ)=−c231(N+h)2(1−τ2)+const.

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

 P = [const−1(N+h)2(1−τ2)]c23(1−τ2)(1−e2τ2)2[h(1−e2τ2)+N(1−e2)]2 (54) = c5(1−τ2)(1−e2τ2)2[h(1−e2τ2)+N(1−e2)]2−c23(1−e2τ2)2(N+h)2[h(1−e2τ2)+N(1−e2)]2 (55) = c51(h+M)2(1−τ2−c23(N+h)2)=c51−τ2G(τ)(1−c23E(τ))=p2. (56)

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

 P1 = c5cos2ϕ1(1−e2sin2ϕ1)2[h(1−e2sin2ϕ1)+N1(1−e2)]2−c23(1−e2sin2ϕ1)2(N1+h)2[h(1−e2sin2ϕ1)+N1(1−e2)]2=p21 (57) = c5cos2ϕ1(1−e2sin2ϕ1)2[h(1−e2sin2ϕ1)+N1(1−e2)]2−(dλ/ds)21(N1+h)2cos4ϕ1(1−e2sin2ϕ1)2[h(1−e2sin2ϕ1)+N1(1−e2)]2. (58)

Solving for yields

 c5 = (59) = (dλds∣1)2E1+(dϕds∣1)2G1=c23(N1+h)2cos2ϕ1+(dϕds∣1)2(M1+h)2.

Compared with the differential version of (22),

 ds2=Edλ2+Gdϕ2;1=E(dλds)2+G(dϕds)2, (60)

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

 p=±1−e2τ2h(1−e2τ2)+N(1−e2) ⎷1−τ2−c23(N+h)2=±√1−τ2−c23(N+h)2h+M(τ)=dτds. (61)

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

 (1−τ2m)[N(τm)+h]2=c23 (62)

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.

### 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),

 dλds=c31(N+h)2(1−τ2)=c3E, (63)

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),

 dτdλ=dτdsdsdλ=dτds/dλds=±(N+h)2(1−τ2)√1−τ2−c23(N+h)2c3(h+M). (64)

Separation of both variables yields

 ±∫1dτc3(h+M)(N+h)2(1−τ2)√1−τ2−c23(N+h)2=λ∣1. (65)

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),

 e4h4(τ2m)4+2e2h2[ρ2e−h2+e2(c23−h2)](τ2m)3+[(ρ2e−h2)2+2e2(2h4−2c23h2−2