**2. The forward expansions of the rectifying, conformal and authalic latitudes**

Cartographers prefer to adopt sphere as a basis of the map projection for convenience since calculation on the ellipsoid are significantly more complex than on the sphere. Formulas for the spherical form of a given map projection may be adapted for use with the ellipsoid by substitution of one of various "auxiliary latitudes" in place of the geodetic latitude. In using them, the ellipsoidal earth is, in effect, transformed to a sphere under certain restraints such as conformality or equal area, and the sphere is then projected onto a plane (Snyder, 1987). If the proper auxiliary latitudes are chosen, the sphere may have either true areas, true distances in certain directions, or conformality, relative to the ellipsoid. Spherical map projection formulas may then be used for the ellipsoid solely with the substitution of the appropriate auxiliary latitudes.

The rectifying, conformal and authalic latitudes are often used as auxiliary ones in cartography. The direct transformations form geodetic latitude to these auxiliary ones are expressed as transcendental functions or non-integrable ones. Adams (1921), Yang (1989, 2000) had derived forward expansions of these auxiliary latitudes form geodetic one through complicated formulation. Due to historical condition limitation, the derivation processes were done by hand and orders of these expansions could not be very high. Due to these reasons, the forward expansions for these auxiliary latitudes are reformulated by means of Mathematica. Readers will see that new expansions are expressed in a power series of the eccentricity of the reference ellipsoid *e* and extended up to tenth-order terms of *e* . The expansion processes become much easier under the system Mathematica.

#### **2.1. The forward expansion of the rectifying latitude**

The meridian arc from the equator *B* 0 to *B* is

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

where *X* is the meridian arc; *B* is the geodetic latitude; *a* is the semi-major axis of the reference ellipsoid;

(1) is an elliptic integral of the second kind and there is no analytical solution. Expanding the integrand by binomial theorem and itntegrating it item by item yield:

$$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

2 Cartography – A Tool for Spatial Analysis

**latitudes** 

appropriate auxiliary latitudes.

basis of cartography to a certain extent.

results indicate that the derivation efficiency can be significantly improved and formulas impossible to be obtained by hand can be easily derived with the help of Mathematica, which renovates the traditional analysis methods and enriches the mathematical theory

The main contents and research results presented in this chapter are organized as follows. In Section II, we discuss the direct transformations from geodetic latitude to three kinds of auxiliary latitudes often used in cartography, and the direct transformations from these auxiliary latitudes to geodetic latitude are studied in Section III. In Section IV, the direct expansions of transformations between meridian arc, isometric latitude, and authalic functions are derived. In Section V, we discuss the non-iterative expressions of the forward and inverse Gauss projections by complex numbers. Finally in Section VI, we make a brief summary. It is

assumed that the readers are somewhat conversant with Mathematica and its syntax.

**2. The forward expansions of the rectifying, conformal and authalic** 

Cartographers prefer to adopt sphere as a basis of the map projection for convenience since calculation on the ellipsoid are significantly more complex than on the sphere. Formulas for the spherical form of a given map projection may be adapted for use with the ellipsoid by substitution of one of various "auxiliary latitudes" in place of the geodetic latitude. In using them, the ellipsoidal earth is, in effect, transformed to a sphere under certain restraints such as conformality or equal area, and the sphere is then projected onto a plane (Snyder, 1987). If the proper auxiliary latitudes are chosen, the sphere may have either true areas, true distances in certain directions, or conformality, relative to the ellipsoid. Spherical map projection formulas may then be used for the ellipsoid solely with the substitution of the

The rectifying, conformal and authalic latitudes are often used as auxiliary ones in cartography. The direct transformations form geodetic latitude to these auxiliary ones are expressed as transcendental functions or non-integrable ones. Adams (1921), Yang (1989, 2000) had derived forward expansions of these auxiliary latitudes form geodetic one through complicated formulation. Due to historical condition limitation, the derivation processes were done by hand and orders of these expansions could not be very high. Due to these reasons, the forward expansions for these auxiliary latitudes are reformulated by means of Mathematica. Readers will see that new expansions are expressed in a power series of the eccentricity of the reference ellipsoid *e* and extended up to tenth-order terms of

> 2 223/2 <sup>0</sup> (1 ) (1 sin ) *B*

*X a e e B dB* (1)

*e* . The expansion processes become much easier under the system Mathematica.

**2.1. The forward expansion of the rectifying latitude** 

The meridian arc from the equator *B* 0 to *B* is

$$\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}{4096}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} \tag{3}$$

The rectifying latitude is defined as

$$\psi = \frac{X}{X(\frac{\pi}{2})} \cdot \frac{\pi}{2} \tag{4}$$

Inserting (2) into (4) yields

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

where

$$\begin{cases} \begin{aligned} \alpha\_2 &= K\_2 \, / K\_0\\ \alpha\_4 &= K\_4 \, / K\_0\\ \alpha\_6 &= K\_6 \, / K\_0\\ \alpha\_8 &= K\_8 \, / K\_0\\ \alpha\_{10} &= K\_{10} \, / K\_0 \end{aligned} \tag{6}$$

Yang (1989, 2000) gave an expansion similar to (5) but expanded up to sin8*B* . For simplicity and computing efficiency, it is better to expand (6) into a power series of the eccentricity. This process is easily done by means of Mathematica. As a result, (6) becomes:

$$\begin{cases} \alpha\_2 = -\frac{3}{8}e^2 - \frac{3}{16}e^4 - \frac{111}{1024}e^6 - \frac{141}{2048}e^8 - \frac{1533}{32768}e^{10} \\ \alpha\_4 = \frac{15}{256}e^4 + \frac{15}{256}e^6 + \frac{405}{8192}e^8 + \frac{165}{4096}e^{10} \\ \alpha\_6 = -\frac{35}{3072}e^6 - \frac{35}{2048}e^8 - \frac{4935}{262144}e^{10} \\ \alpha\_8 = \frac{315}{131072}e^8 + \frac{315}{65536}e^{10} \\ \alpha\_{10} = -\frac{693}{1310720}e^{10} \end{cases} \tag{7}$$

Mathematical Analysis in Cartography by Means of Computer Algebra System 5

 *BBBBB* sin 2 sin 4 sin6 sin8 (13)

> 

246 8 10

*eee e e*

4 6 8 10

*ee e e*

1 5 3 281 7 2 24 32 5760 240 5 7 697 93 48 80 11520 2240 13 461 1693 480 13440 53760

6 8 10

*eee*

8 10

*e e*

10

*e*

1237 131 161280 10080 367 161280

 , <sup>4</sup> , <sup>6</sup> , <sup>8</sup> 

Through the tedious expansion process, Yang (1989, 2000) gave a power series of the

in (13) derived by hand are only expanded to eighth-order terms of *e* . They are not accurate

In fact, Mathematica works perfectly in solving derivatives of any complicated functions. By means of Mathematica, the new derived forward expansion expanded to tenth-order terms

246810

2

 

 

 

**Table 1.** The comparison of coefficients of the forward expansion of conformal latitude derived by

Table 1 shows that the eighth order terms of *e* in coefficients given by Yang(1989, 2000) are

From the knowledge of mapping projection theory, the area of a section of a lune with a

cos sin 1 1 sin (1 ) (1 ) ln

*<sup>B</sup> B B <sup>e</sup> <sup>B</sup> Fa e dB a e*

(1 sin ) 2(1 sin ) 4 1 sin

*e B e B e eB* (15)

0 2 22 2 2

4

6

8

10

 

*BBBBB B* sin 2 sin 4 sin6 sin8 sin10 (14)

 

Coefficients derived by Yang(1989, 2000) Coefficients derived by the author

 

as

 

Due to that (11) is a very complicated transcendental function, the coefficients <sup>2</sup>

as expected and there are some mistakes in the eighth-order terms of *e* .

 

246 8

**2.3. The forward expansion of the authalic latitude** 

2 2 2 2

*eee e*

46 8

*ee e*

5 7 689 48 80 17920 13 1363 480 53760

1 5 3 1399 2 24 32 53760

6 8

*e e*

8

width of a unit interval of longitude *F* is

where *F* is also called authalic latitude function.

*e*

677 17520

Yang (1989, 2000) and the author

The derived coefficients in (13) and (14) are listed in Table 1 for comparison.

as the conventional usage in mathematical cartography.

2468

eccentricity *e* for the conformal latitude

of *e* reads

2

 

erroneous.

4

6

8

#### **2.2 The forward expansion of the conformal latitude**

Omitting the derivation process, the explicit expression for the isometric latitude *q* is

$$\begin{aligned} q &= \int\_0^B \frac{1 - e^2}{(1 - e^2 \sin^2 B)\cos B} dB = \ln \left| \tan \left( \frac{\pi}{4} + \frac{B}{2} \right) \left( \frac{1 - e \sin B}{1 + e \sin B} \right)^{e/2} \right| \\\ q &= \arctan \ln(\sin B) \cdot e \cdot \arctan \ln(e \sin B) \end{aligned} \tag{8}$$

For the sphere, putting *e* 0 and rewriting the geodetic latitude as the conformal one , (8) becomes

$$q = \ln\left[\tan(\frac{\pi}{4} + \frac{\rho}{2})\right] = \arctan\ln(\sin\rho) \tag{9}$$

Comparing (9) with (8) leads to

$$\tan(\frac{\pi}{4} + \frac{\rho}{2}) = \tan(\frac{\pi}{4} + \frac{B}{2})(\frac{1 - e\sin B}{1 + e\sin B})^{e/2} \tag{10}$$

Therefore, it holds

$$\varphi = 2 \arctan \left[ \tan(\frac{\pi}{4} + \frac{B}{2}) (\frac{1 - e \sin B}{1 + e \sin B})^{e/2} \right] - \frac{\pi}{2} \tag{11}$$

Since the eccentricity is small, the conformal latitude is close to the geodetic one. Though (11) is an analytical solution of , (11) is usually expanded into a power series of the eccentricity

$$\begin{aligned} \varphi(B, e) &= \varphi(B, 0) + \left. \frac{\partial \varphi}{\partial e} \right|\_{e=0} e + \frac{1}{2!} \frac{\partial^2 \varphi}{\partial e^2} \bigg|\_{e=0} e^2 + \frac{1}{3!} \frac{\partial^3 \varphi}{\partial e^3} \bigg|\_{e=0} e^3 + \dotsb \frac{1}{9!} \frac{\partial^9 \varphi}{\partial e^9} \bigg|\_{e=0} \\\ e^9 + \frac{1}{10!} \frac{\partial^{10} \varphi}{\partial e^{10}} \bigg|\_{e=0} e^{10} + \dotsb \end{aligned} \tag{12}$$

as the conventional usage in mathematical cartography.

4 Cartography – A Tool for Spatial Analysis

2

 

 

4

6

8

**2.2 The forward expansion of the conformal latitude** 

arctan h(sin )- arctan h( sin )

ln tan( ) arctan h(sin ) 4 2 *<sup>q</sup>*

1 sin /2 tan( ) tan( )( ) 4 2 4 2 1 sin

1 sin /2 2arctan tan( )( ) 4 2 1 sin 2

10 9 10 10 0

*e e e*

 

 

*e*

1 10! 10

 

0 2 2

becomes

Comparing (9) with (8) leads to

(11) is an analytical solution of

Therefore, it holds

eccentricity

2 4 6 8 10

/2 2

(8)

*B eB <sup>e</sup> e B*

*B eB <sup>e</sup> e B*

(9)

(11)

, (11) is usually expanded into a power series of the

23 9 2 3 23 9 0 00 0

 

*e ee e*

*e ee e*

1 1 sin ln tan (1 sin )cos 4 2 1 sin

*e BB e B*

(7)

, (8)

(10)

 

(12)

*ee e e e*

46 8 10

*ee e e*

3 3 111 141 1533 8 16 1024 2048 32768

6 8 10

*ee e*

8 10

Omitting the derivation process, the explicit expression for the isometric latitude *q* is

*Be e B*

 

 

For the sphere, putting *e* 0 and rewriting the geodetic latitude as the conformal one

 

11 1 ( , ) ( ,0) 2! 3! 9!

*Be B e e e*

Since the eccentricity is small, the conformal latitude is close to the geodetic one. Though

*<sup>e</sup> <sup>B</sup> e B <sup>e</sup> <sup>B</sup> q dB*

*e e*

15 15 405 165 256 256 8192 4096 35 35 4935 3072 2048 262144

10

*e*

315 315 131072 65536 693 1310720

Through the tedious expansion process, Yang (1989, 2000) gave a power series of the eccentricity *e* for the conformal latitude as

$$\boldsymbol{\sigma} = \boldsymbol{B} + \boldsymbol{\beta}\_2 \sin 2\boldsymbol{B} + \boldsymbol{\beta}\_4 \sin 4\boldsymbol{B} + \boldsymbol{\beta}\_6 \sin 6\boldsymbol{B} + \boldsymbol{\beta}\_8 \sin 8\boldsymbol{B} \tag{13}$$

Due to that (11) is a very complicated transcendental function, the coefficients <sup>2</sup> , <sup>4</sup> , <sup>6</sup> , <sup>8</sup> in (13) derived by hand are only expanded to eighth-order terms of *e* . They are not accurate as expected and there are some mistakes in the eighth-order terms of *e* .

In fact, Mathematica works perfectly in solving derivatives of any complicated functions. By means of Mathematica, the new derived forward expansion expanded to tenth-order terms of *e* reads

$$\mathbf{p} = \mathbf{B} + \beta\_2 \sin 2\mathbf{B} + \beta\_4 \sin 4\mathbf{B} + \beta\_6 \sin 6\mathbf{B} + \beta\_8 \sin 8\mathbf{B} + \beta\_{10} \sin 10\mathbf{B} \tag{14}$$

The derived coefficients in (13) and (14) are listed in Table 1 for comparison.


**Table 1.** The comparison of coefficients of the forward expansion of conformal latitude derived by Yang (1989, 2000) and the author

Table 1 shows that the eighth order terms of *e* in coefficients given by Yang(1989, 2000) are erroneous.

#### **2.3. The forward expansion of the authalic latitude**

From the knowledge of mapping projection theory, the area of a section of a lune with a width of a unit interval of longitude *F* is

$$F = a^2(1 - e^2) \int\_0^B \frac{\cos B}{\left(1 - e^2 \sin^2 B\right)^2} dB = a^2(1 - e^2) \left[ \frac{\sin B}{2(1 - e^2 \sin^2 B)} + \frac{1}{4e} \ln \frac{1 + e \sin B}{1 - e \sin B} \right] \tag{15}$$

where *F* is also called authalic latitude function.

Denote

$$A = \frac{1}{2(1 - e^2)} + \frac{1}{4e} \ln \frac{1 + e}{1 - e} \tag{16}$$

Mathematical Analysis in Cartography by Means of Computer Algebra System 7

Table 2 shows that the eighth-orders terms of *e* in coefficients given by Yang(1989, 2000) are

In order to validate the exactness and reliability of the forward expansions of rectifying, conformal and authalic latitudes derived by the author, one has examined their accuracies choosing the CGCS2000 (China Geodetic Coordinate System 2000) reference ellipsoid with parameters *a* 6378137m , *f* 1 / 298.257222101 (Chen, 2008; Yang, 2009), where *f* is the flattening. The accuracies of the forward expansions derived by Yang (1989, 2000) are also examined for comparison. The results show that the accuracy of the forward expansion of rectifying latitude derived by Yang (1989, 2000) is higher than 10-5″, while the accuracy of the forward expansion (5) derived by the author is higher than 10-7″. The accuracies of the forward expansion of conformal and authalic latitudes derived by Yang (1989, 2000) are higher than 10-4″, while the accuracies of the forward expansions derived by the author are higher than 10-8″ . The accuracies of forward expansions derived by the author are improved by 2~4 orders of magnitude compared to forward expansions derived by Yang (1989, 2000).

**3. The inverse expansions of rectifying, conformal and authalic latitudes** 

The inverse expansions of these auxiliary latitudes are much more difficult to derive than their forward ones. In this case, the differential equations are usually expressed as implicit functions of the geodetic latitude. There are neither any analytical solutions nor obvious expansions. For the inverse cases, to find geodetic latitude from auxiliary ones, one usually adopts iterative methods based on the forward expansions or an approximate series form. Yang (1989, 2000) had given the direct expansions of the inverse transformation by means of Lagrange series method, but their coefficients are expressed as polynomials of coefficients of the forward expansions, which are not convenient for practical use. Adams (1921) 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. Due to these reasons, new inverse expansions are derived using the power series method by means of Mathematica. Their coefficients are uniformly expressed as a power series of the eccentricity and extended up to

The processes to derive the inverse expansions using the power series method are as

1. To obtain their various order derivatives in terms of the chain rule of implicit

3. To integrate these series item by item and yield the final inverse expansions.

**3.1. The inverse expansions using the power series method** 

2. To compute the coefficients of their power series expansions;

erroneous.

tenth-order terms of *e* .

differentation;

follows:

**2.4. Accuracies of the forward expansions** 

Suppose that there is an imaginary sphere with a radius

$$R = a\sqrt{(1 - e^2)A} \tag{17}$$

and whose area from the spherical equator 0 to spherical latitude with a width of a unit interval of longitude is equal to the ellipsoidal area *F* , it holds

$$R^2 \sin \theta = a^2 (1 - e^2) A \sin \theta = F \tag{18}$$

Therefore, it yields

$$\mathcal{G} = \arcsin\left|\frac{1}{A}\left(\frac{\sin B}{2(1 - e^2 \sin^2 B)} + \frac{1}{4e} \ln \frac{1 + e \sin B}{1 - e \sin B}\right)\right|\tag{19}$$

where is authalic latitude. Yang(1989, 2000) gave its series expansion as

$$\mathcal{B} = B + \underline{\gamma}\_2 \sin 2B + \underline{\gamma}\_4 \sin 4B + \underline{\gamma}\_6 \sin 6B + \underline{\gamma}\_8 \sin 8B \tag{20}$$

(19) is a complicated transcendental function. It is almost impossible to derive its eighthorder derivate by hand. There are some mistakes in the high order terms of coefficients <sup>2</sup> , 4 , <sup>6</sup> , <sup>8</sup> .The new derived forward expansion expanded to tenth-order terms of *e* by means of Mathematica reads

$$\mathcal{B} = B + \mathcal{\chi}\_2 \sin 2B + \mathcal{\chi}\_4 \sin 4B + \mathcal{\chi}\_6 \sin 6B + \mathcal{\chi}\_8 \sin 8B + \mathcal{\chi}\_{10} \sin 10B \tag{21}$$

The derived coefficients in (20) and (21) are listed in Table 2 for comparison.


**Table 2.** The comparison of coefficients of the forward expansion of conformal latitude derived by Yang (1989, 2000) and the author

Table 2 shows that the eighth-orders terms of *e* in coefficients given by Yang(1989, 2000) are erroneous.

#### **2.4. Accuracies of the forward expansions**

6 Cartography – A Tool for Spatial Analysis

2

and whose area from the spherical equator

2 2

Suppose that there is an imaginary sphere with a radius

unit interval of longitude is equal to the ellipsoidal area *F* , it holds

2 22 *R a eA F* sin (1 ) sin

2468

 

The derived coefficients in (20) and (21) are listed in Table 2 for comparison.

1 11 ln 2(1 ) 4 1 *<sup>e</sup> <sup>A</sup>*

1 sin 1 1 sin arcsin ln

(19) is a complicated transcendental function. It is almost impossible to derive its eighthorder derivate by hand. There are some mistakes in the high order terms of coefficients <sup>2</sup>

246810

2

4

6

8

10

 

Coefficients derived by Yang(1989, 2000) Coefficients derived by the author

 **Table 2.** The comparison of coefficients of the forward expansion of conformal latitude derived by

is authalic latitude. Yang(1989, 2000) gave its series expansion as

 

2(1 sin ) 4 1 sin

*A e B e eB*

 

.The new derived forward expansion expanded to tenth-order terms of *e* by

 

*BBBBB B* sin 2 sin 4 sin6 sin8 sin10 (21)

*e e e*

0 to spherical latitude

 

*B e B*

 *BBBBB* sin 2 sin 4 sin6 sin8 (20)

> 

246 8 10

*eee e e*

46 8 10

*ee e e*

*ee e*

1 31 59 42811 605399 3 180 560 604800 11975040 17 61 76969 215431 360 1260 1814400 5987520 383 3347 1751791 45360 259200 119750400

68 10

8 10

*e e*

6007 201293 3628800 59875200 5839 171

 <sup>10</sup> 07200 *e*

(16)

<sup>2</sup> *R a eA* (1 ) (17)

(18)

with a width of a

(19)

,

Denote

Therefore, it yields

means of Mathematica reads

246 8

*eee e*

46 8

*ee e*

1 31 59 126853 3 180 560 518400 17 61 3622447 360 1260 94089600 383 6688039 43560 658627200

6 8

*e e*

8

*e*

Yang (1989, 2000) and the author

27787 23522400

where 

4 , <sup>6</sup> , <sup>8</sup> 

2

 

4

6

8

In order to validate the exactness and reliability of the forward expansions of rectifying, conformal and authalic latitudes derived by the author, one has examined their accuracies choosing the CGCS2000 (China Geodetic Coordinate System 2000) reference ellipsoid with parameters *a* 6378137m , *f* 1 / 298.257222101 (Chen, 2008; Yang, 2009), where *f* is the flattening. The accuracies of the forward expansions derived by Yang (1989, 2000) are also examined for comparison. The results show that the accuracy of the forward expansion of rectifying latitude derived by Yang (1989, 2000) is higher than 10-5″, while the accuracy of the forward expansion (5) derived by the author is higher than 10-7″. The accuracies of the forward expansion of conformal and authalic latitudes derived by Yang (1989, 2000) are higher than 10-4″, while the accuracies of the forward expansions derived by the author are higher than 10-8″ . The accuracies of forward expansions derived by the author are improved by 2~4 orders of magnitude compared to forward expansions derived by Yang (1989, 2000).

#### **3. The inverse expansions of rectifying, conformal and authalic latitudes**

The inverse expansions of these auxiliary latitudes are much more difficult to derive than their forward ones. In this case, the differential equations are usually expressed as implicit functions of the geodetic latitude. There are neither any analytical solutions nor obvious expansions. For the inverse cases, to find geodetic latitude from auxiliary ones, one usually adopts iterative methods based on the forward expansions or an approximate series form. Yang (1989, 2000) had given the direct expansions of the inverse transformation by means of Lagrange series method, but their coefficients are expressed as polynomials of coefficients of the forward expansions, which are not convenient for practical use. Adams (1921) 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. Due to these reasons, new inverse expansions are derived using the power series method by means of Mathematica. Their coefficients are uniformly expressed as a power series of the eccentricity and extended up to tenth-order terms of *e* .

#### **3.1. The inverse expansions using the power series method**

The processes to derive the inverse expansions using the power series method are as follows:


One can image that these procedures are quite complicated. Mathematica output shows that the expression of the sixth order derivative is up to 40 pages long! Therefore, it is unimaginable to derive the so long expression by hand. These procedures, however, will become much easier and be significantly simplified by means of Mathematica. As a result, the more simple and accurate expansions yield.

#### *3.1.1. The inverse expansion of the rectifying Latitude*

Differentiation to the both sides of (1) yields

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

Mathematical Analysis in Cartography by Means of Computer Algebra System 9

 

 

(31)

(32)

(33)

2468 10

(30)

246810

24 6 8 10

1 sin sin sin sin sin *dB AAAAA*

4 6 8 10

Multiplying *K*0 and integrating at the both sides of (30) give the inverse expansion of

 sin 2 sin 4 sin6 sin8 sin10 

46 8 10

3 3 213 255 20861 8 16 2048 4096 524288

*ae e e e e*

6 8 10

21 21 533 197 256 256 8192 4096 151 151 5019 6144 4096 131072 1097 1097 131072 65536 8011 2621440

*aee e e*

8 10

10

*a ee*

*d*

*a e*

*aee e*

2 4 6 8 10

2

(34)

(35)

2 2 1 cos (1 sin )cos *d e dB e BB*

2 2 2 (1 sin )cos (1 )cos *dB e B B*

*e*

<sup>246810</sup> *Ba a a a a*

3 27 729 4329 381645 2 8 128 512 32768 21 621 11987 757215 8 64 512 16384 151 775 621445 32 32 8192

6 8 10

*Ae e e e*

*A ee e e e*

8 10

*A ee e*

10

1097 57607 128 1024 8011 512

*Ae e*

0

where

2

 

 

2

 

 

 

Differentiating the both sides of (10) yields

Therefore, it holds

4

6

8

10

*3.1.2. The inverse expansion of the conformal latitude* 

rectifying latitude as

where

4

6

8

10

 

*A e*

*K d*

From (4) and (2), one knows

$$X = \mathfrak{a} (1 - e^2) K\_0 \mathfrak{y} \tag{23}$$

$$\text{Inserting (23) into (22) yields }$$

$$\frac{d\mathbf{B}}{K\_0 d\mu} = (1 - e^2 \sin^2 B)^{3/2} \tag{24}$$

To expand (24) into a power series of sin, we introduce the following new variable

$$t = \sin \varphi \tag{25}$$

therefore

$$\frac{d\nu}{dt} = \frac{1}{\cos \nu} \tag{26}$$

and then denote

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

Making use of the chain rule of implicit differentiation

$$f\_t' = \frac{df}{d\mathbf{B}} \frac{d\mathbf{B}}{d\boldsymbol{\nu}} \frac{d\boldsymbol{\nu}}{d\mathbf{t}} + \frac{df}{d\boldsymbol{\nu}} \frac{d\boldsymbol{\nu}}{d\mathbf{t}}, \quad f\_t'' = \frac{df'}{d\mathbf{B}} \frac{d\mathbf{B}}{d\boldsymbol{\nu}} \frac{d\boldsymbol{\nu}}{d\mathbf{t}} + \frac{df'}{d\boldsymbol{\nu}} \frac{d\boldsymbol{\nu}}{d\mathbf{t}}, \quad \cdots \tag{28}$$

It is easy to expand (27) into a power series of sin

$$f(t) = f(0) + f\_t'(0)t + \frac{1}{2!}f\_t'(0)t^2 + \frac{1}{3!}f\_t'(0)t^3 + \dots + \frac{1}{10!}f\_t^{(10)}(0)t^{10} + \dotsb \tag{29}$$

Omitting the detailed procedure, one arrives at

Mathematical Analysis in Cartography by Means of Computer Algebra System 9

$$\frac{d\mathcal{B}}{K\_0 d\boldsymbol{\wp}} = 1 + A\_2 \sin^2 \boldsymbol{\wp} + A\_4 \sin^4 \boldsymbol{\wp} + A\_6 \sin^6 \boldsymbol{\wp} + A\_8 \sin^8 \boldsymbol{\wp} + A\_{10} \sin^{10} \boldsymbol{\wp} \tag{30}$$

where

8 Cartography – A Tool for Spatial Analysis

From (4) and (2), one knows

Inserting (23) into (22) yields

therefore

and then denote

To expand (24) into a power series of sin

the more simple and accurate expansions yield.

Differentiation to the both sides of (1) yields

*3.1.1. The inverse expansion of the rectifying Latitude* 

2 2 3/2

2 2 3/2

Making use of the chain rule of implicit differentiation

It is easy to expand (27) into a power series of sin

Omitting the detailed procedure, one arrives at

0

*d dt* 

0

*K d*

*K d*

One can image that these procedures are quite complicated. Mathematica output shows that the expression of the sixth order derivative is up to 40 pages long! Therefore, it is unimaginable to derive the so long expression by hand. These procedures, however, will become much easier and be significantly simplified by means of Mathematica. As a result,

> 2 2 2 3/2 (1 ) (1 sin )

2 <sup>0</sup> *X a eK* (1 )

(1 sin ) *dB e B*

*t* sin

( ) (1 sin ) *dB f t e B*

, , *t t df df df df dB d d dB d d f f dB d dt d dt dB d dt d dt*

1 cos

(28)

11 1 23 (10) 10 ( ) (0) (0) (0) (0) (0) 2! 3! 10! *tt t t ft f f t f t f t f t* (29)

(22)

(24)

, we introduce the following new variable

(23)

(25)

(26)

(27)

*dX a e dB e B*

$$\begin{cases} A\_2 = -\frac{3}{2}e^2 - \frac{27}{8}e^4 - \frac{729}{128}e^6 - \frac{4329}{512}e^8 - \frac{381645}{32768}e^{10} \\ A\_4 = \frac{21}{8}e^4 + \frac{621}{64}e^6 + \frac{11987}{512}e^8 + \frac{757215}{16384}e^{10} \\ A\_6 = -\frac{151}{32}e^6 - \frac{775}{32}e^8 - \frac{621445}{8192}e^{10} \\ A\_8 = \frac{1097}{128}e^8 + \frac{57607}{1024}e^{10} \\ A\_{10} = -\frac{8011}{512}e^{10} \end{cases} \tag{31}$$

Multiplying *K*0 and integrating at the both sides of (30) give the inverse expansion of rectifying latitude 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{32}$$

where

$$\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{33}$$

#### *3.1.2. The inverse expansion of the conformal latitude*

Differentiating the both sides of (10) yields

$$\frac{d\phi}{\cos\phi} = \frac{1 - e^2}{(1 - e^2 \sin^2 B)\cos B} dB \tag{34}$$

Therefore, it holds

$$\frac{d\mathcal{B}}{d\rho} = \frac{(1 - e^2 \sin^2 \mathcal{B}) \cos \mathcal{B}}{(1 - e^2) \cos \rho} \tag{35}$$

For the same reason, we introduce the following new variable

$$t = \sin \varphi \tag{36}$$

Mathematical Analysis in Cartography by Means of Computer Algebra System 11

(46)

(45)

(42)

(43)

(44)

. Using the chain rule of implicit function

 

(47)

(48)

 

0 2 22

2 22 (1 sin ) cos cos

cos

2 4 6 8 10

46 8 10

4 41 4108 58427 28547 3 15 945 9450 3465 92 6574 223469 2768558 45 945 14175 93555 3044 28901 21018157 945 1890 467775

*C ee e e e*

*Ce e e e*

68 10

8 10

To get the inverse expansion of the authalic latitude, one integrates (46) and arrives at

<sup>246810</sup> *Bc c c c c*

 sin 2 sin 4 sin6 sin8 sin10 

*C ee e*

10

24236 2086784 4725 66825 768272 93555

*Ce e*

2 22 (1 sin ) cos ( )

*dB A e B f t d B*

2468 10

<sup>246810</sup> ( ) sin sin sin sin sin *dB ft A C C C C C <sup>d</sup>*

(1 sin ) *<sup>B</sup> <sup>B</sup> <sup>A</sup> dB*

*e B*

cos sin

*dB A e B d B*

For the same reason, we introduce the folllowing new variable

*t* sin

(45) can be expanded into a power series of sin

differentiation, one similarly arrives at

2

 

  4

6

8

10

*C e*

and then denote

where

Differentiating the both sides of (42) yields

and then denote

$$f(t) = \frac{dB}{d\rho} = \frac{(1 - e^2 \sin^2 B)\cos B}{(1 - e^2)\cos\rho} \tag{37}$$

Using the same procedure as described in the former section, one arrives at

$$\frac{dB}{d\rho} = \frac{1}{1 - e^2} + B\_2 \sin^2 \phi + B\_4 \sin^4 \phi + B\_6 \sin^6 \phi + B\_8 \sin^8 \phi + B\_{10} \sin^{10} \phi \tag{38}$$

where

$$\begin{cases} B\_2 = -2e^2 - \frac{11}{2}e^4 - \frac{21}{2}e^6 - 17e^8 - 25e^{10} \\ B\_4 = \frac{14}{3}e^4 + \frac{62}{3}e^6 + \frac{1369}{24}e^8 + \frac{3005}{24}e^{10} \\ B\_6 = -\frac{56}{5}e^6 - \frac{614}{9}e^8 - \frac{4909}{20}e^{10} \\ B\_8 = \frac{8558}{315}e^8 + \frac{7367}{35}e^{10} \\ B\_{10} = -\frac{4174}{63}e^{10} \end{cases} \tag{39}$$

Integrating the both sides of (38) gives the inverse expansion of conformal latitude as

$$B = \phi + b\_2 \sin 2\phi + b\_4 \sin 4\phi + b\_6 \sin 6\phi + b\_8 \sin 8\phi + b\_{10} \sin 10\phi \tag{40}$$

where

$$\begin{cases} b\_2 = \frac{1}{2}e^2 + \frac{5}{24}e^4 + \frac{1}{12}e^6 + \frac{13}{360}e^8 + \frac{3}{160}e^{10} \\ b\_4 = \frac{7}{48}e^4 + \frac{29}{240}e^6 + \frac{811}{11520}e^8 + \frac{81}{2240}e^{10} \\ b\_6 = \frac{7}{120}e^6 + \frac{81}{1120}e^8 + \frac{3029}{53760}e^{10} \\ b\_8 = \frac{4279}{161280}e^8 + \frac{883}{20160}e^{10} \\ b\_{10} = \frac{2087}{161280}e^{10} \end{cases} \tag{41}$$

*3.1.3 The inverse expansion of the authalic latitude* 

Inserting (18) into (15) yields

Mathematical Analysis in Cartography by Means of Computer Algebra System 11

$$A\sin\mathcal{G} = \int\_0^B \frac{\cos B}{\left(1 - e^2 \sin^2 B\right)^2} dB \tag{42}$$

Differentiating the both sides of (42) yields

$$\frac{d\mathcal{B}}{d\mathcal{B}} = \frac{A(1 - e^2 \sin^2 B)^2 \cos \mathcal{B}}{\cos B} \tag{43}$$

For the same reason, we introduce the folllowing new variable

$$t = \sin \mathcal{G}\tag{44}$$

and then denote

10 Cartography – A Tool for Spatial Analysis

and then denote

where

where

For the same reason, we introduce the following new variable

2 2 2

(38)

(37)

(1 sin )cos ( ) (1 )cos *dB e B B f t d*

2 246810 <sup>1</sup> sin sin sin sin sin

2 4 6 8 10

11 21 <sup>2</sup> 17 25 2 2 14 62 1369 3005 3 3 24 24 56 614 4909 5 9 20

*Be e eee*

*Bee e e*

4 6 8 10

6 8 10

8 10

Integrating the both sides of (38) gives the inverse expansion of conformal latitude as

<sup>246810</sup> *Bb b b b b*

 sin 2 sin 4 sin6 sin8 sin10 

> 1 5 1 13 3 2 24 12 360 160 7 29 811 81 48 240 11520 2240 7 81 3029 120 1120 53760 4279 883 161280 20160

*be e e e e*

*be e e e*

2 4 6 8 10

4 6 8 10

6 8 10

8 10

10

*b ee*

*be e e*

2087 161280

*b e*

10

*B ee e*

8558 7367 315 35 4174 63

*Bee*

*e*

Using the same procedure as described in the former section, one arrives at

2

 

 

4

6

8

10

2

 

 

 

*3.1.3 The inverse expansion of the authalic latitude* 

Inserting (18) into (15) yields

4

6

8

10

 

*B e*

2468 10

*dB BBBBB*

(36)

 

(39)

(40)

(41)

 

sin *t*

1

*d e*

$$f(t) = \frac{dB}{d\mathcal{B}} = \frac{A(1 - e^2 \sin^2 B)^2 \cos \mathcal{B}}{\cos B} \tag{45}$$

(45) can be expanded into a power series of sin . Using the chain rule of implicit function differentiation, one similarly arrives at

$$f(t) = \frac{d\mathcal{B}}{d\mathcal{B}} = A + \mathcal{C}\_2 \sin^2 \mathcal{B} + \mathcal{C}\_4 \sin^4 \mathcal{B} + \mathcal{C}\_6 \sin^6 \mathcal{B} + \mathcal{C}\_8 \sin^8 \mathcal{B} + \mathcal{C}\_{10} \sin^{10} \mathcal{B} \tag{46}$$

where

$$\begin{cases} \mathbf{C}\_{2} = -\frac{4}{3}e^{2} - \frac{41}{15}e^{4} - \frac{4108}{945}e^{6} - \frac{58427}{9450}e^{8} - \frac{28547}{3465}e^{10} \\ \mathbf{C}\_{4} = \frac{92}{45}e^{4} + \frac{6574}{945}e^{6} + \frac{223469}{14175}e^{8} + \frac{2768558}{93555}e^{10} \\ \mathbf{C}\_{6} = -\frac{3044}{945}e^{6} - \frac{28901}{1890}e^{8} - \frac{21018157}{467775}e^{10} \\ \mathbf{C}\_{8} = \frac{24236}{4725}e^{8} + \frac{2086784}{66825}e^{10} \\ \mathbf{C}\_{10} = -\frac{768272}{93555}e^{10} \end{cases} \tag{47}$$

To get the inverse expansion of the authalic latitude, one integrates (46) and arrives at

$$B = \mathcal{Y} + c\_2 \sin 2\mathcal{Y} + c\_4 \sin 4\mathcal{Y} + c\_6 \sin 6\mathcal{Y} + c\_8 \sin 8\mathcal{Y} + c\_{10} \sin 10\mathcal{Y} \tag{48}$$

where

$$\begin{cases} c\_2 = \frac{1}{3}e^2 + \frac{31}{180}e^4 + \frac{517}{5040}e^6 + \frac{120389}{181400}e^8 + \frac{1362253}{29937600}e^{10} \\ c\_4 = \frac{23}{360}e^4 + \frac{251}{3780}e^6 + \frac{102287}{1814400}e^8 + \frac{450739}{997920}e^{10} \\ c\_6 = \frac{761}{45360}e^6 + \frac{47561}{1814400}e^8 + \frac{434501}{14968800}e^{10} \\ c\_8 = \frac{6059}{1209600}e^8 + \frac{625511}{59875200}e^{10} \\ c\_{10} = \frac{48017}{29937600}e^{10} \end{cases} \tag{49}$$

Mathematical Analysis in Cartography by Means of Computer Algebra System 13

1. To apply the Lagrange theorem to a trigonometric series;

**3.4. Accuracies of the inverse expansions** 

**latitude and authalic latitude function** 

latitude (14), one obtains the conformal latitude

latitude *q* . The whole formulas are as follows:

been examined choosing the CGCS2000 reference ellipsoid.

eccentricity.

by Yang (1989, 2000).

Mathematica.

**isometric latitude** 

formulas, too.

2. To write the inverse expansions of the rectifying, conformal and authalic latitude; 3. To express the coefficients of the inverse expansions as a power series of the

The detailed derivation processes are given by Li (2010). Confined to the length of the chapter, they are also omitted. Comparing the results derived by this method with those in 3.1 and 3.2, one will find that they are all consistent with each other even though they are also formulated in different ways. This fact substantiates the correctness of the derived

The accuracies of the inverse expansions derived by Yang (1989, 2000) and the author has

The results show that the accuracy of the inverse expansion of rectifying latitude is higher than 10-5″, while the accuracy of the inverse expansion (32) derived by the author is higher than 10-7″. The accuracies of the inverse expansion of conformal and authalic latitudes derived by Yang (1989, 2000) are higher than 10-4″, while the accuracies of the inverse expansions derived by the author are higher than 10-8″. The accuracies of inverse expansions derived by the author are improved by 2~4 orders of magnitude compared to those derived

**4. The direct expansions of transformations between meridian arc, isometric** 

**4.1. The direct expansions of transformations between meridian arc and** 

Inserting the known meridian arc *X* into (23) yields the rectifying latitude

*4.1.1. The direct expansion of the transformation from meridian arc to isometric latitude* 

inverse expansion of the rectifying latitude (32) and the forward expansion of the conformal

. Inserting it into (9) yields the isometric

. Using the

The meridian arc, isometric latitude and authalic latitude function are functions of rectifying, conformal and authalic latitudes correspondingly. The transformations between the three variables are indirectly realized by computing the geodetic latitude in the past literatures such as Yang (1989, 2000), Snyder (1987). The computation processes are tedious and time-consuming. In order to simplify the computation processes and improve the computation efficiency, the direct expansions of transformations between meridian arc, isometric latitude and authalic latitude function are comprehensively derived by means of

#### **3.2. The inverse expansions using the Hermite interpolation method**

In mathematical analysis, interpolation with functional values and their derivative values is called Hermite interpolation. The processes to derive the inverse expansions using this method are as follows:


The detailed derivation processes are given by Li (2008, 2010). Confined to the length of the chapter, they are omitted. Comparing the results derived by this method with those in 3.1, one will find that they are consistent with each other even though they are formulated in different ways. This fact substantiates the correctness of the derived formulas.

#### **3.3. The inverse expansions using the Lagrange's theorem method**

We wish to investigate the inversion of an equation such as

$$y = \mathfrak{x} + f(\mathfrak{x})\tag{50}$$

with *fx x* ( ) and *y x* . The Lagrange's theorem states that in a suitable domain the solution of (50) is

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

The proof of this theorem is given by Whittaker (1902) and Peter (2008).

The processes to derive the inverse expansions using the Lagrange series method are as follows:


The detailed derivation processes are given by Li (2010). Confined to the length of the chapter, they are also omitted. Comparing the results derived by this method with those in 3.1 and 3.2, one will find that they are all consistent with each other even though they are also formulated in different ways. This fact substantiates the correctness of the derived formulas, too.

#### **3.4. Accuracies of the inverse expansions**

12 Cartography – A Tool for Spatial Analysis

2

 

4

6

8

10

arcs with coefficients to be determined;

*c*

method are as follows:

coefficients.

solution of (50) is

follows:

2 4 6 8 10

(49)

46 8 10

1 31 517 120389 1362253 3 180 5040 181400 29937600 23 251 102287 450739 360 3780 1814400 997920 761 47561 434501 45360 1814400 14968800

*ce e e e e*

68 10

In mathematical analysis, interpolation with functional values and their derivative values is called Hermite interpolation. The processes to derive the inverse expansions using this

1. To suppose the inverse expansions are expressed in a series of the sines of the multiple

3. To solve linear equations according to interpolation constraints and obtain the

The detailed derivation processes are given by Li (2008, 2010). Confined to the length of the chapter, they are omitted. Comparing the results derived by this method with those in 3.1, one will find that they are consistent with each other even though they are formulated in

with *fx x* ( ) and *y x* . The Lagrange's theorem states that in a suitable domain the

1

The proof of this theorem is given by Whittaker (1902) and Peter (2008).

*<sup>n</sup> <sup>n</sup>*

The processes to derive the inverse expansions using the Lagrange series method are as

1 1

( 1) [ ( )] ! *n n <sup>n</sup>*

*d x y f y n dy* 

*y x fx* ( ) (50)

(51)

2. To compute the functional values and their derivative values at specific points;

different ways. This fact substantiates the correctness of the derived formulas.

**3.3. The inverse expansions using the Lagrange's theorem method** 

We wish to investigate the inversion of an equation such as

8 10

*ce e e e*

*ce e e*

6059 625511 1209600 59875200

*e*

**3.2. The inverse expansions using the Hermite interpolation method** 

*ce e*

48017

 <sup>10</sup> 29937600

where

The accuracies of the inverse expansions derived by Yang (1989, 2000) and the author has been examined choosing the CGCS2000 reference ellipsoid.

The results show that the accuracy of the inverse expansion of rectifying latitude is higher than 10-5″, while the accuracy of the inverse expansion (32) derived by the author is higher than 10-7″. The accuracies of the inverse expansion of conformal and authalic latitudes derived by Yang (1989, 2000) are higher than 10-4″, while the accuracies of the inverse expansions derived by the author are higher than 10-8″. The accuracies of inverse expansions derived by the author are improved by 2~4 orders of magnitude compared to those derived by Yang (1989, 2000).

### **4. The direct expansions of transformations between meridian arc, isometric latitude and authalic latitude function**

The meridian arc, isometric latitude and authalic latitude function are functions of rectifying, conformal and authalic latitudes correspondingly. The transformations between the three variables are indirectly realized by computing the geodetic latitude in the past literatures such as Yang (1989, 2000), Snyder (1987). The computation processes are tedious and time-consuming. In order to simplify the computation processes and improve the computation efficiency, the direct expansions of transformations between meridian arc, isometric latitude and authalic latitude function are comprehensively derived by means of Mathematica.

### **4.1. The direct expansions of transformations between meridian arc and isometric latitude**

#### *4.1.1. The direct expansion of the transformation from meridian arc to isometric latitude*

Inserting the known meridian arc *X* into (23) yields the rectifying latitude . Using the inverse expansion of the rectifying latitude (32) and the forward expansion of the conformal latitude (14), one obtains the conformal latitude . Inserting it into (9) yields the isometric latitude *q* . The whole formulas are as follows:

$$\begin{cases} \mathscr{Y} = \frac{X}{a(1-e^2)K\_0} \\ B = \mathscr{Y} + a\_2 \sin 2\mathscr{Y} + a\_4 \sin 4\mathscr{Y} + a\_6 \sin 6\mathscr{Y} + a\_8 \sin 8\mathscr{Y} + a\_{10} \sin 10\mathscr{Y} \\ \mathscr{O} = B + \mathscr{O}\_2 \sin 2B + \mathscr{O}\_4 \sin 4B + \mathscr{O}\_6 \sin 6B + \mathscr{O}\_8 \sin 8B + \mathscr{O}\_{10} \sin 10B \\ q = \arctan(\sin \varphi) \end{cases} \tag{52}$$

Mathematical Analysis in Cartography by Means of Computer Algebra System 15

 

 

> 

(58)

(59)

 

> 

(56)

(57)

0 2 <sup>4</sup> <sup>6</sup> <sup>8</sup> <sup>10</sup>

*X aj j j j j j*

sin 2 sin 4 sin6 sin8 sin10

24 6 8 10

2 4 6 8 10

4 6 8 10

*e*

6 8 10

175087 165150720

**4.2. The direct expansions of transformations between meridian arc and authalic** 

The whole formulas for the transformation from meridian arc to authalic latitude function

2 4 6 8 10 246810

Expanding *F* as a power series of the eccentricity *e* at *e* 0 by means of Mathematica yields the direct expansion of the transformation from meridian arc to authalic latitude

sin sin 3 sin 5 sin7 sin9 sin11

 

13 5 7 9 11

 

 

*BBBBB B*

*Ba a a a a*

 

 

> 

sin 2 sin 4 sin6 sin8 sin10 sin 2 sin 4 sin6 sin8 sin10

 

*4.2.1. The direct expansion of the transformation from meridian arc to authalic latitude* 

1 3 5 175 441 <sup>1</sup> 4 64 256 16384 65536 1 1 9 901 16381 8 96 1024 184320 5898240 13 17 311 18931 768 5120 737280 20643840 61 899 14977 15360 430080 27525120

*je e e e e*

*je e e e*

*je e e*

8 10

10

*j ee e e e*

arcsin(tanh )

0

 

2

4

6

8

 

10

49561 41287680

34729 82575360

*j e*

*j e*

2 0

(1 )

*X a eK*

2

2 0

(1 )

*X a eK*

*F R*

2

*F a*

sin

where

**latitude function** 

*function* 

function

where

are as follows:

*q*

Since the coefficients 2*<sup>i</sup> a* , 2*<sup>i</sup>* ( 1,2, 5 *i* ) are expressed in a power series of the eccentricity, one could expand *q* as a power series of the eccentricity *e* at *e* 0 in order to obtain the direct expansion of the transformation from *X* to *q* . It is hardly completed by hand, but could be easily realized by means of Mathematica. Omitting the main operations by means of Mathematica yields the direct expansion of the transformation from meridian arc to isometric latitude

$$\begin{cases} \boldsymbol{\nu} = \frac{\boldsymbol{X}}{a(1 - e^2)\boldsymbol{K}\_0} \\ \boldsymbol{q} = \operatorname{arctanh}(\sin \boldsymbol{\nu}) + \boldsymbol{\xi}\_1 \sin \boldsymbol{\nu} + \boldsymbol{\xi}\_3 \sin 3\boldsymbol{\nu} + \boldsymbol{\xi}\_5 \sin 5\boldsymbol{\nu} + \boldsymbol{\xi}\_7 \sin 7\boldsymbol{\nu} + \boldsymbol{\xi}\_9 \sin 9\boldsymbol{\nu} \end{cases} \tag{53}$$

where

$$\begin{cases} \dot{\xi}\_1 = -\frac{1}{4}e^2 - \frac{1}{64}e^4 + \frac{1}{3072}e^6 + \frac{33}{16384}e^8 + \frac{2363}{1310720}e^{10} \\ \dot{\xi}\_3 = -\frac{1}{96}e^4 - \frac{13}{3072}e^6 - \frac{13}{8192}e^8 - \frac{1057}{1966080}e^{10} \\ \dot{\xi}\_5 = -\frac{11}{7680}e^6 - \frac{29}{24576}e^8 - \frac{2897}{3932160}e^{10} \\ \dot{\xi}\_7 = -\frac{25}{86016}e^8 - \frac{727}{1966080}e^{10} \\ \dot{\xi}\_9 = -\frac{53}{737280}e^{10} \end{cases} \tag{54}$$

#### *4.1.2. The direct expansion of the transformation from isometric latitude to meridian arc*

The whole formulas for the transformation from isometric latitude to meridian arc are as follows:

$$\begin{cases} \varphi = \arcsin(\tanh q) \\ B = \varphi + b\_2 \sin 2\varphi + b\_4 \sin 4\varphi + b\_6 \sin 6\varphi + b\_8 \sin 8\varphi + b\_{10} \sin 10\varphi \\ \nu = B + a\_2 \sin 2B + a\_4 \sin 4B + a\_6 \sin 6B + a\_8 \sin 8B + a\_{10} \sin 10B \\ X = a(1 - e^2)K\_0 \nu \end{cases} \tag{55}$$

Expanding *X* as a power series of the eccentricity *e* at *e* 0 by means of Mathematica yields the direct expansion of the transformation from isometric latitude to meridian arc

#### Mathematical Analysis in Cartography by Means of Computer Algebra System 15

$$\begin{cases} \rho = \arcsin(\tanh q) \\ X = a \left( j\_0 \rho + j\_2 \sin 2\rho + j\_4 \sin 4\rho + j\_6 \sin 6\rho + j\_8 \sin 8\rho + j\_{10} \sin 10\rho \right) \end{cases} \tag{56}$$

where

14 Cartography – A Tool for Spatial Analysis

*q*

Since the coefficients 2*<sup>i</sup> a* , 2*<sup>i</sup>*

isometric latitude

where

follows:

*q*

 

2 0

1

 

 

3

5

7

9

2 0

(1 )

*X a eK*

arcsin(tanh )

*q*

 

  

(1 )

*X a eK* 2 0

arctanh(sin )

(1 )

*X a eK*

> 2 4 6 8 10 2 4 6 8 10

 

one could expand *q* as a power series of the eccentricity *e* at *e* 0 in order to obtain the direct expansion of the transformation from *X* to *q* . It is hardly completed by hand, but could be easily realized by means of Mathematica. Omitting the main operations by means of Mathematica yields the direct expansion of the transformation from meridian arc to

*BBBBB B*

*Ba a a a a*

 

 

sin 2 sin 4 sin6 sin8 sin10 sin 2 sin 4 sin6 sin8 sin10

 

13 5 7 9

2 4 6 8 10

*ee e e e*

4 6 8 10

*eee e*

1 1 1 33 2363 4 64 3072 16384 1310720

6 8 10

*ee e*

*4.1.2. The direct expansion of the transformation from isometric latitude to meridian arc* 

The whole formulas for the transformation from isometric latitude to meridian arc are as

246810 246810

 

Expanding *X* as a power series of the eccentricity *e* at *e* 0 by means of Mathematica yields the direct expansion of the transformation from isometric latitude to meridian arc

*Bb b b b b*

sin 2 sin 4 sin6 sin8 sin10 sin 2 sin 4 sin6 sin8 sin10

*BBBBB B*

 

1 13 13 1057 96 3072 8192 1966080 11 29 2897 7680 24576 3932160

8 10

*e e*

10

*e*

53 737280

25 727 86016 1966080

 

arctanh(sin ) sin sin 3 sin 5 sin7 sin9

 

( 1,2, 5 *i* ) are expressed in a power series of the eccentricity,

  (52)

 

> 

> >

   

(54)

(55)

(53)

$$\begin{cases} j\_0 = 1 - \frac{1}{4}e^2 - \frac{3}{64}e^4 - \frac{5}{256}e^6 - \frac{175}{16384}e^8 - \frac{441}{65536}e^{10} \\ j\_2 = \frac{1}{8}e^2 - \frac{1}{96}e^4 - \frac{9}{1024}e^6 - \frac{901}{184320}e^8 - \frac{16381}{5898240}e^{10} \\ j\_4 = \frac{13}{768}e^4 + \frac{17}{5120}e^6 - \frac{311}{737280}e^8 - \frac{18931}{20643840}e^{10} \\ j\_6 = \frac{61}{15360}e^6 + \frac{899}{430080}e^8 + \frac{14977}{27525210}e^{10} \\ j\_8 = \frac{49561}{41287680}e^8 + \frac{175087}{165150720}e^{10} \\ j\_{10} = \frac{34729}{82575360}e^{10} \end{cases} \tag{57}$$

#### **4.2. The direct expansions of transformations between meridian arc and authalic latitude function**

*4.2.1. The direct expansion of the transformation from meridian arc to authalic latitude function* 

The whole formulas for the transformation from meridian arc to authalic latitude function are as follows:

$$\begin{cases} \mathcal{Y} = \frac{X}{a(1 - e^2)K\_0} \\ B = \mathcal{Y} + a\_2 \sin 2\mathcal{Y} + a\_4 \sin 4\mathcal{Y} + a\_6 \sin 6\mathcal{Y} + a\_8 \sin 8\mathcal{Y} + a\_{10} \sin 10\mathcal{Y} \\ \mathcal{G} = B + \mathcal{Y}\_2 \sin 2B + \mathcal{Y}\_4 \sin 4B + \mathcal{Y}\_6 \sin 6B + \mathcal{Y}\_8 \sin 8B + \mathcal{Y}\_{10} \sin 10B \\ F = R^2 \sin \mathcal{G} \end{cases} \tag{58}$$

Expanding *F* as a power series of the eccentricity *e* at *e* 0 by means of Mathematica yields the direct expansion of the transformation from meridian arc to authalic latitude function

$$\begin{cases} \boldsymbol{\nu} = \frac{\boldsymbol{X}}{a(1 - e^2)\boldsymbol{K}\_0} \\ \boldsymbol{F} = a^2 \left( \boldsymbol{\varepsilon}\_1 \sin \boldsymbol{\nu} + \boldsymbol{\varepsilon}\_3 \sin 3\boldsymbol{\nu} + \boldsymbol{\varepsilon}\_5 \sin 5\boldsymbol{\nu} + \boldsymbol{\varepsilon}\_7 \sin 7\boldsymbol{\nu} + \boldsymbol{\varepsilon}\_9 \sin 9\boldsymbol{\nu} + \boldsymbol{\varepsilon}\_{11} \sin 11\boldsymbol{\nu} \right) \end{cases} \tag{59}$$

where

$$\begin{cases} \varepsilon\_1 = 1 - \frac{5}{16}e^2 - \frac{17}{256}e^4 - \frac{121}{4096}e^6 - \frac{137}{8192}e^8 - \frac{1407}{131072}e^{10} \\ \varepsilon\_3 = \frac{1}{48}e^2 + \frac{1}{384}e^4 - \frac{103}{196608}e^8 - \frac{1775}{3145728}e^{10} \\ \varepsilon\_5 = \frac{3}{1280}e^4 + \frac{43}{30720}e^6 + \frac{17}{24576}e^8 + \frac{467}{1572864}e^{10} \\ \varepsilon\_7 = \frac{37}{86016}e^6 + \frac{5}{10752}e^8 + \frac{563}{1572864}e^{10} \\ \varepsilon\_9 = \frac{59}{589824}e^8 + \frac{1853}{11796480}e^{10} \\ \varepsilon\_{11} = \frac{1543}{57671680}e^{10} \end{cases} \tag{60}$$

Mathematical Analysis in Cartography by Means of Computer Algebra System 17

 

> 

(66)

(64)

(65)

 

> 

**4.3. The direct expansions of transformations between isometric latitude and** 

*4.3.1. The direct expansion of the transformation from isometric latitude to authalic* 

The whole formulas for the transformation from isometric latitude to authalic latitude

246810 246810

 

Expanding *F* as a power series of the eccentricity *e* at *e* 0 by means of Mathematica yields the direct expansion of the transformation from isometric latitude to authalic latitude

1 3 5 7 9 11

268 10

*ee e e*

*eee e*

1 7 1 43 12 960 192 13440 1 1 1 13 60 192 20160 11520 31 7 41 6720 2304 40320

1 1 7 113 7 <sup>1</sup> 4 12 192 5760 560

6 8 10

*ee e*

8 10

*4.3.2. The direct expansion of the transformation from authalic latitude function to* 

The whole formulas for the transformation from authalic latitude function to isometric

*e e*

10

*e*

41 1 26880 640 167 295680

4 6 8 10

 

sin sin 3 sin 5 sin7 sin9 sin11

24 6 8 10

*ee e e e*

 

*Bb b b b b*

<sup>2</sup>

 

sin 2 sin 4 sin6 sin8 sin10 sin 2 sin 4 sin6 sin8 sin10

 

*BBBBB B*

**authalic latitude function** 

*latitude function* 

function

where

*isometric latitude* 

latitude are as follows:

function are as follows:

2

arcsin(tanh )

 *q*

1

 

3

5

7

9

11

*F R*

*F a* 

 

sin

arcsin(tanh )

*q*

 

#### *4.2.2. The direct expansion of the transformation from authalic latitude function to meridian arc*

The whole formulas for the transformation from authalic latitude function to meridian arc are as follows:

$$\begin{cases} \mathcal{B} = \arcsin(\frac{F}{R^2}) \\ B = \mathcal{B} + c\_2 \sin 2\mathcal{B} + c\_4 \sin 4\mathcal{B} + c\_6 \sin 6\mathcal{B} + c\_8 \sin 8\mathcal{B} + c\_{10} \sin 10\mathcal{B} \\ \mathcal{Y} = B + \alpha\_2 \sin 2B + \alpha\_4 \sin 4B + \alpha\_6 \sin 6B + \alpha\_8 \sin 8B + \alpha\_{10} \sin 10B \\ X = a(1 - e^2)K\_0\mathcal{Y} \end{cases} \tag{61}$$

Expanding *X* as a power series of the eccentricity *e* at *e* 0 by means of Mathematica yields the direct expansion of the transformation from authalic latitude function to meridian arc

$$\begin{cases} \mathcal{G} = \arcsin(\frac{F}{R^2})\\ X = a \left(k\_0 \mathcal{G} + k\_2 \sin 2\mathcal{G} + k\_4 \sin 4\mathcal{G} + k\_6 \sin 6\mathcal{G} + k\_8 \sin 8\mathcal{G} + k\_{10} \sin 10\mathcal{G}\right) \end{cases} \tag{62}$$

where

$$\begin{aligned} k\_{k} &= 1 - \frac{1}{4}e^{2} - \frac{3}{64}e^{4} - \frac{5}{256}e^{6} - \frac{175}{16384}e^{8} - \frac{441}{65536}e^{10} \\ k\_{2} &= -\frac{1}{24}e^{2} - \frac{7}{1440}e^{4} - \frac{61}{107520}e^{6} + \frac{2719}{8294400}e^{8} + \frac{30578453}{61312204800}e^{10} \\ k\_{4} &= -\frac{29}{11520}e^{4} - \frac{1411}{967680}e^{6} - \frac{180269}{232243200}e^{8} - \frac{4110829}{10218700800}e^{10} \\ k\_{6} &= -\frac{1003}{2903040}e^{6} - \frac{341}{921600}e^{8} - \frac{36598301}{1226244096600}e^{10} \\ k\_{8} &= -\frac{40457}{619315200}e^{8} - \frac{3602683}{35035545600}e^{10} \\ k\_{10} &= -\frac{1800439}{122624409600}e^{10} \end{aligned} \tag{63}$$

#### **4.3. The direct expansions of transformations between isometric latitude and authalic latitude function**

*4.3.1. The direct expansion of the transformation from isometric latitude to authalic latitude function* 

The whole formulas for the transformation from isometric latitude to authalic latitude function are as follows:

$$\begin{cases} \rho = \arcsin(\tanh q) \\ B = \rho + b\_2 \sin 2\rho + b\_4 \sin 4\rho + b\_6 \sin 6\rho + b\_8 \sin 8\rho + b\_{10} \sin 10\rho \\ \theta = B + \mathcal{Y}\_2 \sin 2B + \mathcal{Y}\_4 \sin 4B + \mathcal{Y}\_6 \sin 6B + \mathcal{Y}\_8 \sin 8B + \mathcal{Y}\_{10} \sin 10B \\ F = R^2 \sin \mathcal{Y} \end{cases} \tag{64}$$

Expanding *F* as a power series of the eccentricity *e* at *e* 0 by means of Mathematica yields the direct expansion of the transformation from isometric latitude to authalic latitude function

$$\begin{cases} \phi = \arcsin(\tanh \eta) \\ F = a^2 \left( \eta\_1 \sin \phi + \eta\_3 \sin 3\phi + \eta\_5 \sin 5\phi + \eta\_7 \sin 7\phi + \eta\_9 \sin 9\phi + \eta\_{11} \sin 11\phi \right) \end{cases} \tag{65}$$

where

16 Cartography – A Tool for Spatial Analysis

1

 

3

5

7

9

11

*meridian arc* 

as follows:

where

2 4 6 8 10

(60)

(61)

(62)

(63)

 

> 

 

*ee e e e*

2 4 8 10

*ee e e*

1 1 103 1775 48 384 196608 3145728 3 43 17 467 1280 30720 24576 1572864

37 5 563 86016 10752 1572864

*e e*

*4.2.2. The direct expansion of the transformation from authalic latitude function to* 

The whole formulas for the transformation from authalic latitude function to meridian arc are

2 4 6 8 10 246810

 

Expanding *X* as a power series of the eccentricity *e* at *e* 0 by means of Mathematica yields the direct expansion of the transformation from authalic latitude function to meridian arc

*Bc c c c c*

 

sin 2 sin 4 sin6 sin8 sin10 sin 2 sin 4 sin6 sin8 sin10

*BBBBB B*

sin 2 sin 4 sin6 sin8 sin10

24 6 8 10

1 7 61 2719 30578453 24 1440 107520 8294400 61312204800 29 1411 180269 4110829 11520 967680 232243200 10218700800

600 122624409600

*ke e e e e*

*k ee e e*

46 8 10

*e e*

36598301

0 2 4 6 8 10

24 6 8 10

*X ak k k k k k*

1 3 5 175 441 <sup>1</sup> 4 64 256 16384 65536

8 10

8 10

*k ee e e e*

10

40457 3602683 619315200 35035545600

*ke e*

 

5 17 121 137 1407 <sup>1</sup> 16 256 4096 8192 131072

68 10

*ee e*

8

1543 57671680

*e*

2

 

6

1800439 122624409600

*k e*

1003 341 2903040 921

*F R*

2 0

2

*F R*

(1 )

*X a eK*

arcsin( )

*k e*

arcsin( )

0

 

2

4

6

 

8

10

 

 

59 1853 589824 11796480

<sup>10</sup>

10

4 6 8 10

*eee e*

$$\begin{aligned} \eta\_1 &= 1 - \frac{1}{4}e^2 - \frac{1}{12}e^4 - \frac{7}{192}e^6 - \frac{113}{5760}e^8 - \frac{7}{560}e^{10} \\ \eta\_3 &= \frac{1}{12}e^2 - \frac{7}{960}e^6 - \frac{1}{192}e^8 - \frac{43}{13440}e^{10} \\ \eta\_5 &= \frac{1}{60}e^4 + \frac{1}{192}e^6 + \frac{1}{20160}e^8 - \frac{13}{11520}e^{10} \\ \eta\_7 &= \frac{31}{6720}e^6 + \frac{7}{2304}e^8 + \frac{41}{40320}e^{10} \\ \eta\_9 &= \frac{41}{26880}e^8 + \frac{1}{640}e^{10} \\ \eta\_{11} &= \frac{167}{295680}e^{10} \end{aligned} \tag{66}$$

#### *4.3.2. The direct expansion of the transformation from authalic latitude function to isometric latitude*

The whole formulas for the transformation from authalic latitude function to isometric latitude are as follows:

$$\begin{cases} \mathcal{B} = \arcsin(\frac{F}{R^2}) \\ \mathcal{B} = \mathcal{B} + c\_2 \sin 2\mathcal{B} + c\_4 \sin 4\mathcal{B} + c\_6 \sin 6\mathcal{B} + c\_8 \sin 8\mathcal{B} + c\_{10} \sin 10\mathcal{B} \\ \mathcal{q} = \mathcal{B} + \mathcal{J}\_2 \sin 2\mathcal{B} + \mathcal{J}\_4 \sin 4\mathcal{B} + \mathcal{J}\_6 \sin 6\mathcal{B} + \mathcal{J}\_8 \sin 8\mathcal{B} + \mathcal{J}\_{10} \sin 10\mathcal{B} \\ \mathcal{q} = \text{arctanh}(\sin \boldsymbol{\rho}) \end{cases} \tag{67}$$

Mathematical Analysis in Cartography by Means of Computer Algebra System 19

N

x

S

O

y

higher than 7 2 10 km . The accuracy of the indirect expansion of the transformation from authalic latitude function to isometric latitude is higher than 10-2″, while the accuracy of the direct expansion (67) is higher than 10-6″. The accuracies of the direct expansions derived by the author are improved by 2~6 orders of magnitude compared to the indirect ones derived

Gauss projection plays a fundamental role in ellipsoidal geodesy, surveying, map projection and geographical information system (GIS). It has found its wide application in those areas.

**Figure 1.** Gauss Projection, where *x* and *y* are the vertical and horizontal axes after the projection

② the central meridian mapped as a straight line (usually chosen as a vertical axis of *x* )

Traditional expressions of the forward and inverse Gauss projections are real functions in a power series of longitude difference. Though real functions are easy to understand for most people, they make Gauss projection expressions very tedious. Mathematically speaking, there is natural relationship between the conformal mapping and analytical complex functions which automatically meet the differential equation of the conformal mapping, or the "Cauchy-Riemann Equations". Complex functions, a powerful mathematical method, play a very special and key role in the conformal mapping. Bowring (1990) and Klotz (1993) have discussed Gauss projection by complex numbers. But the expressions they derived require iterations, which makes the computation process very fussy. In terms of the direct expansions of transformations between meridian arc and isometric latitude given in section

As shown in Figure 1, Gauss projection has to meet the following three constraints:

equator

**5. The non-iterative expressions of the forward and inverse Gauss** 

by Yang (1989, 2000).

① conformal mapping;

after projection;

**projections by complex numbers** 

respectively, *O* is the origin of the projection coordinates.

O

N

S

③ scale being true along the central meridian.

Expanding *q* as a power series of the eccentricity at *e* 0 by means of Mathematica yields the direct expansion of the transformation form authalic latitude function to isometric latitude

$$\begin{cases} \mathcal{G} = \arcsin(\frac{F}{R^2})\\ q = \operatorname{arctanh}(\sin \mathcal{G}) + l\_1 \sin \mathcal{G} + l\_3 \sin 3\mathcal{G} + l\_5 \sin 5\mathcal{G} + l\_7 \sin 7\mathcal{G} + l\_9 \sin 9\mathcal{G} \end{cases} \tag{68}$$

where

$$\begin{cases} l\_1 = -\frac{1}{3}e^2 - \frac{1}{30}e^4 - \frac{11}{1890}e^6 - \frac{107}{32400}e^8 + \frac{1513}{1663200}e^{10} \\ l\_3 = -\frac{1}{90}e^4 - \frac{61}{11340}e^6 - \frac{2321}{907200}e^8 - \frac{1021}{831600}e^{10} \\ l\_5 = -\frac{1}{756}e^6 - \frac{5}{4032}e^8 - \frac{151}{166320}e^{10} \\ l\_7 = -\frac{71}{302400}e^8 - \frac{41}{123200}e^{10} \\ l\_9 = -\frac{61}{1197504}e^{10} \end{cases} \tag{69}$$

#### **4.4 Accuracies of the direct expansions**

The accuracies of the indirect and direct expansions given by Yang(1989, 2000) derived by the author has been examined choosing the CGCS2000 reference ellipsoid.

The results show that the accuracy of the indirect expansion of the transformation from meridian arc to isometric latitude is higher than 10-3″, while the accuracy of the direct expansion (53) is higher than 10-7″. The accuracy of the indirect expansion of the transformation from isometric latitude to meridian arc is higher than 10-2 m, while the accuracy of the direct expansion (56) is higher than 10-7 m. The accuracy of the indirect expansion of the transformation from meridian arc to authalic latitude function is higher than 0.1 <sup>2</sup> km , while the accuracy of the direct expansion (59) is higher than 5 2 10 km . The accuracy of the indirect expansion of the transformation from authalic latitude function to meridian arc is higher than 10-2 m, while the direct expansion (62) is higher than 10-4 m. The accuracy of the indirect expansion of the transformation from isometric latitude to authalic latitude function is higher than 0.1 <sup>2</sup> km , while the accuracy of the direct expansion (65) is higher than 7 2 10 km . The accuracy of the indirect expansion of the transformation from authalic latitude function to isometric latitude is higher than 10-2″, while the accuracy of the direct expansion (67) is higher than 10-6″. The accuracies of the direct expansions derived by the author are improved by 2~6 orders of magnitude compared to the indirect ones derived by Yang (1989, 2000).
