Mathematical Analysis in Cartography by Means of Computer Algebra System

1. Introduction

Theory of map projections is a branch of cartography studying the ways of projecting the curved surface of the earth and other heavenly bodies into the plane, and it is often called mathematical cartography. There are many fussy symbolic problems to be dealt with in map projections, such as power series expansions of elliptical functions, high order differential of transcendental functions, elliptical integrals and the operation of complex numbers. Many famous cartographers such as Adams (1921), Snyder (1987), Yang (1989, 2000) have made great efforts to solve these problems. Due to historical condition limitation, there were no advanced computer algebra systems at that time, so they had to dispose of these problems by hand, which had often required a paper and a pen. Some derivations and computations were however long and labor intensive such that one gave up midway. Briefly reviewing the existing methods, one will find that these problems were not perfectly and ideally solved yet. Formulas derived by hand often have quite complex and prolix forms, and their orders could not be very high. The most serious problem is that some higher terms of the formulas are erroneous because of the adopted approximate disposal.

With the advent of computers, the paper and pen approach is slowly being replaced by software developed to undertake symbolic derivations tasks. Specially, where symbolic rather than numerical solutions are desired, this software normally comes in handy. Software which performs symbolic computations is called computer algebra system. Nowadays, computer algebra systems like Maple, Mathcad, and Mathematica are widely used by scientists and engineers in different fields (Awang, 2005; Bian, 2006). By means of computer algebra system Mathematica, we have already solved many complicated mathematical problems in special fields of cartography in the past few years. Our research 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 basis of cartography to a certain extent.

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 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){\displaystyle {\int }_0^B{(1-e^2{\sin }^{2}B)}^{-3/2}}dB$$E1

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

where

$$\{\begin{array}{l}{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{array}$$E3

The rectifying latitude $\psi$is defined as

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

Inserting (2) into (4) yields

$$\psi =B+{\alpha }_2\sin 2B+{\alpha }_4\sin 4B+{\alpha }_6\sin 6B+{\alpha }_8\sin 8B+{\alpha }_{10}\sin 10B$$E5

where

$$\{\begin{array}{c}{\alpha }_2=K_2/K_0\\ {\alpha }_4=K_4/K_0\\ {\alpha }_6=K_6/K_0\\ \begin{array}{l}{\alpha }_8=K_8/K_0\\ {\alpha }_{10}=K_{10}/K_0\end{array}\end{array}$$E6

Yang (1989, 2000) gave an expansion similar to (5) but expanded $\psi$up to$$\sin 8B$$. 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{array}{c}{\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}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}\\ {\alpha }_6=-\frac{35}{3072}{e}^6-\frac{35}{2048}{e}^8-\frac{4935}{262144}{e}^{10}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}\\ \begin{array}{l}{\alpha }_8=\frac{315}{131072}{e}^8+\frac{315}{65536}{e}^{10}\begin{array}{cccc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}& \end{array}& & & \end{array}\\ {\alpha }_{10}=-\frac{693}{1310720}{e}^{10}\end{array}\end{array}$$E7

2.2. The forward expansion of the conformal latitude

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

$$\begin{array}{l}q={\displaystyle {\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]\\ =\arctan \text{h(}\sin B\text{)-}e\cdot \arctan \text{h}(e\sin B)\end{array}$$E8

For the sphere, putting$e=0$and rewriting the geodetic latitude as the conformal one$\phi$, (8) becomes

$$q=\ln \left[\tan (\frac{\pi }{4}+\frac{\phi }{2})\right]=\arctan \text{h(sin}\phi \text{)}$$E9

Comparing (9) with (8) leads to

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

Therefore, it holds

$$\phi =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}$$E11

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

$$\begin{array}{l}\phi (B,e)=\phi (B,0)+{\frac{\partial \phi }{\partial e}|}_{e=0}e+\frac{1}{2!}{\frac{{\partial }^2\phi }{\partial e^2}|}_{e=0}{e}^2+\frac{1}{3!}{\frac{{\partial }^3\phi }{\partial e^3}|}_{e=0}{e}^3+\cdots \frac{1}{9!}{\frac{{\partial }^9\phi }{\partial e^9}|}_{e=0}\\ e^9+\frac{1}{10!}{\frac{{\partial }^{10}\phi }{\partial e^{10}}|}_{e=0}{e}^{10}+\cdots \end{array}$$E12

as the conventional usage in mathematical cartography.

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

$$\phi =B+{\beta }_2\sin 2B+{\beta }_4\sin 4B+{\beta }_6\sin 6B+{\beta }_8\sin 8B$$E13

Due to that (11) is a very complicated transcendental function, the coefficients${\beta }_2$,${\beta }_4$,${\beta }_6$,${\beta }_8$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

$$\phi =B+{\beta }_2\sin 2B+{\beta }_4\sin 4B+{\beta }_6\sin 6B+{\beta }_8\sin 8B+{\beta }_{10}\sin 10B$$E14

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

Coefficients derived by Yang(1989, 2000)

Coefficients derived by the author

$$\{\begin{array}{c}{\beta }_2=-\frac{1}{2}{e}^2-\frac{5}{24}{e}^4-\frac{3}{32}{e}^6-\frac{1399}{53760}{e}^8\\ {\beta }_4=\frac{5}{48}{e}^4+\frac{7}{80}{e}^6+\frac{689}{17920}{e}^8\begin{array}{ccc}& & \end{array}\\ {\beta }_6=-\frac{13}{480}{e}^6-\frac{1363}{53760}{e}^8\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}& \end{array}\\ {\beta }_8=\frac{677}{17520}{e}^8\begin{array}{cccc}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}& & & \end{array}\end{array}$$

$$\{\begin{array}{c}{\beta }_2=-\frac{1}{2}{e}^2-\frac{5}{24}{e}^4-\frac{3}{32}{e}^6-\frac{281}{5760}{e}^8-\frac{7}{240}{e}^{10}\\ {\beta }_4=\frac{5}{48}{e}^4+\frac{7}{80}{e}^6+\frac{697}{11520}{e}^8+\frac{93}{2240}{e}^{10}\begin{array}{cc}& \end{array}\\ {\beta }_6=-\frac{13}{480}{e}^6-\frac{461}{13440}{e}^8-\frac{1693}{53760}{e}^{10}\begin{array}{ccc}& & \end{array}\\ \begin{array}{l}{\beta }_8=\frac{1237}{161280}{e}^8+\frac{131}{10080}{e}^{10}\begin{array}{cccc}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}& & & \end{array}\\ {\beta }_{10}=-\frac{367}{161280}{e}^{10}\end{array}\end{array}$$

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){\displaystyle {\int }_0^B\frac{\cos B}{{(1-e^2{\sin }^{2}B)}^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]$$E15

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

Suppose that there is an imaginary sphere with a radius

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

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

$$R^2\sin \vartheta =a^2(1-e^2)A\sin \vartheta =F$$E18

Therefore, it yields

$$\vartheta =\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]$$E19

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

$$\vartheta =B+{\gamma }_2\sin 2B+{\gamma }_4\sin 4B+{\gamma }_6\sin 6B+{\gamma }_8\sin 8B$$E20

(19) is a complicated transcendental function. It is almost impossible to derive its eighth-order derivate by hand. There are some mistakes in the high order terms of coefficients${\gamma }_2$,${\gamma }_4$,${\gamma }_6$,${\gamma }_8$.The new derived forward expansion expanded to tenth-order terms of $e$by means of Mathematica reads

$$\vartheta =B+{\gamma }_2\sin 2B+{\gamma }_4\sin 4B+{\gamma }_6\sin 6B+{\gamma }_8\sin 8B+{\gamma }_{10}\sin 10B$$E21

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

Coefficients derived by Yang(1989, 2000)

Coefficients derived by the author

$$\{\begin{array}{c}{\gamma }_2=-\frac{1}{3}{e}^2-\frac{31}{180}{e}^4-\frac{59}{560}{e}^6-\frac{126853}{518400}{e}^8\begin{array}{cc}& \end{array}\\ {\gamma }_4=\frac{17}{360}{e}^4+\frac{61}{1260}{e}^6+\frac{3622447}{94089600}{e}^8\begin{array}{ccc}& & \end{array}\\ {\gamma }_6=-\frac{383}{43560}{e}^6-\frac{6688039}{658627200}{e}^8\begin{array}{cccc}\begin{array}{cc}& \end{array}& & & \end{array}\\ {\gamma }_8=-\frac{27787}{23522400}{e}^8\begin{array}{cccc}\begin{array}{cc}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}& \end{array}& & & \end{array}\end{array}$$

$$\{\begin{array}{c}{\gamma }_2=-\frac{1}{3}{e}^2-\frac{31}{180}{e}^4-\frac{59}{560}{e}^6-\frac{42811}{604800}{e}^8-\frac{605399}{11975040}{e}^{10}\begin{array}{c}\end{array}\\ {\gamma }_4=\frac{17}{360}{e}^4+\frac{61}{1260}{e}^6+\frac{76969}{1814400}{e}^8+\frac{215431}{5987520}{e}^{10}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}\\ {\gamma }_6=-\frac{383}{45360}{e}^6-\frac{3347}{259200}{e}^8-\frac{1751791}{119750400}{e}^{10}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}\\ \begin{array}{l}{\gamma }_8=\frac{6007}{3628800}{e}^8+\frac{201293}{59875200}{e}^{10}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cc}& \end{array}& & & \end{array}& & & \end{array}\\ {\gamma }_{10}=-\frac{5839}{17107200}{e}^{10}\end{array}\end{array}$$

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

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=6378137\text{m}$,$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:

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

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

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

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)}{{(1-e^2{\sin }^{2}B)}^{3/2}}$$E22

From (4) and (2), one knows

$$X=a(1-e^2)K_0\psi$$E23

Inserting (23) into (22) yields

$$\frac{dB}{K_{0}d\psi }={(1-e^2{\sin }^{2}B)}^{3/2}$$E24

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

$$t=\sin \psi$$E25

therefore

$$\frac{d\psi }{dt}=\frac{1}{\cos \psi }$$E26

and then denote

$$f(t)=\frac{dB}{K_{0}d\psi }={(1-e^2{\sin }^{2}B)}^{3/2}$$E27

Making use of the chain rule of implicit differentiation

$${f^{\prime }}_t=\frac{df}{dB}\frac{dB}{d\psi }\frac{d\psi }{dt}+\frac{df}{d\psi }\frac{d\psi }{dt},\begin{array}{c}\end{array}{f^{″}}_t=\frac{d{f}^{\prime }}{dB}\frac{dB}{d\psi }\frac{d\psi }{dt}+\frac{d{f}^{\prime }}{d\psi }\frac{d\psi }{dt},\begin{array}{c}\end{array}\cdots \cdots$$E28

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

$$f(t)=f(0)+{f^{\prime }}_t(0)t+\frac{1}{2!}{f^{″}}_t(0)t^2+\frac{1}{3!}{f^{‴}}_t(0)t^3+\cdots +\frac{1}{10!}{f}_t{}^{(10)}(0)t^{10}+\cdots$$E29

Omitting the detailed procedure, one arrives at

$$\frac{dB}{K_{0}d\psi }=1+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$$E30

where

$$\{\begin{array}{c}{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}\begin{array}{cc}& \end{array}\\ A_4=\frac{21}{8}{e}^4+\frac{621}{64}{e}^6+\frac{11987}{512}{e}^8+\frac{757215}{16384}{e}^{10}\begin{array}{ccc}& & \begin{array}{cc}& \end{array}\end{array}\\ A_6=-\frac{151}{32}{e}^6-\frac{775}{32}{e}^8-\frac{621445}{8192}{e}^{10}\begin{array}{cc}& \end{array}\begin{array}{cccc}\begin{array}{cc}& \end{array}& & & \end{array}\\ \begin{array}{l}{A}_8=\frac{1097}{128}{e}^8+\frac{57607}{1024}{e}^{10}\begin{array}{cc}& \end{array}\begin{array}{cccc}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}& & & \begin{array}{cc}& \end{array}\end{array}\\ A_{10}=-\frac{8011}{512}{e}^{10}\end{array}\end{array}$$E31

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

where

$$\{\begin{array}{c}{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}\begin{array}{c}\end{array}\\ a_4=\frac{21}{256}{e}^4+\frac{21}{256}{e}^6+\frac{533}{8192}{e}^8+\frac{197}{4096}{e}^{10}\begin{array}{c}\end{array}\begin{array}{ccc}& & \begin{array}{c}\end{array}\end{array}\\ a_6=\frac{151}{6144}{e}^6+\frac{151}{4096}{e}^8+\frac{5019}{131072}{e}^{10}\begin{array}{cccc}\begin{array}{cc}& \end{array}& & & \end{array}\\ \begin{array}{l}{a}_8=\frac{1097}{131072}{e}^8+\frac{1097}{65536}{e}^{10}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}\begin{array}{cc}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}& \end{array}\\ a_{10}=\frac{8011}{2621440}{e}^{10}\end{array}\end{array}$$E33

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

Therefore, it holds

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

For the same reason, we introduce the following new variable

$$t=\sin \phi$$E36

and then denote

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

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

$$\frac{dB}{d\phi }=\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$$E38

where

$$\{\begin{array}{c}{B}_2=-2e^2-\frac{11}{2}{e}^4-\frac{21}{2}{e}^6-17e^8-25e^{10}\begin{array}{cc}& \end{array}\\ B_4=\frac{14}{3}{e}^4+\frac{62}{3}{e}^6+\frac{1369}{24}{e}^8+\frac{3005}{24}{e}^{10}\begin{array}{cc}& \end{array}\\ B_6=-\frac{56}{5}{e}^6-\frac{614}{9}{e}^8-\frac{4909}{20}{e}^{10}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}& \end{array}\\ \begin{array}{l}{B}_8=\frac{8558}{315}{e}^8+\frac{7367}{35}{e}^{10}\begin{array}{c}\end{array}\begin{array}{ccc}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}& & \end{array}\\ B_{10}=-\frac{4174}{63}{e}^{10}\end{array}\end{array}$$E39

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

where

$$\{\begin{array}{c}{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}\begin{array}{c}\end{array}\\ b_6=\frac{7}{120}{e}^6+\frac{81}{1120}{e}^8+\frac{3029}{53760}{e}^{10}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}\\ \begin{array}{l}{b}_8=\frac{4279}{161280}{e}^8+\frac{883}{20160}{e}^{10}\begin{array}{cccc}\begin{array}{cc}& \end{array}& & & \end{array}\\ b_{10}=\frac{2087}{161280}{e}^{10}\end{array}\end{array}$$E41

3.1.3. The inverse expansion of the authalic latitude

Inserting (18) into (15) yields

$$A\sin \vartheta ={\displaystyle {\int }_0^B\frac{\cos B}{{(1-e^2{\sin }^{2}B)}^2}dB}$$E42

Differentiating the both sides of (42) yields

$$\frac{dB}{d\vartheta }=\frac{A{(1-e^2{\sin }^{2}B)}^2\cos \vartheta }{\cos B}$$E43

For the same reason, we introduce the folllowing new variable

$$t=\sin \vartheta$$E44

and then denote

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

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

$$f(t)=\frac{dB}{d\vartheta }=A+C_2{\sin }^2\vartheta +C_4{\sin }^4\vartheta +C_6{\sin }^6\vartheta +C_8{\sin }^8\vartheta +C_{10}{\sin }^{10}\vartheta$$E46

where

$$\{\begin{array}{c}{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}\begin{array}{c}\end{array}\\ C_4=\frac{92}{45}{e}^4+\frac{6574}{945}{e}^6+\frac{223469}{14175}{e}^8+\frac{2768558}{93555}{e}^{10}\begin{array}{cc}& \end{array}\\ C_6=-\frac{3044}{945}{e}^6-\frac{28901}{1890}{e}^8-\frac{21018157}{467775}{e}^{10}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}\\ \begin{array}{l}{C}_8=\frac{24236}{4725}{e}^8+\frac{2086784}{66825}{e}^{10}\begin{array}{ccc}\begin{array}{ccc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}& \end{array}& & \end{array}& & \end{array}\\ C_{10}=-\frac{768272}{93555}{e}^{10}\end{array}\end{array}$$E47

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

$$B=\vartheta +c_2\sin 2\vartheta +c_4\sin 4\vartheta +c_6\sin 6\vartheta +c_8\sin 8\vartheta +c_{10}\sin 10\vartheta$$E48

where

$$\{\begin{array}{c}{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}\begin{array}{ccc}& & \end{array}\\ c_6=\frac{761}{45360}{e}^6+\frac{47561}{1814400}{e}^8+\frac{434501}{14968800}{e}^{10}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}\\ \begin{array}{l}{c}_8=\frac{6059}{1209600}{e}^8+\frac{625511}{59875200}{e}^{10}\begin{array}{c}\end{array}\begin{array}{ccc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}& \end{array}& \end{array}& & \end{array}\\ c_{10}=\frac{48017}{29937600}{e}^{10}\end{array}\end{array}$$E49

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:

1. To suppose the inverse expansions are expressed in a series of the sines of the multiple arcs with coefficients to be determined;

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

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

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=x+f(x)$$E50

with $\left|f(x)\right|\ll \left|x\right|$and$y\approx x$. The Lagrange’s theorem states that in a suitable domain the solution of (50) is

$$x=y+{\displaystyle \sum _{n=1}^{\infty }\frac{{(-1)}^n}{n!}}\frac{d^{n-1}}{d{y}^{n-1}}{[f(y)]}^n$$E51

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:

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

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 eccentricity.

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

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$\psi$. Using the inverse expansion of the rectifying latitude (32) and the forward expansion of the conformal latitude (14), one obtains the conformal latitude$\phi$. Inserting it into (9) yields the isometric latitude$q$. The whole formulas are as follows:

$$\{\begin{array}{l}\psi =\frac{X}{a(1-e^2)K_0}\\ 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 \\ \phi =B+{\beta }_2\sin 2B+{\beta }_4\sin 4B+{\beta }_6\sin 6B+{\beta }_8\sin 8B+{\beta }_{10}\sin 10B\\ q=\mathrm{arctanh}(\sin \phi )\end{array}$$E52

Since the coefficients$a_{2i}$, ${\beta }_{2i}$($i=1,2,\cdots 5$) 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{array}{l}\psi =\frac{X}{a(1-e^2)K_0}\\ q=\mathrm{arctanh}(\sin \psi )+{\xi }_1\sin \psi +{\xi }_3\sin 3\psi +{\xi }_5\sin 5\psi +{\xi }_7\sin 7\psi +{\xi }_9\sin 9\psi \end{array}$$E53

where

$$\{\begin{array}{l}{\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}\\ {\xi }_3=-\frac{1}{96}{e}^4-\frac{13}{3072}{e}^6-\frac{13}{8192}{e}^8-\frac{1057}{1966080}{e}^{10}\\ {\xi }_5=-\frac{11}{7680}{e}^6-\frac{29}{24576}{e}^8-\frac{2897}{3932160}{e}^{10}\\ {\xi }_7=-\frac{25}{86016}{e}^8-\frac{727}{1966080}{e}^{10}\\ {\xi }_9=-\frac{53}{737280}{e}^{10}\end{array}$$E54

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{array}{l}\phi =\arcsin (\tanh q)\\ 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 \\ \psi =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\psi \end{array}$$E55

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

$$\{\begin{array}{l}\phi =\arcsin (\tanh q)\\ X=a\left(j_0\phi +j_2\sin 2\phi +j_4\sin 4\phi +j_6\sin 6\phi +j_8\sin 8\phi +j_{10}\sin 10\phi \right)\end{array}$$E56

where

$$\{\begin{array}{l}{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}{27525120}{e}^{10}\\ j_8=\frac{49561}{41287680}{e}^8+\frac{175087}{165150720}{e}^{10}\\ j_{10}=\frac{34729}{82575360}{e}^{10}\end{array}$$E57

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{array}{l}\psi =\frac{X}{a(1-e^2)K_0}\\ 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 \\ \vartheta =B+{\gamma }_2\sin 2B+{\gamma }_4\sin 4B+{\gamma }_6\sin 6B+{\gamma }_8\sin 8B+{\gamma }_{10}\sin 10B\\ F=R^2\sin \vartheta \end{array}$$E58

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{array}{l}\psi =\frac{X}{a(1-e^2)K_0}\\ F=a^2\left({\epsilon }_1\sin \psi +{\epsilon }_3\sin 3\psi +{\epsilon }_5\sin 5\psi +{\epsilon }_7\sin 7\psi +{\epsilon }_9\sin 9\psi +{\epsilon }_{11}\sin 11\psi \right)\end{array}$$E59

where

$$\{\begin{array}{l}{\epsilon }_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}\\ {\epsilon }_3=\frac{1}{48}{e}^2+\frac{1}{384}{e}^4-\frac{103}{196608}{e}^8-\frac{1775}{3145728}{e}^{10}\\ {\epsilon }_5=\frac{3}{1280}{e}^4+\frac{43}{30720}{e}^6+\frac{17}{24576}{e}^8+\frac{467}{1572864}{e}^{10}\\ {\epsilon }_7=\frac{37}{86016}{e}^6+\frac{5}{10752}{e}^8+\frac{563}{1572864}{e}^{10}\\ {\epsilon }_9=\frac{59}{589824}{e}^8+\frac{1853}{11796480}{e}^{10}\\ {\epsilon }_{11}=\frac{1543}{57671680}{e}^{10}\end{array}$$E60

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{array}{l}\vartheta =\arcsin (\frac{F}{R^2})\\ B=\vartheta +c_2\sin 2\vartheta +c_4\sin 4\vartheta +c_6\sin 6\vartheta +c_8\sin 8\vartheta +c_{10}\sin 10\vartheta \\ \psi =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\psi \end{array}$$E61

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{array}{l}\vartheta =\arcsin (\frac{F}{R^2})\\ X=a\left(k_0\vartheta +k_2\sin 2\vartheta +k_4\sin 4\vartheta +k_6\sin 6\vartheta +k_8\sin 8\vartheta +k_{10}\sin 10\vartheta \right)\end{array}$$E62

where

$$\{\begin{array}{l}{k}_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}\\ 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}{122624409600}{e}^{10}\\ k_8=-\frac{40457}{619315200}{e}^8-\frac{3602683}{35035545600}{e}^{10}\\ k_{10}=-\frac{1800439}{122624409600}{e}^{10}\end{array}$$E63

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{array}{l}\phi =\arcsin (\tanh q)\\ 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 \\ \vartheta =B+{\gamma }_2\sin 2B+{\gamma }_4\sin 4B+{\gamma }_6\sin 6B+{\gamma }_8\sin 8B+{\gamma }_{10}\sin 10B\\ F=R^2\sin \vartheta \end{array}$$E64

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{array}{l}\varphi =\arcsin (\tanh q)\\ F=a^2\left({\eta }_1\sin \varphi +{\eta }_3\sin 3\varphi +{\eta }_5\sin 5\varphi +{\eta }_7\sin 7\varphi +{\eta }_9\sin 9\varphi +{\eta }_{11}\sin 11\varphi \right)\end{array}$$E65

where

$$\{\begin{array}{l}{\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{array}$$E66

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{array}{l}\vartheta =\arcsin (\frac{F}{R^2})\\ B=\vartheta +c_2\sin 2\vartheta +c_4\sin 4\vartheta +c_6\sin 6\vartheta +c_8\sin 8\vartheta +c_{10}\sin 10\vartheta \\ \phi =B+{\beta }_2\sin 2B+{\beta }_4\sin 4B+{\beta }_{\text{6}}\sin 6B+{\beta }_{\text{8}}\sin 8B+{\beta }_{\text{10}}\sin 10B\\ q=\mathrm{arctanh}(\sin \phi )\end{array}$$E67

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{array}{l}\vartheta =\arcsin (\frac{F}{R^2})\\ q=arctanh(\sin \vartheta )+l_1\sin \vartheta +l_3\sin 3\vartheta +l_5\sin 5\vartheta +l_7\sin 7\vartheta +l_9\sin 9\vartheta \end{array}$$E68

where

$$\{\begin{array}{l}{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{array}$$E69

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${\text{km}}^{\text{2}}$, while the accuracy of the direct expansion (59) is higher than${10}^{-5}{\text{km}}^{\text{2}}$. 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${\text{km}}^{\text{2}}$, while the accuracy of the direct expansion (65) is higher than${10}^{-7}{\text{km}}^{\text{2}}$. 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).

5. The non-iterative expressions of the forward and inverse Gauss projections by complex numbers

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 respectively,$O$is the origin of the projection coordinates.

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

① conformal mapping;

② the central meridian mapped as a straight line (usually chosen as a vertical axis of$x$) after projection;

③ scale being true along the central meridian.

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 Ⅳ, the non-iterative expressions of the forward and inverse Gauss projections by complex numbers are derived.

5.1. The non-iterative expressions of the forward Gauss projection by complex numbers

Let $w$be complex numbers consisting of isometric latitude $q$and longitude difference $l$before projection, $z$be complex numbers consisting of corresponding coordinates$x$,$y$after projection.

$$\{\begin{array}{l}w=q+il\\ z=x+iy\end{array}$$E70

where

$$i=\sqrt{-1}$$E71

.

In terms of complex functions theory, analytical functions meet conformal mapping naturally. Therefore, to meet the conformal mapping constraint, the forward Gauss projection should be in the following form

$$z=x+iy=f(w)=f(q+il)$$E72

where $f$is an arbitrary analytical function in the complex numbers domain. According to the second constraint, when$l=0$, imaginary part disappears and only real part exists, (71) becomes

$$x=f(q)$$E73

(72) shows that the central meridian is a straight line after the projection when$l=0$.

Finally, from the third constraint, “scale is true along the central meridian”, one knows that $x$in (72) should be nothing else but the meridian arc$X$, and (72) is essentially consist with the direct expansion of the transformation from isometric latitude to meridian arc (56). Substituting $X$in (56) with $x$gives the explicit form of (72)

$$\{\begin{array}{l}\phi =\arcsin (\tanh q)\\ x=a\left(j_0\phi +j_2\sin 2\phi +j_4\sin 4\phi +j_6\sin 6\phi +j_8\sin 8\phi +j_{10}\sin 10\phi \right)\end{array}$$E74

(73) defines the functional relationship between meridian arc and isometric latitude. If one extends the definition of $q$in a real number variable to a complex numbers variable, or substitutes $q$with$w=q+il$, the original real number conformal latitude $\phi$will be automatically extended as a complex numbers variable. We denote the corresponding complex numbers latitude as$\Phi$, and insert it into (73). Rewriting a real variable $x$at the left-hand of the second equation in (73) as a complex numbers variable$z=x+iy$, one arrives at

$$\{\begin{array}{l}\Phi =\arcsin (\tanh w)\\ z=x+iy=a\left(j_0\phi +j_2\sin 2\phi +j_4\sin 4\phi +j_6\sin 6\phi +j_8\sin 8\phi +j_{10}\sin 10\phi \right)\end{array}$$E75

(74) is the solution of the forward Gauss projection by complex numbers. Its correctness can be explained as follows:

The two equations in (74) are all elementary complex functions. Because elementary functions in their basic interval are all analytical ones in the complex numbers domain, the mapping defined by (74) form $w=q+il$to $z=x+iy$meets the conformal mapping constraint. When$l=0$, the imaginary part disappears and (74) restores to (73). Therefore, (74) meets the second and third constraints of Gauss projection when$l=0$. Hence, it is clear that (74) is the solution of the forward Gauss projection indeed.

5.2. The non-iterative expressions of the inverse Gauss projection by complex numbers

In principle, the inverse Gauss projection can be iteratively solved in terms of the forward Gauss projection (74). In order to eliminate the iteration, one more practical approach is proposed based on the direct expansion of the transformation from meridian arc to isometric latitude (53).

In order to meet the conformal mapping constraint, the inverse Gauss projection should be in the following form

$$w=q+il=f^{-1}(z)=f^{-1}(x+iy)$$E76

where $f^{-1}$is the inverse function of$f$. According to the second constraint, when$l=0$, imaginary part disappears and only real part exists, (75) becomes

$$q=f^{-1}(x)$$E77

Finally, from the third constraint, one knows that $x$in (76) should be the meridian arc$X$, and (76) is essentially consist with the direct expansion of the transformation from meridian arc to isometric latitude as (53) shows. Substituting $X$in (53) with $x$gives the explicit form of (76)

$$\{\begin{array}{l}\psi =\frac{x}{a(1-e^2)K_0}\\ q=arctanh(\sin \psi )+{\xi }_1\sin \psi +{\xi }_3\sin 3\psi +{\xi }_5\sin 5\psi +{\xi }_7\sin 7\psi +{\xi }_9\sin 9\psi \end{array}$$E78

If one extends the definition of $x$in a real number variable to a complex numbers variable, or substitutes $x$with$z=x+iy$, the original real number rectifying latitude $\psi$will be automatically extended as a complex numbers variable. We denote the corresponding complex number latitude as$\Psi$, and insert it into (77). Rewriting a real variable $q$at the left-hand of the second equation in (77) as a complex numbers variable$w=q+il$, one arrives at

$$\{\begin{array}{l}\Psi =\frac{x+iy}{a(1-e^2)K_0}\\ w=q+il=arctanh(\sin \Psi )+{\xi }_1\sin \Psi +{\xi }_3\sin 3\Psi +{\xi }_5\sin 5\Psi +{\xi }_7\sin 7\Psi +{\xi }_9\sin 9\Psi \end{array}$$E79

Therefore, the isometric latitude $q$and longitude $l$is known. Inserting $q$into (78) yields the conformal latitude

$$\phi =\arcsin (\tanh q)$$E80

Then one can compute the geodetic latitude through the inverse expansion of the conformal latitude (40).

(77) is the solution of the inverse Gauss projection by complex numbers. Its correctness can be explained as follows:

The two equations in (78) are all elementary complex functions, so the mapping defined by (78) form $z=x+iy$to $w=q+il$meets the conformal mapping constraint. When$l=0$, the imaginary part disappears and (78) restores to (77). Therefore, (78) meets the second and third constraints of Gauss projection when$l=0$. Hence, it is clear that (78) is the solution of the inverse Gauss projection indeed.

6. Conclusions

Some typical mathematical problems in map projections are solved by means of computer algebra system which has powerful function of symbolical operation. The main contents and research results presented in this chapter are as follows:

1. Forward expansions of rectifying, conformal and authalic latitudes are derived, and some mistakes once made in the high orders of traditional forward formulas are pointed out and corrected. Inverse expansions of rectifying, conformal and authalic latitudes are derived using power series expansion, Hermite interpolation and Language’s theorem methods respectively. These expansions are expressed in a series of the sines of the multiple arcs. Their coefficients are expressed in a power series of the first eccentricity of the reference ellipsoid and extended up to its tenth-order terms. The accuracies of these expansions are analyzed through numerical examples. The results show that the accuracies of these expansions derived by means of computer algebra system are improved by 2~4 orders of magnitude compared to the formulas derived by hand.

2. Direct expansions of transformations between meridian arc, isometric latitude and authalic latitude function are derived. Their coefficients are expressed in a power series of the first eccentricity of the reference ellipsoid, and extended up to its tenth-order terms. Numerical examples show that the accuracies of these direct expansions are improved by 2~6 orders of magnitude compared to the traditional indirect formulas.

3. Gauss projection is discussed in terms of complex numbers theory. The non-iterative expressions of the forward and inverse Gauss projections by complex numbers are derived based on the direct expansions of transformations between meridian arc and isometric latitude, which enriches the theory of conformal projection. In USA, Universal Transverse Mercator Projection (or UTM) is usually implemented. Mathematically speaking, there is no essential difference between UTM and Gauss projections. The only difference is that the scale factor of UTM is 0.9996 rather than 1. With a slight modification, the non-iterative expressions of the forward and inverse Gauss projections can be extended to UTM projection accordingly.

Acknowledgement

This work was financially supported by 973 Program (2012CB719902), National Natural Science Foundation of China (No. 41071295 and 40904018), and Key Laboratory of Surveying and Mapping Technology on Island and Reef, State Bureau of Surveying, Mapping and Geoinformation, China (No.2010B04).