[Proj] Bessel and the geodesic problem

[Charles Karney]

I was reading Bessel's paper on geodesics

  F.W. Bessel,
  Ueber die Berechnung der geographischen Laengen und Breiten aus
  geodaetischen Vermessungen,
  Astron. Nachr. 4, 241-254 (1826),
  http://adsabs.harvard.edu/abs/1826AN......4..241B

(brushing off some very rusty German I learned in Vienna, long long
ago), and was struck by two things:

(1) How much of the problem he tackles in this one paper

  * calculus of variations for a rotationally symmetric figure
  * introduction of the auxiliary sphere;
  * specialization to ellipsoid with the differential equations for
    distance and longitude;
  * identification of integrals as elliptic integrals;
  * series expansions of the integrals;
  * an assessment of the number of terms necessary to provide accuracy
    to 0.001" and 0.0001" (3cm and 3mm);
  * tabulation of the coefficients of the expansions valid for arbitrary
    eccentricity (given to 8 places for the leading terms);
  * a worked example (Seeberg to Dunkirk, a distance of 300818 toises =
    586 km, with the result accurate to 0.001").

Essentially, he cleans up the direct geodesic problem, with the minor
caveat that he has to pull a hack to express the longitude integral in
terms of a single parameter so that the coefficients can be tabulated
(he pushes the dependence on the second parameter off into a higher
order term).

(2) However, the real eye-opener was that Bessel uses a change of
variable so that half the terms in his expansions drop out!  (Having to
evaluate a series several hundred times in order to compile an
exhaustive table is powerful motivation to look for such tricks.)

This change of variable essentially exposes a symmetry between prolate
and oblate ellipsoids and this symmetry causes the coefficients of
sin(2*k*sigma)*eps^m to vanish if k+m is odd.

The peculiar thing is that this technique seems to have been promptly
forgotten.  Certainly most of the 20th century literature (Rainsford,
Vicenty, Bowring, Rapp) doesn't use it.

Here's a comparison of the expansions for the distance integral for the
geodesic problem (DISPLAY WITH A FIXED WIDTH FONT!).

The integral is question is

  s/b = integral(sqrt(1 + u2 * sin(sigma)^2, sigma)

where u2 = e^2/(1-e^2) * cos(alpha0)^2 and e^2 = f * (2-f), s =
distance, b = minor radius, f = flattening, sigma = arc length on
auxiliary sphere.

The expansion [Rainsford, Rapp, et al.] of s/b as a Taylor series in u2
up to O(f^4) is:

                          2       3         4
                 u2   3 u2    5 u2    175 u2
   sigma * ( 1 + -- - ----- + ----- - ------- + . . . )
                 4     64      256     16384

                             2        3        4
                      u2   u2    15 u2    35 u2
 - sin(2 * sigma) * ( -- - --- + ------ - ------ + . . . )
                      8    32     1024     4096

                        2       3        4
                      u2    3 u2    35 u2
 - sin(4 * sigma) * ( --- - ----- + ------ + . . . )
                      256   1024    16384

                        3        4
                      u2     5 u2
 - sin(6 * sigma) * ( ---- - ----- + . . . )
                      3072   12288

                          4
                      5 u2
 - sin(8 * sigma) * ( ------ + . . . )
                      131072


Bessel (1926) substitutes

  u2 = 4 * eps / (1 - eps)^2

and expands (1-eps) * s/b as a Taylor series in eps up to O(f^4) to
give:

                     2      4
                  eps    eps
   sigma * ( 1 + ---- + ---- + . . . )
                   4      64

                               3
                      eps   eps
 - sin(2 * sigma) * ( --- - ---- + . . . )
                       2     16

                         2      4
                      eps    eps
 - sin(4 * sigma) * ( ---- - ---- + . . . )
                       16     64

                         3
                      eps
 - sin(6 * sigma) * ( ---- + . . . )
                       48

                           4
                      5 eps
 - sin(8 * sigma) * ( ------ + . . . )
                       512

With Bessel's substitution, there are half the number of terms and the
convergence is faster.  (Bessel writes the coefficients so that the
pattern is clear.  See the expressions for A, B, and C on p. 246 of his
paper.)

A final note: The notation in Bessel's paper is thoroughly modern.  A
colleague reminded me that vector notation had to wait another 75 years
for Gibbs.  Aside from this, the minor pecularities are

  * RR, ee, for e^2, R^2 (but the exponent notation is used elsewhere);

  * variable notation for trig functions (e.g., tgt and tang for tan);

  * "cos m^2" for (cos m)^2 (now printed as cos^2 m);

  * log10(0.2) = 9.301, instead of the 20th century notation of bar1.301
    (and of course there is no 21th century notation for logs as a means
    of hand calculation).

You've got to love that 19th century mathematics!

-- 
Charles Karney <ckarney at sarnoff.com>
Sarnoff Corporation, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662