Numerical problems in the computation of ellipsoidal harmonics

The correct use of ellipsoidal coordinates and related ellipsoidal harmonic functions can provide a representation of linearized Geodetic Boundary Value Problems (GBVP) much closer to the exact ones than what is usually done in spherical approximation: this becomes important in the present age, since terms of the type e(2)N, possibly amounting to several dozens of centimetres, are nowadays observable. Although the theory of ellipsoidal harmonics has been introduced into geodesy by several authors to treat gravity global models, the numerical computation of ellipsoidal harmonics of high degree and order seems to be more critical than it has been recognized. In particular, exact recursive relations display a quite unstable behaviour, no matter what normalization constants are used; it is only through particular representation of hypergeometric functions that it is possible to find a sound method for numerical manipulation. Also the asymptotic approximations, exploiting the smallness of the eccentricity, e(2), are analysed in relation to their critical behaviour for particular values of degree and order; it is shown that a limit layer theory can provide a simpler, better, and stable approximation of the exact values of ellipsoidal harmonics.