2. The forward and inverse expansions of the meridian arc in geometric geodesy

The forward and inverse problem of the meridian arc is one of the fundamental problems in geometric geodesy, which seems to be a solved one. Briefly reviewing the existing methods, however, one will find that the inverse problem was not perfectly and ideally solved yet. This situation is due to the complexity of the problem itself and the lack of advanced computer algebra systems. Yang had given the direct expansions of the inverse transformation by means of the Lagrange series method, but their coefficients are expressed as polynomials of coefficients of the forward expansions, which are not convenient for practical use, see [6, 7]. Adams expressed the coefficients of inverse expansions as a power series of the eccentricity e by hand, but expanded them up to eighth-order terms of e at most, see [1]. Due to these reasons, the forward and inverse expansions of the meridian arc are discussed by means of Mathematica in the following sections. Their coefficients are uniformly expressed as a power series of the eccentricity and extended up to tenth-order terms of e.

#### 2.1. The forward expansion of the meridian arc

The meridian arc from the equator where the latitude is from B ¼ 0 to B is

$$X = a \left( 1 - e^2 \right) \int\_0^B \left( 1 - e^2 \sin^2 B \right)^{-3/2} dB \tag{1}$$

where X is the meridian arc; Bis the geodetic latitude; a is the semi-major axis of the reference ellipsoid; e is the first eccentricity of the reference ellipsoid.

Expanding the integrand in Eq. (1) and integrating it item by item using Mathematica as follows:

$$\mathbf{S} = \mathbf{S} \mathbf{series} \left[ \left( \mathbf{1} - \mathbf{e}^2 \star \mathbf{S} \mathbf{in} \left[ \mathbf{B} \right]^2 \right)^{-3/2}, \left\{ \mathbf{e}, \begin{array}{c} \mathbf{0}, \ \mathbf{0} \end{array} \mathbf{10} \right\} \right]$$

$$\frac{1}{12} + \frac{3}{2} \cdot \mathfrak{sl} \dot{\mathfrak{su}} \cdot \mathfrak{v} \Big|^{2} \bullet ^{2} + \frac{15}{8} \cdot \mathfrak{sl} \dot{\mathfrak{su}} \cdot \mathfrak{v} \Big|^{4} \bullet ^{4} + \frac{35}{16} \cdot \mathfrak{sl} \dot{\mathfrak{su}} \cdot \mathfrak{v} \Big|^{6} \bullet ^{6} + \frac{315}{128} \cdot \mathfrak{sl} \dot{\mathfrak{su}} \cdot \mathfrak{v} \Big|^{6} \bullet ^{6} + \frac{693}{256} \cdot \mathfrak{sl} \dot{\mathfrak{su}} \cdot \mathfrak{v} \Big|^{20} \bullet ^{10} + \mathfrak{o} \circ \mathfrak{o} \Big|^{110} \bullet ^{120} + \mathfrak{o} \circ \mathfrak{so} \Big|^{211} \bullet ^{220} + \mathfrak{o} \circ \mathfrak{so} \Big|^{222} \bullet ^{233} + \mathfrak{o} \circ \mathfrak{so} \Big|^{234} \bullet ^{24} + \mathfrak{o} \circ \mathfrak{so} \Big|^{243} \bullet ^{25} + \mathfrak{o} \circ \mathfrak{so} \Big|^{25} \bullet ^{26} + \mathfrak{o} \circ \mathfrak{so} \Big|^{27} \bullet ^{28} + \mathfrak{o} \circ \mathfrak{so} \Big|^{28} \bullet ^{29} + \mathfrak{o} \circ \mathfrak{so} \Big|^{29} \bullet ^{29} + \mathfrak{o} \circ \mathfrak{so} \Big| ^{210} \bullet ^{20} + \mathfrak{o} \circ \mathfrak{so} \Big| ^{211} \bullet ^{20} + \mathfrak{o} \circ \mathfrak{so} \Big| ^{212} \bullet ^{20} + \mathfrak{o} \circ \mathfrak{so} \Big| ^$$

$$\begin{aligned} \text{10 } &+\frac{3}{4} \left( \text{3 } - \text{Go} \text{ [B] } \text{ 3\'in [B] } \right) \text{ } \bullet \text{ }^2 + \frac{125}{246} \left( \text{12 } \text{в} - \text{@\\_@\\_[12 \text{B}] } + \text{@\\_[14 \text{B}] } \right) \text{ } \bullet \text{ }^4 - \\ \hline \text{35 } &(-6 \text{ } \text{B} + 4 \text{B} \text{ \sin \{2 \text{B} \} - \text{ } \text{B \sin \{4 \text{B} \} + \text{@\\_[15 \text{B}] } \text{ } \bullet \text{ }^4 - \\ \hline \text{103 } &(4 \text{B} \text{ } \text{B} - 6 \text{ } \text{B} \text{ \sin \{2 \text{B} \} + 16 \text{B} \text{ \sin \{4 \text{B} \} - 32 \text{ } \text{\sin \{6 \text{B} \} + 3 \text{ } \text{\sin \{6 \text{B} \} \} \text{ } \bullet \text{ }^4 + \frac{1}{2 \text{B} \text{ } \text{2$$

Then one arrives at

$$X = a(1 - e^2)(K\_0 B + K\_2 \sin 2B + K\_4 \sin 4B + K\_6 \sin 6B + K\_8 \sin 8B + K\_{10} \sin 10B) \tag{2}$$

where

transformation. Many geodesists have made great efforts to solve these problems, see [1–8]. Due to historical condition limitation, they mainly disposed of these problems by hand, which were not perfectly solved yet. Traditional algorithms derived by hand mainly have the following problems: (1) The expressions are complex and lengthy, which makes the computation process very complicated and time-consuming. (2) Some approximate disposal is adopted, which influences the computation accuracy. (3) Some formulae are numerical and only apply

In computational mathematics, computer algebra, also called symbolic computation, is a scientific area that refers to the study and development of algorithms and software for manipulating mathematical expressions and other mathematical objects. Software applications that perform symbolic calculations are called computer algebra systems, which are more popular today. Computer algebra systems, like Mathematica, Maple and Mathcad, are powerful tools to solve some mathematical derivation in geodesy, see [9–11]. By means of computer algebra analysis method and computer algebra system Mathematica, we have already solved many complicated mathematical problems in special fields of geodesy in the past few years; see [12–15].

The main contents and research results presented in this chapter are organized as follows: In Section 2, we discuss the forward and inverse expansions of the meridian arc often used in geometric geodesy. In Section 3, the nonsingular expressions of singular integration in physical geodesy are derived. In Section 4, we discuss series expansions of direct transformations between three anomalies in satellite geodesy. Finally in Section 5, we make a brief summary.

2. The forward and inverse expansions of the meridian arc in geometric

The forward and inverse problem of the meridian arc is one of the fundamental problems in geometric geodesy, which seems to be a solved one. Briefly reviewing the existing methods, however, one will find that the inverse problem was not perfectly and ideally solved yet. This situation is due to the complexity of the problem itself and the lack of advanced computer algebra systems. Yang had given the direct expansions of the inverse transformation by means of the Lagrange series method, but their coefficients are expressed as polynomials of coefficients of the forward expansions, which are not convenient for practical use, see [6, 7]. Adams expressed the coefficients of inverse expansions as a power series of the eccentricity e by hand, but expanded them up to eighth-order terms of e at most, see [1]. Due to these reasons, the forward and inverse expansions of the meridian arc are discussed by means of Mathematica in the following sections. Their coefficients are uniformly expressed as a power series of the

eccentricity and extended up to tenth-order terms of e.

The meridian arc from the equator where the latitude is from B ¼ 0 to B is

<sup>2</sup> � � <sup>ð</sup><sup>B</sup>

0

1 � e

<sup>2</sup> sin <sup>2</sup> B � ��3=<sup>2</sup>

dB (1)

X ¼ a 1 � e

2.1. The forward expansion of the meridian arc

geodesy

to a specific reference ellipsoid, which are not convenient to be generalized.

68 Trends in Geomatics - An Earth Science Perspective

$$\begin{cases} K\_0 = 1 + \frac{3}{4}e^2 + \frac{45}{64}e^4 + \frac{175}{256}e^6 + \frac{11025}{16384}e^8 + \frac{43659}{65536}e^{10} \\ K\_2 = -\frac{3}{8}e^2 - \frac{15}{32}e^4 - \frac{525}{1024}e^6 - \frac{2205}{4996}e^8 - \frac{72765}{131072}e^{10} \\ K\_4 = \frac{15}{256}e^4 + \frac{105}{1024}e^6 + \frac{2205}{16384}e^8 + \frac{10395}{65536}e^{10} \\ K\_6 = -\frac{35}{3072}e^6 - \frac{105}{4096}e^8 - \frac{10395}{262144}e^{10} \\ K\_8 = \frac{315}{131072}e^8 + \frac{3465}{524288}e^{10} \\ K\_{10} = -\frac{693}{1310720}e^{10} \end{cases} (3)$$

Eqs. (2) and (3) can also be derived using the binomial theorem by hand which consumes much more time, see [6–9]. The denominator 693 of the last coefficient K<sup>10</sup> is mistaken as 639 in Ref. [9].

#### 2.2. The inverse expansion of the meridian arc using the Hermite interpolation method

Differentiation to the both sides of Eq. (1) yields:

$$\frac{dX}{dB} = \frac{a\left(1 - e^2\right)}{\left(1 - e^2 \sin^2 B\right)^{3/2}}\tag{4}$$

<sup>B</sup>ð Þ<sup>5</sup> ð Þ¼ <sup>0</sup> <sup>12</sup><sup>e</sup>

0

BBBBBB@

The solution to Eq. (13) is

0

BBBBBB@

Suppose that

Eq. (16) is

a2 a4 a6 a8 a<sup>10</sup> 1

0

BBBBBB@

<sup>a</sup><sup>2</sup> <sup>¼</sup> <sup>3</sup> 8 e 2 þ 3 16 e 4 þ 213 <sup>2048</sup> <sup>e</sup> 6 þ 255 <sup>4096</sup> <sup>e</sup> 8 þ

8

>>>>>>>>>>>>>><

>>>>>>>>>>>>>>:

<sup>a</sup><sup>4</sup> <sup>¼</sup> <sup>21</sup> <sup>256</sup> <sup>e</sup> 4 þ 21 <sup>256</sup> <sup>e</sup> 6 þ 533 <sup>8192</sup> <sup>e</sup> 8 þ 197 <sup>4096</sup> <sup>e</sup> 10

<sup>a</sup><sup>6</sup> <sup>¼</sup> <sup>151</sup> <sup>6144</sup> <sup>e</sup> 6 þ 151 <sup>4096</sup> <sup>e</sup> 8 þ

<sup>a</sup><sup>8</sup> <sup>¼</sup> <sup>1097</sup> <sup>131072</sup> <sup>e</sup>

<sup>a</sup><sup>10</sup> <sup>¼</sup> <sup>8011</sup> <sup>2621440</sup> <sup>e</sup>

CCCCCCA ¼ <sup>2</sup> <sup>þ</sup> <sup>90</sup><sup>e</sup>

2 4 6 8 10 �2 4 �6 8 �10 �8 �64 �216 �512 �1000 8 �64 �216 �512 1000 32 1024 7776 32768 100000

4 þ 4455 <sup>16</sup> <sup>e</sup> 6 þ

correspondingly, one arrives at a set of linear equations for the unknown coefficients

2 4 6 8 10 �2 4 �6 8 �10 �8 �64 �216 �512 �1000 8 �64 �216 �512 1000 32 1024 7776 32768 100000

Omitting the main operations by means of Mathematica, one arrives at

8 þ

10

2.3. The inverse expansion of the meridian arc using Lagrange's theorem method

with j j f xð Þ << j j x and y ≈ x. Lagrange's theorem states that in a suitable domain the solution of

1097 <sup>65536</sup> <sup>e</sup> 10

Making use of the five interpolation constraints in Eqs. (8–12) and differentiating Eq. (7)

20145 <sup>32</sup> <sup>e</sup> 8 þ

Mathematical Analysis of Some Typical Problems in Geodesy by Means of Computer Algebra

1

0

BBBBBB@

a2 a4 a6 a8 a<sup>10</sup> 1

CCCCCCA ¼

1

�<sup>1</sup> <sup>B</sup><sup>0</sup>

0

BBBBBBBBB@

20861 <sup>524288</sup> <sup>e</sup>

y ¼ x þ f xð Þ (16)

CCCCCCA

CCCCCCA

5019 <sup>131072</sup> <sup>e</sup>

10

4924935 <sup>4096</sup> <sup>e</sup>

> B0 ð Þ� 0 1 <sup>B</sup><sup>0</sup> <sup>π</sup> 2 � � � <sup>1</sup>

0

BBBBBBBBB@

B‴ð Þ0 <sup>B</sup>‴ <sup>π</sup> 2 � �

<sup>B</sup>ð Þ<sup>5</sup> ð Þ<sup>0</sup>

ð Þ� 0 1 <sup>B</sup><sup>0</sup> <sup>π</sup> 2 � � � <sup>1</sup>

> B‴ð Þ0 <sup>B</sup>‴ <sup>π</sup> 2 � �

<sup>B</sup>ð Þ<sup>5</sup> ð Þ<sup>0</sup>

10

<sup>10</sup> (12)

1

http://dx.doi.org/10.5772/intechopen.81586

CCCCCCCCCA

1

CCCCCCCCCA

(13)

71

(14)

(15)

Define ψ as

$$
\psi = \frac{X}{a(1 - \varepsilon^2)K\_0} \tag{5}
$$

Substituting Eq. (5) into Eq. (4) yields:

$$\frac{dB}{K\_0 d\psi} = \left(1 - e^2 \sin^2 B\right)^{3/2} \tag{6}$$

Suppose that the inverse solution of Eq. (6) has the following form:

$$B = \psi + a\_2 \sin 2\psi + a\_4 \sin 4\psi + a\_6 \sin 6\psi + a\_8 \sin 8\psi + a\_{10} \sin 10\psi \tag{7}$$

Eq. (7) has five coefficients to be determined. Once these coefficients are known, the inverse problem can be solved.

We consider that the values of differentiation Eq. (6) at the beginning and end points can be treated as interpolation constraints. It is generally known that

$$B'(0) = K\_0 \tag{8}$$

$$B'\left(\frac{\pi}{2}\right) = K\_0 \left(1 - e^2\right)^{3/2} \tag{9}$$

The further derivation of Eq. (6) with respect to ψ yields B00ð Þ ψ . Unfortunately, B00ð Þ ψ is equal to zero at <sup>ψ</sup> <sup>¼</sup> 0, <sup>ψ</sup> <sup>¼</sup> <sup>π</sup> <sup>2</sup>. Hence, one differentiates Eq. (6) twice and it yields B‴ð Þ ψ . Omitting the derivative procedure by means of Mathematica, one arrives at

$$B''(0) = -3\,e^2 - \frac{27}{4}e^4 - \frac{729}{64}e^6 - \frac{4329}{256}e^8 - \frac{381645}{16384}e^{10} \tag{10}$$

$$B''\left(\frac{\pi}{2}\right) = 3\left.e^2 - \frac{15}{4}e^4 + \frac{57}{64}e^6 + \frac{3}{256}e^8 - \frac{51}{16384}e^{10}\right.\tag{11}$$

The further derivation of <sup>B</sup>‴ð Þ <sup>ψ</sup> with respect to <sup>ψ</sup> yields <sup>B</sup>ð Þ<sup>4</sup> ð Þ <sup>ψ</sup> , but <sup>B</sup>ð Þ<sup>4</sup> ð Þ <sup>ψ</sup> is equal to zero at <sup>ψ</sup> <sup>¼</sup> 0, <sup>ψ</sup> <sup>¼</sup> <sup>π</sup> <sup>2</sup>. Hence, one differentiates <sup>B</sup>‴ð Þ <sup>ψ</sup> twice and it yieldsBð Þ<sup>5</sup> ð Þ <sup>ψ</sup> . Then one arrives at

Mathematical Analysis of Some Typical Problems in Geodesy by Means of Computer Algebra http://dx.doi.org/10.5772/intechopen.81586 71

$$B^{(5)}(0) = 12e^2 + 90e^4 + \frac{4455}{16}e^6 + \frac{20145}{32}e^8 + \frac{4924935}{4096}e^{10} \tag{12}$$

Making use of the five interpolation constraints in Eqs. (8–12) and differentiating Eq. (7) correspondingly, one arrives at a set of linear equations for the unknown coefficients

$$
\begin{pmatrix} 2 & 4 & 6 & 8 & 10 \\ -2 & 4 & -6 & 8 & -10 \\ -8 & -64 & -216 & -512 & -1000 \\ 8 & -64 & -216 & -512 & 1000 \\ 32 & 1024 & 7776 & 32768 & 100000 \end{pmatrix} \begin{pmatrix} a\_2 \\ a\_4 \\ a\_6 \\ a\_8 \\ a\_{10} \end{pmatrix} = \begin{pmatrix} B'(0) - 1 \\ B'\left(\frac{\pi}{2}\right) - 1 \\ B'\left(0\right) \\ B'\left(\frac{\pi}{2}\right) \\ B\left(\frac{\pi}{2}\right) \\ B^{(5)}(0) \end{pmatrix} \tag{13}
$$

The solution to Eq. (13) is

Eqs. (2) and (3) can also be derived using the binomial theorem by hand which consumes much more time, see [6–9]. The denominator 693 of the last coefficient K<sup>10</sup> is mistaken as 639 in

2.2. The inverse expansion of the meridian arc using the Hermite interpolation method

dB <sup>¼</sup> <sup>a</sup> <sup>1</sup> � <sup>e</sup><sup>2</sup>

<sup>ψ</sup> <sup>¼</sup> <sup>X</sup>

a 1 � e<sup>2</sup> ð ÞK<sup>0</sup>

<sup>2</sup> sin <sup>2</sup>

Eq. (7) has five coefficients to be determined. Once these coefficients are known, the inverse

We consider that the values of differentiation Eq. (6) at the beginning and end points can be

¼ K<sup>0</sup> 1 � e

The further derivation of Eq. (6) with respect to ψ yields B00ð Þ ψ . Unfortunately, B00ð Þ ψ is equal to

<sup>4</sup> � <sup>729</sup> <sup>64</sup> <sup>e</sup>

The further derivation of <sup>B</sup>‴ð Þ <sup>ψ</sup> with respect to <sup>ψ</sup> yields <sup>B</sup>ð Þ<sup>4</sup> ð Þ <sup>ψ</sup> , but <sup>B</sup>ð Þ<sup>4</sup> ð Þ <sup>ψ</sup> is equal to zero at

<sup>2</sup>. Hence, one differentiates <sup>B</sup>‴ð Þ <sup>ψ</sup> twice and it yieldsBð Þ<sup>5</sup> ð Þ <sup>ψ</sup> . Then one arrives at

B0

<sup>B</sup><sup>0</sup> <sup>π</sup> 2 

<sup>2</sup> � <sup>27</sup> 4 e

> <sup>2</sup> � <sup>15</sup> 4 e 4 þ 57 64 e 6 þ 3 <sup>256</sup> <sup>e</sup>

B ¼ ψ þ a<sup>2</sup> sin 2ψ þ a<sup>4</sup> sin 4ψ þ a<sup>6</sup> sin 6ψ þ a<sup>8</sup> sin 8ψ þ a<sup>10</sup> sin 10ψ (7)

<sup>2</sup>. Hence, one differentiates Eq. (6) twice and it yields B‴ð Þ ψ . Omitting the

<sup>8</sup> � <sup>381645</sup> <sup>16384</sup> <sup>e</sup>

<sup>8</sup> � <sup>51</sup> <sup>16384</sup> <sup>e</sup>

<sup>6</sup> � <sup>4329</sup> <sup>256</sup> <sup>e</sup>

<sup>1</sup> � <sup>e</sup><sup>2</sup> sin ð Þ <sup>2</sup><sup>B</sup> <sup>3</sup>=<sup>2</sup> (4)

B <sup>3</sup>=<sup>2</sup> (6)

ð Þ¼ 0 K<sup>0</sup> (8)

<sup>2</sup> <sup>3</sup>=<sup>2</sup> (9)

<sup>10</sup> (10)

<sup>10</sup> (11)

(5)

dX

dB

Suppose that the inverse solution of Eq. (6) has the following form:

treated as interpolation constraints. It is generally known that

derivative procedure by means of Mathematica, one arrives at

¼ 3 e

B‴ð Þ¼� 0 3 e

<sup>B</sup>‴ <sup>π</sup> 2  <sup>K</sup>0d<sup>ψ</sup> <sup>¼</sup> <sup>1</sup> � <sup>e</sup>

Differentiation to the both sides of Eq. (1) yields:

70 Trends in Geomatics - An Earth Science Perspective

Substituting Eq. (5) into Eq. (4) yields:

problem can be solved.

zero at <sup>ψ</sup> <sup>¼</sup> 0, <sup>ψ</sup> <sup>¼</sup> <sup>π</sup>

<sup>ψ</sup> <sup>¼</sup> 0, <sup>ψ</sup> <sup>¼</sup> <sup>π</sup>

Ref. [9].

Define ψ as

$$
\begin{pmatrix} a\_2 \\ a\_4 \\ a\_6 \\ a\_8 \\ a\_{10} \end{pmatrix} = \begin{pmatrix} 2 & 4 & 6 & 8 & 10 \\ -2 & 4 & -6 & 8 & -10 \\ -8 & -64 & -216 & -512 & -1000 \\ 8 & -64 & -216 & -512 & 1000 \\ 32 & 1024 & 7776 & 32768 & 100000 \end{pmatrix}^{-1} \begin{pmatrix} B'(0) - 1 \\ B'\left(\frac{\pi}{2}\right) - 1 \\ B'\left(\frac{\pi}{2}\right) \\ B'\left(\frac{\pi}{2}\right) \\ B'\left(\frac{\pi}{2}\right) \\ B^{(5)}(0) \end{pmatrix} \tag{14}
$$

Omitting the main operations by means of Mathematica, one arrives at

$$\begin{cases} a\_2 = \frac{3}{8}e^2 + \frac{3}{16}e^4 + \frac{213}{2048}e^6 + \frac{255}{4096}e^8 + \frac{20861}{524288}e^{10} \\ a\_4 = \frac{21}{256}e^4 + \frac{21}{256}e^6 + \frac{533}{8192}e^8 + \frac{197}{4096}e^{10} \\ a\_6 = \frac{151}{6144}e^6 + \frac{151}{4096}e^8 + \frac{5019}{131072}e^{10} \\ a\_8 = \frac{1097}{131072}e^8 + \frac{1097}{65536}e^{10} \\ a\_{10} = \frac{8011}{2621440}e^{10} \end{cases} \tag{15}$$

#### 2.3. The inverse expansion of the meridian arc using Lagrange's theorem method

Suppose that

$$y = \mathbf{x} + f(\mathbf{x})\tag{16}$$

with j j f xð Þ << j j x and y ≈ x. Lagrange's theorem states that in a suitable domain the solution of Eq. (16) is

$$\infty = y + \sum\_{n=1}^{\infty} \frac{(-1)^n}{n!} \frac{d^{n-1}}{dy^{n-1}} [f(y)]^n \tag{17}$$

a<sup>2</sup> ¼ �α<sup>2</sup> � α2α<sup>4</sup> � α4α<sup>6</sup> þ

<sup>a</sup><sup>6</sup> ¼ �α<sup>6</sup> <sup>þ</sup> <sup>3</sup>α2α<sup>4</sup> � <sup>3</sup>α2α<sup>8</sup> � <sup>3</sup>

<sup>a</sup><sup>10</sup> ¼ �α<sup>10</sup> <sup>þ</sup> <sup>5</sup>α2α<sup>8</sup> <sup>þ</sup> <sup>5</sup>α4α<sup>6</sup> � <sup>25</sup>

<sup>a</sup><sup>2</sup> <sup>¼</sup> <sup>3</sup> 8 e 2 þ 3 16 e 4 þ 213 <sup>2048</sup> <sup>e</sup> 6 þ 255 <sup>4096</sup> <sup>e</sup> 8 þ

8

>>>>>>>>>>>>>>>>><

>>>>>>>>>>>>>>>>>:

correctness of the derived formula.

The differences ΔBi ¼ Bi � B0, ΔB<sup>0</sup>

Table 1.

<sup>a</sup><sup>4</sup> <sup>¼</sup> <sup>21</sup> <sup>256</sup> <sup>e</sup> 4 þ 21 <sup>256</sup> <sup>e</sup> 6 þ 533 <sup>8192</sup> <sup>e</sup> 8 þ 197 <sup>4096</sup> <sup>e</sup> 10

<sup>a</sup><sup>6</sup> <sup>¼</sup> <sup>151</sup> <sup>6144</sup> <sup>e</sup> 6 þ 151 <sup>4096</sup> <sup>e</sup> 8 þ

<sup>a</sup><sup>8</sup> <sup>¼</sup> <sup>1097</sup> <sup>131072</sup> <sup>e</sup>

<sup>a</sup><sup>10</sup> <sup>¼</sup> <sup>8011</sup> <sup>2621440</sup> <sup>e</sup>

2.4. The accuracy of the inverse expansion of the meridian arc

derived by Yang (see [6, 7]) is also examined for comparison.

<sup>i</sup> ¼ B<sup>0</sup>

8 þ

10

1097 <sup>65536</sup> <sup>e</sup> 10

The results in Eq. (24) are consistent with the coefficients in Eq. (15), which substantiates the

In order to validate the exactness of the inverse expansions of meridian arc derived by the author, one has examined its accuracy choosing the CGCS2000 reference ellipsoid with parametersa ¼ 6378137, e ¼ 0:08181919104281579. The accuracy of the inverse expansions

One makes use of Eq. (1) and Eq. (5) to obtain the theoretical value ψ0at given geodetic latitude B0. Then one makes use of the inverse expansions derived by Yang (see [6, 7]) to obtain the computation value B1. Substituting ψ0into Eq. (22), one arrives at the computation value B<sup>0</sup>

expansions derived by Yang (see [6, 7]) and the author respectively. These errors are listed in

means of Mathematica. Omitting the main procedure, one arrives at

<sup>a</sup><sup>4</sup> ¼ �α<sup>4</sup> <sup>þ</sup> <sup>α</sup><sup>2</sup>

8

>>>>>>>>>>>>>>>>>><

>>>>>>>>>>>>>>>>>>:

<sup>a</sup><sup>8</sup> ¼ �α<sup>8</sup> <sup>þ</sup> <sup>2</sup>α<sup>2</sup>

1 2 α3 <sup>2</sup> <sup>þ</sup> <sup>α</sup>2α<sup>2</sup>

> <sup>2</sup>α<sup>4</sup> � <sup>4</sup> 3 α4 2

> > <sup>2</sup>α<sup>4</sup> þ 8 3 α4 2

> > > <sup>2</sup> <sup>α</sup>2α<sup>2</sup>

The coefficients in Eq. (23) are also easily expressed in a power series of the eccentricity by

<sup>4</sup> � <sup>25</sup> 2 α2 <sup>2</sup>α<sup>6</sup> þ

5019 <sup>131072</sup> <sup>e</sup>

10

2 α3 <sup>2</sup> þ 9 2 α2α<sup>2</sup>

<sup>2</sup> � <sup>2</sup>α2α<sup>6</sup> <sup>þ</sup> <sup>4</sup>α<sup>2</sup>

<sup>4</sup> <sup>þ</sup> <sup>4</sup>α2α<sup>6</sup> � <sup>8</sup>α<sup>2</sup>

<sup>4</sup> � <sup>1</sup> 2 α2 <sup>2</sup>α<sup>6</sup> þ 1 3 α3 <sup>2</sup>α<sup>4</sup> � <sup>1</sup> <sup>12</sup> <sup>α</sup><sup>5</sup> 2

<sup>4</sup> <sup>þ</sup> <sup>9</sup>α<sup>2</sup>

Mathematical Analysis of Some Typical Problems in Geodesy by Means of Computer Algebra

<sup>2</sup>α<sup>6</sup> � <sup>27</sup> 2 α3 <sup>2</sup>α<sup>4</sup> þ

> 125 <sup>6</sup> <sup>α</sup><sup>3</sup>

20861 <sup>524288</sup> <sup>e</sup>

<sup>i</sup> � B0ð Þ i ¼ 1; 2 indicate the accuracies of the inverse

10

http://dx.doi.org/10.5772/intechopen.81586

(23)

73

(24)

1.

<sup>2</sup>α<sup>4</sup> � <sup>125</sup> <sup>24</sup> <sup>α</sup><sup>5</sup> 2

Supposef xð Þis defined by the following finite trigonometric series:

$$f(\mathbf{x}) = \alpha \sin 2\mathbf{x} + \beta \sin 4\mathbf{x} + \gamma \sin 6\mathbf{x} + \delta \sin 8\mathbf{x} + \varepsilon \sin 10\mathbf{x} \tag{18}$$

where the coefficients <sup>α</sup> <sup>¼</sup> O e<sup>2</sup> � �, <sup>β</sup> <sup>¼</sup> O e<sup>4</sup> � �, <sup>γ</sup> <sup>¼</sup> O e<sup>6</sup> � �, <sup>δ</sup> <sup>¼</sup> O e<sup>8</sup> � �, <sup>ε</sup> <sup>¼</sup> O e<sup>10</sup> � � are small enough for the condition j j f xð Þ << j j x . In deriving the inversion we shall truncate the infinite Lagrange expansion at eighth-order terms of e and drop higher powers. Inserting Eq. (18) into Eq. (17), one arrives at

$$\ln x = y - f(y) + \frac{1}{2!} \frac{d}{dy} [f(y)]^2 - \frac{1}{3!} \frac{d^2}{dy^2} [f(y)]^3 + \frac{1}{4!} \frac{d^3}{dy^3} [f(y)]^4 - \frac{1}{5!} \frac{d^4}{dy^4} [f(y)]^5 \tag{19}$$

One should make use of several trigonometric identities to calculate the derivatives, substituting them into Eq. (19) and grouping terms according to the trigonometric functions. It is a difficult and time-consuming work to do by hand, but could be easily realized by means of Mathematica. Omitting the main procedure, one arrives at

$$\mathbf{x} = y + d\_2 \sin 2y + d\_4 \sin 4y + d\_6 \sin 6y + d\_8 \sin 8y + d\_{10} \sin 10y \tag{20}$$

where

$$\begin{cases} d\_2 = -\alpha - a\beta - \beta\gamma + \frac{1}{2}\alpha^3 + a\beta^2 - \frac{1}{2}\alpha^2\gamma + \frac{1}{3}\alpha^3\beta - \frac{1}{12}\alpha^5 \\\\ d\_4 = -\beta + \alpha^2 - 2\alpha\gamma + 4\alpha^2\beta - \frac{4}{3}\alpha^4 \\\\ d\_6 = -\gamma + 3\alpha\beta - 3\alpha\delta - \frac{3}{2}\alpha^3 + \frac{9}{2}\alpha\beta^2 + 9\alpha^2\gamma - \frac{27}{2}\alpha^3\beta + \frac{27}{8}\alpha^5 \\\\ d\_8 = -\delta + 2\beta^2 + 4\alpha\gamma - 8\alpha^2\beta + \frac{8}{3}\alpha^4 \\\\ d\_{10} = -\epsilon + 5\alpha\delta + 5\beta\gamma - \frac{25}{2}\alpha\beta^2 - \frac{25}{2}\alpha^2\gamma + \frac{125}{6}\alpha^3\beta - \frac{125}{24}\alpha^5 \end{cases} \tag{21}$$

Substituting x for Band y for ψ in Eq. (16), the coefficients α, β, γ, δ, ε are consistent with α2, α4, α6, α<sup>8</sup> in Eq. (3). According to Eq. (20) and denoting a2, a4, a6, a8, a<sup>10</sup> as the new coefficients, the inverse expansion of the meridian arc can be written as

$$B = \psi + a\_2 \sin 2\psi + a\_4 \sin 4\psi + a\_6 \sin 6\psi + a\_8 \sin 8\psi + a\_{10} \sin 10\psi \tag{22}$$

where

Mathematical Analysis of Some Typical Problems in Geodesy by Means of Computer Algebra http://dx.doi.org/10.5772/intechopen.81586 73

$$\begin{cases} a\_2 = -a\_2 - a\_2 a\_4 - \alpha\_4 a\_6 + \frac{1}{2} a\_2^3 + \alpha\_2 a\_4^2 - \frac{1}{2} a\_2^2 a\_6 + \frac{1}{3} a\_2^3 a\_4 - \frac{1}{12} a\_2^5 \\\\ a\_4 = -a\_4 + a\_2^2 - 2a\_2 a\_6 + 4a\_2^2 a\_4 - \frac{4}{3} a\_2^4 \\\\ a\_6 = -a\_6 + 3a\_2 a\_4 - 3a\_2 a\_8 - \frac{3}{2} a\_2^3 + \frac{9}{2} a\_2 a\_4^2 + 9a\_2^2 a\_6 - \frac{27}{2} a\_2^3 a\_4 + \frac{27}{8} a\_2^5 \\\\ a\_8 = -a\_8 + 2a\_4^2 + 4a\_2 a\_6 - 8a\_2^2 a\_4 + \frac{8}{3} a\_2^4 \\\\ a\_{10} = -a\_{10} + 5a\_2 a\_8 + 5a\_4 a\_6 - \frac{25}{2} a\_2 a\_4^2 - \frac{25}{2} a\_2^2 a\_6 + \frac{125}{6} a\_2^3 a\_4 - \frac{125}{24} a\_2^5 \end{cases} \tag{23}$$

The coefficients in Eq. (23) are also easily expressed in a power series of the eccentricity by means of Mathematica. Omitting the main procedure, one arrives at

$$\begin{cases} a\_2 = \frac{3}{8}e^2 + \frac{3}{16}e^4 + \frac{213}{2048}e^6 + \frac{255}{4096}e^8 + \frac{20861}{524288}e^{10} \\\\ a\_4 = \frac{21}{256}e^4 + \frac{21}{256}e^6 + \frac{533}{8192}e^8 + \frac{197}{4096}e^{10} \\\\ a\_6 = \frac{151}{6144}e^6 + \frac{151}{4096}e^8 + \frac{5019}{131072}e^{10} \\\\ a\_8 = \frac{1097}{131072}e^8 + \frac{1097}{65536}e^{10} \\\\ a\_{10} = \frac{8011}{2621440}e^{10} \end{cases} \tag{24}$$

The results in Eq. (24) are consistent with the coefficients in Eq. (15), which substantiates the correctness of the derived formula.

#### 2.4. The accuracy of the inverse expansion of the meridian arc

<sup>x</sup> <sup>¼</sup> <sup>y</sup> <sup>þ</sup>X<sup>∞</sup>

Supposef xð Þis defined by the following finite trigonometric series:

Eq. (17), one arrives at

where

where

<sup>x</sup> <sup>¼</sup> <sup>y</sup> � f yð Þþ <sup>1</sup>

72 Trends in Geomatics - An Earth Science Perspective

8

>>>>>>>>>>>>>>>>>><

>>>>>>>>>>>>>>>>>>:

2! d

Mathematica. Omitting the main procedure, one arrives at

d<sup>2</sup> ¼ �α � αβ � βγ þ

<sup>d</sup><sup>4</sup> ¼ �<sup>β</sup> <sup>þ</sup> <sup>α</sup><sup>2</sup> � <sup>2</sup>αγ <sup>þ</sup> <sup>4</sup>α<sup>2</sup><sup>β</sup> � <sup>4</sup>

<sup>d</sup><sup>8</sup> ¼ �<sup>δ</sup> <sup>þ</sup> <sup>2</sup>β<sup>2</sup> <sup>þ</sup> <sup>4</sup>αγ � <sup>8</sup>α<sup>2</sup><sup>β</sup> <sup>þ</sup>

<sup>d</sup><sup>10</sup> ¼ �<sup>ε</sup> <sup>þ</sup> <sup>5</sup>αδ <sup>þ</sup> <sup>5</sup>βγ � <sup>25</sup>

inverse expansion of the meridian arc can be written as

<sup>d</sup><sup>6</sup> ¼ �<sup>γ</sup> <sup>þ</sup> <sup>3</sup>αβ � <sup>3</sup>αδ � <sup>3</sup>

dy ½ � f yð Þ <sup>2</sup> � <sup>1</sup>

3! d2

1 2

> 2 <sup>α</sup><sup>3</sup> <sup>þ</sup> 9 2

n¼1

ð Þ �<sup>1</sup> <sup>n</sup> n!

where the coefficients <sup>α</sup> <sup>¼</sup> O e<sup>2</sup> � �, <sup>β</sup> <sup>¼</sup> O e<sup>4</sup> � �, <sup>γ</sup> <sup>¼</sup> O e<sup>6</sup> � �, <sup>δ</sup> <sup>¼</sup> O e<sup>8</sup> � �, <sup>ε</sup> <sup>¼</sup> O e<sup>10</sup> � � are small enough for the condition j j f xð Þ << j j x . In deriving the inversion we shall truncate the infinite Lagrange expansion at eighth-order terms of e and drop higher powers. Inserting Eq. (18) into

dy<sup>2</sup> ½ � f yð Þ <sup>3</sup> <sup>þ</sup>

One should make use of several trigonometric identities to calculate the derivatives, substituting them into Eq. (19) and grouping terms according to the trigonometric functions. It is a difficult and time-consuming work to do by hand, but could be easily realized by means of

<sup>α</sup><sup>3</sup> <sup>þ</sup> αβ<sup>2</sup> � <sup>1</sup>

3 α4

> 8 3 α4

<sup>2</sup> αβ<sup>2</sup> � <sup>25</sup>

Substituting x for Band y for ψ in Eq. (16), the coefficients α, β, γ, δ, ε are consistent with α2, α4, α6, α<sup>8</sup> in Eq. (3). According to Eq. (20) and denoting a2, a4, a6, a8, a<sup>10</sup> as the new coefficients, the

d<sup>n</sup>�<sup>1</sup>

f xð Þ¼ α sin 2x þ β sin 4x þ γ sin 6x þ δ sin 8x þ ε sin 10x (18)

1 4! d3

x ¼ y þ d<sup>2</sup> sin 2y þ d<sup>4</sup> sin 4y þ d<sup>6</sup> sin 6y þ d<sup>8</sup> sin 8y þ d<sup>10</sup> sin 10y (20)

2 α2 γ þ 1 3 α3 <sup>β</sup> � <sup>1</sup> <sup>12</sup> <sup>α</sup><sup>5</sup>

αβ<sup>2</sup> <sup>þ</sup> <sup>9</sup>α<sup>2</sup>

2 α2 γ þ 125 <sup>6</sup> <sup>α</sup><sup>3</sup>

B ¼ ψ þ a<sup>2</sup> sin 2ψ þ a<sup>4</sup> sin 4ψ þ a<sup>6</sup> sin 6ψ þ a<sup>8</sup> sin 8ψ þ a<sup>10</sup> sin 10ψ (22)

<sup>γ</sup> � <sup>27</sup> 2 α3 β þ 27 8 α5

> <sup>β</sup> � <sup>125</sup> <sup>24</sup> <sup>α</sup><sup>5</sup>

dy<sup>3</sup> ½ � f yð Þ <sup>4</sup> � <sup>1</sup>

5! d4

dy<sup>4</sup> ½ � f yð Þ <sup>5</sup> (19)

(21)

dyn�<sup>1</sup> ½ � f yð Þ <sup>n</sup> (17)

In order to validate the exactness of the inverse expansions of meridian arc derived by the author, one has examined its accuracy choosing the CGCS2000 reference ellipsoid with parametersa ¼ 6378137, e ¼ 0:08181919104281579. The accuracy of the inverse expansions derived by Yang (see [6, 7]) is also examined for comparison.

One makes use of Eq. (1) and Eq. (5) to obtain the theoretical value ψ0at given geodetic latitude B0. Then one makes use of the inverse expansions derived by Yang (see [6, 7]) to obtain the computation value B1. Substituting ψ0into Eq. (22), one arrives at the computation value B<sup>0</sup> 1. The differences ΔBi ¼ Bi � B0, ΔB<sup>0</sup> <sup>i</sup> ¼ B<sup>0</sup> <sup>i</sup> � B0ð Þ i ¼ 1; 2 indicate the accuracies of the inverse expansions derived by Yang (see [6, 7]) and the author respectively. These errors are listed in Table 1.


Table 1. Errors of the inverse expansions of meridian arc.

From Table 1, one could find that the accuracy of the inverse expansion of meridian arc derived by Yang (see [6, 7]) is higher than 10�500, and the accuracy of the inverse expansion Eq. (22) derived by the author is higher than 10�700. The accuracy is improved by 2 orders of magnitude by means of computer algebra.

### 3. The singular integration in physical geodesy

Singular integrals associated with the reciprocal distance usually exist in the computations of physical geodesy and geophysics. For example, the integral expressions of height anomaly, deflections of the vertical and vertical gradient of gravity anomaly can be written in planar approximation as

$$\mathcal{L} = \frac{1}{2\pi\gamma} \iiint \frac{\Delta\mathbf{g}}{r} d\mathbf{x} dy \tag{25}$$

As shown in Figure 1, let the innermost area be the rectangular σ∈ ½ � �a < x < a; �b < y < b (a > 0, b > 0) due to the convergence of meridian, and the gravity anomaly is expressed as

> x a � �<sup>i</sup> X 2

Inserting the innermost area into Eqs. (25)–(28), and let the contributions be Δζ, Δξ, Δηand ΔL,

Δg xð Þ ; y

xΔg xð Þ ; y

yΔg xð Þ ; y

Δg xð Þ� ; y Δgð Þ 0; 0

<sup>u</sup> <sup>¼</sup> <sup>x</sup> a <sup>v</sup> <sup>¼</sup> <sup>y</sup> b

8 ><

>:

<sup>r</sup> dxdy (30)

<sup>r</sup><sup>3</sup> dxdy (31)

<sup>r</sup><sup>3</sup> dxdy (32)

<sup>r</sup><sup>3</sup> dxdy (33)

j¼0 αij y b � �<sup>j</sup>

Mathematical Analysis of Some Typical Problems in Geodesy by Means of Computer Algebra

http://dx.doi.org/10.5772/intechopen.81586

(29)

75

(34)

2

i¼0

<sup>2</sup>πγ ðð σ

<sup>2</sup>πγ ðð σ

<sup>Δ</sup>g xð Þ¼ ; <sup>y</sup> <sup>X</sup>

<sup>Δ</sup><sup>ζ</sup> <sup>¼</sup> <sup>1</sup> <sup>2</sup>πγ ðð σ

<sup>Δ</sup><sup>ξ</sup> ¼ � <sup>1</sup>

<sup>Δ</sup><sup>η</sup> ¼ � <sup>1</sup>

σ

<sup>Δ</sup><sup>L</sup> <sup>¼</sup> <sup>1</sup> 2π ðð

The following transformation is introduced for σ

double quadratic polynomial:

Figure 1. Integrals in the rectangular area.

one arrives at

$$\xi = -\frac{1}{2\pi\gamma} \iiint \frac{\mathbf{x} \Delta \mathbf{g}}{r^3} d\mathbf{x} dy \tag{26}$$

$$\eta = -\frac{1}{2\pi\eta} \iiint \frac{y\Delta\mathbf{g}}{r^3} d\mathbf{x} dy \tag{27}$$

$$L = \frac{1}{2\pi} \iiint \frac{\Delta \mathbf{g} - \Delta \mathbf{g}\_0}{r^3} d\mathbf{x} dy\tag{28}$$

where <sup>r</sup> <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>y</sup><sup>2</sup> <sup>p</sup> , <sup>Δ</sup>g<sup>0</sup> is the gravity anomaly at the computation point. When <sup>r</sup> ! 0, the above integrals become singular and need special treatment at the computation point. The past treatments are with respect to template computations, which regards the innermost area as a circle, see [3, 16]. But the gravity anomalies are given in the rectangular grids(such as 2<sup>0</sup> � 2<sup>0</sup> ), if the approximate disposal is used, some significant error may be introduced. Sunkel and Wang expressed the gravity anomalies block by block as an interpolation polynomial and derived the analytic values of the integrals, see [17, 18]. However, the integrals of the rational functions are very complicated, especially when related interpolation polynomials contain many terms. One only can give the analytic values of the corresponding linear approximation. In this chapter, we use the nonsingular integration transformations proposed by Bian (see [19]) to compute the above integrals precisely.

Mathematical Analysis of Some Typical Problems in Geodesy by Means of Computer Algebra http://dx.doi.org/10.5772/intechopen.81586 75

Figure 1. Integrals in the rectangular area.

From Table 1, one could find that the accuracy of the inverse expansion of meridian arc derived by Yang (see [6, 7]) is higher than 10�500, and the accuracy of the inverse expansion Eq. (22) derived by the author is higher than 10�700. The accuracy is improved by 2 orders of

<sup>Δ</sup>B1<sup>=</sup> <sup>00</sup> ð Þ �3:<sup>2</sup> � <sup>10</sup>�<sup>7</sup> �1:<sup>6</sup> � <sup>10</sup>�<sup>6</sup> �1:<sup>8</sup> � <sup>10</sup>�<sup>6</sup> �1:<sup>4</sup> � <sup>10</sup>�<sup>6</sup>

<sup>1</sup><sup>=</sup> <sup>00</sup> ð Þ <sup>2</sup>:<sup>7</sup> � <sup>10</sup>�<sup>9</sup> <sup>8</sup>:<sup>6</sup> � <sup>10</sup>�<sup>9</sup> <sup>1</sup>:<sup>3</sup> � <sup>10</sup>�<sup>8</sup> <sup>1</sup>:<sup>7</sup> � <sup>10</sup>�<sup>8</sup>

<sup>B</sup>0<sup>=</sup> <sup>∘</sup> ð Þ <sup>20</sup> <sup>40</sup> <sup>60</sup> <sup>80</sup>

Singular integrals associated with the reciprocal distance usually exist in the computations of physical geodesy and geophysics. For example, the integral expressions of height anomaly, deflections of the vertical and vertical gradient of gravity anomaly can be written in planar

ðð Δg

ðð xΔg

ðð yΔg

ðð <sup>Δ</sup><sup>g</sup> � <sup>Δ</sup>g<sup>0</sup>

above integrals become singular and need special treatment at the computation point. The past treatments are with respect to template computations, which regards the innermost area as a circle, see [3, 16]. But the gravity anomalies are given in the rectangular grids(such as 2<sup>0</sup> � 2<sup>0</sup>

the approximate disposal is used, some significant error may be introduced. Sunkel and Wang expressed the gravity anomalies block by block as an interpolation polynomial and derived the analytic values of the integrals, see [17, 18]. However, the integrals of the rational functions are very complicated, especially when related interpolation polynomials contain many terms. One only can give the analytic values of the corresponding linear approximation. In this chapter, we use the nonsingular integration transformations proposed by Bian (see [19]) to compute the

<sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>y</sup><sup>2</sup> <sup>p</sup> , <sup>Δ</sup>g<sup>0</sup> is the gravity anomaly at the computation point. When <sup>r</sup> ! 0, the

<sup>r</sup> dxdy (25)

<sup>r</sup><sup>3</sup> dxdy (26)

<sup>r</sup><sup>3</sup> dxdy (27)

<sup>r</sup><sup>3</sup> dxdy (28)

), if

<sup>ζ</sup> <sup>¼</sup> <sup>1</sup> 2πγ

<sup>ξ</sup> ¼ � <sup>1</sup> 2πγ

<sup>η</sup> ¼ � <sup>1</sup> 2πγ

<sup>L</sup> <sup>¼</sup> <sup>1</sup> 2π

magnitude by means of computer algebra.

Table 1. Errors of the inverse expansions of meridian arc.

74 Trends in Geomatics - An Earth Science Perspective

approximation as

ΔB<sup>0</sup>

where <sup>r</sup> <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

above integrals precisely.

3. The singular integration in physical geodesy

As shown in Figure 1, let the innermost area be the rectangular σ∈ ½ � �a < x < a; �b < y < b (a > 0, b > 0) due to the convergence of meridian, and the gravity anomaly is expressed as double quadratic polynomial:

$$\Delta \mathbf{g}(\mathbf{x}, \mathbf{y}) = \sum\_{i=0}^{2} \left(\frac{\mathbf{x}}{a}\right)^{i} \sum\_{j=0}^{2} a\_{ij} \left(\frac{\mathbf{y}}{b}\right)^{j} \tag{29}$$

Inserting the innermost area into Eqs. (25)–(28), and let the contributions be Δζ, Δξ, Δηand ΔL, one arrives at

$$
\Delta \zeta = \frac{1}{2\pi \gamma} \iiint\_{\sigma} \frac{\Delta \mathbf{g}(\mathbf{x}, y)}{r} d\mathbf{x} dy \tag{30}
$$

$$
\Delta\xi = -\frac{1}{2\pi\gamma} \iiint\_{\mathcal{\sigma}} \frac{\mathbf{x} \Delta\mathbf{g}(\mathbf{x}, y)}{r^3} d\mathbf{x} dy \tag{31}
$$

$$
\Delta \eta = -\frac{1}{2\pi \gamma} \iiint\_{\sigma} \frac{y \Delta g(\mathbf{x}, y)}{r^3} d\mathbf{x} dy \tag{32}
$$

$$
\Delta L = \frac{1}{2\pi} \iint\limits\_{\boldsymbol{\phi}} \frac{\Delta \mathbf{g}(\mathbf{x}, \boldsymbol{y}) - \Delta \mathbf{g}(\mathbf{0}, \mathbf{0})}{r^3} d\mathbf{x} d\mathbf{y} \tag{33}
$$

The following transformation is introduced for σ

$$\begin{cases} u = \frac{\chi}{a} \\ v = \frac{\mathcal{Y}}{b} \end{cases} \tag{34}$$

Using the properties of the integration for even/odd functions and exploiting the symmetry of the integration area, one arrives at

$$
\Delta\zeta = \frac{4ab}{2\pi\gamma'} \int\_0^1 \int\_0^1 \frac{a\_{00} + a\_{20}u^2 + a\_{02}v^2 + a\_{22}u^2v^2}{\sqrt{a^2u^2 + b^2v^2}} du dv \tag{35}
$$

$$
\Delta\xi = \frac{-4a^2b}{2\pi\gamma} \int\_0^1 \int\_0^1 \frac{u^2 \left(\alpha\_{10} + \alpha\_{12}v^2\right)}{\left(a^2u^2 + b^2v^2\right)^{3/2}} du dv \tag{36}
$$

<sup>Δ</sup><sup>L</sup> <sup>¼</sup> <sup>4</sup>ab 2π ð1 0

> þ 4ab 2π ð1 0

<sup>α</sup><sup>20</sup> <sup>þ</sup> <sup>α</sup>02k<sup>2</sup> <sup>þ</sup>

<sup>α</sup><sup>00</sup> <sup>þ</sup> <sup>α</sup>20λ<sup>2</sup> <sup>þ</sup>

1 3 α22k 2

1 3 α22λ<sup>2</sup> � � dλ

Now we can see that the denominators are greater than zero after transformation Eq. (39) and Eq. (40), and the singularities have been eliminated. The integrals in x and y directions are converted to the integrals of the powers of k andλ. This basically changes the double integrals to single variable integrals, which could easily be calculated in Mathematica as follows:

<sup>a</sup><sup>2</sup> <sup>þ</sup> <sup>b</sup><sup>2</sup> k <sup>2</sup> � �<sup>3</sup>=<sup>2</sup>

Mathematical Analysis of Some Typical Problems in Geodesy by Means of Computer Algebra

<sup>b</sup><sup>2</sup> <sup>þ</sup> <sup>a</sup><sup>2</sup>λ<sup>2</sup> � �<sup>3</sup>=<sup>2</sup>

http://dx.doi.org/10.5772/intechopen.81586

(44)

77

ð45Þ

ð46Þ

� � dk

$$
\Delta \eta = \frac{-4ab^2}{2\pi \gamma'} \int\_0^1 \int\_0^1 \frac{v^2 \left(\alpha\_{01} + \alpha\_{21} v^2\right)}{\left(a^2 u^2 + b^2 v^2\right)^{3/2}} du dv \tag{37}
$$

$$
\Delta L = \frac{4ab}{2\pi} \int\_0^1 \int\_0^1 \frac{\alpha\_{20}u^2 + \alpha\_{02}v^2 + \alpha\_{22}u^2v^2}{\left(a^2u^2 + b^2v^2\right)^{3/2}} du dv \tag{38}
$$

Drawing a line from the origin to the upper right corner, it divides the upper right quadrant. into σ<sup>1</sup> ∈½ � 0 < u < 1; 0 < v < u , σ<sup>2</sup> ∈½ � 0 < u < v; 0 < u < 1 .

The following nonsingular integration transformation is introduced for σ<sup>1</sup>

$$\begin{cases} u = v \\ k = \frac{v}{u} \end{cases} \tag{39}$$

The following nonsingular integration transformation is introduced for σ<sup>2</sup>

$$\begin{cases} v = v \\ \lambda = \frac{u}{v} \end{cases} \tag{40}$$

Insertingv ¼ ku(oru ¼ λv)into Eqs. (35)–(37),one arrives at

<sup>Δ</sup><sup>ζ</sup> <sup>¼</sup> <sup>4</sup>ab <sup>2</sup>πγ <sup>ð</sup><sup>1</sup> 0 α<sup>00</sup> þ 1 3 α<sup>20</sup> þ 1 3 <sup>α</sup>02k<sup>2</sup> <sup>þ</sup> 1 5 α22k<sup>2</sup> � � dk ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>a</sup><sup>2</sup> <sup>þ</sup> <sup>b</sup><sup>2</sup> k <sup>2</sup> p þ 4ab <sup>2</sup>πγ <sup>ð</sup><sup>1</sup> 0 α<sup>00</sup> þ 1 3 α<sup>02</sup> þ 1 3 <sup>α</sup>20λ<sup>2</sup> <sup>þ</sup> 1 5 α22λ<sup>2</sup> � � dλ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>b</sup><sup>2</sup> <sup>þ</sup> <sup>a</sup><sup>2</sup>λ<sup>2</sup> <sup>p</sup> (41) <sup>Δ</sup><sup>ξ</sup> <sup>¼</sup> �4a<sup>2</sup><sup>b</sup> <sup>2</sup>πγ <sup>ð</sup><sup>1</sup> 0 α<sup>10</sup> þ 1 3 α12k<sup>2</sup> � � dk <sup>a</sup><sup>2</sup> <sup>þ</sup> <sup>b</sup><sup>2</sup> k<sup>2</sup> � �<sup>3</sup>=<sup>2</sup> � <sup>4</sup>a<sup>2</sup><sup>b</sup> <sup>2</sup>πγ <sup>ð</sup><sup>1</sup> 0 α<sup>10</sup> þ 1 3 <sup>α</sup><sup>12</sup> � � <sup>λ</sup><sup>2</sup> dλ <sup>b</sup><sup>2</sup> <sup>þ</sup> <sup>a</sup><sup>2</sup>λ<sup>2</sup> � �<sup>3</sup>=<sup>2</sup> (42) <sup>Δ</sup><sup>η</sup> <sup>¼</sup> �4ab<sup>2</sup> <sup>2</sup>πγ <sup>ð</sup><sup>1</sup> 0 α<sup>01</sup> þ 1 3 <sup>α</sup><sup>21</sup> � � <sup>k</sup> 2 dk <sup>a</sup><sup>2</sup> <sup>þ</sup> <sup>b</sup><sup>2</sup> k <sup>2</sup> � �<sup>3</sup>=<sup>2</sup> � <sup>4</sup>ab<sup>2</sup> <sup>2</sup>πγ <sup>ð</sup><sup>1</sup> 0 α<sup>01</sup> þ 1 3 α21λ<sup>2</sup> � � dλ <sup>b</sup><sup>2</sup> <sup>þ</sup> <sup>a</sup><sup>2</sup>λ<sup>2</sup> � �<sup>3</sup>=<sup>2</sup> (43)

Mathematical Analysis of Some Typical Problems in Geodesy by Means of Computer Algebra http://dx.doi.org/10.5772/intechopen.81586 77

$$\begin{split} \Delta L &= \frac{4ab}{2\pi} \Big|\_{0}^{1} \left( \alpha\_{20} + \alpha\_{02} k^{2} + \frac{1}{3} \alpha\_{22} k^{2} \right) \frac{dk}{\left( a^{2} + b^{2} k^{2} \right)^{3/2}} \\ &+ \frac{4ab}{2\pi} \Big|\_{0}^{1} \left( \alpha\_{00} + \alpha\_{20} \lambda^{2} + \frac{1}{3} \alpha\_{22} \lambda^{2} \right) \frac{d\lambda}{\left( b^{2} + a^{2} \lambda^{2} \right)^{3/2}} \end{split} \tag{44}$$

Now we can see that the denominators are greater than zero after transformation Eq. (39) and Eq. (40), and the singularities have been eliminated. The integrals in x and y directions are converted to the integrals of the powers of k andλ. This basically changes the double integrals to single variable integrals, which could easily be calculated in Mathematica as follows:

Using the properties of the integration for even/odd functions and exploiting the symmetry of

<sup>α</sup><sup>00</sup> <sup>þ</sup> <sup>α</sup>20u<sup>2</sup> <sup>þ</sup> <sup>α</sup>02v<sup>2</sup> <sup>þ</sup> <sup>α</sup>22u<sup>2</sup>v<sup>2</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>a</sup><sup>2</sup>u<sup>2</sup> <sup>þ</sup> <sup>b</sup><sup>2</sup>

<sup>u</sup><sup>2</sup> <sup>α</sup><sup>10</sup> <sup>þ</sup> <sup>α</sup>12v<sup>2</sup> � �

<sup>v</sup><sup>2</sup> <sup>α</sup><sup>01</sup> <sup>þ</sup> <sup>α</sup>21v<sup>2</sup> � �

<sup>a</sup><sup>2</sup>u<sup>2</sup> <sup>þ</sup> <sup>b</sup><sup>2</sup>

<sup>a</sup><sup>2</sup>u<sup>2</sup> <sup>þ</sup> <sup>b</sup><sup>2</sup>

<sup>α</sup>20u<sup>2</sup> <sup>þ</sup> <sup>α</sup>02v<sup>2</sup> <sup>þ</sup> <sup>α</sup>22u<sup>2</sup>v<sup>2</sup> <sup>a</sup><sup>2</sup>u<sup>2</sup> <sup>þ</sup> <sup>b</sup><sup>2</sup>

Drawing a line from the origin to the upper right corner, it divides the upper right quadrant.

u ¼ v <sup>k</sup> <sup>¼</sup> <sup>v</sup> u

v ¼ v <sup>λ</sup> <sup>¼</sup> <sup>u</sup> v

<sup>α</sup>02k<sup>2</sup> <sup>þ</sup>

<sup>α</sup>20λ<sup>2</sup> <sup>þ</sup>

� � dλ

� � dk

1 5 α22k<sup>2</sup>

> 1 5 α22λ<sup>2</sup>

<sup>a</sup><sup>2</sup> <sup>þ</sup> <sup>b</sup><sup>2</sup> k <sup>2</sup> � �<sup>3</sup>=<sup>2</sup>

dλ <sup>b</sup><sup>2</sup> <sup>þ</sup> <sup>a</sup><sup>2</sup>λ<sup>2</sup> � �<sup>3</sup>=<sup>2</sup>

> 2 dk

<sup>b</sup><sup>2</sup> <sup>þ</sup> <sup>a</sup><sup>2</sup>λ<sup>2</sup> � �<sup>3</sup>=<sup>2</sup>

<sup>a</sup><sup>2</sup> <sup>þ</sup> <sup>b</sup><sup>2</sup> k <sup>2</sup> � �<sup>3</sup>=<sup>2</sup>

(

(

v2

<sup>p</sup> dudv (35)

<sup>v</sup><sup>2</sup> � �<sup>3</sup>=<sup>2</sup> dudv (36)

<sup>v</sup><sup>2</sup> � �<sup>3</sup>=<sup>2</sup> dudv (37)

<sup>v</sup><sup>2</sup> � �<sup>3</sup>=<sup>2</sup> dudv (38)

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>a</sup><sup>2</sup> <sup>þ</sup> <sup>b</sup><sup>2</sup> k <sup>2</sup> p

> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>b</sup><sup>2</sup> <sup>þ</sup> <sup>a</sup><sup>2</sup>λ<sup>2</sup> <sup>p</sup>

(39)

(40)

(41)

(42)

(43)

the integration area, one arrives at

76 Trends in Geomatics - An Earth Science Perspective

<sup>Δ</sup><sup>ζ</sup> <sup>¼</sup> <sup>4</sup>ab 2πγ ð1 0 ð1 0

<sup>Δ</sup><sup>ξ</sup> <sup>¼</sup> �4a<sup>2</sup><sup>b</sup> 2πγ

<sup>Δ</sup><sup>η</sup> <sup>¼</sup> �4ab<sup>2</sup> 2πγ

<sup>Δ</sup><sup>L</sup> <sup>¼</sup> <sup>4</sup>ab 2π ð1 0 ð1 0

into σ<sup>1</sup> ∈½ � 0 < u < 1; 0 < v < u , σ<sup>2</sup> ∈½ � 0 < u < v; 0 < u < 1 .

Insertingv ¼ ku(oru ¼ λv)into Eqs. (35)–(37),one arrives at

ð1 0

> ð1 0

<sup>Δ</sup><sup>ξ</sup> <sup>¼</sup> �4a<sup>2</sup><sup>b</sup> 2πγ

> � <sup>4</sup>a<sup>2</sup><sup>b</sup> 2πγ

<sup>Δ</sup><sup>η</sup> <sup>¼</sup> �4ab<sup>2</sup> 2πγ

> � <sup>4</sup>ab<sup>2</sup> 2πγ

α<sup>00</sup> þ 1 3 α<sup>20</sup> þ 1 3

α<sup>00</sup> þ 1 3 α<sup>02</sup> þ 1 3

> ð1 0

ð1 0

> ð1 0

ð1 0

α<sup>10</sup> þ 1 3 α12k<sup>2</sup> � � dk

α<sup>10</sup> þ 1 3 α<sup>12</sup> � � λ<sup>2</sup>

α<sup>01</sup> þ 1 3 α<sup>21</sup> � � k

α<sup>01</sup> þ 1 3 α21λ<sup>2</sup> � � dλ

<sup>Δ</sup><sup>ζ</sup> <sup>¼</sup> <sup>4</sup>ab 2πγ

> þ 4ab 2πγ

The following nonsingular integration transformation is introduced for σ<sup>1</sup>

The following nonsingular integration transformation is introduced for σ<sup>2</sup>

ð1 0 ð1 0

ð1 0 ð1 0

ð45Þ

ð46Þ

$$\begin{aligned} \text{a. } &= \begin{bmatrix} \text{a. } & \text{a. } & \text{a. } \\ \hline 2\pi \text{ } & \text{a. } & \text{a. } \\ \hline 2\pi \text{ } & \text{a. } & \text{a. } \\ \hline 2\pi \text{ } & \text{a. } & \text{a. } \\ \hline 2\pi \text{ } & \text{a. } & \text{a. } \\ \end{bmatrix} \begin{aligned} \text{a. } &= \begin{bmatrix} \text{a. } & \text{a. } & \text{a. } \\ \hline & \text{a. } & \text{a. } \\ \hline & \text{a. } & \text{a. } \\ \hline & \text{a. } & \text{a. } \\ \end{bmatrix} \begin{aligned} \text{a. } &= \begin{bmatrix} \text{a. } & \text{a. } \\ \hline & \text{a. } \\ \hline & \text{a. } \\ \end{bmatrix} \begin{aligned} \text{a. } &= \begin{bmatrix} \text{a. } & \text{a. } \\ \hline & \text{a. } \\ \end{bmatrix} \end{aligned} \\ &= \begin{bmatrix} \text{a. } & \text{a. } \\ \hline & \text{a. } & \text{a. } \\ \end{bmatrix} \begin{aligned} \text{a. } &= \begin{bmatrix} \text{a. } & \text{a. } \\ \hline & \text{a. } \\ \end{bmatrix} \end{aligned} \\ &= \begin{bmatrix} \text{a. } & \text{a. } \\ \hline & \text{a. } & \text{a. } \\ \end{bmatrix} \begin{aligned} \text{a. } & \text{a. } \\ & \text{a. } \\ \end{aligned} \\ &= \begin{bmatrix} \text{a. } & \text{a. } \\ \end{bmatrix} \end{aligned}$$

ð52Þ

79

http://dx.doi.org/10.5772/intechopen.81586

One could find that it is greatly fussy to complete these integrals by hand, which could be

Mathematical Analysis of Some Typical Problems in Geodesy by Means of Computer Algebra

The determination of satellite orbit is one of the fundamental problems in satellite geodesy. A

Eccentric, mean and true anomalies are used to describe the movement of satellites. Their transformations are often to be dealt with in satellite ephemeris computation and orbit determination of the spacecraft. In Figure 2, E is the Eccentric anomaly, υis the true anomaly. In order to realize the direct transformations between these anomalies, the series expansions of their transformations are derived using the power series method with the help of computer algebra system Mathematica. Their coefficients are expressed in a power series of the orbital

4. The series expansions of direct transformations between three

graphical representation of the Keplerian orbit is given in Figure 2, see [20].

eccentricity e and extended up to eighth-order terms of the orbital eccentricity.

Let the mean anomaly be M. Mcan be expressed by E as follows:

Differentiating the both sides of Eq. (53) yields

4.1. The series expansions of the direct transformation between eccentric and mean

M ¼ E � e sin E (53)

easily realized using some commands of computer algebra system.

anomalies in satellite geodesy

anomalies

Figure 2. Keplerian orbit.

ð48Þ

In case of a square grid with a unit length, Eqs. (45)–(48) can be simplified in Mathematica as.

$$\frac{60 \text{ \~Ar\text{\footnotesizeSink}\{1\} \text{ \kern-1.72\text{\footnotesize-h\text{\footnotesize$$

$$\begin{array}{c} \text{-6 \texttt{Arc\\$} \| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,\| \,$$

$$\frac{-6\operatorname{Arc\ $}\sinh\{1\}\operatorname{\mathfrak{a}\_{\mathsf{I},\mathsf{L}}\star\operatorname{2}\left(\sqrt{2}\ -2\operatorname{Arc\$ }\sinh\{1\}\right)\operatorname{\mathfrak{a}\_{\mathsf{I},\mathsf{L}}}}{3\operatorname{\mathfrak{a}\,\mathsf{Y}}}\qquad(51)$$

ð52Þ

One could find that it is greatly fussy to complete these integrals by hand, which could be easily realized using some commands of computer algebra system.
