About the geometry of the Earth geodetic reference surfaces

Journal of Geometry and Physics 120 (2017) 192–207

Contents lists available at ScienceDirect

Journal of Geometry and Physics journal homepage: www.elsevier.com/locate/geomphys

About the geometry of the Earth geodetic reference surfaces Ladislav Husár a , Peter Švaral b , Juraj Janák a, * a

Department of Theoretical Geodesy, Faculty of Civil Engineering, Slovak University of Technology, Radlinského 11, 81005 Bratislava, Slovakia b Cortizo Slovakia, a.s., Železničný rad 29, 968 01 Nová Baňa, Slovakia

article

info

Article history: Received 7 November 2015 Received in revised form 5 March 2017 Accepted 22 May 2017

Keywords: Metric tensor Curvilinear coordinates Triaxial ellipsoid Geodesy

a b s t r a c t The paper focuses on the comparison of metrics of three most common reference surfaces of the Earth used in geodesy (excluding the plane which also belongs to reference surfaces used in geodesy when dealing with small areas): a sphere, an ellipsoid of revolution and a triaxial ellipsoid. The two latter surfaces are treated in a more detailed way. First, the mathematical form of the metric tensors using three types of coordinates is derived and the lengths of meridian and parallel arcs between the two types of ellipsoids are compared. Three kinds of parallels, according to the type of latitude, can be defined on a triaxial ellipsoid. We show that two types of parallels are spatial curves and one is represented by ellipses. The differences of curvature of both kinds of ellipsoid are analysed using the normal curvature radii. Priority of the chosen triaxial ellipsoid is documented by its better fit with respect to the high-degree geoid model EIGEN6c4 computed up to degree and order 2160. © 2017 Elsevier B.V. All rights reserved.

0. Introduction Geodesy uses several reference surfaces which approximate the physical surface of the Earth for solving problems. These surfaces generalize its true form by suitable approximation while the level of approximation depends on the nature of the problem. Such surfaces are a sphere, ellipsoid of revolution, triaxial ellipsoid, spheroid and geoid. This paper deals with the basic geometric aspects of these surfaces focusing on differences in their metrics expressed by the metric tensors in various coordinates. It focuses mainly on differences between the metrics of ellipsoid of revolution and triaxial ellipsoid. The historical overview of triaxial ellipsoid as a reference surface of the Earth is well described in [1]. The basic equations concerning on this topic were published by Burša and Pěč [2] and Feltens [3]. This paper extends and further develops their work in several aspects. From the practical point of view the paper analyses several quantitative metrics and curvature indicators and also qualitative differences between the shapes of the coordinate lines, changes in geographic grid orthogonality or differences between the geographic coordinates defined on particular surfaces. Based on the obtained results the applications in which it would be useful to replace the ellipsoid of revolution by triaxial ellipsoid are suggested in Conclusion. The possible fields of application are geodesy in general (cadastre, mapping, engineering surveying, cartography), global geophysics, satellite geodesy or planetary sciences. The lack of symmetry of the triaxial ellipsoid is still a strong reason to stick with the biaxial ellipsoid which is mathematically more simple and easier to use. On the other hand, with the progress of computational tools this problem can be manageable.

*

Corresponding author. E-mail address: [email protected] (J. Janák).

http://dx.doi.org/10.1016/j.geomphys.2017.05.016 0393-0440/© 2017 Elsevier B.V. All rights reserved.

L. Husár et al. / Journal of Geometry and Physics 120 (2017) 192–207

193

1. Metric tensors of Earth’s reference surfaces The metric tensor is a quantity that carries the basic information about the metric characteristics of a particular surface. On a surface we can express a point location using a pair of general curvilinear coordinates uα (α = 1, 2) and also a triplet spatial Cartesian coordinates xi (i = 1, 2, 3). Conventionally, indices in 2D space are designated using the small case Greek alphabet while in 3D space by the Latin alphabet. On each reference surface we use the specific pair of geographic coordinates, longitude and latitude uα (α = 1, 2), u1 – longitude and u2 – latitude to establish a right-handed system of the basis vectors on a surface b1 , b2 or b1 , b2 respectively. We distinguish the spherical coordinates u1 = Λs , u2 = Φs , the geodetic coordinates u1 = λ, u2 = ϕ , the geocentric coordinates u1 = λ, u2 = Φc and the reduced coordinates u1 = λr , u2 = β . The Einstein summation convention is used throughout the paper. 1.1. Sphere The general equation of the sphere fs = fs (Λs , Φs ) with the radius R can be written as follows x1 = R cos Λs cos Φs , fs :

x21

+

x22

+

x23

=R ,

x2 = R sin Λs cos Φs ,

2

R=

√

x21 + x22 + x23 ,

Λs = arctan

(

x2

)

x1

,

(1)

⎞

⎛ x3 = R sin Φs ,

Φs = arctan ⎝ √

x3 x21

+ x22

⎠.

Relations (1) between the Cartesian coordinates xi and curvilinear coordinates Λs , Φs determine the form of the covariant basis bα and the contravariant basis bβ [4, p. 135]

⎞ ∂ x1 ⎜ ∂ Φs ⎟ ( ) ⎜ ⎟ −R cos Λs sin Φs ⎜ ∂ x2 ⎟ ⎜ ⎟ b2 = ⎜ ⎟ = −R sin Λs sin Φs , ⎜ ∂ Φs ⎟ R cos Φs ⎝ ∂ x3 ⎠

⎞ ∂ x1 ⎜ ∂ Λs ⎟ ( ) ⎜ ⎟ −R sin Λs cos Φs ⎜ ∂ x2 ⎟ ⎜ ⎟ b1 = ⎜ ⎟ = R cos Λs cos Φs , ⎜ ∂ Λs ⎟ 0 ⎝ ∂ x3 ⎠

⎛

⎛

∂ Λs

(2a)

∂ Φs

⎞ −x1 x3 √ ⎜ 2 2 ⎟ ⎛ − cos Λs sin Φs ⎞ ∂ Φs ⎜ R x1 + x22 ⎟ ⎜ ∂ x1 ⎟ ⎜ ⎟ ⎟ R ⎜ ⎟ ⎜ −x2 x3 ⎟ ⎜ ⎟ ⎜ ⎜ ⎟ ⎟ ⎜ ∂ Φ ⎜ − sin Λs sin Φs ⎟ s⎟ √ 2 ⎜ ⎜ ⎟ b =⎜ = = ⎜ ⎟. 2⎟ ⎟ ⎜ 2 2 ⎟ R ⎜ ∂ x2 ⎟ ⎜ R x1 + x2 ⎟ ⎜ ⎝ ⎠ ⎟ ⎝ ∂ Φs ⎠ ⎜ √ cos Φs ⎜ ⎟ 2 2 ⎝ x1 + x2 ⎠ R ∂x ⎛

⎛ −x ⎞ ⎛ ⎞ ∂ Λs 2 − sin Λs ⎜ ∂ x1 ⎟ ⎜ ⎟ ⎜ x2 + x22 ⎟ ⎜ R cos Φs ⎟ ⎜ ∂ Λs ⎟ ⎜ 1 ⎟ ⎜ ⎟ 1 ⎜ ⎟ = ⎜ x1 ⎟ = ⎜ cos Λs ⎟ , b =⎜ ⎜ ⎜ ⎟ ⎟ ⎟ ⎜ ∂ x2 ⎟ ⎝ x21 + x22 ⎠ ⎝ R cos Φs ⎠ ⎝ ∂ Λs ⎠ ⎛

⎞

0

0

∂ x3

⎛

⎞

(2b)

3

R2 The basis vectors bα , by definition (2a), are tangential to the coordinate lines on a sphere. The basis vectors bβ , from definition (2b), are perpendicular to the coordinate surfaces. Components of the covariant and contravariant metric tensors of a surface gαβ and gαβ respectively, are defined as inner products of the basis vectors, see Spiegel [4], so that on a sphere after applying Eqs. (2a) and (2b) hold

⎛ ( gαβ = bα · bβ =

R cos Φs 0 2

2

0 R2

)

,

1

2 ⎜ 2 gαβ = bα · bβ = ⎝ R cos Φs 0

⎞ 0

. 1⎠

⎟

(3)

R2

From Eq. (3) one can see that gαβ gαβ = E ⇒ gαβ = (gαβ )−1 , where E is a unit matrix. This results from the property of the β β conjugate bases bα , bβ , which fulfil the reciprocity conditions bα · bβ = δα , where δα is the Kronecker delta. 1.2. Ellipsoid of revolution On an ellipsoid of revolution we use three types of curvilinear coordinates: geodetic coordinates u1 = λ, u2 = ϕ , reduced coordinates u1 = λ, u2 = β and geocentric coordinates u1 = λ, u2 = Φc . We demonstrate the advantage of using the

194

L. Husár et al. / Journal of Geometry and Physics 120 (2017) 192–207

reduced coordinates. Note that these types of coordinates differ only in latitude. Mutual relations among particular latitudes are shown in Fig. 1 and described by Eq. (4), see e.g. [5] tan Φc tan β = √ , 1 − e2

tan β =

√

1 − e2 tan ϕ,

tan Φc = (1 − e2 ) tan ϕ,

(4)

where e = (a2 − c 2 )/a2 is the first numerical eccentricity of ellipsoid of revolution with the major semi-axis a and the minor semi-axis c. The mathematical expression of the surface of the ellipsoid of revolution fre in different coordinates fre (λ, ϕ ), fre (λ, β ), fre (λ, Φc ) results from Eqs. (4) and is given by relations

√

Coordinate type

fre :

x21 + x22 a2

+

x23 c2

= 1,

Geodetic - G

Reduced - R

x1 (λ, ϕ ) = N cos λ cos ϕ,

x1 (λ, β ) = a cos λ cos β,

x2 (λ, ϕ ) = N sin λ cos ϕ,

x2 (λ, β ) = a sin λ cos β,

x3 (λ, ϕ ) = N(1 − e2 ) sin ϕ,

x3 (λ, β ) = c sin β,

a=

√

c2 + E2,

(5)

Geocentric - C x1 (λ, Φc ) = ρ cos λ cos Φc x2 (λ, Φc ) = ρ sin λ cos Φc x3 (λ, Φc ) = ρ sin Φc , where E = ae is the linear eccentricity and N is the prime vertical radius of curvature. The covariant basis vectors bα on an ellipsoid of revolution result from definition (2a) and Eqs. (5) coordinate type G R C ⎞ ∂ x1 ⎜ ∂ u1 ⎟ ( ) ) ( ) ( ⎜ ⎟ −N sin λ cos ϕ −a sin λ cos β −ρ sin λ cos Φc ⎜ ∂ x2 ⎟ b1 = ⎜ ⎟ = N cos λ cos ϕ = a cos λ cos β = ρ cos λ cos Φc , ⎜ ∂ u1 ⎟ 0 0 0 ⎝ ∂x ⎠ 3 ∂ u1

⎛

coordinate type

G

R

(6a)

C

( ) ⎞ ⎛ ⎞ ∂N ∂ x1 cos λ cos ϕ − N sin ϕ ⎜ ⎟ ∂ϕ ⎜ ∂ u2 ⎟ ⎜ ( ) ( ) ⎟ ⎜ ⎟ ⎟ ⎜ −a cos λ sin β ∂N ⎜ ∂ x2 ⎟ ⎜ ⎟ b2 = ⎜ cos ϕ − N sin ϕ ⎟ = −a sin λ sin β ⎟ = ⎜ sin λ ⎜ ∂ u2 ⎟ ⎜ ⎟ ∂ϕ c cos β ⎝ ∂x ⎠ ⎜ )⎟ ( 3 ⎝ ⎠ ∂N 2 (1 − e ) sin ϕ + N cos ϕ ∂ u2 ∂ϕ ( )⎞ ⎛ ∂ρ cos λ cos Φc − ρ sin Φc ⎜ ⎟ ∂ Φc ⎜ ( )⎟ ⎜ ⎟ ∂ρ ⎟ =⎜ ⎜ sin λ ∂ Φ cos Φc − ρ sin Φc ⎟ . c ⎜ ⎟ ⎝ ⎠ ∂ρ sin Φc + ρ cos Φc ∂ Φc ⎛

(6b)

√

Partial derivatives in Eq. (6b) can be derived by√ differentiation of the known expressions of N (ϕ) = a/ 1 − e2 sin2 ϕ [5] √ and the geocentric radius vector ρ (Φc ) = ρ = a 1 − e2 / 1 − e2 cos2 Φc , which can be obtained by simple modification of Eqs. (5). We get

∂N Ne2 sin ϕ cos ϕ = , ∂ϕ 1 − e2 sin2 ϕ

∂ρ −ρ e2 sin Φc cos Φc = . ∂ Φc 1 − e2 cos2 Φc

(7)

From Eqs. (5) and (6) follow that the most convenient coordinates for an ellipsoid of revolution are reduced coordinates as they do not contain the terms depended on N or ρ which substantially simplify the basis vectors and also the covariant

L. Husár et al. / Journal of Geometry and Physics 120 (2017) 192–207

195

Fig. 1. Different latitudes defined on an ellipsoid of revolution.

metric tensor gαβ that can be obtained as inner products of the basis vectors b1 , b2 from Eqs. (6)

⎛ G

gαβ

N 2 cos2 ϕ,

=⎝

gRαβ =

(

⎞

0 N2 +

a2 cos2 β, 0,

∂N ∂ϕ

)2

)2

( 2 2 N cos ϕ, ⎠ ) = 2

∂N + N cos ϕ + sin ϕ e4 − 2e ∂ϕ ⎛ 2 2 ⎞ ) ρ cos Φc , 0 )2 ( 0 C ⎠, , gαβ = ⎝ ∂ρ c 2 + E 2 sin2 β 0, ρ2 + ∂ Φc

0,

(

(

(

0

0 M2

)

1

(8)

where M = a(1 − e2 )/(1 − e2 sin2 ϕ )3/2 , see e.g. [5], is the meridian radius of curvature. The upper index in the metric tensor indicates the type of coordinates. The contravariant form of the metric tensor gαβ can be obtained simply by inversion of gαβ . Obviously, the simplest form of the metric tensor of the ellipsoid of revolution is given in terms of reduced coordinates. The zero off-diagonal elements of the metric tensors in Eq. (6b) confirm that all three types of curvilinear coordinates uα are orthogonal. 1.3. Triaxial ellipsoid When dealing with a more complicated surface we use geocentric polar coordinates ρ, λ, Φc . The mathematical expression of the triaxial ellipsoid in Cartesian coordinates is then given by the equation x21

x22

x23

= 1. (9a) c2 Substituting the Cartesian coordinates xi by the geocentric coordinates ρ, λ, Φc from Eq. (5) we get after simple manipulation the parametric expression of the geocentric radius-vector of the triaxial ellipsoid in form ρ3 = ρ3 (λ, Φc ) (Burša and Pěč, 1988) f3e :

a2

+

b2

+

√

ρ3 = √

a (1 − e2 )(1 − e21 ) 1 − e2 cos2 Φc − e21 sin2 Φc − e21 (1 − e2 )cos2 Φc cos2 (λ0 )

,

λ0 = λ − λa .

(9b)

In Eqs. (9) e is the numerical eccentricity e2 = (a2 − c 2 )/a2 , e1 is the equatorial numerical eccentricity e21 = (a2 − b2 )/a2 , λa is the geographical longitude of the major semi-axis a of the equatorial ellipse, b is the minor semi-axis of the equatorial ellipse, c is the polar semi-axis of the triaxial ellipsoid and index 3 in the geocentric radius-vector is used in order to distinguish it from the radius-vector of the ellipsoid of revolution ρ . We can also define reduced latitude β for the triaxial ellipsoid, similarly as β for the ellipsoid of revolution, and because b ̸ = a also reduced longitude λr as functions of geocentric coordinates λ0 , Φc . The ratio of the reduced longitude λr to the geocentric longitude λ0 is the same as the ratio of β to Φc , expressed by the first relationship of Eq. (4). For determination of λr just replacement of the numerical eccentricity e by e1 is needed. The geometry of the equatorial ellipse of a triaxial ellipsoid is shown in Fig. 2. G 1 Simplification of the term g 2 2 = g11 (du1 )2 + g22 (du2 )2 = 22 = M in gαβ also results from comparison with the surface length element: ds

(N cos ϕ )2 dλ2 + M 2 dϕ 2 .

196

L. Husár et al. / Journal of Geometry and Physics 120 (2017) 192–207

Fig. 2. Geometry of the equatorial ellipse of a triaxial ellipsoid.

The reduced latitude β on the triaxial ellipsoid is again given by √ ( the first)relationship of Eq. (4), however, the parameters

of the meridian ellipse are changing with the longitude: e2λ0 =

a2λ0 − c 2 /a2λ0 and aλ0 = ρ3 (λ0 , Φc = 0) = a 1 − e21 /

√

1 − e21 cos2 λ0 , so for the reduced coordinates hold tan λ0

tan λr = √

,

1 − e21

√

a 1 − e21 tan Φc

tan Φc

tan β = √

1 − e2λ0

= √ c

1 − e21 cos2 λ0

.

(10)

The advantage of using reduced coordinates is again evident from the simplification of the radius-vector ρ3 (λr , β ), the parametric expression of the surface xi = xi (uα ) = xi (λr , β )

ρ3 =

√

(x01 )2 + (x02 )2 + (x03 )2 ,

x01 = a cos λr cos β,

x02 = b sin λr cos β,

x03 = c sin β,

(11)

as well as the covariant basis vectors and the metric tensor b3e 1

( ) −a sin λr cos β ∂ x0i = = b cos λr cos β , ∂λr 0

gαβ = a

2

( ) −a cos λr sin β ∂ x0i = = −b sin λr sin β , ∂β c cos β

(12)

(cos β )2 (sin λr )2 + 1 − e21 (cos λr )2 ,

e21 sin λr cos λr sin β cos β

e21 sin λr cos λr sin β cos β,

(sin β )2 (cos λr )2 + 1 − e21 (sin λr )2 + 1 − e2 (cos β )2

( 3e

b3e 2

[

(

)

]

[

(

)

]

(

) )

.

The problems of the metric tensor of irregular surfaces were studied by Burša, see Burša and Pěč, (1988), in case of the geoid. He expressed the metric tensor in geocentric coordinates ρ3 , λ, ϑ = 90 − Φc as ρ3 = ρ3 (ϑ, λ). If we use geocentric latitude Φc instead of polar distance ϑ and apply Eqs. (6) to the modified form of Eqs. (5), where we substitute ρ by ρ3 = ρ3 (Φc , λ) from Eq. (9b), we get the basis vectors b3e α on a triaxial ellipsoid in geocentric coordinates

) ⎞ ∂ρ3 cos λ − ρ3 sin λ cos Φc ⎜ ∂λ ⎟ ⎜( ⎟ ) ⎜ ⎟ ⎜ ∂ρ3 sin λ + ρ cos λ cos Φ ⎟ , b3e = 1 3 c⎟ ⎜ ⎜ ∂λ ⎟ ⎝ ⎠ ∂ρ3 sin Φc ∂λ ⎛(

) ⎞ ∂ρ3 cos Φc − ρ3 sin Φc cos λ ⎟ ⎜ ∂ Φc ⎜( ⎟ ) ⎟ ⎜ ⎜ ∂ρ3 cos Φc − ρ3 sin Φc sin λ ⎟ . b3e = 2 ⎜ ∂Φ ⎟ c ⎟ ⎜ ⎠ ⎝ ∂ρ3 sin Φc + ρ3 cos Φc ∂ Φc ⎛(

(13)

The metric tensor of a triaxial ellipsoid in geocentric coordinates can be obtained as inner products of the basis vectors given by Eq. (13)

⎛ ⎜ρ ⎜ ⎝

g3e αβ = ⎜

2 2 3 cos Φc

( +

∂ρ3 ∂λ

∂ρ3 ∂ρ3 , ∂λ ∂ Φc

)2

,

⎞ ∂ρ3 ∂ρ3 ∂λ ∂ Φc ⎟ ⎟ ( )2 ⎟ . ⎠ ∂ρ 3 ρ32 + ∂ Φc

(14)

L. Husár et al. / Journal of Geometry and Physics 120 (2017) 192–207

197

Eq. (14) can be seen as an alternative expression of Eq. (12) but in contrast to Eq. (12) it can also be applied to an irregular parametric surface expressed as ρ = ρ (Φc , λ), i.e. also to the geoid. The partial derivatives in Eq. (14) we get from Eq. (9b)

−ρ33 e21 cos2 Φc sin λ0 cos λ0 ∂ρ3 = , λ0 = λ − λa , ∂λ a2 (1 − e21 ) { [ ]} −ρ33 sin Φc cos Φc e2 − e21 1 − (1 − e2 )cos2 λ0 ∂ρ3 = . ∂ Φc a2 (1 − e2 )(1 − e21 )

(15)

The metric tensor of a triaxial ellipsoid in geocentric coordinates is thus given by Eqs. (14), (9b) and (15) which are substantially more complicated than the analogous expression, Eq. (12), in reduced coordinates. The metric tensor gαβ characterizes the metrics of the above mentioned reference surfaces of the Earth. Thus it is present in formulae defining the geometric invariants on surfaces such as the angle ω between two curves p, q on the surface with tangent vectors p˙ = pα bα , q˙ = qβ bβ , the length s of the surface curve or the area A of certain part of the surface, see Budinský and Kepr [6] cos ω = (√

gαβ pα qβ gγ δ pγ pδ

√

gεη qε qη

), s =

t2

∫

√

gαβ duα duβ ,

A=

∫ ∫ √

g11 g22 − (g12 )2 du1 du2 .

(16)

t1

Eqs. (16) are invariant with respect to the choice of the surface parameters uα , thus we use them in the next section to demonstrate the measure of difference between the ellipsoid of revolution and the )triaxial ellipsoid. Integration ( ) boundaries ( in second Eq. (16) are defined by a parameter t of the curve y (t ) = x u1 (t ) , u2 (t ) on the surface x u1 , u2 . 2. Differences in metrics of the ellipsoid of revolution and the triaxial ellipsoid The measure of difference of two types of ellipsoids will be demonstrated using differentials of several chosen metric quantities. Numerical parameters of the ellipsoids used for our comparison are shown in Table 1. As a reference surface we chose the generally used geocentric ellipsoid of revolution WGS84 with two basic geometric parameters a, 1/f , see Moritz [7]. The parameters of the triaxial ellipsoid a, 1/f , 1/f1 , λa were taken from Burša and Pěč [2] derived as the average of several geopotential models (GEM 7, 8, 10B, GRIM 2, SE II). The symbols f , f1 in Table 1 represent the flattening of the meridian ellipse in geographical longitude λa and flattening of the equatorial ellipse respectively. 2.1. Differences between the lengths of meridian and parallel sections We assessed the differences of the ellipsoids shown in Table 1 using the length differences of meridian and parallel arcs, by different curvature of geographical parallels and non-orthogonality of coordinate frame. The triaxial ellipsoid is symmetric with respect to equator and also with respect to two meridian planes corresponding to geographical longitudes λa and λb . These symmetrical properties allow us to restrict our comparison to longitudes λ ∈ ⟨λa , λa + 90◦ ⟩ ⇒ λ0 ∈ ⟨0◦ , 90◦ ⟩ and northern latitudes Φc ∈ ⟨0◦ , 90◦ ⟩. The constant grid step ∆λ = 1◦ and ∆Φc = 1◦ was chosen for numerical computation. All integrals were solved numerically in Mathcad 14. The meridian arcs mψ on the ellipsoid of revolution represent one type of coordinate lines. The lengths of these arcs were computed using the simplified form of the surface integral in Eq. (16) mψ = sψ (λconst ) =

ψ

∫

√ 0

˜ gψ˜ ψ˜ dψ,

(17)

where ψ denotes an arbitrary type of latitude. The lengths of meridian arcs are not dependent on λ and for a particular type of latitude ψ have been computed in accordance with g22 , see Eq. (8), so we get βi

∫

R

c2

mβi =

+

E 2 sin2

β dβ =

0

mC(Φc )i = mGϕi =

√

βi

∫

√

c 2 + a2 e2 sin2 β dβ,

(18)

0 (Φc )i

∫

ρ2 +

0

ϕi

∫

√

(

∂ρ ∂ Φc

)2

√

dΦc = a 1 − e2

(Φc )i

∫ 0

√

1 1 − e2 cos2 Φc

+

1

e4 sin2 Φc cos2 Φc

1 − e2 cos2 Φc (1 − e2 cos2 Φc )2

dΦc ,

M(ϕ ) dϕ

0

ϕi

∫ = 0

 ]2 [ ]2 [  ae2 sin ϕ cos2 ϕ a sin ϕ ae2 sin2 ϕ cos ϕ a cos ϕ √ − + + (1 − e2 )2 dϕ. 3 1 3 1 (1 − e2 sin2 ϕ ) 2 (1 − e2 sin2 ϕ ) 2 (1 − e2 sin2 ϕ ) 2 (1 − e2 sin2 ϕ ) 2

Considering the fact that Eqs. (18) represent the length of the curve which is a geometrical quantity invariant with respect to the choice of the coordinate system, the equality mRβi = mC(Φc ) = mGϕi must hold. Eqs. (18) demonstrate the advantage i of reduced coordinates in simplification of mathematical formulae as compared to geocentric or geodetic coordinates. In

198

L. Husár et al. / Journal of Geometry and Physics 120 (2017) 192–207

Table 1 Basic parameters of chosen ellipsoids.

WGS 84 Triaxial ellipsoid

a

1/f = a/(a − c)

1/f1 = a/(a − b)

λa

6 378 137 m 6 378 172 m

298.257 297.776

92 000

345.1◦

Fig. 3. Differences between lengths of meridian arcs on the ellipsoid of revolution WGS84 and the triaxial ellipsoid defined in Table 1. The arcs were integrated (a) over the interval Φc ∈ (0◦ , Φc ) (b) over the 1◦ interval. Units of the vertical axes are metres. The numbers at the latitude axis represent the variable upper integration limit of compared arcs.

order to confirm the correctness of Eqs. (18), the equality of arc lengths expressed in all types of latitude has been tested numerically. The meridian arcs on a triaxial ellipsoid may correspond to different ellipses with the common minor semi-axis c and a variable major semi-axis aλ0 ∈ ⟨a, b⟩ which depends on the longitude λ0 , see Eqs. (10). The final expression for meridian arcs on a triaxial ellipsoid m3 can be obtained after substitution of ρ3 (aλ , λ0 ) and ∂ρ3 /∂ Φc from Eqs. (9b) and (15) into the general expression (Φc )i

∫

√

m3λj ,(Φc )i =

ρ32 +

0

(

∂ρ3 ∂ Φc

)2

dΦc .

(19)

The elliptic integral me 3 for particular meridian ellipses with semi-axes aλ , c, see Eq. (20), was used for verification of Eq. (19) (Φc )i

∫ me 3λj ,(Φc )i = 0

√

a2λ cos2 Φc + c 2 sin2 Φc dΦc

(20)

3e as well as the modified integral (19) expressed in reduced coordinates λr , β using the g22 (λr , β ) term from Eqs. (12) and integration limits from 0 to βi = f [(Φc )i ] given by second expression of Eqs. (10). The differences in lengths of meridian arcs of the ellipsoid of revolution and the triaxial ellipsoid, defined by parameters in Table 1, were computed as Eqs. (19)–(18) for a range of geographical longitudes from 0◦ to 90◦ with the step of 1◦ . For each longitude, we gradually extended the arc length from 1◦ to 90◦ with the step ∆Φc = 1◦ (Fig. 3). The extreme differences between the meridian arcs reached the maximum value of 27.77 m in longitude λ0 = 0◦ , integrated over the entire quadrant Φc ∈ ⟨0◦ , 90◦ ⟩, and minimum value of −26.75 m in longitude λ0 = 90◦ , integrated over the shorter interval Φc ∈ ⟨0◦ , 84◦ ⟩. Similarly, the differences of one-degree meridian arcs were also computed which varied from 0.61 m in (Φc = 0◦ , λ0 = 0◦ ) to −0.60 m in (Φc = 0◦ , λ0 = 90◦ ). These differences relative to the Earth size are ±5.5 · 10−6 which corresponds to 5.5 mm per 1 km. The meridian arc differences computed over one quadrant are plotted in Fig. 3. In comparison of parallel arcs the situation is more complicated due to three different parallel types defined on a triaxial ellipsoid according to particular types of the latitudes. Thus we distinguish reduced parallels (βi = const i ), geocentric parallels (Φc(i) = const i ), and geodetic parallels (ϕi = const i ). This fact is documented in Fig. 4 where there are several meridian ellipses with common minor semi-axis c and the major semi-axes aλ0 ∈ ⟨a, b⟩. The ellipses are rotated into the meridian plane λ = λa and the point P, defining a certain parallel by geocentric latitude Φc(i) = consti = Φc(P) , is chosen on the meridian λa . The projected part of the reduced parallel PR, depicted in Fig. 4 as a horizontal dashed line, is in reality a part of the ellipse in accordance with the definition of the reduced latitude. From the definition of the Cartesian coordinate x3(i) = c sin βi follows that if the reduced latitude on a reduced parallel is by definition

L. Husár et al. / Journal of Geometry and Physics 120 (2017) 192–207

199

Fig. 4. Projection of geocentric (PG), geodetic (Pg), and reduced (PR) parallel arcs containing the point P situated in meridian plane λ = λa .

constant, the x3 coordinate will also be constant. The projected part of the geocentric parallel is a spatial curve PG with the constant latitude Φc which is formed as an intersection of the triaxial ellipsoid surface with the conical surface having vertex O and vertex angle π − 2Φc . The geodetic parallel defined by ϕi = consti is a spatial curve Pg with constant latitude ϕ . The fact that geocentric and geodetic parallels on a triaxial ellipsoid are not planar curves is caused by a variable numerical eccentricity of particular meridian ellipses. The term parallel in these cases does not represent reality because spatial curves are not parallel with the equatorial plane. A more convenient term for these curves would be the latitude isoline. Lengths of reduced parallel arcs p were computed as parts of elliptic arcs using the following modification of the second expression in Eqs. (18) pβi ,j = Ai

√

1−

e2i

λ0 j

∫ 0

√

1 1 − e2i cos2 λ0

+

1

e4i sin2 λ0 cos2 λ0

1 − e2i cos2 λ0 (1 − e2i cos2 λ0 )2

dλ0 ,

(21)

with semi-axes Ai , Bi and numerical eccentricities ei of the parallel ellipses in the ith reduced latitude βi : Ai = a cos βi , Bi = b cos βi , e2i = (A2i − B2i )/A2i . For the numerical computations√we used reduced latitudes βi generated on the zero meridian

λ0 = 0 from the first expression of Eqs. (4): tg βi = tg (Φc )i / 1 − e20 , e20 = e2 = (a2 − c 2 )/a2 where (Φc )i ∈ ⟨0◦ , 90◦ ⟩ with the interval 1◦ . 3e The lengths of geocentric parallel arcs are defined using the first diagonal component of the metric tensor g11 from Eq. (14) by formula λ0j

∫ pΦc(i),j =

√ 2 ρ3(i) cos2 (Φc )i +

0

(

∂ρ3 ∂λ

)2

dλ,

(22)

where ρ3 and ∂ρ3 /∂λ can be obtained from Eqs. (9b) and (15). Lengths of geodetic parallel arcs have not been computed due to the complexity of the function ρ3 (ϕ ). However, we partially describe their spatial behaviour by variation of the Cartesian coordinate x3 . Lengths of reduced and geocentric parallel arcs on a triaxial ellipsoid were compared with the lengths of circular parallel arcs pi,j on the ellipsoid of revolution WGS84 pi,j = ρ (λ, Φc(i) ) cos Φc(i) λ0j .

(23)

We computed and analysed the three types of differences

∆pβi,j = pβi ,j − pi,j , ∆pΦi,j = p(Φc )i ,j − pi,j , ∆pΦ βi,j = p(Φc )i ,j − pβi ,j .

(24)

The maximum difference ∆pβmax = 17.60 m in reduced parallel arc was detected on the equator over the interval ∆λ0 ∈ ⟨0◦ , 45◦ ⟩, its minimum ∆pβmin = −20.66 m at latitude β = 55◦ over the interval ∆λ0 ∈ ⟨0◦ , 90◦ ⟩. The range of the ∆pΦ was from zero to ∆pΦmax = 17, 60 m on the equator over the interval ∆λ0 ∈ ⟨0◦ , 45◦ ⟩. The difference ∆pΦ β reached its maximum value ∆pΦ βmax = 20, 96 m at latitude Φc = 55◦ over the interval ∆λ0 ∈ ⟨0◦ , 90◦ ⟩. Graphical representations of the differences defined by Eqs. (24) are shown in Fig. 5. As the lengths of parallel arcs decrease with growing latitude due to convergence of meridians towards the poles, Fig. 5 is distorted. Therefore we decided to compute also the relative differences ∆re β = ∆pβi,j /pi,j and ∆re pΦi,j = pΦ c(i) i ,j /pi,j . The extreme values of the relative differences were similar to the relative meridian differences, i.e. ±5.5 · 10−6 . Deviations of geocentric (C) and geodetic (G) parallels on a triaxial ellipsoid from the flatness manifest as a variation of the Cartesian coordinate x3 which is constant in case of the parallels on an ellipsoid of revolution. We computed the differences

200

L. Husár et al. / Journal of Geometry and Physics 120 (2017) 192–207

Fig. 5. Length differences of parallel arcs between the triaxial ellipsoid and the ellipsoid of revolution WGS84 computed over the interval ∆λ0 ∈ (0◦ , λ0 ). (a) Reduced parallels, (b) geocentric parallels, (c) reduced versus geocentric on a triaxial ellipsoid. Units: metres.

of x3 -coordinates between the jth point and the first point of particular parallel with the longitude step of 1◦

∆xC3(j,i) = xC3(j,i) − xC3(1,i) , ∆xG3(j,i) = xG3(j,i) − xG3(1,i) ,

xC3(j,i) = ρ3 (λ0j , Φci ) sin(Φci ), ϕ ϕ xG3(j,i) = ρ3 (λ0j , Φc i ) sin(Φc i )j .

(25)

In each latitude the 90 parallel arcs were compared from the shortest arc (λ0(2) = 1◦ ) to the longest one (λ0(91) = 90◦ ), with the extension step of 1◦ . The first point was selected at the geographic longitude λ0(1) = 0◦ . The parallels were calculated in 1◦ step of the geocentric latitude Φc . In Eqs. (25) Φci is the ith geocentric latitude measured along the meridian λ0 = 0◦ and ρ3 (λ0j , Φci ) is the radius vector ϕ of the jth point on ith parallel. Similarly (Φc i )j are the geocentric latitudes of the jth point on ith geodetic parallel with the ϕ constant geodetic latitude ϕi = consti and ρ3 (λ0j , Φc i ) is the corresponding geocentric radius vector computed according ϕ to Eq. (9b). For all geocentric latitudes on the 0th meridian (Φc i )j=1 the corresponding geodetic latitudes ϕi were generated ϕ ◦ according to simplified Eq. (26) for λ01 = 0 : tan ϕi = tan Φci /(1 − e2 ). Consequently the geocentric latitudes (Φc i )j of all 90 points on geodetic parallel have been computed from Eq. (26) ϕ

tan(Φc i )j = Qj tan ϕi ,

Qj =

√

(1 − e2 )2 cos2 λ0j + (1 − e22 )2 sin2 λ0j .

(26)

Thus we established that all three kinds of parallels pass through the identical point (λ01 , Φci ). However, with the growing longitude λ0j both the geocentric and the geodetic parallel gradually straggle from the plane defined by the reduced parallel xi3 = consti . For derivation of Eq. (26) we used the dual expression of the x3 coordinate, see Fig. 8 x3 = ρ3 sin Φc = ⏐P1⏐ sin ϕ,

⏐

⏐

(27)

L. Husár et al. / Journal of Geometry and Physics 120 (2017) 192–207

201

Fig. 6. Spatiality of geocentric (∆xC3 /m) and geodetic (∆xG3 /m) parallels.

where we used the coordinates of the point P [x1 , x2 ,⏐x3 ]⏐and of the point 1 x1 e2 , x2 e22 , 0 (see Fig. 8 and Eq. (32)). After the simple manipulation we get the length of the vector ⏐P1⏐ in the following form

[

]

√ ⏐ ⏐ ⏐P1⏐ = ρ3 cos Φc (1 − e2 )2 cos2 λ0 + (1 − e2 )2 sin2 λ0 + tan2 Φc .

(28)

2

Eq. (26) results from the substitution of Eq. (28) into Eq. (27) after the simple arrangement. The deviations of xG3 and xC3 coordinates from the constant value are symmetric to the planes parallel to the equator, see Fig. 6, and the extreme values ∆xG3(max) = +26.594 m and ∆xC3(min) = −26.594 m are located at Φc = 35◦ , λ0 = 90◦ . Spatiality of the geocentric parallels on a triaxial ellipsoid causes the non-orthogonality of the geocentric coordinate lines. Geocentric parallels are not perpendicular to the meridians. We computed the deviations ∆ω of the angle between geocentric parallels and meridians from the right angle according to the formulae

∆ωΦc(i), λ0j =

π 2

− ωΦc(i) ,λ0j ,

) ( g12 cos ωΦc(i), λ0j = √ √ , g11 g22

(29)

where the components of the metric tensor gαβ were evaluated according to Eq. (14). As the values of ∆ω are too small, we decided to present them in milli-arc-seconds (mas = 0.001′′ ). The relation between ∆ω and the geographical location is shown in Fig. 7. The maximum value ∆ω = 5.8 mas appears at Φc = 35◦ , λ0 = 44◦ . 2.2. Differences in curvature of the ellipsoid of revolution and the triaxial ellipsoid A suitable quantity characterizing the discrepancies between the surfaces of an ellipsoid of revolution and a triaxial ellipsoid is normal curvature (curvature of the normal section of the surface in a particular azimuth). Let us assign the extreme values of the normal curvature as κmax , κmin , the mean curvature (their average) as κS and the normal curvature of the section slewed by specific angle θ as κθ . The quantity κθ is defined by the Euler formula [6]

κs =

1 2

(κmax + κmin ) ,

Rs =

1

κs

,

κθ = κmax cos2 θ + κmin sin2 θ,

Rθ =

1

κθ

.

(30)

On an ellipsoid of revolution, the angle θ is identical with the azimuth and the extreme normal curvatures are the meridian and prime vertical curvatures κmax = κM , κmin = κN , to which the meridian and prime vertical radii of curvature M, N correspond, see Eqs. (6b) and (8). To compare the curvatures of two surfaces more demonstratively we chose the metric quantity, the curvature radius Rµ , (µ = s, θ ). The situation on a triaxial ellipsoid is more complicated because although the meridian planes are perpendicular to the equator (they all contain the x3 axis) they cannot be considered as normal planes because they do not contain the ellipsoidal normal. The ellipsoidal normals are deflected from the meridian planes as a consequence of the inequality of the equatorial semi-axes b ̸ = a, see Fig. 8. Ellipsoidal normal n passing through the point P is deflected from the meridian plane (which contains the axis x3 and intersects the equatorial plane x1 x2 in a trace me ) and intersects the equator in point 1, the plane x1 x3 in point 2 and plane

202

L. Husár et al. / Journal of Geometry and Physics 120 (2017) 192–207

Fig. 7. Deviations from orthogonality of geocentric coordinate lines presented by deviations ∆ω/mas from the right angle.

Fig. 8. Normal n on a triaxial ellipsoid.

x2 x3 in point 3, see Fig. 8. The direction of the normal is given by the gradient of the f3e function (9a)

( n = ∇ f3e =

∂ f3e ∂ f3e ∂ f3e , , ∂ x1 ∂ x2 ∂ x3

)T

( ( x x x )T 2 1 2 3 = 2 2 , 2 , 2 = 2 x1 , a

b

c

a

x2

, 2

x3

)T

1 − e1 1 − e2

,

(31)

where symbol T represents the transposition. The coordinates of the points 1, 2, 3 on a normal n, see Fig. 8, were determined by Feltens [3]

)T

x(1) = x1 e2 , x2 e22 , 0 ,

(

2

e′ =

(a2 − c 2 )

,

2

e′ 1 =

(

2

x(2) = x1 e21 , 0, −x3 e′ 2 (a2 − b2 )

,

e22 =

)T

,

(b2 − c 2 )

(

2

x(4) = 0, −x2 e′ 1 , −x3 e′

,

2

e′ 2 =

2

)T

,

(32)

(b2 − c 2 )

c2 b2 b2 c2 and were used to test the normal computed by Eq. (31) by expressing the normal vector n by three couples of points P1 : n = x − x(1) , P2 : n = x − x(2) and P3 : n = x − x(3) where x, x(1) − x(4) are the geocentric radius vectors of

L. Husár et al. / Journal of Geometry and Physics 120 (2017) 192–207

203

the points P, 1, 2 and 3. As the vectors P1, P2 and P3 are collinear and differ only in absolute value for the test purpose they were transformed to unit vectors. The points 1, 2, 3 can also serve for a quantification of the deviation of the normal from the meridian plane. For this purpose we used the angle ε , see Fig. 8, between the meridian plane passing through the point P and the so called normal meridian plane ν0 which is the only normal plane of P perpendicular to equatorial plane and therefore suitable for the azimuth definition. Note that the geodetic azimuth α on a triaxial ellipsoid is defined on a local horizontal plane of the point P as an angle from the normal meridian plane parallel to the axis x3 to the normal plane containing the target, in accordance with the astronomical azimuth definition [8, p. 19]. Both, meridian plane and normal meridian plane ν0 are perpendicular to the equatorial plane and intersect in the line PP0 . The vectors me , ne defined as x1 x2 0

( ) me =

) ρ3 cos λ cos Φc = ρ3 sin λ cos Φc , (

n1 n2 0

( ) ne =

0

⎛

⎞

x1 2 ⎜ x2 ⎟ = 2⎝ ⎠ 1 − e21 a 0

(33)

also form an angle ε . When only the latitude is changing, the moving point still remains on the same meridian and the vector me remains constant, which results from the slope of the tangent x2 /x1 = tan λ0 which is independent on the latitude. Similarly the vector ne does not depend on latitude because its slope of the tangent, see Eq. (33), is only a function of longitude: n2 /n1 = tan λ0 /(1 − e21 ). The angle ε is thus solely a function of λ0 and can be determined either from the inner product of the vectors me , ne or from the difference of the slopes of tangents

ε = arccos

(

me · n e

|me | |ne |

)

,

ε = arctan

(

tan λ0 1 − e21

)

− λ0 .

(34)

From Eq. (34) follows that ε gets the maximum value of 2.24′′ at latitude 45◦ which represents the extreme deviation of the normal meridian plane from the meridian on the triaxial ellipsoid defined in Table 1. As the azimuth on both ellipsoids is defined from different meridians we chose the angle θ , known from the Euler formula Eq. (30), to be the convenient quantity for a comparison of normal curvature. The angle θ is measured from the normal section with the maximal curvature. At first we determined the main directions with extreme curvatures, defined by the contravariant coordinates du γ , (γ = 1, 2) of the tangent vector to the normal section lying in the tangential plane to the ellipsoid, using the condition for the main direction [6, p. 213]

⏐ ⏐(du1 )2 ⏐ ⏐ g22 ⏐ ⏐ h22

du1 du2 g12 h12

⏐

(du2 )2 ⏐⏐ g11 ⏐⏐ = 0, h11 ⏐

(35)

which after some manipulation gives us the solution of unknown roots σ1,2 = (du1 /du2 )1,2 of a quadratic equation dependent only on the first and second surface tensors gαβ , hαβ

√ (g11 h22 − g22 h11 )2 − 4 (g11 h12 − g12 h11 ) (g12 h22 − g22 h12 ) σ1,2 = . (36) 2 (g11 h12 − g12 h11 ) The extreme (main) normal curvatures κmax , κmin occur in main directions and the final expression for their computation is obtained by introducing the main directions σ1,2 into a general formula for the normal curvature [6, p. 213] − (g11 h22 − g22 h11 ) ±

κmax,min =

h11 (σ1,2 )2 + 2h12 σ1,2 + h22 g11 (σ1,2 )2 + 2g12 σ1,2 + g22

.

(37)

A development of the normal curvature κθ depending on the direction θ at selected points on the ellipsoid was determined by the Euler formula, Eq. (30). Fig. 9 shows the difference of the curvature radii which correspond to the mean curvatures W ◦ ◦ ∆Rs = R3e s − Rs on the triaxial and WGS84 ellipsoids. The differences were computed by Eq. (30) on a regular grid 1 × 1 ◦ ◦ for northern latitudes and longitudes λ0 ∈ ⟨0 , 90 ⟩. Extreme differences occurred on the equator: maximum 69.303 m at longitude λ0 = 90◦ , minimum −68.655 m at longitude λ0 = 0◦ while towards the pole they inclined to zero. The differences of curvature radii which correspond to the normal curvatures κθ depending on the angle θ were computed on three chosen meridians λ0 = 0◦ (∆Ra), 45◦ (∆R45 ) and 90◦ (∆Rb) at a latitude interval ⟨0◦ , 90◦ ⟩, see Fig. 10. The values vary from −103.655 m on the equator and a longitude λ0 = 90◦ on the first vertical to 104.329 m on the equator and a longitude λ0 = 0◦ also on the first vertical. 3. Contribution of a triaxial ellipsoid to a refinement of the Earth reference body Despite of relatively small differences in several above mentioned metric quantities between a biaxial and a triaxial ellipsoid we could expect the more significant contribution of a triaxial ellipsoid in the overall fit of the shape of the planet Earth. This should manifest itself with the better conformity of the triaxial ellipsoid with the geoid, representing the shape of the real Earth. The comparison has been made on a regular geographical grid 1◦ ×1◦ covering almost entire Earth, except the polar areas above the 80◦ and below the −80◦ of geocentric latitude. For the evaluation of the geoid undulation ∆hWGS with respect to the biaxial ellipsoid WGS84 we used the Graflab software [9] and EIGEN6c4 global geoid model up to spherical

204

L. Husár et al. / Journal of Geometry and Physics 120 (2017) 192–207

Fig. 9. Differences of the curvature radii ∆Rs /m corresponding to the mean curvatures on the triaxial and WGS84 ellipsoids.

Table 2 Basic statistics of geoidal heights computed using EIGEN6c4 global geopotential model with respect to the biaxial ellipsoid WGS84 (∆hWGS ), triaxial ellipsoid of Burša (∆h3 ) and differences between the geocentric radii of both ellipsoids ∆ρ3WGS on a nearly global regular grid 1◦ ×1◦ . Min Max Mean Range

∆hWGS (m)

∆h3 (m)

∆ρ3WGS (m)

−106.171 84.937 −0.371 191.109

−72.431 67.489 −0.701 139.921

−34.328 35.000 0.330 69.328

harmonic degree 2160 [10]. The differences between the geocentric radius vectors of the triaxial and biaxial ellipsoids ∆ρ3WGS = ρ3 − ρ have been computed using Eqs. (5) and (9b) and are depicted in Fig. 11. The geoid heights with respect to the triaxial ellipsoid ∆h3 have been computed by subtracting the difference ∆ρ3WGS from the geoid height ∆hWGS . The geoid heights with respect to both ellipsoids, WGS84 and the triaxial, are shown in Fig. 12. The basic statistical parameters of ∆hWGS , ∆h3 and ∆ρ3WGS are shown in Table 2. The range of the geoid undulation using the triaxial instead of biaxial ellipsoid decreased of about 51 metres which corresponds to 26%, see Table 2. The minimum located in the Indian Ocean in absolute value decreased of 33.740 m and the maximum located in Western Pacific Ocean decreased of 17.448 m.

4. Conclusion

The aim of this paper was to quantify the differences between the most common reference surfaces of the Earth body. In order to express the differences in metrics and curvature we used the tensor quantities, mainly the metric tensor and the second surface tensor, which were expressed in various curvilinear coordinates. Suitability of certain type of coordinates was demonstrated by simplification of the above mentioned tensors and thereby simplification of mathematical operations on the surfaces. From the length and angular differences between a biaxial and triaxial ellipsoid follow that for certain applications (cadastre, mapping, engineering surveying, cartography) they are negligible. In these applications the better fit of the triaxial ellipsoid would not be evident and its use would not be adequate due to computational cost. This confirms the justification of a biaxial ellipsoid as a general conventional reference body of the Earth. The extreme values of differences between the length of meridians and also parallels are 5.5 mm per 1 km. Very small equatorial flattening of the Earth body, less than

L. Husár et al. / Journal of Geometry and Physics 120 (2017) 192–207

205

Fig. 10. Differences of the curvature radii of normal sections depending on the angle θ (a) on a meridian λ0 = 0◦ , (b) on a meridian λ0 = 90◦ , (c) on a meridian λ0 = 45◦ and (d) is a detailed view of (c). The value of geocentric latitude is assigned to each curve.

70 m in absolute value, can be revealed only by accurate satellite measurements. This equatorial flattening causes the nonorthogonality of parallels and meridians which can reach values up to several mas. However, the triaxial ellipsoid reduces the geoid heights and its variability for more than 26%. Thus for some specific scientific tasks the triaxial ellipsoid can be useful, mainly in satellite geodesy, physical geodesy and space geodesy. In satellite geodesy, the triaxiality is responsible for perturbations in the precession of the node of the satellite orbit [1]. The equatorial flattening of the Earth also reveals some insight to solid Earth geophysics connected to deep Earth-mass structure. The main utilization of a triaxial ellipsoid is in planetary geodesy for modelling of extraterrestrial objects (planets, their natural satellites and asteroids) with the striking shape anomalies [11–15] or for a refinement of the motion equations and dynamic properties of these bodies [16,17].

Acknowledgements This paper was created with the support of Grant No. 1/0642/13 of the VEGA Grant Agency of the Slovak Republic and with the support of the Ministry of Education, Science and Sport of the Slovak Republic within the Research and Development Operational Programme for the project ‘‘University Science Park of STU Bratislava,’’ ITMS 26240220084, co-founded by the European Regional Development Fund.

206

L. Husár et al. / Journal of Geometry and Physics 120 (2017) 192–207

40 m 35 m 30 m 25 m

50

20 m 15 m 10 m

0

5m 0m -5 m

-50

-10 m -15 m -20 m

0

50

100

150

200

250

300

350

-25 m -30 m -35 m

Fig. 11. Differences ∆ρ3WGS = ρ3 − ρ between the biaxial and triaxial ellipsoids.

50

0

-50

0

50

100

150

200

250

300

350

0

50

100

150

200

250

300

350

50

0

90 m 80 m 70 m 60 m 50 m 40 m 30 m 20 m 10 m 0m -10 m -20 m -30 m -40 m -50 m -60 m -70 m -80 m -90 m -100 m -110 m

-50

Fig. 12. Geoid heights with respect to the biaxial ellipsoid WGS84 (top) and the Burša triaxial ellipsoid (bottom).

References [1] W.A. Heiskanen, Is the earth a triaxial ellipsoid? J. Geophys. Res. 67 (1) (1962) 321–327. [2] M. Burša, K. Pěč, Gravity Field and Dynamics of the Earth, Springer-Verlag, Berlin, 1993. (Translation from the Czech original: Tíhové pole a dynamika Zem˘e, Academia, Prague, 1988). [3] J. Feltens, Vector method to compute the cartesian (X , Y , Z ) to geodetic (Φ , λ, h) transformation on a triaxial elipsoid, J. Geod. 83 (2009) 129–137. [4] M.R. Spiegel, Theory and Problems of Vector Analysis and an Introduction to Tensor Analysis SI (Metric), McGraw-Hill International Book Company, London: McGraw-Hill Book Company (UK), Dűsseldorf: McGraw-Hill Book Company GmbH, McGraw-Hill Book Company Australia Pty, Limited, Johannesburg: McGraw-Hill Book Company (SA) (Pty) Limited, New York, 1992.

L. Husár et al. / Journal of Geometry and Physics 120 (2017) 192–207

207

[5] W. Torge, Geodesy, second ed., De Gruyter, Berlin, 1991. [6] B. Budinský, B. Kepr, Basics of Differential Geometry with Engineering Applications. (Základy Diferenciální geometrie s technickými aplikacemi), SNTL, Prague, 1970 (in Czech). [7] H. Moritz, Geodetic reference system 1980, in: Bulletin Geodesique, The Geodesist’s Handbook, vol. 62, No. 3, Anne 1988. [8] I.I. Mueller, Spherical and Practical Astronomy As Applied to Geodesy, Frederick Ungar Publishing Co., Inc., New York, 1969. [9] B. Bucha, J. Janák, A MATLAB-based graphical user interface program for computing functionals of the geopotential up to ultra-high degrees and orders, Comput. Geosci. 56 (2013) 186–196. [10] Ch. Förste, S.L. Bruinsma, O. Abrikosov, J.M. Lemoine, T. Schaller, H.J. Götze, J. Ebbing, J.C. Marty, F. Flechtner, G. Balmino, R. Biancale, EIGEN-6C4 The latest combined global gravity field model including GOCE data up to degree and order 2190 of GFZ Potsdam and GRGS Toulouse, 2015. http://dx.doi.org/105880/icgem.20151. [11] W.R. Wang, F. Li, J.J. Liu, X. Ren, X.D. Zou, L.L. Mu, J.G. Yan, J.L. Zou, H.B. Zhang, Ch. Lu, J.Z. Liu, W. Zuo, Y. Su, W.B. Wen, W. Bian, M. Wang, Ch.L. Li, Z.Y. Ouyang, Triaxial ellipsoid models of the moon based on the laser altimetry data of chang’e-1, Sci. China Earth Sci. 53 (11) (2010) 1594–1601. [12] Y. Ohba, M. Abe, S. Hasegawa, M. Ishiguro, T. Kwiatkowski, F. Colas, B. Dermawan, A. Fujiwara, Pole orientation and triaxial ellipsoid shape of (25143) 1998 SF36, a target asteroid of the MUSES- C mission, Earth Planets Space 55 (2003) 341–347. [13] K.B. Bugaevsky, L.M. Shingareva, Triaxial ellipsoid as a reference-surface on mars mapping, in: The Fifth International Conference on Mars, July 19–24, 1999, Pasadena, California, abstract no. 6219. [14] J.D. Drummond, A. Conrad, J.W. Merline, B. Carry, C.R. Chapman, H.A. Weaver, P.M. Tamblyn, J.C. Christou, C. Dumas, The triaxial ellipsoid dimensions, rotational pole, and bulk density of esa rosetta target asteroid (21) lutetia, A&A 523 (2010) A93. [15] J.D. Drummond, J.W. Merline, A. Conrad, C. Dumas, P. Tamblyn, J. Christou, B. Carry, C. Chapman, The Triaxial Ellipsoid Diameters and Rotational Pole of Asteroid (9) Metis from AO at Gemini and Keck, American Astronomical Society, DPS meeting #44, #30209, 2012. [16] O. Vasilkova, Three-dimensional periodic motion in the vicinity of the equilibrium points of an asteroid, A&A 430 (2005) 713–723. http://dx.doi.org/ 10.1051/0004-6361:20034414. [17] J. Bellerose, J.S. Daniel, Restricted Full Three-Body Problem: Application to Binary System 1999 KW4, J. Guid. Control Dyn. 31 (1) (2008).