A Simple Perturbation Algorithm for Inverting the Cartesian to Geodetic Transformation

A singularity-free perturbation solution is presented for inverting the Cartesian to Geodetic transformation. Geocentric latitude is used to model the satellite ground track position vector. A natural geometric perturbation variable is identified as the ratio of the major and minor Earth ellipse radii minus one. A rapidly converging perturbation solution is developed by expanding the satellite height above theEarth and the geocentric latitude as a perturbation power series in the geometric perturbation variable.The solution avoids the classical problem encountered of having to deal with highly nonlinear solutions for quartic equations. Simulation results are presented that compare the solution accuracy and algorithm performance for applications spanning the LEO-to-GEO range of missions.


Introduction
A frequent calculation for satellites in low Earth orbit (LEO) to geosynchronous Earth orbit (GEO) involves inverting transformations between 3D satellite Cartesian Earth centred coordinates and geodetic coordinates.The geodetic coordinates consist of   ,   , and ℎ, which denote the geodetic longitude of the satellite subpoint , the geodetic latitude of the satellite, and the height of the satellite above the reference Earth elliptical surface along the surface normal from the geodetic ellipsoid to the satellite position.Referring to Figure 1, the transformation from geodetic coordinates to Cartesian (  ,   ,   ) coordinates is given by [1]   = ( (  ) + ℎ) cos   cos   ,   = ( (  ) + ℎ) cos   sin   ,   = ( (  ) (1 −  2 ) + ℎ) sin   , (1) where (  ) = /√ 1 −  2 sin 2   denotes the ellipsoid radius of curvature in the prime vertical plane defined by vectors n (ellipsoid outward normal) and τ (local east), ℎ is assumed to lie along n,  denotes the semimajor axis,  denotes the semiminor axis (Figure 2), and  denotes the eccentricity of the Earth's reference ellipsoid.The solution for   = tan −1 (  /  ) is obtained by elementary methods.
Because of the fundamental problem of nonlinearity, successive approximation strategies are required for   and ℎ.A two-step process is introduced to solve for   .First, one inverts for the geocentric latitude,   , Figure 2; second, using standard trigonometric identities,   is recovered.A successive approximation strategy is developed by introducing a naturally available geometric perturbation variable, which is defined as  = / − 1 ∼ 0.0034.Rapidly convergent approximations are obtained for   and ℎ in the  −  plane by developing power series in the expansion variable .Elementary vector methods are introduced for inverting for the satellite height, ℎ.The resulting analytic perturbation solutions are remarkably simple and computationally efficient.
Many methods have been proposed for implementing the inverse of the transformation presented in (1).The nonlinear Cartesian-to-Geodetic transformation problem is challenging, as geometrical singularities plague many solution strategies.The solution for the geodetic longitude,  however, is elementary and noniterative.The most common problem encountered is the need for handling sensitive quartic polynomial solutions [2][3][4][5].The analytic complexity of the problem arises because the geodetic latitude and satellite height solution algorithms are coupled and highly nonlinear.Three classes of methods have been proposed: (i) closed-form solutions for cubic and quartic polynomials, (ii) perturbation methods, and (iii) successive approximation algorithms.The closed-form class of solution algorithms typically introduces sequences of trigonometric transformations that exploit identities to simplify the governing equation.Important examples of this approach include tha following: (i) the very well-known solution in [6], where the reduced latitude is iterated in Newton's method; (ii) a closed-form solution for a high-order algebraic equation [7]; (iii) introducing the geodetic height of the satellite to develop an elliptic integral-based arc-length solution [8]; (iv) development of an approximate closed-form solution [9]; and (v) introduction of complicated algebraic transformations to develop a series solution [10].Not only are the proposed closed-form solutions highly accurate, but they are also computationally expensive to perform.

Many iterative techniques have been proposed
. Early examples of this approach include the work of [1] which influenced the GPS-based need for the geodetic transformation methods developed by [11][12][13].Several innovative problem formulations have been proposed, including the work of [3,5,14,15].Unfortunately, geometric singularities plague many of these iterative strategies.To avoid troublesome singularities, several authors have investigated vector methods, including the work of [4,16].In [17], an elegant optimization-based strategy is presented.Accelerated convergence techniques are considered by [18] who has presented a third-order version of Newton's method (known as Halley's method).Recently, [19] has presented a very fast singularity-free second-order perturbation solution that introduces an artificial perturbation variable to transform the classical quartic solution problem into a singularity-free noniterative quadratic equation problem.In a more recent addition to iterative methods [20], the projection of a point on the reference ellipsoid is used to solve a system of nonlinear equations using second-and thirdorder Newton's method.The results presented by the authors show millimeter accuracy in height and 10 −8 degree accuracy in latitude with the third-order approach.
The main contribution of this paper is the presentation of a noniterative series-based solution algorithm that effectively provides a closed-form solution for the Cartesian-to-Geodetic transformation throughout the LEO-to-GEO range of applications.

Mathematical Formulation
The problem is formulated by introducing a local coordinate system that tracks the local - axis motion of the satellite.In the local coordinate system, a simplified perturbation solution is developed in the  −  plane by defining a vector constraint of the form where r = (  , ) denotes the satellite position vector, with (3) Clearly, the equations are highly nonlinear.To begin the simplification process, one replaces  in (3) with which exploits the natural parameter for the problem and transforms (3) into An approximate solution is recovered by assuming that the geocentric latitude and satellite height are expanded in the power series representations as follows: Introducing ( 6) into (5) and collecting terms in powers of  yield the cascade of necessary conditions as follows: Table 1: List of polynomial coefficients.

Coefficient Expression
Simple algebraic manipulations yield the ten coefficients appearing in (7) as the polynomials presented in Table 1.
These analytic results are very compact for a fourth-order perturbation expansion.The conversion from the geocentric to the geodetic latitude is given by

Numerical Results
The perturbation expansion method is used to carry out the coordinate transformation for several cases of LEO-to-GEO orbits.Using the WGS84, the forward transformation is carried out first, then the perturbation solution is applied, and the results are compared with the original values, which represent exact values for the inverse solution.For the sake of demonstration a longitude angle of 30 ∘ is utilized.The geodetic latitude,   , is swept for angles from −90 to 90 degrees and the height is swept from 200 KM (LEO) to 35,000 KM (GEO).First, the expansion is carried to second order, and the errors in latitude and height are plotted as functions of the true latitudes and heights as shown in Figures 3 and 4, respectively.The expansion is then carried out to third order, and the errors in latitude and height are plotted in Figures 5 and 6, respectively.Finally, the fourth-order expansion is utilized, and results are shown in Figures 7 and  8, respectively.The improvement of accuracy is quite obvious as the order of expansion is increased.A two-order-of-magnitude  improvement is achieved by adding the third-order terms to each of the coordinates.Another two-order-of-magnitude improvement is achieved with the fourth-order terms.In height, millimeter accuracy is achieved at the fourth-order expansion level.This shows the fast convergence nature and the accuracy of the perturbation solution.These results demonstrate that higher-order approximations do not provide additional useful information for the inversion process.

Conclusion
Earth-Centered Earth-Fixed (ECEF) to geodetic coordinate transformation has been examined with several numerical and analytical approaches throughout the literature.A noniterative expansion-based approach inspired by the Earth's perturbed geometry is introduced in this work, where the expansion parameter is nothing but the ratio of the Earth semimajor axis and semi-minor axis subtracted from 1.The expansion is carried out to second, third, and fourth orders.
A numerical example is introduced to compare the accuracies at each order of expansion.Accuracies showed significant improvements as the order of expansion is increased, and the at fourth order, millimeter accuracy is achieved in height and 10 −11 degree error in latitude.Those errors at such low orders of the expansion are proof of the effectiveness of the method and its potential in solving such a highly nonlinear transformation noniteratively.The method can be further streamlined for timing studies, but in general it is a clean straightforward approach to the coordinate transformation problem that utilizes a physical perturbation parameter and that proved to be very accurate and efficient.