# 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 φ1 and φ2. 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  ε  =  e2/(1 − e2)  =  (a2b2)/b2.

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 + εcos2φ)−3/2  by the binomial theorem, taking as many terms as are required to achieve the desired accuracy. The powers  cos2mφ  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  cm  are power series in  ε. To reduce the effect of rounding errors when  φ2 − φ1  is small, it may be better to rewrite  sin2mφ2 − sin2mφ1  as  2sinm2 − φ1)cosm2 + φ1). The coefficients  cm  are given below up to  ε8  (equivalent to e16), 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 φ1 to φ2 can be found by applying the method to φ1 and φ2, and taking the difference.

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

 Km = ∫ φ cos2mφdφ, 0

we find easily that  2mKm = (2m − 1)Km − 1 + sinφcos2m − 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  (= a2/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.

```const MAX_POWER_EPSILON = 8; // arbitrary: higher if machine precision warrants it var coeff : array [0..MAX_POWER_EPSILON] of extended; // Delphi "extended" is floating-point with 64-bit precision (19 sig. fig.) ————————————————————————————————— Compute the coefficients of phi, sin(phi)*cos(phi), sin(phi)*cos^3(phi), etc. This procedure need only be called once for a given ellipsoid (a, b). procedure ComputeCoefficients( a, b : extended); var epsilon : extended; j, k : integer; denom, mult : extended; begin epsilon := (a - b)*(a + b) / (b*b); // Compute sums shown inside large () in Web text j := MAX_POWER_EPSILON; denom := 4.0*j*j; coeff[j] := (denom - 1.0)/denom; for j := MAX_POWER_EPSILON - 1 downto 0 do begin for k := MAX_POWER_EPSILON downto j + 1 do coeff[k] := -epsilon * coeff[k]; if (j > 0) then begin denom := 4.0*j*j; coeff[j] := (coeff[j+1] + 1.0) * (denom - 1.0)/denom; end else coeff[0] := coeff[1] + 1.0; end; // Apply factors 1/2, 1.3/2.4, 1.3.5/2.4.6, etc. and dimensions of ellipsoid. mult := a*a/b; coeff[0] := mult*coeff[0]; coeff[1] := mult*coeff[1]; for j := 2 to MAX_POWER_EPSILON do begin mult := mult*(2*j - 3)/(2*j - 2); coeff[j] := mult*coeff[j]; end; end; ————————————————————————————————— Apply the above coefficients to compute arc length on meridian from equator to passed-in latitude phi. function ComputeArcLengthFromEquator( phi_degrees : extended) : extended; var phi, c, c2, total : extended; j : integer; begin phi := (PI/180.0) * phi_degrees; // degrees to radians c := Cos( phi); c2 := c*c; total := coeff[MAX_POWER_EPSILON]; for j := MAX_POWER_EPSILON - 1 downto 1 do total := total*c2 + coeff[j]; result := phi*coeff[0] + Sin(phi)*c*total; end; ```

## Reference

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