Length of meridian arc on an ellipsoid

Suppose the Earth is modelled as an ellipsoid with semi-major axis a, semi-minor axis b, and eccentricity e. A standard problem is to find the length of a meridian arc between two geographical latitudes, say φ₁ and φ₂. Three methods are suggested below. For more on this topic see e.g. the lecture notes by Deakin and Hunter [1].

The radius of curvature in meridian at latitude φ is given by

ρ = a(1 − e2)(1 − e2sin2φ)−3/2 = a(1 + ε)1/2(1 + εcos2φ)−3/2

where  ε  =  e²/(1 − e²)  =  (a² − b²)/b².

We work here in terms of ε, which is assumed to be small (ε ≈ 1/150 for the Earth). The arc length is given by

a(1 + ε)1/2 ∫ φ2 (1 + εcos2φ)−3/2dφ. φ1

Method 1: Numerical integration

In theory this method is subject to rounding errors, due to the addition of many terms. In practice it seems to work well enough, provided one uses an accurate integration formula such as Weddle’s rule (rather than the better known Simpson’s and three-eighths rules). In Weddle’s rule the number of intervals is a multiple of 6, and the integral over a block of 6 intervals is estimated as

∫ 6h y (x)dx = 3h  [ y (0) + 5y (h) + y (2h) + 6y (3h) + y (4h) + 5y (5h) + y (6h)]. 0 10

The number of intervals can be repeatedly doubled until the value stays constant within the desired order of accuracy. This method is slow, but is useful as a check on other methods.

Method 2: integration of the power series

Here we expand the integrand  (1 + εcos²φ)^−3/2  by the binomial theorem, taking as many terms as are required to achieve the desired accuracy. The powers  cos^2mφ  are integrated by rewriting them in terms of multiple angles:

cos2mφ = 21 − 2m [ 1 ( 2m )  +  ( 2m ) cos2φ +  ( 2m ) cos4φ +  ( 2m ) cos6φ + … + cos2mφ ] 2 m m − 1 m − 2 m − 3

We may as well multiply the result by the binomial expansion of (1 + ε)^1/2, since this does not make it any more complicated. If this is done, the length of the meridian arc is

a [c0(φ2 − φ1) + c1(sin2φ2 − sin2φ1) + c2(sin4φ2 − sin4φ1) + c3(sin6φ2 − sin6φ1) + … ],

where the  c_m  are power series in  ε. To reduce the effect of rounding errors when  φ₂ − φ₁  is small, it may be better to rewrite  sin2mφ₂ − sin2mφ₁  as  2sinm(φ₂ − φ₁)cosm(φ₂ + φ₁). The coefficients  c_m  are given below up to  ε⁸  (equivalent to e¹⁶), which should be more than enough for most purposes.

c0 = 1 −  1 ε +  13 ε2 −  45 ε3 +  2 577 ε4 −  9 417 ε5 +  139 613 ε6 −  522 821 ε7 +  126 287 705 ε8 − … 4 64 256 16 384 65 536 1 048 576 4 194 304 1 073 741 824

c1 = −  3 ε +  9 ε2 −  237 ε3 +  819 ε4 −  23 325 ε5 +  84 711 ε6 −  4 993 233 ε7 +  18 593 103 ε8 − … 8 32 1 024 4 096 131 072 524 288 33 554 432 134 217 728

c2 = 15 ε2 −  75 ε3 +  1 245 ε4 −  4 905 ε5 +  607 125 ε6 −  2 332 785 ε7 +  35 785 995 ε8 − … 256 1 024 16 384 65 536 8 388 608 33 554 432 536 870 912

c3 = −  35 ε3 +  245 ε4 −  19 985 ε5 +  90 475 ε6 −  3 093 755 ε7 +  12 812 765 ε8 − … 3 072 12 288 786 432 3 145 728 100 663 296 402 653 184

c4 = 315 ε4 −  2 835 ε5 +  68 859 ε6 −  354 627 ε7 +  26 776 197 ε8 − … 131 072 524 288 8 388 608 33 554 432 2 147 483 648

c5 = −  693 ε5 +  7 623 ε6 −  430 353 ε7 +  498 267 ε8 − … 1 310 720 5 242 880 167 772 160 134 217 728

c6 = 1 001 ε6 −  13 013 ε7 +  418 847 ε8 − … 8 388 608 33 554 432 536 870 912

c7 = −  6 435 ε7 +  96 525 ε8 − … 234 881 024 939 524 096

c8 = 109 395 ε8 − … 17 179 869 184

Method 3: Another power series method

The following method may be of interest in that it can be extended by computer to any desired degree of accuracy, without the need to program-in complicated coefficients like those in Method 2. It finds the arc length along the meridian from the equator to a given latitude φ. The arc length from φ₁ to φ₂ can be found by applying the method to φ₁ and φ₂, and taking the difference.

We expand  (1 + εcos²φ)^−3/2  by the binomial theorem, as in Method 2, but now integrate the  cos^2mφ  as follows. Setting

Km =  ∫ φ cos2mφdφ, 0

we find easily that  2mK_m = (2m − 1)K_{m − 1} + sinφcos^{2m − 1}φ, whence

K0 = φ,      K1 =  1 φ +  1 sinφcosφ,      K2 =  3 φ +  3 sinφcosφ +  1 sinφcos3φ,      etc. 2 2 8 8 4

A routine calculation shows that the integral, which has to be multiplied by a(1 + ε)^1/2  (= a²/b)  to get the arc length, is

φ ( 1  −  3 ε 4  +  3 · 15 ε2 4 16  −  3 · 15 · 35 ε3 4 16 36  +  3 · 15 · 35 · 63 ε4 4 16 36 64  −  3 · 15 · 35 · 63 · 99 ε5 + … ) 4 16 36 64 100 + sinφ cosφ (  −  3 ε 4  +  3 · 15 ε2 4 16  −  3 · 15 · 35 ε3 4 16 36  +  3 · 15 · 35 · 63 ε4 4 16 36 64  −  3 · 15 · 35 · 63 · 99 ε5 + … ) 4 16 36 64 100 +  1 sinφ cos3φ 2 ( 15 ε2 16  −  15 · 35 ε3 16 36  +  15 · 35 · 63 ε4 16 36 64  −  15 · 35 · 63 · 99 ε5 + … ) 16 36 64 100 +  1 · 3 sinφ cos5φ 2 4 (  −  35 ε3 36  +  35 · 63 ε4 36 64  −  35 · 63 · 99 ε5 + … ) 36 64 100 +  1 · 3 · 5 sinφ cos7φ 2 4 6 ( 63 ε4 64  −  63 · 99 ε5 + … ) 64 100 +  1 · 3 · 5 · 7 sinφ cos9φ 2 4 6 8 (  −  99 ε5 + … ) 100 + …

It is clear how the pattern can be continued indefinitely. The following tested piece of code, written in Delphi Pascal, shows how the method can be applied to give any degree of precision up to the limit of the language used.

Reference

[1] R.E. Deakin and M.N. Hunter, Geometric Geodesy Part A (PDF lecture notes, March 2006)