# A Simple Perturbation Algorithm for Inverting the Cartesian to Geodetic Transformation

Authors: James D. Turner; Tarek Elgohary
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

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 the Earth 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.


---

## Body

## 1. 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λg,ϕg, and h, which denote the geodetic longitude of the satellite subpoint g, 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 (xs,ys,zs) coordinates is given by [1]
(1)xs=(N(ϕg)+h)cosϕgcosλg,ys=(N(ϕg)+h)cosϕgsinλg,zs=(N(ϕg)(1-e2)+h)sinϕg,
where N(ϕg)=a/1-e2sin2ϕg denotes the ellipsoid radius of curvature in the prime vertical plane defined by vectors n^ (ellipsoid outward normal) and τ^ (local east), h is assumed to lie along n^, a denotes the semimajor axis, b denotes the semiminor axis (Figure 2), and e denotes the eccentricity of the Earth’s reference ellipsoid. The solution for λg=tan-1(ys/xs) is obtained by elementary methods.Figure 1
Geodetic and Cartesian coordinates.Figure 2
Geodetic versus geocentric latitude.Because of the fundamental problem of nonlinearity, successive approximation strategies are required forϕg and h. A two-step process is introduced to solve for ϕg. First, one inverts for the geocentric latitude, ϕc, Figure 2; second, using standard trigonometric identities, ϕg is recovered. A successive approximation strategy is developed by introducing a naturally available geometric perturbation variable, which is defined as p=a/b-1~0.0034. Rapidly convergent approximations are obtained for ϕc and h in the τ-z plane by developing power series in the expansion variable p. Elementary vector methods are introduced for inverting for the satellite height, h. 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–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–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 third-order Newton’s method. The results presented by the authors show millimeter accuracy in height and 10−8degree 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 aclosed-form solution for the Cartesian-to-Geodetic transformation throughout the LEO-to-GEO range of applications.


## 2. Mathematical Formulation

The problem is formulated by introducing a local coordinate system that tracks the localx-y axis motion of the satellite. In the local coordinate system, a simplified perturbation solution is developed in the τ-z plane by defining a vector constraint of the form
(2)r-rg-hn^=0,
where r=(rxy,z) denotes the satellite position vector, with rxy=x2+y2, rg=(acos(ϕc),bsin(ϕc)) denotes the satellite ground track point, ϕc denotes the geocentric latitude, h denotes the height of the satellite above the Earth’s surface, and n^=(cos(ϕc)/a,sin(ϕc)/b)/((cos(ϕc)/a)2+(sin(ϕc)/b)2) denotes the unit vector that is normal to the Earth’s surface and points at the satellite. Expanding (2) provides two necessary conditions:
(3)rxy-acos(ϕc)-hcos(ϕc)a((cos(ϕc)/a)2+(sin(ϕc)/b)2)=0,z-bsin(ϕc)-hsin(ϕc)b((cos(ϕc)/a)2+(sin(ϕc)/b)2)=0.Clearly, the equations are highly nonlinear. To begin the simplification process, one replacesa in (3) with
(4)a=b(1+p),p≈0.00314,
which exploits the natural parameter for the problem and transforms (3) into
(5)rxy-b(1+p)cos(ϕc)-hcos(ϕc)((cos(ϕc))2+((1+p)sin(ϕc))2)=0,z-bsin(ϕc)-hsin(ϕc)((cos(ϕc)/(1+p))2+(sin(ϕc))2)=0.An approximate solution is recovered by assuming that the geocentric latitude and satellite height are expanded in the power series representations as follows:(6)ϕc=ϕ0+pϕ1+p2ϕ2+⋯,h=h0+ph1+p2h2+⋯.Introducing (6) into (5) and collecting terms in powers of p yield the cascade of necessary conditions as follows:
(7)O(p0):h0=rxy2+z2-b,ϕ0=2tan-1(rxy2+z2z-rxyz),O(p1):h1=-bcos(ϕ0)2,ϕ1=b-h02(b+h0)sin(2ϕ0),O(p2):h2=b(3b-h0)8(b+h0)sin(2ϕ0)2,ϕ2=h02-4bh0+3b28(b+h0)2sin(4ϕ0)+sin(2ϕ0)4,O(p3):h3=bsin2(2ϕ0)8(b+h0)2[(h0-3b)2cos2(ϕ0)+4b(h0-b)],ϕ3=sin(2ϕ0)6(b+h0)3[C1cos4(ϕ0)+C2cos2(ϕ0)+C3],O(p4):h4=bsin2(2ϕ0)32(b+h0)3[C4cos4(ϕ0)+C5cos2(ϕ0)+C6],ϕ4=sin(2ϕ0)4(b+h0)4[C7cos6(ϕ0)+C8cos4(ϕ0)+C9cos2(ϕ0)+C10].Simple algebraic manipulations yield the ten coefficients appearing in (7) as the polynomials presented in Table 1.Table 1
List of polynomial coefficients.
Coefficient
Expression
C
1
-
4
h
0
3
+
37
b
3
-
66
b
2
h
0
+
33
b
h
0
2
C
2
h
0
3
-
31
b
3
+
75
b
2
h
0
-
33
b
h
0
2
C
3
3
b
(
b
2
+
3
h
0
2
-
6
b
h
0
)
C
4
139
b
3
-
5
h
0
3
+
49
b
h
0
2
-
143
b
2
h
0
C
5
h
0
3
-
127
b
3
+
163
b
2
h
0
-
45
b
h
0
2
C
6
4
b
(
-
10
b
h
0
+
h
0
2
+
5
b
2
)
C
7
4
h
0
4
+
118
b
4
+
198
b
2
h
0
2
-
266
b
3
h
0
-
54
b
h
0
3
C
8
-
155
b
4
-
2
h
0
4
+
67
b
h
0
3
+
421
b
3
h
0
-
315
b
2
h
0
2
C
9
49
b
4
-
15
b
h
0
3
-
185
b
3
h
0
+
135
b
3
h
0
2
C
10
2
b
2
(
10
b
h
0
-
5
h
0
2
-
b
2
)These analytic results are very compact for a fourth-order perturbation expansion. The conversion from the geocentric to the geodetic latitude is given by(8)ϕg=tan-1(cosϕca/bsinϕc).


## 3.  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, ϕg, 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.Figure 3
Errors in latitude, second-order expansion.Figure 4
Errors in height, second-order expansion.Figure 5
Errors in latitude, third-order expansion.Figure 6
Errors in height, third-order expansion.Figure 7
Errors in latitude, fourth-order expansion.Figure 8
Errors in height, fourth-order expansion.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.


## 4. 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.



---

*Source: 712729-2013-07-10.xml*