D geodesy

D.1 Convention on coordinate systems

D.1.1 Resting coordinates

Here you choose a base $\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3$, which is at the center of the earth and does not rotate with the earth. $\mathbf{e}_1$ and $\mathbf{e}_2$ lie in the equatorial plane, $\mathbf{e}_3$ points towards the North Pole. The wind field $\mathbf{v}$ is given as

\[ \begin{align} \mathbf{v} = u_1\mathbf{e}_1 + u_2\mathbf{e}_2 + u_3\mathbf{e}_3 \end{align} \]

written. The wind field is not the same as the particle movement measured in this system, as the rotation of the Earth is also superimposed.

D.1.2 Global coordinates

Here you choose a base $\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z$, which is at the center of the earth just like the base of the resting coordinates, but also rotates. $\mathbf{e}_x$ points in the direction of the intersection between the prime meridian and the equator, $\mathbf{e}_y$ points in the direction of the intersection between ninetieth longitude and the equator and $\mathbf{e}_z$ points in the direction of the North Pole. The wind field is written as

\[ \begin{align} \mathbf{v} = u_x\mathbf{e}_x + u_y\mathbf{e}_y + u_z\mathbf{e}_z. \end{align} \]

D.1.3 Geographic coordinates

Alternatively, you can also use geographical coordinates $\left(r, \varphi, \lambda\right)$, where $r\geq 0$ is the distance from the center of the Earth, $-\pi/2\leq \varphi\leq\pi/2$ is the angle that a point makes with the equatorial plane (the geographical latitude assuming a spherical Earth) and $0\leq\lambda<2\pi$ is the angle a point with the xz plane includes (the geographical longitude). This can be transformed into ordinary spherical coordinates $\left(r, \theta, \phi\right)$, with $\theta + \varphi = \pi/2$

\[ \begin{align} \sin\left(\varphi\right) &= \sin\left(\pi/2 - \theta\right) = \cos\left(\theta\right), \tag{D.3}\label{eq:kugel_zu_kugel_met_1}\\ \cos\left(\varphi\right) &= \cos\left(\pi/2 - \theta\right) = \sin\left(\theta\right), \tag{D.4}\label{eq:kugel_zu_kugel_met_2}\\ \tan\left(\theta\right) &= \frac{\sin\left(\theta\right)}{\cos\left(\theta\right)} = \frac{\cos\left(\varphi\right)}{\sin\left(\varphi\right)} = \frac{1}{\tan\left(\varphi\right)}, \tag{D.5}\label{eq:kugel_zu_kugel_met_3}\\ v_\theta &= -v_y, \tag{D.6}\label{eq:kugel_zu_kugel_met_4}\\ \frac{\partial}{\partial\theta} &= -\frac{\partial}{\partial\varphi}\tag{D.7}\label{eq:kugel_zu_kugel_met_5}. \end{align} \]

Locally, the use of the spherical coordinate system is very confusing if you want to note vectors. Therefore, a Cartesian, right-handed coordinate system is used locally, which stands on the sphere and whose x-axis points to the east, y-axis to north and z-axis to the top. The respective basis vectors are called $\mathbf{i}, \mathbf{j}, \mathbf{k}$. The wind field $\mathbf{v}$ is written as

\[ \begin{align} \mathbf{v} = u\mathbf{i} + v\mathbf{j} + w\mathbf{k}. \end{align} \]

D.2 Spherical geometry

D.2.1 Geodesics

Here we want to determine the shortest connecting line between two points $\left(\varphi_i, \lambda_i\right)$ on a sphere with $i = 1, 2$ along the surface of the sphere. The length $L$ of any connecting line is calculated as:

\[ \begin{align} L = \int ds = \int_{0}^{1}\left|\frac{d\mathbf{r}}{d\tau}\right|d\tau, \end{align} \]

where $\tau\in\left[0, 1\right]$ is a parameter that allows one to move along the connecting line and $\mathbf{r}\left(\tau\right)$ sets the trajectory. In particular, $\mathbf{r}\left(0\right) = \mathbf{r}\left(\varphi_1, \lambda_1\right)$ and $\mathbf{r}\left(1\right) = \mathbf{r}\left(\varphi_2, \lambda_2\right)$.

First start from the sphere, here to determine the functions $\varphi\left(\tau\right), \lambda\left(\tau\right)$ you can first start from the unit sphere and ignore the radius. Then apply

\[ \begin{align} \mathbf{r}\left(0\right) &= \left(\begin{array}{c} \cos\left(\varphi_1\right)\cos\left(\lambda_1\right)\\ \cos\left(\varphi_1\right)\sin\left(\lambda_1\right)\\ \sin\left(\varphi_1\right) \end{array}\right),\\ \mathbf{r}\left(1\right) &= \left(\begin{array}{c} \cos\left(\varphi_2\right)\cos\left(\lambda_2\right)\\ \cos\left(\varphi_2\right)\sin\left(\lambda_2\right)\\ \sin\left(\varphi_2\right) \end{array}\right). \end{align} \]

Define a function

\[ \begin{align} \mathbf{r}'\left(\tau'\right)&\coloneqq\tau'\mathbf{r}\left(1\right) + \left(1 - \tau'\right)\mathbf{r}\left(0\right). \end{align} \]

Then applies

\[ \begin{align} \varphi\left(\tau'\right) = \arcsin\left(\tau'\sin\left(\varphi_2\right) + \left(1 - \tau'\right)\sin\left(\varphi_1\right)\right). \end{align} \]

Furthermore applies

\[ \begin{align} \lambda\left(\tau'\right) &= \arctan2\big[\tau'\cos\left(\varphi_2\right)\sin\left(\lambda_2\right) + \left(1 - \tau'\right)\cos\left(\varphi_1\right)\sin\left(\lambda_1\right), \nonumber\\ & \tau'\cos\left(\varphi_2\right)\cos\left(\lambda_2\right) + \left(1 - \tau'\right)\cos\left(\varphi_1\right)\cos\left(\lambda_1\right)\big]\left( + 2\pi\text{ im Fall }\arctan2<0\right). \end{align} \]

For the distance $d$ between two points $\left(r_a, \varphi_a, \lambda_a\right)$ and $\left(r_b, \varphi_b, \lambda_b\right)$ applies

\[ \begin{align} d &= \sqrt{\left(x_b - x_a\right)^2 + \left(y_b - y_a\right)^2 + \left(z_b - z_a\right)^2}\nonumber\\ &= \sqrt{\left(r_b\cos\left(\varphi_b\right)\cos\left(\lambda_b\right) - r_a\cos\left(\varphi_a\right)\cos\left(\lambda_a\right)\right)^2 + \left(r_b\cos\left(\varphi_b\right)\sin\left(\lambda_b\right) - r_a\cos\left(\varphi_a\right)\sin\left(\lambda_a\right)\right)^2}\nonumber\\ & \newoverline{ + \left(r_b\sin\left(\varphi_b\right) - r_a\sin\left(\varphi_a\right)\right)^2}\nonumber\\ &= \sqrt{r_b^2\cos^2\left(\varphi_b\right) + r_a^2\cos^2\left(\varphi_a\right) + r_b^2\sin^2\left(\varphi_b\right) + r_a^2\sin^2\left(\varphi_a\right) - 2r_ar_b\cos\left(\varphi_b\right)\cos\left(\lambda_b\right)\cos\left(\varphi_a\right)\cos\left(\lambda_a\right)}\nonumber\\ & \newoverline{ - 2r_ar_b\cos\left(\varphi_b\right)\sin\left(\lambda_b\right)\cos\left(\varphi_a\right)\sin\left(\lambda_a\right) - 2r_br_a\sin\left(\varphi_b\right)\sin\left(\varphi_a\right)}\nonumber\\ &= \sqrt{r_a^2 + r_b^2 - 2r_ar_b\cos\left(\varphi_b\right)\cos\left(\varphi_a\right)\left(\cos\left(\lambda_b\right)\cos\left(\lambda_a\right) + \sin\left(\lambda_b\right)\sin\left(\lambda_a\right)\right) - 2r_br_a\sin\left(\varphi_b\right)\sin\left(\varphi_a\right)}\nonumber\\ &= \sqrt{r_a^2 + r_b^2 - 2r_ar_b\left(\cos\left(\varphi_a\right)\cos\left(\varphi_b\right)\cos\left(\lambda_a - \lambda_b\right) + \sin\left(\varphi_a\right)\sin\left(\varphi_b\right)\right)} \end{align} \]

In the case $r_a = r_b = r$ applies

The distance between two points on a sphere can be easily calculated from this. Let a sphere with radius $r$ be given. Then the distance $\Delta$ between two points on the surface, measured along the surface, is given by

\[ \begin{align} \Delta = r\theta, \end{align} \]

where $\theta$ is the angle between the two points, for this angle applies

\[ \begin{align} \sin\left(\frac{\theta}{2}\right) = \frac{d}{2r}. \end{align} \]

So the equation applies because of $\frac{\theta}{2}\in\left[0, \frac{\pi}{2}\right]$ for the distance $\Delta$

Now $\tau$ has to be converted into $\tau'$. It applies

\[ \begin{align} \tau &= \frac{1}{2} + \frac{1}{\theta}\arctan\left(\frac{\tau' - \frac{1}{2}}{\sqrt{r^2 - \frac{d^2}{4}}}d\right) = \frac{1}{2} + \frac{1}{\theta}\arctan\left(\frac{\tau' - \frac{1}{2}}{\sqrt{\frac{r^2}{d^2} - \frac{1}{4}}}\right)\nonumber\\ \Leftrightarrow\tau - \frac{1}{2} &= \frac{1}{\theta}\arctan\left(\frac{\tau' - \frac{1}{2}}{\sqrt{\frac{r^2}{d^2} - \frac{1}{4}}}\right)\nonumber\\ \Leftrightarrow\theta\left(\tau - \frac{1}{2}\right) &= \arctan\left(\frac{\tau' - \frac{1}{2}}{\sqrt{\frac{r^2}{d^2} - \frac{1}{4}}}\right)\nonumber\\ \Leftrightarrow\tan\left[\theta\left(\tau - \frac{1}{2}\right)\right] &= \frac{\tau' - \frac{1}{2}}{\sqrt{\frac{r^2}{d^2} - \frac{1}{4}}}\nonumber\\ \Leftrightarrow\sqrt{\frac{r^2}{d^2} - \frac{1}{4}}\tan\left[\theta\left(\tau - \frac{1}{2}\right)\right] &= \tau' - \frac{1}{2}\nonumber. \end{align} \]

From this it follows

D.2.2 Vertical surfaces

For vertical surfaces $A_v$ with the base length $L$ and the same designation for the radii applies

\[ \begin{align} A_v = \int_{r_0}^{r_1}L\frac{r}{r_0}dr = \frac{L}{2r_0}\left(r_1^2 - r_0^2\right). \end{align} \]

D.2.3 Volumes

Volumes $V$ with the base area $A$, the inner radius $r_0$ and the outer radius $r_1$ are calculated via

\[ \begin{align} V = \int_{r_0}^{r_1}A\frac{r^2}{r_0^2}dr = \frac{A}{r_0^2}\frac{1}{3}\left(r_1^3 - r_0^3\right). \end{align} \]

D.2.4 Spherical triangles

D.2.5 Spherical Voronoi lattices

D.3 Earth's gravity field

D.3.1 Properties of ellipses and ellipsoids

An ellipse in the xz plane is described by

Here, $a > 0$ is the semi-major axis (maximum extent of the ellipse in the x-direction) and $0 < c < a$ is the semi-minor axis (maximum extent of the ellipse in the z-direction). The eccentricity becomes

defined. This gives you

\[ \begin{align} c = a\left(1 - \epsilon\right). \end{align} \]

Follow in polar coordinates with

\[ \begin{align} x = r\cos\left(\varphi\right), & {} & z = r\sin\left(\varphi\right) \end{align} \]

the equations

\[ \begin{align} \frac{r^2}{a^2}\cos^2\left(\varphi\right) + \frac{r^2}{c^2}\sin^2\left(\varphi\right) &= 1\nonumber\\ \Rightarrow r\left(\varphi\right) &= \left[\frac{\cos^2\left(\varphi\right)}{a^2} + \frac{\sin^2\left(\varphi\right)}{c^2}\right]^{-1/2} = a\left[\cos^2\left(\varphi\right) + \frac{a^2}{c^2}\sin^2\left(\varphi\right)\right]^{-1/2}\nonumber\\ \Rightarrow r\left(\varphi\right) &= a\left[1 + \left(\frac{a^2}{c^2} - 1\right)\sin^2\left(\varphi\right)\right]^{-1/2}. \end{align} \]

With

\[ \begin{align} \frac{a^2}{c^2} - 1 &= \frac{1}{\left(1 - \epsilon\right)^2} - 1 = \frac{1 - \left(1 - \epsilon\right)^2}{\left(1 - \epsilon\right)^2} = \frac{2\epsilon - \epsilon^2}{\left(1 - \epsilon\right)^2}\tag{D.29}\label{eq:ellips_prop_1} \end{align} \]

you can do this

rewrite. It follows from equation (D.29)

\[ \begin{align} a^2 - c^2 &= c^2\frac{2\epsilon - \epsilon^2}{\left(1 - \epsilon\right)^2} = a^2\epsilon\left(2 - \epsilon\right)\tag{D.31}\label{eq:ellips_prop_2}. \end{align} \]

A ellipsoid of rotation or simply ellipsoid is created by letting the ellipse just examined rotate around the z-axis. Its volume $V$ is calculated as

\[ \begin{align} V &= \int_{ - c}^{c}\int_{0}^{a\sqrt{1 - \frac{z^2}{c^2}}}2\pi xdxdz = 2\pi\int_{ - c}^{c}\left[\frac{1}{2}x^2\right]_0^{a\sqrt{1 - \frac{z^2}{c^2}}}dz\nonumber\\ &= 2\pi\int_{ - c}^{c}\frac{1}{2}a^2\left(1 - \frac{z^2}{c^2}\right)dz = \pi\int_{ - c}^{c}a^2 - \frac{a^2}{c^2}z^2dz = \pi 2a^2c - \pi\frac{a^2}{c^2}\frac{2}{3}c^3\nonumber\\ &= \pi 2a^2c - \frac{2}{3}\pi a^2c\nonumber \end{align} \]

D.3.2 Gravity field of an ellipsoid

The ellipsoidal coordinates will be based on the Earth's geopotential $\phi_g$, for which applies

\[ \begin{align} \phi_g = \phi_z + \phi_0, \end{align} \]

Here $\phi_z$ is the centrifugal potential and $\phi_0$ is the gravitational potential of the earth. The following applies to the fleeing potential

\[ \begin{align} \phi_z = -\frac{1}{2}\left(\omega^2\left(x^2 + y^2\right)\right), \tag{D.34}\label{eq:zentrifugal_potential} \end{align} \]

because gradient formation results in

\[ \begin{align} - \nabla\phi_z = \omega^2\left(\begin{array}{c} x\\ y\\ 0 \end{array}\right) = -\left(\begin{array}{c} 0\\ 0\\ \omega \end{array}\right)\times\left(\begin{array}{c} - \omega y\\ \omega x\\ 0 \end{array}\right) = - \omegabi\times\left(\omegabi\times\mathbf{r}\right). \end{align} \]

The context applies

\[ \begin{align} \mathbf{a} = -\nabla\phi_g \end{align} \]

between the gravitational potential $\phi_g$ and the gravitational acceleration $\mathbf{a}$. Newton's law of gravity can then be written as

\[ \begin{align} \nabla\cdot\mathbf{a} = -4\pi G\rho \end{align} \]

with $\rho$ as the mass distribution, because in the case of a point mass this results in $M$ at the origin using Gauss's theorem

\[ \begin{align} - 4\pi G M = a4\pi r^2\Leftrightarrow a = -G\frac{M}{r^2}. \end{align} \]

In terms of the gravitational potential this becomes

\[ \begin{align} \Delta\phi_g = 4\pi G\rho \end{align} \]

or

\[ \begin{align} \Delta\phi_g = 0\tag{D.40}\label{eq:laplace_dgl} \end{align} \]

outside the earth. This is called a Laplace differential equation. This has to be solved in order to determine the gravitational potential of the earth. In general you can use an integral

\[ \begin{align} \phi_0\left(\mathbf{r}\right) = -G\rho\int_{E}^{}\frac{1}{\left|\mathbf{r} - \mathbf{r}'\right|}d^3r' \end{align} \]

set up for this, here $E$ is the earth, but this cannot be solved analytically in the case of an ellipsoidal set $E$. Therefore, $\phi_g$ is developed according to certain functions. The approach used is $\mu \coloneqq\cos\left(\theta\right)$, a separation approach

\[ \begin{align} \phi_g\left(r, \mu, \phi\right) = \frac{U\left(r\right)}{r}P\left(\mu\right)Q\left(\phi\right).\tag{D.42}\label{eq:ansatz_poisson} \end{align} \]

made. If you substitute this into $\Delta\phi_g = 0$, you get with Eq. (B.91) and the statement

\[ \begin{align} \frac{1}{r}\frac{\partial^2}{\partial r^2}\left(r f\right) = \frac{1}{r}\frac{\partial}{\partial r}\left(f + rf'\right) = \frac{2}{r}f' + f'' = \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial f}{\partial r}\right), \end{align} \]

here $f = f\left(r\right)$ and $f' \coloneqq \frac{\partial f}{\partial r}$

\[ \begin{align} PQ\frac{U''}{r} + UQ\frac{1}{r^3}\frac{\partial}{\partial\mu}\left(1 - \mu^2\right)P' + \frac{UP}{r^3\left(1 - \mu^2\right)}Q'' &= 0\nonumber\\ \Leftrightarrow r^2\left(1 - \mu^2\right)\frac{1}{v}U'' + \frac{1}{P}\left(1 - \mu^2\right)\frac{\partial}{\partial\mu}\left(1 - \mu^2\right)P' &= -\frac{Q''}{Q}.\tag{D.44}\label{eq:deriv_legendre_dgl} \end{align} \]

The left side of Eq. (D.44) does not depend on $\phi$, while the right-hand side does not depend on $\mu$ and $r$, therefore both are equal to a separation constant $m^2$. This therefore applies to $Q$

\[ \begin{align} Q'' = -m^2Q, \end{align} \]

what is solved by

\[ \begin{align} Q\left(\phi\right) = \exp\left(im\phi\right). \end{align} \]

Since $Q\left(\phi\right) = Q\left(\phi + 2\pi\right)$, $m\in \mathbb{Z}$. Now divide Eq. (D.44) by $1 - \mu^2$ and sets the right side equal to $m^2$, you then get

\[ \begin{align} r^2\frac{U''}{v} = -\frac{1}{P}\frac{\partial}{\partial\mu}\left(1 - \mu^2\right)P'\left(\mu\right) + \frac{m^2}{1 - \mu^2}. \end{align} \]

Both sides are equal to a new separation constant $\lambda$. You get two differential equations, one for $P$ and one for $U$, which reads for $P$

\[ \begin{align} \frac{d}{d\mu}\left(1 - \mu^2\right)P'\left(\mu\right) = P\left(-\lambda + \frac{m^2}{1 - \mu^2}\right)\tag{D.48}\label{eq:legendre_vorform} \end{align} \]

and that for $U$ is

\[ \begin{align} r^2U'' - \lambda U = 0.\tag{D.49}\label{eq:dgl_u} \end{align} \]

Cylindrical symmetry is sufficient for the geopotential, i.e. $m = 0$. Then you get

\[ \begin{align} \frac{d}{d\mu}\left(1 - \mu^2\right)P'\left(\mu\right) = -\lambda P.\tag{D.50}\label{eq:legendre_dgl} \end{align} \]

If you put in Eq. (D.49) $\lambda = n\left(n + 1\right)$, you get

\[ \begin{align} r^2U'' - n\left(n + 1\right)U = 0. \end{align} \]

$U$ depends on $n$, write $U_n$. It is found that

\[ \begin{align} U_n\left(r\right) = r^{n + 1}\tag{D.52}\label{eq:radial_kugel_dgl} \end{align} \]

and

\[ \begin{align} U_n\left(r\right) = r^{-n} \end{align} \]

solutions are and thus the general solution of the ordinary second-order DGL Eq. (D.52) is given by

\[ \begin{align} U_n\left(r\right) = a_nr^{n + 1} + \frac{b_n}{r^n}. \end{align} \]

Now solutions are of the form

\[ \begin{align} \phi_g^{(n)}\left(r, \mu\right) = \frac{U_n\left(r\right)}{r}P_n\left(\mu\right) \end{align} \]

found, since the Laplace equation (D.40) is linear in $\phi_g$, any linear combinations of these are again solutions and one obtains as a general solution

\[ \begin{align} \phi_g\left(r, \mu\right) = \sum_{n = 0}^{\infty}\left(a_nr^n + \frac{b_n}{r^{n + 1}}\right)P_n\left(\mu\right).\tag{D.56}\label{eq:allg_laplace_loesung} \end{align} \]

The specific solution for the rotating ellipsoid case is determined by the boundary conditions. This is:

\[ \begin{align} \phi_g\left(x, z\right) = 0 \end{align} \]

for

\[ \begin{align} \frac{x^2}{a^2} + \frac{z^2}{c^2} = 1, \end{align} \]

Here $a$ is the semi-major and $c$ the semi-minor axis of the Earth. At this point the sum in Eq. (D.56) is terminated at $n = 2$ and only the gravitational component is considered:

\[ \begin{align} \phi_0 &= \left(a_0 + \frac{b_0}{r}\right)P_0\left(\sin\left(\chi\right)\right) + \left(a_1r + \frac{b_1}{r^2}\right)P_1\left(\sin\left(\chi\right)\right) + \left(a_2r^2 + \frac{b_2}{r^3}\right)P_2\left(\sin\left(\chi\right)\right) \end{align} \]

$\chi$ is the geocentric latitude (as opposed to the geodetic latitude). We have $\chi + \theta = \pi/2$. The first three Legendre polynomials are

\[ \begin{align} P_0\left(\sin\left(\chi\right)\right) &= 1,\\ P_1\left(\sin\left(\chi\right)\right) &= \frac{1}{2}2\sin\left(\chi\right) = \sin\left(\chi\right),\\ P_2\left(\sin\left(\chi\right)\right) &= \frac{1}{8}\left(12\sin^2\left(\chi\right) - 4\right) = \frac{3}{2}\sin^2\left(\chi\right) - \frac{1}{2}. \end{align} \]

So the expression for $\phi_g$ is

\[ \begin{align} \phi_0 = a_0 + \frac{b_0}{r} + \left(a_1r + \frac{b_1}{r^2}\right)\sin\left(\chi\right) + \left(\frac{a_2}{2}r^2 + \frac{b_2^2}{2r^3}\right)\left(3\sin^2\left(\chi\right) - 1\right). \end{align} \]

The goal is now to calculate the coefficients $a_0, a_1, a_2, b_0, b_1, b_2$. However, because of the usual boundary condition $\lim\limits_{r\to\infty}\phi_0\left(r\right) = 0$ applies

\[ \begin{align} a_0 = a_1 = a_2 = 0, \end{align} \]

furthermore, the solution should become $-G\frac{M}{r}$ at an infinite distance, since the ellipsoidal mass distribution looks like a point mass from a distance, so it is expected

\[ \begin{align} b_0 = -GM. \end{align} \]

The coefficients $b_1$ and $b_2$ are still free:

\[ \begin{align} \phi_0 = -\frac{GM}{r} + \frac{b_1}{r^2}\sin\left(\chi\right) + \frac{b_2}{2r^3}\left(3\sin^2\left(\chi\right) - 1\right). \end{align} \]

To determine this, one calculates the gravitational potential on the z-axis ($z = r> c$). There we have $\chi = \pi/2$

\[ \begin{align} \phi_0\left(r\right) = -\frac{GM}{r} + \frac{b_1\left(\epsilon\right)}{r^2} + \frac{b_2\left(\epsilon\right)}{r^3}.\tag{D.67}\label{eq:grav_z_dritte_ordnung} \end{align} \]

It was assumed that the $b_i$ depend on the eccentricity $\epsilon$. So one assumes a homogeneous earth of density $\rho$ and integrates

\[ \begin{align} \phi_0\left(r\right) &= -G\rho\int_{z = -c}^{c}\int_{0}^{a\sqrt{1 - \frac{z^2}{c^2}}}\frac{1}{\sqrt{\left(r - z\right)^2 + x^2}}2\pi xdxdz = -2\pi G\rho\int_{z = -c}^{c}\left[\sqrt{\left(r - z\right)^2 + x^2}\right]_0^{a\sqrt{1 - \frac{z^2}{c^2}}}dz\nonumber\\ &= -2\pi G\rho\int_{ - c}^c\sqrt{\left(z - r\right)^2 + a^2\left(1 - \frac{z^2}{c^2}\right)} - r + zdz = -2\pi G\rho\int_{ - c}^{c}\sqrt{r^2 + z^2\left(1 - \frac{a^2}{c^2}\right) - 2zr + a^2} - r + zdz\nonumber\\ &= 2\pi G\rho r2c - 2\pi G\rho\int_{ - c}^{c}\sqrt{\left(1 - \frac{a^2}{c^2}\right)z^2 - 2zr + a^2 + r^2}dz\nonumber\\ &= 4\pi G\rho cr - 2\pi G\rho\int_{ - c}^{c}\sqrt{Az^2 + Bz + C}dz = 2\pi G\rho\left(2cr - \int_{ - c}^{c}\sqrt{Az^2 + Bz + C}dz\right) \end{align} \]

with

\[ \begin{align} A &= 1 - \frac{a^2}{c^2},\\ B &= -2r,\\ C &= a^2 + r^2. \end{align} \]

This formula can first be checked for a simple limiting case, namely the case of a ball $a = c$. In this case you get

\[ \begin{align} \phi_0\left(r\right) &= 2\pi G\rho\left(2ar - \int_{ - c}^{c}\sqrt{r^2 - 2zr + a^2}dz\right) = 2\pi G\rho\left(2ar - \left[-\frac{1}{3r}\left(r^2 - 2zr + a^2\right)^{3/2}\right]_{ - c}^c\right)\nonumber\\ &= 2\pi G\rho\left(2ar + \left[\frac{1}{3r}\left(r^2 - 2ar + a^2\right)^{3/2} - \frac{1}{3r}\left(r^2 + 2ar + a^2\right)\right]\right)\nonumber\\ &= 2\pi G\rho\left(2ar + \frac{1}{3r}\left[\left(r - a\right)^{3} - \left(r + a\right)^{3}\right]\right)\nonumber\\ &= 2\pi G\rho\left(2ar + \frac{1}{3r}\left[\left(r^2 - a^2 - 2ar\right)\left(r - a\right) - \left(r^2 + a^2 + 2ar\right)\left(r + a\right)\right]\right)\nonumber\\ &= 2\pi G\rho(2ar + \frac{1}{3r}[r^3 + r^2a + a^2r - a^3 - 2ar^2 + 2a^2r - r^3 - r^2a - a^2r - a^3 - 2ar^2 - 2a^2r])\nonumber\\ &= 2\pi G\rho\left(2ar + \frac{1}{3r}\left[-6r^2a - 2a^3\right]\right) = 2\pi G\rho\left(2ar - 2ra - \frac{2}{3}\frac{a^3}{r}\right) = -\frac{4\pi G\rho a^3}{3r}\nonumber\\ &= - GM\frac{1}{r}. \end{align} \]

The gravitational potential of a homogeneous solid sphere outside the sphere is equal to that of a point mass at the location of the sphere's center of gravity. According to [11] Eq. (3.3.37)

\[ \begin{align} \int\sqrt{Ax^2 + Bx + C}dx &= \frac{2Ax + B}{4A}\sqrt{Ax^2 + Bx + C} + \frac{4AC - B^2}{8A}\int\frac{dx}{\sqrt{Ax^2 + Bx + C}}. \end{align} \]

Here is

\[ \begin{align} \int\frac{dx}{\sqrt{Ax^2 + Bx + C}} &= -\frac{1}{\sqrt{ - A}}\arcsin\frac{2Ax + B}{\sqrt{B^2 - 4AC}}. \end{align} \]

The requirements for this are

\[ \begin{align} A&<0,\\ B^2>4AC\Leftrightarrow 4r^2 &> \left(1 - \frac{a^2}{c^2}\right)\left(a^2 + r^2\right)\nonumber\\ \Leftrightarrow r^2&>a^2 + r^2 - \frac{a^4}{c^2} - a^2\frac{r^2}{c^2}, \nonumber\\ |2Ax + B|<\sqrt{B^2 - 4AC}\Leftrightarrow \left|2\left(1 - \frac{a^2}{c^2}\right)x - 2r\right|&<\sqrt{4r^2 + 4\left(a^2 + r^2\right)\left(\frac{a^2}{c^2} - 1\right)}\nonumber\\ \Leftrightarrow \left|x\left(1 - \frac{a^2}{c^2}\right) - r\right|&<\sqrt{r^2 - a^2 - r^2 + \frac{a^4}{c^2} + r^2\frac{a^2}{c^2}}\nonumber\\ \Leftarrow c\left(\frac{a^2}{c^2} - 1\right) + r&

For $r = c$ the last two expressions become equal. The derivative of the left expression with respect to $r$ is

\[ \begin{align} \frac{d}{dr}\left(c\left(\frac{a^2}{c^2} - 1\right) + r\right) = 1, \end{align} \]

is the derivative of the right-hand expression with respect to $r$

\[ \begin{align} \frac{d}{dr}a\sqrt{\frac{a^2}{c^2} - 1 + \frac{r^2}{c^2}} = a\frac{2r}{c^2}\frac{1}{2}\frac{1}{\sqrt{\frac{a^2}{c^2} - 1 + \frac{r^2}{c^2}}} = \frac{ra}{c\sqrt{a^2 - c^2 + r^2}}. \end{align} \]

This expression is equal to one for $r = c$ and $a/c>1$ for $r\to\infty$, the derivative

\[ \begin{align} \frac{d}{dr}\frac{ra}{c\sqrt{a^2 - c^2 + r^2}}&\propto& ac\sqrt{a^2 - c^2 + r^2} - racr\frac{1}{\sqrt{a^2 - c^2 + r^2}}\nonumber\\ &\propto&ac\left(a^2 - c^2 + r^2\right) - r^2ac > 0 \end{align} \]

is positive everywhere and therefore is

\[ \begin{align} \frac{d}{dr}a\sqrt{\frac{a^2}{c^2} - 1 + \frac{r^2}{c^2}}>1 \end{align} \]

for $r>c$. Thus the inequality (16.242) for $r> c$ is satisfied and the integral formula

\[ \begin{align} \int\sqrt{Ax^2 + Bx + C}dx &= \frac{2Ax + B}{4A}\sqrt{Ax^2 + Bx + C} - \frac{4AC - B^2}{8A}\frac{1}{\sqrt{ - A}}\arcsin\frac{2Ax + B}{\sqrt{B^2 - 4AC}} \end{align} \]

can be applied. So it applies

\[ \begin{align} & \int_{ - c}^{c}\sqrt{Ax^2 + Bx + C}dx\nonumber\\ &= \frac{2Ac + B}{4A}\cdot\sqrt{Ac^2 + Bc + C} - \frac{4AC - B^2}{8A}\frac{1}{\sqrt{ - A}}\arcsin\left(\frac{2Ac + B}{\sqrt{B^2 - 4AC}}\right) - \frac{B - 2Ac}{4A}\cdot\sqrt{Ac^2 - Bc + C}\nonumber\\ & + \frac{4AC - B^2}{8A}\frac{1}{\sqrt{ - A}}\arcsin\left(\frac{ - 2Ac + B}{\sqrt{B^2 - 4AC}}\right) = \frac{2Ac + B}{4A}\sqrt{Ac^2 + Bc + C} + \frac{4AC - B^2}{8A\sqrt{ - A}}\cdot\nonumber\\ & \cdot\left[\arcsin\left(\frac{B - 2Ac}{\sqrt{B^2 - 4AC}}\right) - \arcsin\left(\frac{2Ac + B}{\sqrt{B^2 - 4AC}}\right)\right] - \frac{B - 2Ac}{4A}\sqrt{Ac^2 - Bc + C}. \end{align} \]

Here you now insert $A, B$ and $C$ and get

\[ \begin{align} \sqrt{Ac^2 + Bc + C} &= \sqrt{c^2 - a^2 - 2rc + a^2 + r^2} = \sqrt{c^2 + r^2 - 2rc} = r - c,\\ \sqrt{Ac^2 - Bc + C} &= \sqrt{c^2 - a^2 + 2rc + a^2 + r^2} = r + c,\\ B - 2Ac &= -2r - 2c + 2\frac{a^2}{c} = 2\left(\frac{a^2}{c} - r - c\right),\\ B + 2Ac &= -2r + 2c - 2\frac{a^2}{c} = 2\left(c - r - \frac{a^2}{c}\right),\\ B^2 - 4AC &= 4r^2 - 4\left(1 - \frac{a^2}{c^2}\right)\left(a^2 + r^2\right) = -4a^2 + 4\frac{a^4}{c^2} + 4\frac{a^2r^2}{c^2} = 4\frac{a^2}{c^2}\left(a^2 + r^2 - c^2\right),\\ & \frac{2Ac + B}{4A}\sqrt{Ac^2 + Bc + C} - \frac{B - 2Ac}{4A}\sqrt{Ac^2 - Bc + C}\nonumber\\ &= \frac{2\left(c - r - \frac{a^2}{c}\right)}{4 - 4\frac{a^2}{c^2}}\left(r - c\right) - \frac{2\left(\frac{a^2}{c} - r - c\right)}{4 - 4\frac{a^2}{c^2}}\left(r + c\right)\nonumber\\ &= \frac{4cr - 2a^2r/c}{2 - 2a^2/c^2} = \frac{2c^3r - a^2rc}{c^2 - a^2} = rc\frac{a^2 - 2c^2}{a^2 - c^2},\\ \frac{4AC - B^2}{8A\sqrt{ - A}} &= -\frac{4\frac{a^2}{c^2}\left(a^2 + r^2 - c^2\right)}{8\left(1 - \frac{a^2}{c^2}\right)\sqrt{\frac{a^2}{c^2} - 1}} = -\frac{a^2\left(a^2 + r^2 - c^2\right)}{2\left(c^2 - a^2\right)\sqrt{\frac{a^2}{c^2} - 1}}\nonumber\\ &= -\frac{a^2 + r^2 - c^2}{2\left(\frac{c^2}{a^2} - 1\right)\sqrt{\frac{a^2}{c^2} - 1}} = a^2c\frac{a^2 + r^2 - c^2}{2\left(a^2 - c^2\right)\sqrt{a^2 - c^2}},\\ \frac{B - 2Ac}{\sqrt{B^2 - 4AC}} &= \frac{2\left(\frac{a^2}{c} - r - c\right)}{2\frac{a}{c}\sqrt{a^2 + r^2 - c^2}} = \frac{a^2 - rc - c^2}{a\sqrt{a^2 + r^2 - c^2}},\\ \frac{2Ac + B}{\sqrt{B^2 - 4AC}} &= \frac{2\left(c - r - \frac{a^2}{c}\right)}{2\frac{a}{c}\sqrt{a^2 + r^2 - c^2}} = \frac{c^2 - rc - a^2}{a\sqrt{a^2 + r^2 - c^2}} \end{align} \]

the equation

\[ \begin{align} \int_{ - c}^c\dotsc dx &= a^2c\frac{a^2 + r^2 - c^2}{2\left(a^2 - c^2\right)^{3/2}}\cdot\bigg[\arcsin\left(\frac{a^2 - rc - c^2}{a\sqrt{a^2 + r^2 - c^2}}\right) - \arcsin\left(\frac{c^2 - rc - a^2}{a\sqrt{a^2 + r^2 - c^2}}\right)\bigg]\nonumber\\ & + rc\frac{a^2 - 2c^2}{a^2 - c^2}. \end{align} \]

So for the gravitational potential on the z-axis it follows

\[ \begin{align} 4cr - 2cr\frac{a^2 - 2c^2}{a^2 - c^2} = rc\frac{4a^2 - 4c^2 - 2a^2 + 4c^2}{a^2 - c^2} = 2rc\frac{a^2}{a^2 - c^2} \end{align} \]

the bill

\[ \begin{align} \phi_0\left(r\right) &= \pi G\rho\bigg(4cr - 2rc\frac{a^2 - 2c^2}{a^2 - c^2} - a^2c\frac{a^2 + r^2 - c^2}{\left(a^2 - c^2\right)^{3/2}}\cdot\left[\arcsin\left(\frac{a^2 - rc - c^2}{a\sqrt{a^2 + r^2 - c^2}}\right) - \arcsin\left(\frac{c^2 - rc - a^2}{a\sqrt{a^2 + r^2 - c^2}}\right)\right]\bigg)\nonumber \end{align} \] \[ \begin{align} &= \pi G\rho\bigg(rc\frac{2a^2}{a^2 - c^2} - a^2c\frac{a^2 + r^2 - c^2}{\left(a^2 - c^2\right)^{3/2}}\cdot\left[\arcsin\left(\frac{a^2 - rc - c^2}{a\sqrt{a^2 + r^2 - c^2}}\right) - \arcsin\left(\frac{c^2 - rc - a^2}{a\sqrt{a^2 + r^2 - c^2}}\right)\right]\bigg)\nonumber\\ &= \pi G\rho\frac{a^2c}{a^2 - c^2}\bigg(2r - \frac{a^2 + r^2 - c^2}{\sqrt{a^2 - c^2}}\cdot\left[\arcsin\left(\frac{a^2 - rc - c^2}{a\sqrt{a^2 + r^2 - c^2}}\right) - \arcsin\left(\frac{c^2 - rc - a^2}{a\sqrt{a^2 + r^2 - c^2}}\right)\right]\bigg). \end{align} \]

Here you can see the mass of the earth

\[ \begin{align} M = \frac{4}{3}\pi\rho a^2c\Rightarrow \pi\rho a^2c = \frac{3M}{4} \end{align} \]

install:

\[ \begin{align} \phi_0\left(r\right) &= G\frac{M}{a^2 - c^2}\bigg(\frac{3}{2}r - \frac{3}{4}\frac{a^2 + r^2 - c^2}{\sqrt{a^2 - c^2}}\cdot\nonumber\\ & \cdot\left[\arcsin\left(\frac{a^2 - rc - c^2}{a\sqrt{a^2 + r^2 - c^2}}\right) - \arcsin\left(\frac{c^2 - rc - a^2}{a\sqrt{a^2 + r^2 - c^2}}\right)\right]\bigg) \end{align} \]

With Eq. (D.31) follows for the gravitational potential

\[ \begin{align} \phi_0\left(r\right) &= G\frac{M}{a^2\epsilon\left(2 - \epsilon\right)}\bigg(\frac{3}{2}r - \frac{3}{4}\frac{a^2\epsilon\left(2 - \epsilon\right) + r^2}{a\sqrt{\epsilon\left(2 - \epsilon\right)}}\cdot\nonumber\\ & \cdot\left[\arcsin\left(\frac{a^2\epsilon\left(2 - \epsilon\right) - ra\left(1 - \epsilon\right)}{a\sqrt{a^2\epsilon\left(2 - \epsilon\right) + r^2}}\right) - \arcsin\left(-\frac{a^2\epsilon\left(2 - \epsilon\right) + ra\left(1 - \epsilon\right)}{a\sqrt{a^2\epsilon\left(2 - \epsilon\right) + r^2}}\right)\right]\bigg)\nonumber\\ &= G\frac{M}{a\epsilon\left(2 - \epsilon\right)}\bigg(\frac{3}{2}\frac{r}{a} - \frac{3}{4}\frac{\epsilon\left(2 - \epsilon\right) + r^2/a^2}{\sqrt{\epsilon\left(2 - \epsilon\right)}}\cdot\nonumber\\ & \cdot\left[\arcsin\left(\frac{\epsilon\left(2 - \epsilon\right) - \frac{r}{a}\left(1 - \epsilon\right)}{\sqrt{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}}\right) - \arcsin\left(-\frac{\epsilon\left(2 - \epsilon\right) + \frac{r}{a}\left(1 - \epsilon\right)}{\sqrt{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}}\right)\right]\bigg)\nonumber\\ &= \frac{3GM}{2a\epsilon\left(2 - \epsilon\right)}\bigg(\frac{r}{a} - \frac{\epsilon\left(2 - \epsilon\right) + r^2/a^2}{2\sqrt{\epsilon\left(2 - \epsilon\right)}}\cdot\bigg[\arcsin\left(\frac{\epsilon\left(2 - \epsilon\right) - \frac{r}{a}\left(1 - \epsilon\right)}{\sqrt{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}}\right)\nonumber\\ & + \arcsin\left(\frac{\epsilon\left(2 - \epsilon\right) + \frac{r}{a}\left(1 - \epsilon\right)}{\sqrt{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}}\right)\bigg]\bigg).\tag{D.97}\label{eq:grav_pot_ex_z} \end{align} \]

Eq. (D.97) is the basis for the theoretical investigation of the gravity field of an ellipsoid of revolution. It is exact for eccentricities $\epsilon\not = 0$, but is not applicable for $\epsilon = 0$. Unfortunately, the expression does not yet have the form of Equation. (D.67). To achieve this, it is expanded in the small parameter $\epsilon\ll 1$. To do this, you define the abbreviations

\[ \begin{align} f_1\left(\epsilon\right)& \coloneqq&\frac{\epsilon\left(2 - \epsilon\right) - \frac{r}{a}\left(1 - \epsilon\right)}{\sqrt{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}},\\ f_2\left(\epsilon\right)& \coloneqq&\frac{\epsilon\left(2 - \epsilon\right) + \frac{r}{a}\left(1 - \epsilon\right)}{\sqrt{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}}. \end{align} \]

Preparations are already underway

\[ \begin{align} f_1'\left(\epsilon\right) &= \frac{\left(2 - 2\epsilon + \frac{r}{a}\right)\sqrt{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}} - \left[\epsilon\left(2 - \epsilon\right) - \frac{r}{a}\left(1 - \epsilon\right)\right]\frac{1}{2}\left(2 - 2\epsilon\right)\frac{1}{\sqrt{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}}}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}\nonumber\\ &= \frac{\left(2 - 2\epsilon + \frac{r}{a}\right)\left(\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right) - \left[\epsilon\left(2 - \epsilon\right) - \frac{r}{a}\left(1 - \epsilon\right)\right]\left(1 - \epsilon\right)}{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^{3/2}}\nonumber\\ &= \frac{2\left(1 - \epsilon\right)\epsilon\left(2 - \epsilon\right) + 2\left(1 - \epsilon\right)\frac{r^2}{a^2} + \frac{r}{a}\epsilon\left(2 - \epsilon\right) + \frac{r^3}{a^3} - \epsilon\left(2 - \epsilon\right)\left(1 - \epsilon\right) + \frac{r}{a}\left(1 - \epsilon\right)^2}{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^{3/2}}\nonumber\\ &= \frac{\epsilon\left(1 - \epsilon\right)\left(2 - \epsilon\right) + \frac{r^3}{a^3} + \frac{r^2}{a^2}2\left(1 - \epsilon\right) + \frac{r}{a}}{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^{3/2}} \end{align} \]

and

\[ \begin{align} f_2'\left(\epsilon\right) &= \frac{\left(2 - 2\epsilon - \frac{r}{a}\right)\sqrt{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}} - \left[\epsilon\left(2 - \epsilon\right) + \frac{r}{a}\left(1 - \epsilon\right)\right]\frac{1}{2}\left(2 - 2\epsilon\right)\frac{1}{\sqrt{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}}}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}\nonumber\\ &= \frac{\left(2 - 2\epsilon - \frac{r}{a}\right)\left(\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right) - \left[\epsilon\left(2 - \epsilon\right) + \frac{r}{a}\left(1 - \epsilon\right)\right]\left(1 - \epsilon\right)}{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^{3/2}}\nonumber\\ &= \frac{2\left(1 - \epsilon\right)\epsilon\left(2 - \epsilon\right) + 2\left(1 - \epsilon\right)\frac{r^2}{a^2} - \frac{r}{a}\epsilon\left(2 - \epsilon\right) - \frac{r^3}{a^3} - \epsilon\left(2 - \epsilon\right)\left(1 - \epsilon\right) - \frac{r}{a}\left(1 - \epsilon\right)^2}{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^{3/2}}\nonumber\\ &= \frac{\epsilon\left(1 - \epsilon\right)\left(2 - \epsilon\right) - \frac{r^3}{a^3} + \frac{r^2}{a^2}2\left(1 - \epsilon\right) - \frac{r}{a}}{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^{3/2}}. \end{align} \]

The gravitational potential is written as:

\[ \begin{align} \phi_0\left(r\right) = \frac{3GM}{2a}\frac{2\frac{r}{a}\left(\epsilon\left(2 - \epsilon\right)\right)^{1/2} - \left(\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right)\left[\arcsin\left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right]}{2\left(\epsilon\left(2 - \epsilon\right)\right)^{3/2}}. \end{align} \]

This will be with

\[ \begin{align} g\left(\epsilon\right)&\coloneqq2\frac{r}{a}\left(\epsilon\left(2 - \epsilon\right)\right)^{1/2} - \left(\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right)\left[\arcsin\left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right],\\ h\left(\epsilon\right)&\coloneqq2\left(\epsilon\left(2 - \epsilon\right)\right)^{3/2},\\ F\left(\epsilon\right) &\coloneqq\frac{g\left(\epsilon\right)}{h\left(\epsilon\right)} \end{align} \]

to

\[ \begin{align} \phi_0\left(r\right) = \frac{3GM}{2a}\frac{g\left(\epsilon\right)}{h\left(\epsilon\right)} = \frac{3}{2}\frac{GM}{a}F\left(\epsilon\right). \end{align} \]

The goal now is to expand $F\left(\epsilon\right)$ to the first order in $\epsilon$:

\[ \begin{align} \phi_0\left(r\right) = \frac{3}{2}\frac{GM}{a}\left(F\left(0\right) + F'\left(0\right)\epsilon + \mathcal{O}\left(\epsilon^2\right)\right) \end{align} \]

Because $h\left(0\right) = 0$ this is not possible directly, but you have to do the replacements

\[ \begin{align} F\left(0\right) \to \lim\limits_{x\downarrow 0}F\left(x\right), F'\left(0\right) \to \lim\limits_{x\downarrow 0}F'\left(x\right) \end{align} \]

make. It's because

\[ \begin{align} f_1\left(0\right) = -1 = -f_2\left(0\right) \end{align} \]

also

\[ \begin{align} g\left(0\right) = 0, \end{align} \]

So one can apply L'Hospital's rule to find $\lim\limits_{\epsilon\to 0}\frac{g\left(\epsilon\right)}{h\left(\epsilon\right)}$. It applies

\[ \begin{align} g'\left(\epsilon\right) &= 2\frac{r}{a}\frac{1}{2}\left(2 - 2\epsilon\right)\frac{1}{\sqrt{\epsilon\left(2 - \epsilon\right)}} - \left(2 - 2\epsilon\right)\left[\arcsin\left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right], \nonumber\\ & - \left(\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right)\left[f_1'\left(\epsilon\right)\frac{1}{\sqrt{1 - f_1^2\left(\epsilon\right)}} + f_2'\left(\epsilon\right)\frac{1}{\sqrt{1 - f_2^2\left(\epsilon\right)}}\right]\nonumber\\ &= \frac{2r\left(1 - \epsilon\right)}{a\sqrt{\epsilon\left(2 - \epsilon\right)}} - 2\left(1 - \epsilon\right)\left[\arcsin\left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right]\nonumber\\ & - \left(\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right)\left[f_1'\left(\epsilon\right)\frac{1}{\sqrt{1 - f_1^2\left(\epsilon\right)}} + f_2'\left(\epsilon\right)\frac{1}{\sqrt{1 - f_2^2\left(\epsilon\right)}}\right]. \end{align} \]

They apply

\[ \begin{align} f_1^2\left(\epsilon\right) = \frac{\epsilon^2\left(2 - \epsilon\right)^2 + \frac{r^2}{a^2}\left(1 - \epsilon\right)^2 - 2\epsilon\left(2 - \epsilon\right)\frac{r}{a}\left(1 - \epsilon\right)}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}} \end{align} \]

and

\[ \begin{align} f_2^2\left(\epsilon\right) = \frac{\epsilon^2\left(2 - \epsilon\right)^2 + \frac{r^2}{a^2}\left(1 - \epsilon\right)^2 + 2\epsilon\left(2 - \epsilon\right)\frac{r}{a}\left(1 - \epsilon\right)}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}. \end{align} \]

Follow from this

\[ \begin{align} 1 - f_1^2\left(\epsilon\right) &= \frac{\epsilon\left(2 - \epsilon\right)\left(1 - \epsilon\left(2 - \epsilon\right)\right) + \frac{r^2}{a^2}\left(1 - \left(1 - \epsilon\right)^2\right) + 2\epsilon\left(2 - \epsilon\right)\frac{r}{a}\left(1 - \epsilon\right)}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}\nonumber\\ &= \frac{\epsilon\left(2 - \epsilon\right)\left(1 - 2\epsilon + \epsilon^2\right) + \frac{r^2}{a^2}\left(1 - 1 - \epsilon^2 + 2\epsilon\right) + 2\epsilon\left(2 - \epsilon\right)\frac{r}{a}\left(1 - \epsilon\right)}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}\nonumber\\ &= \frac{\epsilon\left(2 - \epsilon\right)\left(1 - \epsilon\right)^2 + \frac{r^2}{a^2}\epsilon\left(2 - \epsilon\right) + 2\epsilon\frac{r}{a}\left(2 - \epsilon\right)\left(1 - \epsilon\right)}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}} = \left(2 - \epsilon\right)\epsilon\frac{\left(1 - \epsilon\right)^2 + \frac{r^2}{a^2} + 2\frac{r}{a}\left(1 - \epsilon\right)}{\left(2 - \epsilon\right)\epsilon + \frac{r^2}{a^2}}\nonumber\\ &= \left(2 - \epsilon\right)\epsilon\frac{\left[\left(1 - \epsilon\right) + \frac{r}{a}\right]^2}{\left(2 - \epsilon\right)\epsilon + \frac{r^2}{a^2}} \end{align} \]

and

\[ \begin{align} 1 - f_2^2\left(\epsilon\right) &= \frac{\epsilon\left(2 - \epsilon\right)\left(1 - \epsilon\left(2 - \epsilon\right)\right) + \frac{r^2}{a^2}\left(1 - \left(1 - \epsilon\right)^2\right) - 2\epsilon\left(2 - \epsilon\right)\frac{r}{a}\left(1 - \epsilon\right)}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}\nonumber\\ &= \frac{\epsilon\left(2 - \epsilon\right)\left(1 - 2\epsilon + \epsilon^2\right) + \frac{r^2}{a^2}\left(1 - 1 - \epsilon^2 + 2\epsilon\right) - 2\epsilon\left(2 - \epsilon\right)\frac{r}{a}\left(1 - \epsilon\right)}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}\nonumber\\ &= \frac{\epsilon\left(2 - \epsilon\right)\left(1 - \epsilon\right)^2 + \frac{r^2}{a^2}\epsilon\left(2 - \epsilon\right) - 2\epsilon\frac{r}{a}\left(2 - \epsilon\right)\left(1 - \epsilon\right)}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}} = \left(2 - \epsilon\right)\epsilon\frac{\left(1 - \epsilon\right)^2 + \frac{r^2}{a^2} - 2\frac{r}{a}\left(1 - \epsilon\right)}{\left(2 - \epsilon\right)\epsilon + \frac{r^2}{a^2}}\nonumber\\ &= \left(2 - \epsilon\right)\epsilon\frac{\left[\left(1 - \epsilon\right) - \frac{r}{a}\right]^2}{\left(2 - \epsilon\right)\epsilon + \frac{r^2}{a^2}}. \end{align} \]

Follow from this

\[ \begin{align} \frac{1}{\sqrt{1 - f_1^2\left(\epsilon\right)}} = \frac{1}{1 - \epsilon + \frac{r}{a}}\sqrt{\frac{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}{\epsilon\left(2 - \epsilon\right)}} \end{align} \]

and

\[ \begin{align} \frac{1}{\sqrt{1 - f_2^2\left(\epsilon\right)}} = \frac{1}{\epsilon + \frac{r}{a} - 1}\sqrt{\frac{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}{\epsilon\left(2 - \epsilon\right)}}. \end{align} \]

It follows

\[ \begin{align} & \left(\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right)\left(f_1'\left(\epsilon\right)\frac{1}{\sqrt{1 - f_1^2\left(\epsilon\right)}} + f_2'\left(\epsilon\right)\frac{1}{\sqrt{1 - f_2^2\left(\epsilon\right)}}\right)\nonumber\\ &= \frac{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^{3/2}}{\sqrt{\epsilon\left(2 - \epsilon\right)}}\cdot\frac{1}{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^{3/2}}\cdot\bigg[\frac{\epsilon\left(1 - \epsilon\right)\left(2 - \epsilon\right) + \frac{r^3}{a^3} + \frac{r^2}{a^2}2\left(1 - \epsilon\right) + \frac{r}{a}}{1 - \epsilon + \frac{r}{a}}\nonumber\\ & + \frac{\epsilon\left(1 - \epsilon\right)\left(2 - \epsilon\right) - \frac{r^3}{a^3} + \frac{r^2}{a^2}2\left(1 - \epsilon\right) - \frac{r}{a}}{\epsilon + \frac{r}{a} - 1}\bigg]\nonumber\\ &= \frac{1}{\sqrt{\epsilon\left(2 - \epsilon\right)}}\bigg[\frac{\epsilon\left(1 - \epsilon\right)\left(2 - \epsilon\right) + \frac{r^3}{a^3} + \frac{r^2}{a^2}2\left(1 - \epsilon\right) + \frac{r}{a}}{1 - \epsilon + \frac{r}{a}} + \frac{\epsilon\left(1 - \epsilon\right)\left(2 - \epsilon\right) - \frac{r^3}{a^3} + \frac{r^2}{a^2}2\left(1 - \epsilon\right) - \frac{r}{a}}{\epsilon + \frac{r}{a} - 1}\bigg]\nonumber\\ &= \frac{1}{\sqrt{\epsilon\left(2 - \epsilon\right)}}\frac{2\frac{r}{a}\epsilon\left(1 - \epsilon\right)\left(2 - \epsilon\right) + 4\left(1 - \epsilon\right)\frac{r^3}{a^3} - 2\left(1 - \epsilon\right)\frac{r^3}{a^3} - 2\frac{r}{a}\left(1 - \epsilon\right)}{\frac{r^2}{a^2} - 1 - \epsilon^2 + 2\epsilon}\nonumber\\ &= \frac{2\frac{r}{a}\left(1 - \epsilon\right)}{\sqrt{\epsilon\left(2 - \epsilon\right)}}\frac{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2} - 1}{\frac{r^2}{a^2} - 1 - \epsilon^2 + 2\epsilon} = \frac{2\frac{r}{a}\left(1 - \epsilon\right)}{\sqrt{\epsilon\left(2 - \epsilon\right)}}. \end{align} \]

So you have

\[ \begin{align} g'\left(\epsilon\right) &= \frac{2r\left(1 - \epsilon\right)}{a\sqrt{\epsilon\left(2 - \epsilon\right)}} - 2\left(1 - \epsilon\right)\left[\arcsin\left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right],\\ - \frac{2\frac{r}{a}\left(1 - \epsilon\right)}{\sqrt{\epsilon\left(2 - \epsilon\right)}} &= -2\left(1 - \epsilon\right)\left[\arcsin\left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right]. \end{align} \]

Still applies

\[ \begin{align} h'\left(\epsilon\right) &= 2\frac{3}{2}\sqrt{\epsilon\left(2 - \epsilon\right)}\left(2 - 2\epsilon\right) = 6\sqrt{\epsilon\left(2 - \epsilon\right)}\left(1 - \epsilon\right). \end{align} \]

So it applies

\[ \begin{align} \frac{g'\left(\epsilon\right)}{h'\left(\epsilon\right)} &= -\frac{1}{3\sqrt{\epsilon\left(2 - \epsilon\right)}}\left[\arcsin \left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right]. \end{align} \]

Here you have to apply L'Hospital's rule again, it applies

\[ \begin{align} & \frac{d}{d\epsilon}\left[\arcsin\left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right] = \frac{1}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}\frac{2\frac{r}{a}\left(1 - \epsilon\right)}{\sqrt{\epsilon\left(2 - \epsilon\right)}} \end{align} \]

as well as

\[ \begin{align} \frac{d}{d\epsilon}\sqrt{\epsilon\left(2 - \epsilon\right)} = \frac{1 - \epsilon}{\sqrt{\epsilon\left(2 - \epsilon\right)}}. \end{align} \]

This follows

\[ \begin{align} \lim\limits_{\epsilon\to 0}\frac{g'\left(\epsilon\right)}{h'\left(\epsilon\right)} = -\frac{1}{3}\frac{2\frac{r}{a}}{\frac{r^2}{a^2}} = - \frac{2a}{3r}. \end{align} \]

and thus

\[ \begin{align} \lim\limits_{\epsilon\to 0}F\left(\epsilon\right) = -\frac{2a}{3r}. \end{align} \]

In the case of a sphere, the correct limiting case results. For the derivation of $F$ applies

\[ \begin{align} F'\left(\epsilon\right) &= \frac{g'\left(\epsilon\right)h\left(\epsilon\right) - g\left(\epsilon\right)h'\left(\epsilon\right)}{h\left(\epsilon\right)^2} = \frac{ - 2\left(1 - \epsilon\right)\left[\arcsin\left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right]2\left(\epsilon\left(2 - \epsilon\right)\right)^{3/2}}{4\left(\epsilon\left(2 - \epsilon\right)\right)^3}\nonumber\\ & - \frac{\left(2\frac{r}{a}\left(\epsilon\left(2 - \epsilon\right)\right)^{1/2} - \left(\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right)\left[\arcsin\left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right]\right)6\left(1 - \epsilon\right)\sqrt{\epsilon\left(2 - \epsilon\right)}}{4\left(\epsilon\left(2 - \epsilon\right)\right)^3}\nonumber\\ &= \frac{\left(1 - \epsilon\right)\left[\arcsin\left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right]}{2\left(\epsilon\left(2 - \epsilon\right)\right)^{3/2}} + \frac{3r^2}{2a^2}\frac{\left(1 - \epsilon\right)\left[\arcsin\left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right]}{\left(\epsilon\left(2 - \epsilon\right)\right)^{5/2}}\nonumber\\ & - 3\frac{r}{a}\frac{1 - \epsilon}{\left(\epsilon\left(2 - \epsilon\right)\right)^2}. \end{align} \]

Due to the calculation rules for function limits, $\lim\limits_{\epsilon\to 0}F'\left(\epsilon\right) = \frac{1}{4\sqrt{2}}\lim\limits_{\epsilon\to 0}G\left(\epsilon\right)$ with

\[ \begin{align} G\left(\epsilon\right)&\coloneqq\frac{\left(2 - \epsilon\right)\left[\arcsin\left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right]}{2\epsilon^{3/2}} + \frac{3r^2}{2a^2}\frac{\arcsin\left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)}{\epsilon^{5/2}} - 3\frac{r\sqrt{2 - \epsilon}}{a\epsilon^2}\nonumber\\ &= \frac{\left(a^2\epsilon\left(2 - \epsilon\right) + 3r^2\right)\left[\arcsin\left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right] - 6ra\sqrt{\epsilon\left(2 - \epsilon\right)}}{2a^2\epsilon^{5/2}}. \end{align} \]

It applies with L'Hospital's rule

\[ \begin{align} \lim_{\epsilon\to 0}G\left(\epsilon\right) &= \lim_{\epsilon\to 0}\frac{2a^2\left(1 - \epsilon\right)\left[\arcsin\left(f_1\left(\epsilon\right)\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right]}{5a^2\epsilon^{3/2}}\nonumber\\ & + \frac{\left(a^2\epsilon\left(2 - \epsilon\right) + 3r^2\right)\frac{1}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}\frac{2\frac{r}{a}\left(1 - \epsilon\right)}{\sqrt{\epsilon\left(2 - \epsilon\right)}}}{5a^2\epsilon^{3/2}} - \frac{6ra\frac{1 - \epsilon}{\sqrt{\epsilon\left(2 - \epsilon\right)}}}{5a^2\epsilon^{3/2}}\nonumber\\ &= \frac{2}{5a^2}\lim_{\epsilon\to 0}\frac{a^2\left[\arcsin f_1\left(\epsilon\right) + \arcsin\left(f_2\left(\epsilon\right)\right)\right]}{\epsilon^{3/2}} + \frac{\left(a^2\epsilon\left(2 - \epsilon\right) + 3r^2\right)\frac{1}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}\frac{\frac{r}{a}}{\sqrt{\epsilon\left(2 - \epsilon\right)}} - 3\frac{ra}{\sqrt{\epsilon\left(2 - \epsilon\right)}}}{\epsilon^{3/2}}\nonumber\\ &= \frac{4}{15a^2}\lim_{\epsilon\to 0}\frac{2a^2}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}\frac{2\frac{r}{a}\left(1 - \epsilon\right)}{\epsilon\sqrt{2 - \epsilon}} + 3\frac{ra\left(1 - \epsilon\right)}{\epsilon^2\left(2 - \epsilon\right)^{3/2}}\nonumber\\ & - \left(a^2\epsilon\left(2 - \epsilon\right) + 3r^2\right)\frac{2\left(1 - \epsilon\right)}{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^2}\frac{\frac{r}{a}}{\epsilon\sqrt{2 - \epsilon}} - \frac{a^2\epsilon\left(2 - \epsilon\right) + 3r^2}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}\frac{\frac{r}{a}\left(1 - \epsilon\right)}{\epsilon^2\left(2 - \epsilon\right)^{3/2}}. \end{align} \]

With

\[ \begin{align} & \frac{2a^2}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}\frac{2\frac{r}{a}\left(1 - \epsilon\right)}{\epsilon\sqrt{2 - \epsilon}} - \left(a^2\epsilon\left(2 - \epsilon\right) + 3r^2\right)\frac{2\left(1 - \epsilon\right)}{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^2}\frac{\frac{r}{a}}{\epsilon\sqrt{2 - \epsilon}}\nonumber\\ &= \frac{2\frac{r}{a}\left(1 - \epsilon\right)a^2}{\epsilon\sqrt{\left(2 - \epsilon\right)}}\frac{1}{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^2}\left[2\epsilon\left(2 - \epsilon\right) + 2\frac{r^2}{a^2} - \epsilon\left(2 - \epsilon\right) - \frac{3r^2}{a^2}\right]\nonumber\\ &= \frac{2\frac{r}{a}\left(1 - \epsilon\right)a^2}{\epsilon\sqrt{2 - \epsilon}}\frac{1}{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^2}\left[-\frac{r^2}{a^2} + \epsilon\left(2 - \epsilon\right)\right]\nonumber\\ &= -\frac{2r^3\left(1 - \epsilon\right)}{a\epsilon\sqrt{2 - \epsilon}}\frac{1}{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^2} + \frac{2\frac{r}{a}\left(1 - \epsilon\right)a^2\sqrt{2 - \epsilon}}{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^2} \end{align} \]

and

\[ \begin{align} & 3\frac{ra\left(1 - \epsilon\right)}{\epsilon^2\left(2 - \epsilon\right)^{3/2}} - \frac{a^2\epsilon\left(2 - \epsilon\right) + 3r^2}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}\frac{\frac{r}{a}\left(1 - \epsilon\right)}{\epsilon^2\left(2 - \epsilon\right)^{3/2}}\nonumber\\ &= ra\frac{1 - \epsilon}{\epsilon^2\left(2 - \epsilon\right)^{3/2}}\left[3 - \frac{\epsilon\left(2 - \epsilon\right) + 3\frac{r^2}{a^2}}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}}\right] = \frac{ra}{\epsilon^2\left(2 - \epsilon\right)^{3/2}}\frac{2\epsilon\left(2 - \epsilon\right)}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}} \end{align} \]

follows

\[ \begin{align} \lim_{\epsilon\to 0}G\left(\epsilon\right) &= \frac{4}{15a^2\sqrt{2}}\lim_{\epsilon\to 0} - \frac{1r^3}{a\epsilon\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^2} + \frac{ra}{\epsilon^2\left(2 - \epsilon\right)}\frac{2\epsilon\left(2 - \epsilon\right)}{\epsilon\left(2 - \epsilon\right) + r^2} + \frac{2\frac{r}{a}\left(1 - \epsilon\right)a^2\left(2 - \epsilon\right)}{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^2}\nonumber\\ &= \frac{4r}{15a^2\sqrt{2}}\lim_{\epsilon\to 0} - \frac{2r^2}{a\epsilon\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^2} + \frac{a}{\epsilon}\frac{2}{\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}} + \frac{2\frac{1}{a}\left(1 - \epsilon\right)a^2\left(2 - \epsilon\right)}{\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^2}\nonumber\\ &= \frac{4}{15r\sqrt{2}}\lim\limits_{\epsilon\to 0} - \frac{2r^2}{a\epsilon\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]} + \frac{2a}{\epsilon} + \frac{2r^2\left(1 - \epsilon\right)\left(2 - \epsilon\right)}{a\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^2}\nonumber\\ &= \frac{4}{15r\sqrt{2}}\lim\limits_{\epsilon\to 0}\frac{ - 2\frac{r^2}{a} + 2a\epsilon\left(2 - \epsilon\right) + 2\frac{r^2}{a}}{\epsilon\left(\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right)} + \frac{2r^2\left(1 - \epsilon\right)\left(2 - \epsilon\right)}{a\left[\epsilon\left(2 - \epsilon\right) + \frac{r^2}{a^2}\right]^2}\nonumber\\ &= \frac{4}{15r\sqrt{2}}\left[\frac{4a}{\frac{r^2}{a^2}} + \frac{4r^2}{a\frac{r^4}{a^4}}\right] = \frac{4}{15r\sqrt{2}}\frac{8a^3}{r^2} = \frac{32}{15\sqrt{2}}\frac{a^3}{r^3}. \end{align} \]

Thus follows

\[ \begin{align} \lim_{\epsilon\to 0} F'\left(\epsilon\right) = \frac{4}{15}\frac{a^3}{r^3}. \end{align} \]

Analogously, the higher order coefficients can be determined in $\epsilon$. So you have

\[ \begin{align} \phi_0\left(r\right) = \frac{3GM}{2a}\left[-\frac{2a}{3r} + \epsilon\frac{4a^3}{15r^3}\right] = -\frac{GM}{r} + \epsilon\frac{2GM}{5r^3}a^2, \end{align} \]

the terms $\mathcal{O}\left(\epsilon^2\right)$ were neglected. It therefore applies in terms of Eq. (D.67)

\[ \begin{align} b_1 &= 0,\\ b_2\left(\epsilon\right) &= \epsilon\frac{2GMa^2}{5}. \end{align} \]

This therefore applies to the gravitational component of the Earth's gravity potential

\[ \begin{align} \phi_0\left(r, \chi\right) &= -\frac{GM}{r} + \frac{b_2}{2r^3}\left(3\sin^2\left(\chi\right) - 1\right) = -\frac{GM}{r} + \epsilon\frac{GMa^2}{5r^3}\left(3\sin^2\left(\chi\right) - 1\right). \end{align} \]

Now the centrifugal force is taken into account. Define

\[ \begin{align} m \coloneqq \frac{\omega^2a^3}{GM}. \end{align} \]

Then the centrifugal component $\phi_z$ is written using Eq. (D.34) as

\[ \begin{align} \phi_z = -\frac{1}{2}\omega^2r^2\cos^2\left(\chi\right) = -\frac{1}{2}\frac{GMm}{a^3}r^2\cos^2\left(\chi\right) = \frac{GMm}{2a^3}r^2\left(\sin^2\left(\chi\right) - 1\right). \end{align} \]

For the gravity potential one obtains

\[ \begin{align} \phi_g\left(r, \chi\right) &= \phi_0\left(r, \chi\right) + \phi_z\left(r, \chi\right) = - \frac{GM}{r} + \epsilon\frac{GMa^2}{5r^3}\left(3\sin^2\left(\chi\right) - 1\right) + \frac{GMm}{2a^3}r^2\left(\sin^2\left(\chi\right) - 1\right)\nonumber\\ &= \frac{GM}{r}\left[\left\{\frac{3a^2\epsilon}{5r^2} + \frac{mr^3}{2a^3}\right\}\sin^2\left(\chi\right) - 1\right] - \frac{GM}{a}\left\{\frac{a^3}{5r^3}\epsilon + \frac{mr^2}{2a^2}\right\}. \end{align} \]

An additional condition is now required that the surface of the ellipsoid of revolution is an equipotential surface, i.e

\[ \begin{align} \phi_g\left(r, \chi\right) = \text{const}. \end{align} \]

for

\[ \begin{align} \frac{x^2 + y^2}{a^2} + \frac{z^2}{c^2} = 1.\tag{D.142}\label{eq:allg_ellipse} \end{align} \]

A necessary condition for this is that the geopotential at the equator corresponds to that at the pole, i.e

\[ \begin{align} - \frac{GM}{a} - \frac{GM}{a}\left[\frac{\epsilon}{5} + \frac{m}{2}\right] &= \frac{GM}{c}\left[\frac{3a^2\epsilon}{5c^2} + \frac{mc^3}{2a^3} - 1\right] - \frac{GM}{a}\left[\frac{\epsilon a^3}{5c^3} + \frac{mc^2}{2a^2}\right]\nonumber\\ \Leftrightarrow 1 + \frac{\epsilon}{5} + \frac{m}{2} &= -\frac{3a^3}{5c^3}\epsilon - \frac{mc^2}{2a^2} + \frac{a}{c} + \frac{\epsilon a^3}{5c^3} + \frac{mc^2}{2a^2}\nonumber\\ \Leftrightarrow 1 + \frac{\epsilon}{5} + \frac{m}{2} &= -\frac{2a^3}{5c^3}\epsilon + \frac{1}{1 - \epsilon}\nonumber\\ \Leftrightarrow - \epsilon + \left(1 - \epsilon\right)\left(\frac{\epsilon}{5} + \frac{m}{2}\right) &= -\left(1 - \epsilon\right)\frac{2a^3}{5c^3}\epsilon = -\epsilon\left(1 - \epsilon\right)\frac{2}{5}\left(1 - \epsilon\right)^3. \end{align} \]

If you neglect all the terms $\mathcal{O}\left(\epsilon^2\right)$ and also $\mathcal{O}\left(\epsilon m\right)$, you get

\[ \begin{align} - \epsilon + \frac{\epsilon}{5} + \frac{m}{2} &= -\epsilon\frac{2}{5}\Leftrightarrow \frac{m}{2} = \epsilon\frac{2}{5}. \end{align} \]

So it applies in the equilibrium case

\[ \begin{align} \epsilon = \epsilon\left(m\right) = \frac{5m}{4}. \end{align} \]

Given a semi-major axis, the semi-minor axis $c$ is fixed:

\[ \begin{align} c = a\left(1 - \epsilon\left(m\right)\right) = a\left(1 - \frac{5m}{4}\right). \end{align} \]

You use it

\[ \begin{align} \frac{3\epsilon}{5} = \epsilon - \frac{2}{5}\epsilon = \epsilon - \frac{m}{2} = \frac{2\epsilon - m}{2} \end{align} \]

to simplify the expression for the gravity potential:

This is the formula that will be used from now on. Since it is first order in $\epsilon$ and $m$, all subsequent approximations are also only made in first order in $\epsilon$ and $m$.

D.3.3 Requirements for non-spherical coordinates

For the development of a meteorological coordinate system, it is practical if the surfaces with the same vertical coordinate are equipotential surfaces in the gravity field (see section 13.2). It is also expected to be advantageous for these surfaces to be ellipsoidal. For the development of a non-spherical coordinate system, it is also important to be able to specify the trajectories perpendicular to all potential surfaces in terms of the basic functions. In addition, the geopotential approximation used must be as consistent as possible with real physics and at the same time can be treated analytically. This means that the correct Clairaut ratio results again, and that $c = a\left(1 - \epsilon\right)$.

The requirements for an approximation of the gravity field are:

1. Possibility of analytically specifying the vertical trajectories
2. ellipsoidal contour surfaces
3. correct Clairaut ratio
4. correct semi-axis ratio

D.3.4 Ellipsoidal approximation of the gravity field

The first question that arises is whether all contour surfaces of Eq. (D.148) are ellipsoids. For the surface of the Earth (without orography) this has already been ensured to the first order in $\epsilon$ and $m$ by $\epsilon = \frac{5}{4}m$. With $r = \sqrt{x^2 + y^2 + z^2}$ one can derive an expression for the distance $r$ from the center of an ellipsoid of revolution with a semi-major axis $A$ and a semi-minor axis $C$ as a function of the geocentric angle $\chi$ at a given eccentricity:

\[ \begin{align} 1 = \frac{r^2\cos^2\left(\chi\right)}{A^2} + \frac{r^2\sin^2\left(\chi\right)}{C^2} &= \frac{r^2}{A^2}\left(\cos^2\left(\chi\right) + A^2\frac{\sin^2\left(\chi\right)}{C^2}\right)\nonumber\\ \Leftrightarrow r\left(\chi\right) &= A\left[\cos^2\left(\chi\right) + A^2\frac{\sin^2\left(\chi\right)}{C^2}\right]^{-1/2}\tag{D.149}\label{eq:defellipse_1} \end{align} \]

With $A = a$, $C = c = a\left(1 - \epsilon\right)$ it follows

Eq. (D.148) is not of this form: If you keep the left side of the equation constant and multiply by $r^3$, you get a fifth degree polynomial in $r$, which can generally only be solved numerically.

So one has to use an ellipsoidal approximation of Eq. Find (D.148). The following discussion is based on [34]. First of all, make it clear how you can designate a specific geopotential surface. The equatorial radius $R$ of the surface is suitable for this. For him applies

\[ \begin{align} \phi_R = -\frac{GM}{R}\left[1 + \left(\frac{2\epsilon - m}{6}\right)\frac{a^2}{R^2} + \frac{m}{2}\frac{R^3}{a^3}\right]. \end{align} \]

If you now want to specify this area more precisely, i.e. specify its distance $r = r\left(\chi\right)$ from the center of the earth at a geocentric latitude $\phi$, you simply write

\[ \begin{align} \frac{1}{R}\left[1 + \left(\frac{2\epsilon - m}{6}\right)\frac{a^2}{R^2} + \frac{m}{2}\frac{R^3}{a^3}\right] &= \frac{1}{r}\left[1 - \left\{\left(\frac{2\epsilon - m}{2}\right)\frac{a^2}{r^2} + \frac{m}{2}\frac{r^3}{a^3}\right\}\sin^2\left(\chi\right)\right]\nonumber\\ & + \frac{1}{a}\left\{\left(\frac{2\epsilon - m}{6}\right)\frac{a^3}{r^3} + \frac{m}{2}\frac{r^2}{a^2}\right\}. \end{align} \]

Now one makes the assumption that this applies

\[ \begin{align} \frac{r}{R} = 1 + \frac{r - R}{R} = 1 + \mathcal{O}\left(\epsilon\right), \end{align} \]

It is valid

\[ \begin{align} \frac{r - R}{R} \lesssim \frac{c - a}{a} = -\epsilon, \end{align} \]

at least within the weather-relevant atmosphere. You can now write first

\[ \begin{align} \phi_g\left(r, \chi\right) &= -\frac{GM}{r} + \frac{GM}{R}\frac{1}{r/R}\left[\left(\frac{2\epsilon - m}{2}\right)\frac{a^2}{R^2}\frac{R^2}{r^2} + \frac{m}{2}\frac{R^3}{a^3}\frac{r^3}{R^3}\right]\sin^2\left(\chi\right)\nonumber\\ & - \frac{GM}{a}\left[\left(\frac{2\epsilon - m}{6}\right)\frac{a^3}{R^3}\frac{R^3}{r^3} + \frac{m}{2}\frac{R^2}{a^2}\frac{r^2}{R^2}\right]\nonumber\\ &= -\frac{GM}{r} + \frac{GM}{R}\left(1 + \frac{r - R}{R}\right)^{-1}\bigg[\left(\frac{2\epsilon - m}{2}\right)\frac{a^2}{R^2}\left(1 + \frac{r - R}{R}\right)^{-2}\nonumber\\ & + \frac{m}{2}\frac{R^3}{a^3}\left(1 + \frac{r - R}{R}\right)^{3}\bigg]\sin^2\left(\chi\right) - \frac{GM}{a}\bigg[\left(\frac{2\epsilon - m}{6}\right)\frac{a^3}{R^3}\left(1 + \frac{r - R}{R}\right)^{-3}\nonumber\\ & + \frac{m}{2}\frac{R^2}{a^2}\left(1 + \frac{r - R}{R}\right)^{2}\bigg].\tag{D.155}\label{eq:schwere_approx_1} \end{align} \]

By neglecting terms $\mathcal{O}\left(\epsilon^2, m^2, \epsilon m\right)$ this becomes initially the case

\[ \begin{align} \phi_g\left(\chi, R\right) &= -\frac{GM}{r\left(\chi, R\right)} + \frac{GM}{R}\left[1 + \left\{\left(\frac{2\epsilon - m}{2}\right)\frac{a^2}{R^2} + \frac{m}{2}\frac{R^3}{a^3}\right\}\sin^2\left(\chi\right)\right]\nonumber\\ & - \frac{GM}{R}\left[1 + \left(\frac{2\epsilon - m}{6}\right)\frac{a^2}{R^2} + \frac{m}{2}\frac{R^3}{a^3}\right] + \mathcal{O}\left(\epsilon^2, m^2, \epsilon m\right) \end{align} \]

and then continue to

\[ \begin{align} \phi_g\left(\chi, R\right) &\approx -\frac{GM}{r\left(\chi, R\right)} + \frac{GM}{R}\left[1 + \left\{\left(2\epsilon - m\right)\frac{a^2}{R^2} + m\frac{R^3}{a^3}\right\}\sin^2\left(\chi\right)\right]^{1/2}\nonumber\\ & - \frac{GM}{R}\left[1 + \left(\frac{2\epsilon - m}{6}\right)\frac{a^2}{R^2} + \frac{m}{2}\frac{R^3}{a^3}\right]. \end{align} \]

The contour surfaces of this are already ellipsoids. To show this, determine the geocentric distance $r$ as a function of latitude given the equatorial radius $R$:

\[ \begin{align} & -\frac{GM}{R} + \frac{GM}{R} - \frac{GM}{R} - \frac{GM}{R}\left[\left(\frac{2\epsilon - m}{6}\right)\frac{a^2}{R^2} + \frac{m}{2}\frac{R^3}{a^3}\right]\nonumber\\ &= -\frac{GM}{r} + \frac{GM}{R}\left[1 + \left\{\left(2\epsilon - m\right)\frac{a^2}{R^2} + m\frac{R^3}{a^3}\right\}\sin^2\left(\chi\right)\right]^{1/2} - \frac{GM}{R} - \frac{GM}{R}\left[\left(\frac{2\epsilon - m}{6}\right)\frac{a^2}{R^2} + \frac{m}{2}\frac{R^3}{a^3}\right]\nonumber\\ \Leftrightarrow 1 &= \frac{R}{r} - \left[1 + \sin^2\left(\chi\right)\left\{\left(2\epsilon - m\right)\frac{a^2}{R^2} + m\frac{R^3}{a^3}\right\}\right]^{1/2} + 1\nonumber\\ \Leftrightarrow r &= R\frac{1}{\left[1 + \sin^2\left(\chi\right)\left\{\left(2\epsilon - m\right)\frac{a^2}{R^2} + m\frac{R^3}{a^3}\right\}\right]^{1/2}}\nonumber\\ \Leftrightarrow r\left(R, \chi\right) &= R\left[\cos^2\left(\chi\right) + \sin^2\left(\chi\right)\left\{\left(2\epsilon - m\right)\frac{a^2}{R^2} + m\frac{R^3}{a^3} + 1\right\}\right]^{-1/2} \end{align} \]

If one sets $R = A$ and $C = A\left\{\left(2\epsilon - m\right)\frac{a^2}{R^2} + m\frac{R^3}{a^3} + 1\right\}^{-1/2}$ in Eq. (D.149) you get exactly this result. Unfortunately, the trajectories orthogonal to the contour sets of this potential cannot be determined analytically, so Eq. (D.155) needs to be further approximated. To do this, the terms $\frac{a^2}{R^2}$ and $\frac{R^3}{a^3}$ are expanded by $R^2 = a^2$: This is justified because of $\frac{R - a}{a}\lesssim 10^{-2}$.

\[ \begin{align} \frac{a^2}{R^2} &\approx 1 - \frac{1}{a^2}\left(R^2 - a^2\right) = 2 - \frac{R^2}{a^2},\\ \frac{R^3}{a^3} = \frac{\left(R^2\right)^{3/2}}{a^3} &\approx 1 + \frac{3}{2}\frac{R^2 - a^2}{a^2} = -\frac{1}{2} + \frac{3}{2}\frac{R^2}{a^2} \end{align} \]

Putting this into Eq. (D.155), you get

\[ \begin{align} \phi_g\left(\chi, R\right) &\approx -\frac{GM}{r\left(\chi, R\right)} + \frac{GM}{R}\left[1 + \left\{\left(\frac{8\epsilon - 5m}{2}\right) + \left(\frac{5m - 4\epsilon}{2}\right)\frac{R^2}{a^2}\right\}\sin^2\left(\chi\right)\right]^{1/2}\nonumber\\ & - \frac{GM}{R}\left[1 + \left(\frac{8\epsilon - 7m}{12}\right) + \left(\frac{11m - 4\epsilon}{12}\right)\frac{R^2}{a^2}\right]. \end{align} \]

This leads to

\[ \begin{align} r\left(\chi, R\right) = R\left[\cos^2\left(\chi\right) + \sin^2\left(\chi\right)\left\{\left(\frac{8\epsilon - 5m}{2}\right) + \left(\frac{5m - 4\epsilon}{2}\right)\frac{R^2}{a^2} + 1\right\}\right]^{-1/2}, \tag{D.162}\label{eq:schwere_approx_2_parameter} \end{align} \]

with Eq. (D.149) will do this

\[ \begin{align} C = R\left\{1 + \left(\frac{8\epsilon - 5m}{2}\right) + \left(\frac{5m - 4\epsilon}{2}\right)\frac{R^2}{a^2}\right\}^{-1/2}. \end{align} \]

However, a self-consistency problem arises here. If you set $R = a$ in Eq. (D.162), you get

\[ \begin{align} r\left(\chi, a\right) = a\left[1 + 2\epsilon\sin^2\left(\chi\right)\right]^{-1/2}. \end{align} \]

This expression describes the earth's surface without orography. You receive

\[ \begin{align} r\left(0, a\right) = a, & {} & r\left(\pi/2, a\right) = a\frac{1}{\sqrt{1 + 2\epsilon}}. \end{align} \]

With

\[ \begin{align} \frac{d}{d\epsilon}\frac{1}{\sqrt{1 + 2\epsilon}} = -\frac{1}{\left(1 + 2\epsilon\right)^{3/2}} \end{align} \]

follows

\[ \begin{align} r\left(\pi/2\right) = a - a\epsilon + \mathcal{O}\left(\epsilon^2\right)\approx a\left(1 - \epsilon\right), \end{align} \]

So not the desired exact result $r\left(\pi/2\right) = a\left(1 - \epsilon\right)$. To solve this problem, define

\[ \begin{align} \newtilde{\epsilon} \coloneqq \frac{a^2 - c^2}{2c^2} = \frac{1 - \left(1 - \epsilon\right)^2}{2\left(1 - \epsilon\right)^2}. \end{align} \]

Define the auxiliary function $f\left(\epsilon\right)$ by

\[ \begin{align} f\left(\epsilon\right) \coloneqq \frac{1 - \left(1 - \epsilon\right)^2}{2\left(1 - \epsilon\right)^2}, \end{align} \]

then applies

\[ \begin{align} f'\left(\epsilon\right) = \frac{4\left(1 - \epsilon\right)^3 - \left[1 - \left(1 - \epsilon\right)^2\right]2\left(1 - \epsilon\right)^2}{4\left(1 - \epsilon\right)^4} \Rightarrow f'\left(0\right) = 1. \end{align} \]

So is

\[ \begin{align} \newtilde{\epsilon} = \epsilon + \mathcal{O}\left(\epsilon^2\right). \end{align} \]

So in first order you can replace $\epsilon$ with $\newtilde{\epsilon}$ everywhere. You receive

\[ \begin{align} \phi_g\left(\chi, R\right) &\approx -\frac{GM}{r\left(\chi, R\right)} + \frac{GM}{R}\left[1 + \left\{\left(\frac{8\newtilde{\epsilon} - 5m}{2}\right) + \left(\frac{5m - 4\newtilde{\epsilon}}{2}\right)\frac{R^2}{a^2}\right\}\sin^2\left(\chi\right)\right]^{1/2}\nonumber\\ & - \frac{GM}{R}\left[1 + \left(\frac{8\newtilde{\epsilon} - 7m}{12}\right) + \left(\frac{11m - 4\newtilde{\epsilon}}{12}\right)\frac{R^2}{a^2}\right].\tag{D.172}\label{eq:schwere_approx} \end{align} \]

With this formula you also get ellipsoidal contour surfaces,

and also

\[ \begin{align} r\left(\chi, a\right) &= a\left[1 + 2\newtilde{\epsilon}\sin^2\left(\chi\right)\right]^{-1/2} = a\left(1 - \epsilon\right)\left[\left(1 - \epsilon\right)^2 + \left(1 - \left(1 - \epsilon\right)^2\right)\sin^2\left(\chi\right)\right]^{-1/2}. \end{align} \]

This is the desired result, see Eq. (D.150).

The Clairaut's theorem is a statement about the connection between eccentricity $\epsilon$, angular velocity $\omega$ and the ratio of gravity $g_a, g_p$ at the equator and pole. At the pole, the geopotential is $\phi_g = \phi_g^{(p)}$ with Eq. (D.148)

\[ \begin{align} \phi_g^{(p)} = -\frac{GM}{r} + \frac{GM}{r^3}a^2\left(\frac{2\epsilon - m}{3}\right), \end{align} \]

applies at the equator

\[ \begin{align} \phi_g^{(a)} = -\frac{GM}{r} - \frac{GM}{a}\left\{\left(\frac{2\epsilon - m}{6}\right)\frac{a^3}{r^3} + \frac{m}{2}\frac{r^2}{a^2}\right\}. \end{align} \]

Therefore apply

\[ \begin{align} \frac{d\phi_g^{(a)}}{dr} &= \frac{GM}{r^2} - \frac{GM}{a}\left\{ - 3\left(\frac{2\epsilon - m}{6}\right)\frac{a^3}{r^4} + m\frac{r}{a^2}\right\},\\ \frac{d\phi_g^{(p)}}{dr} &= \frac{GM}{r^2} - 3\frac{GM}{r^4}a^2\left(\frac{2\epsilon - m}{3}\right). \end{align} \]

This results in $r = a$ or $r = c = a\left(1 - \epsilon\right)$

\[ \begin{align} g_a &= \frac{GM}{a^2}\left[1 + \epsilon - 3\frac{m}{2}\right],\\ g_p &= \frac{GM}{a^2\left(1 - \epsilon\right)^2}\left[1 - \frac{2\epsilon}{\left(1 - \epsilon\right)^2} + \frac{m}{\left(1 - \epsilon\right)^2}\right]. \end{align} \]

It follows

\[ \begin{align} \frac{g_p - g_a}{g_a} &= \frac{\frac{1}{\left(1 - \epsilon\right)^2}\left[1 - \frac{2}{\left(1 - \epsilon\right)^2}\epsilon + \frac{m}{\left(1 - \epsilon\right)^2}\right] - 1 - \epsilon + 3\frac{m}{2}}{1 + \epsilon - 3\frac{m}{2}}\nonumber\\ &\approx \frac{1 - \frac{2}{\left(1 - \epsilon\right)^2}\epsilon + \frac{m}{\left(1 - \epsilon\right)^2} - 1 - \epsilon + 3\frac{m}{2} + 2\epsilon}{1 - \epsilon - 3\frac{m}{2}}, \end{align} \]

where again terms $\mathcal{O}\left(\epsilon^2, \epsilon m, m^2\right)$ were neglected. With

\[ \begin{align} \frac{d}{d\epsilon}\frac{1}{\left(1 - \epsilon\right)^2} = \frac{2\left(1 - \epsilon\right)}{\left(1 - \epsilon\right)^4} = \frac{2}{\left(1 - \epsilon\right)^3} \end{align} \]

follows

\[ \begin{align} \frac{1}{\left(1 - \epsilon\right)^2}\approx 1 + 2\epsilon \end{align} \]

and thus

\[ \begin{align} \frac{g_p - g_a}{g_a} &\approx&\frac{1 - 2\epsilon + m - 1 - \epsilon + 3\frac{m}{2} + 2\epsilon}{1 - \epsilon - \frac{3m}{2}} = \frac{\frac{5m}{2} - \epsilon}{1 - \epsilon - \frac{3m}{2}}. \end{align} \]

If you insert $\epsilon\approx\frac{5m}{4}$ here, you get

\[ \begin{align} \frac{g_p - g_a}{g_a} \approx m\frac{\frac{5}{4}}{1 - m\left(\frac{5}{4} - \frac{3}{2}\right)} = \frac{m}{4}\frac{5}{1 + \frac{m}{4}} = \frac{5m}{4 + m}. \end{align} \]

Define $f\left(m\right) \coloneqq \frac{5m}{4 + m}$, then holds

\[ \begin{align} \frac{df}{dm} = 5\frac{4 + m - m}{\left(4 + m\right)^2}, \end{align} \]

so

\[ \begin{align} \frac{df}{dm}\left(m = 0\right) = \frac{5}{4}. \end{align} \]

This follows

\[ \begin{align} \frac{g_p - g_a}{g_a} \approx 5\frac{m}{4}.\tag{D.188}\label{eq:clairaut_deriv_1} \end{align} \]

So that applies

This is Clairaut's theorem. This should now be done for Eq. (D.172) can be recalculated. All that needs to be shown is that Eq. (D.188) still applies. At the equator applies

\[ \begin{align} \phi_g^{(a)} = -\frac{GM}{r} - \frac{GMr}{a^2}\left[\left(\frac{8\newtilde{\epsilon} - 7m}{12}\right) + \left(\frac{11m - 4\newtilde{\epsilon}}{12}\right)\right]. \end{align} \]

It follows

\[ \begin{align} g_a = \frac{GM}{a^2}\left[1 + \left(\frac{8\newtilde{\epsilon} - 7m}{12}\right) + \left(\frac{11m - 4\newtilde{\epsilon}}{12}\right)\right]. \end{align} \]

At the pole applies

\[ \begin{align} \phi_g\left(\pi/2, R\right) &\approx -\frac{GM}{r\left(\pi/2, R\right)} + \frac{GM}{R}\left[1 + \left\{\left(\frac{8\newtilde{\epsilon} - 5m}{2}\right) + \left(\frac{5m - 4\newtilde{\epsilon}}{2}\right)\frac{R^2}{a^2}\right\}\right]^{1/2}\nonumber\\ & -\frac{GM}{R}\left[1 + \left(\frac{8\newtilde{\epsilon} - 7m}{12}\right) + \left(\frac{11m - 4\newtilde{\epsilon}}{12}\right)\frac{R^2}{a^2}\right]. \end{align} \]

It follows

\[ \begin{align} g_p &= \frac{d\phi_g}{dr} = \frac{\partial R}{\partial r}\frac{d\phi_g}{dR} = \frac{1}{\partial r/\partial R}\frac{d\phi_g}{dR} = \frac{1}{\partial r/\partial R}g_a. \end{align} \]

For $dr/dR$ it follows from Eq. (D.173)

\[ \begin{align} \frac{\partial r}{\partial R} &= \left[1 + \left\{\left(\frac{8\newtilde{\epsilon} - 5m}{2}\right) + \left(\frac{5m - 4\newtilde{\epsilon}}{2}\right)\frac{R^2}{a^2}\right\}\sin^2\left(\chi\right)\right]^{-1/2}\nonumber\\ & - \frac{R^2}{a^2}\left(\frac{5m - 4\newtilde{\epsilon}}{2}\right)\left[1 + \left\{\left(\frac{8\newtilde{\epsilon} - 5m}{2}\right) + \left(\frac{5m - 4\newtilde{\epsilon}}{2}\right)\frac{R^2}{a^2}\right\}\sin^2\left(\chi\right)\right]^{-3/2}. \end{align} \]

At the point $\chi = \pi/2$, $R = a$ this becomes

\[ \begin{align} \frac{\partial r}{\partial R}\left(R = a, \chi = \pi/2\right) &= \left[1 + 2\newtilde{\epsilon}\right]^{-1/2} - \left(\frac{5m - 4\newtilde{\epsilon}}{2}\right)\left[1 + 2\newtilde{\epsilon}\right]^{-3/2}\nonumber\\ &= \left(1 + 2\newtilde{\epsilon}\right)^{-3/2}\left[1 + 4\newtilde{\epsilon} - 5\frac{m}{2}\right]. \end{align} \]

So you get

\[ \begin{align} g_p = \frac{GM}{a^2}\left[1 + \left(\frac{8\newtilde{\epsilon} - 7m}{12}\right) + \left(\frac{11m - 4\newtilde{\epsilon}}{12}\right)\right]\left(1 + 2\newtilde{\epsilon}\right)^{3/2}\frac{1}{1 + 4\newtilde{\epsilon} - 5\frac{m}{2}}. \end{align} \]

This leads to

\[ \begin{align} \frac{g_p - g_a}{g_a} &= \frac{\partial R}{\partial r} - 1 = \frac{\left(1 + 2\newtilde{\epsilon}\right)^{3/2}}{1 + 4\newtilde{\epsilon} - 5\frac{m}{2}} - 1. \end{align} \]

Define

\[ \begin{align} f\left(x, y\right) \coloneqq \frac{\left(1 + 2x\right)^{3/2}}{1 + 4x - 5\frac{y}{2}}, \end{align} \]

then $f\left(0\right) = 1$ and

\[ \begin{align} \frac{\partial f}{\partial x} = \frac{3\sqrt{1 + 2x}\left(1 + 4x - 5\frac{y}{2}\right) - \left(1 + 2x\right)^{3/2}4}{\left(1 + 4x - 5\frac{y}{2}\right)^2}, \end{align} \]

so

\[ \begin{align} \frac{\partial f}{\partial x}\left(0\right) = -1, \end{align} \]

as well as

\[ \begin{align} \frac{\partial f}{\partial y} = \frac{5}{2}\left(1 + 2x\right)^{3/2}\frac{1}{\left(1 + 4x - 5\frac{y}{2}\right)^2}, \end{align} \]

so

\[ \begin{align} \frac{\partial f}{\partial y}\left(0\right) = \frac{5}{2}. \end{align} \]

This follows

so with Eq. (D.172) also satisfies Clairaut's theorem. As a note, it should be said: Eq. (D.172) is first-order exact in $\epsilon$ and $m$ and therefore no better or worse than Eq. (D.148).

D.3.5 Vertical trajectories

Since horizontal derivatives in a suitable geophysical coordinate system refer to equipotential surfaces in the gravity field, vertical derivatives (if the KS is orthogonal) are perpendicular to this and therefore aligned in the direction of the gravity field. So it's important to know this direction. For this purpose, the trajectories perpendicular to the equipotential surfaces are determined in this section.

Eq. (D.173) describes the equipotential surfaces in the gravity field, more precisely their distance $r$ from the center of the Earth for a given equatorial radius $R$ and geocentric latitude $\chi$. This equation is:

\[ \begin{align} r\left(\chi, R\right) = R\left[1 + \sin^2\left(\chi\right)\left\{\mu + \nu \frac{R^2}{a^2}\right\}\right]^{-1/2} \end{align} \]

This was done

\[ \begin{align} \mu \coloneqq \frac{8\newtilde{\epsilon} - 5m}{2}, \nu \coloneqq \frac{5m - 4\newtilde{\epsilon}}{2} \end{align} \]

defined. This can be transformed into:

\[ \begin{align} r\left(\chi, R\right) = \frac{RC\left(R\right)}{\sqrt{R^2\sin^2\left(\chi\right) + C^2\left(R\right)\cos^2\left(\chi\right)}}\tag{D.206}\label{eq:deriv_senkr_darst_pot_flaechen} \end{align} \]

with

\[ \begin{align} C\left(R\right) \coloneqq R\frac{1}{\sqrt{1 + \mu + \nu \frac{R^2}{a^2}}}. \end{align} \]

Fulfill ellipsoids

\[ \begin{align} \frac{r^2\cos^2\left(\chi\right)}{R^2} + \frac{r^2\sin^2\left(\chi\right)}{c^2} = 1 \end{align} \]

with semi-major axis $R$ and semi-minor axis $c$. If you put Eq. (D.206), you get

\[ \begin{align} \frac{1}{R^2}\frac{R^2C^2\cos^2\left(\chi\right)}{R^2\sin^2\left(\chi\right) + C^2\cos^2\left(\chi\right)} + \frac{1}{c^2}\frac{R^2C^2\sin^2\left(\chi\right)}{R^2\sin^2\left(\chi\right) + C^2\cos^2\left(\chi\right)}\hastobe1. \end{align} \]

One sees that this is true for $c = C$. $C$ is therefore the semi-minor axis of the equipotential surfaces.

Since an ellipsoid is rotationally symmetrical, it is sufficient to first determine the vertical trajectories in the meridian plane $\lambda = 0$. Here you can simply look at ellipses and fulfill them

\[ \begin{align} \frac{x^2}{R^2} + \frac{z^2}{C^2} = 1.\tag{D.210}\label{eq:deriv_senkr_ellipse_2d} \end{align} \]

If you insert the definition of $C$ here, it follows

\[ \begin{align} x^2 + \left(1 + \mu + \nu \frac{R^2}{a^2}\right)z^2 = R^2.\tag{D.211}\label{eq:great_senkr_traj_deriv_1} \end{align} \]

If you convert this to $R^2$, you get

\[ \begin{align} R^2 = \frac{x^2 + \left(1 + \mu\right)z^2}{1 - \nu\frac{z^2}{a^2}}.\tag{D.212}\label{eq:great_senkr_traj_deriv_5} \end{align} \]

Furthermore, it follows from this

\[ \begin{align} 1 + \mu + \nu\frac{x^2}{a^2} &= 1 + \mu + \frac{x^2}{z^2} - \frac{x^2}{z^2} + \nu\frac{x^2}{a^2} = R^2\frac{1}{z^2}\left(1 - \nu\frac{z^2}{a^2}\right) - \frac{x^2}{z^2}\left(1 - \nu\frac{z^2}{a^2}\right)\nonumber\\ &= \left(1 - \nu\frac{z^2}{a^2}\right)\frac{R^2 - x^2}{z^2} = \left(1 - \nu\frac{z^2}{a^2}\right)\left(1 + \mu + \nu\frac{R^2}{a^2}\right).\tag{D.213}\label{eq:great_senkr_traj_deriv_4} \end{align} \]

If $R$ is fixed, the coordinate $z$ of an ellipse in a half-space can be understood as a function of $x$. This function $z\left(x\right)$ satisfies Eq. (D.210). Differentiating this equation in terms of $x$ gives us

\[ \begin{align} \frac{dz}{dx} = -\frac{C^2}{R^2}\frac{x}{z}. \end{align} \]

This is a directional derivative along the potential surface. Geometrically it is clear that the local perpendicular to the potential surface applies

\[ \begin{align} \frac{dz}{dx} = \frac{R^2}{C^2}\frac{z}{x}. \end{align} \]

With the definition of $C$ follows

\[ \begin{align} \frac{R^2}{C^2} = 1 + \mu + \nu\frac{R^2}{a^2}, \end{align} \]

therefore applies

\[ \begin{align} \frac{dz}{dx} = \left(1 + \mu + \nu\frac{R^2}{a^2}\right)\frac{z}{x} = \frac{1 + \mu + \nu\frac{x^2}{a^2}}{1 - \nu\frac{z^2}{a^2}}\frac{z}{x}. \end{align} \]

This can be solved using the method of separating the variables presented in Section A.7. $z\left(x\right)$ is bijective, so the requirements are met. First you do the math

\[ \begin{align} \left(1 - \nu\frac{z^2}{a^2}\right)\frac{dz}{dx} &= \frac{z}{x}\left(1 + \mu + \nu \frac{x^2}{a^2}\right)\Leftrightarrow\left(\frac{1}{z} - \nu\frac{z}{a^2}\right)\frac{dz}{dx} = \frac{1 + \mu}{x} + \nu \frac{x}{a^2}. \end{align} \]

So follow

\[ \begin{align} \int_{}^{}\frac{1}{z} - \nu\frac{z}{a^2}dz &= \int_{}^{}\frac{1 + \mu}{x} + \nu \frac{x}{a^2}dx\nonumber\\ \Leftrightarrow\frac{1}{a}\int_{}^{}\frac{a}{z} - \nu\frac{z}{a}d\left(\frac{z}{a}\right) &= \frac{1}{a}\int_{}^{}\left(1 + \mu\right)\frac{a}{x} + \nu \frac{x}{a}d\left(\frac{x}{a}\right)\nonumber\\ \Leftrightarrow \ln\frac{z}{a} - \frac{\nu}{2}\frac{z^2}{a^2} &= D' + \left(1 + \mu\right)\ln\left(\frac{x}{a}\right) + \frac{\nu}{2a^2}x^2\nonumber\\ \Leftrightarrow \frac{z}{a} &= \exp\left( D' + \left(1 + \mu\right)\ln\left(\frac{x}{a}\right) + \frac{\nu}{2a^2}x^2 + \frac{\nu}{2}\frac{z^2}{a^2}\right)\nonumber\\ \Leftrightarrow z &= aD\left(\frac{x}{a}\right)^{1 + \mu}\exp\left[\frac{\nu}{2a^2}\left(x^2 + z^2\right)\right]\tag{D.219}\label{eq:great_senkr_traj_deriv_2} \end{align} \]

with an integration constant $D \coloneqq \exp\left(D'\right)>0$ with $D'\in \mathbb{R}$. For a given $x$, by solving this equation one can find the $z$ which lies on the desired orthogonal trajectory. Which trajectory this is is determined by the number $D$, so $D$ is a function of a width $\varphi$, it is $D = D\left(\varphi\right)$. The exact form of this function $D = D\left(\varphi\right)$ is by convention. Here it is determined by the intersection of the trajectory with the WGS 84. Let $\left(x_S, z_S\right)^T$ be this intersection. Then applies

\[ \begin{align} \left|z_S\right| &= aD\left(\frac{x_S}{a}\right)^{1 + \mu}\exp\left[\frac{\nu}{2a^2}\left(x_S^2 + z_S^2\right)\right]\nonumber\\ \Rightarrow D &= \frac{\left|z_S\right|}{a}\left(\frac{a}{x_S}\right)^{1 + \mu}\exp\left[-\frac{\nu}{2a^2}\left(x_S^2 + z_S^2\right)\right].\tag{D.220}\label{eq:great_senkr_traj_deriv_3} \end{align} \]

The following equation thus determines the vertical trajectories sought:

D.4 Great-Koordinaten

The Great coordinates are described in [34]. They are suitable for expressing differential operators in ellipsoidal geometries. They are also orthogonal.

First, all lengths are scaled to the semi-major axis of the Earth $a$:

\[ \begin{align} \frac{x}{a} \to x, & {} & \frac{y}{a} \to y, & {} & \frac{z}{a} \to z, & {} & \frac{R}{a} \to R\tag{D.222}\label{eq:great_skal} \end{align} \]

Eq. (D.221) is then

D.4.1 Vertical coordinate

In spherical coordinates, the radial coordinate $r$ denotes a spherical surface. In Great coordinates, the radial coordinate $R$ denotes a specific ellipsoid, $R$ has already been introduced as the equatorial radius of the corresponding ellipsoid. This is the vertical coordinate in the Great system.

D.4.2 Zonal coordinate

The zonal coordinate $\lambda$ of a location $\mathbf{r}$ is the angle by which one must rotate the xz plane about the z axis until $\mathbf{r}$ lies in this plane. This is analogous to the azimuth angle in ordinary spherical coordinates.

D.4.3 Meridional coordinate

To ensure consistency with measurement data, the geographical latitude $\varphi$ is chosen as the meridional coordinate. The geographical latitude of a point $\mathbf{r}$ outside the Earth can be found as follows:

• Find the vertical trajectory on which $\mathbf{r}$ lies.
• Follow this vertical trajectory to its intersection $\mathbf{r}_s$ with the WGS 84.
• Find the perpendicular direction $\mathbf{s}$ to the WGS 84 at the point $\mathbf{r}_s$ and form the straight line $\mathbf{g} = \mathbf{r}_s + \gamma\mathbf{s}$ with $\gamma\in \mathbb{R}$.
• Trace the line $\mathbf{g}$ to its intersection $\mathbf{r}_0$ with the xy-plane of the global coordinates. The angle that $\mathbf{r}_0$ makes with the xy plane is the latitude $\varphi$.

D.4.4 Metric tensor

In order to be able to express the governing equations in Great coordinates, the differential operators must be expressed by partial derivatives with respect to $\varphi, \lambda, R$. To do this, some preparatory work must first be done. So let a point $\mathbf{r} = \left(x, y, z\right)^T$ be given in global coordinates, $\mathbf{r}$ has the Great coordinates $\left(R, \varphi, \lambda\right)$. The transformation of all differential operators is based on the metric factors $h_\lambda, h_\varphi, h_R>0$, defined by

\[ \begin{align} h_\lambda^2 &= \left(\frac{\partial x}{\partial\lambda}\right)^2 + \left(\frac{\partial y}{\partial\lambda}\right)^2 + \left(\frac{\partial z}{\partial\lambda}\right)^2,\\ h_\varphi^2 &= \left(\frac{\partial x}{\partial\varphi}\right)^2 + \left(\frac{\partial y}{\partial\varphi}\right)^2 + \left(\frac{\partial z}{\partial\varphi}\right)^2,\\ h_R^2 &= \left(\frac{\partial x}{\partial R}\right)^2 + \left(\frac{\partial y}{\partial R}\right)^2 + \left(\frac{\partial z}{\partial R}\right)^2. \end{align} \]

To determine $h_R$ and $h_\varphi$, an observation in the meridional plane $\lambda = 0$ is sufficient. If you apply the rescaling equations (D.222) - (D.222) to the equations (D.211) and (D.219), you get

\[ \begin{align} x^2 + \left(1 + \mu + \nu R^2\right)z^2 &= R^2, \tag{D.227}\label{eq:great_deriv_1}\\ z &= D\left(\varphi\right)x^{1 + \mu}\exp\left[\frac{\nu}{2}\left(x^2 + z^2\right)\right].\tag{D.228}\label{eq:great_deriv_2} \end{align} \]

If you take this logarithm, you get

\[ \begin{align} \ln\left(z^2\right) &= \ln\left(R^2 - x^2\right) - \ln\left(1 + \mu + \nu R^2\right), \tag{D.229}\label{eq:great_deriv_3}\\ \ln\left(z\right) &= \ln\left[D\left(\varphi\right)\right] + \left(1 + \mu\right)\ln\left(x\right) + \frac{\nu}{2}\left(x^2 + z^2\right)\tag{D.230}\label{eq:great_deriv_4}. \end{align} \]

$D\left(\varphi\right)$ results from Eq. (D.220), which reads rescaled

\[ \begin{align} D\left(\varphi\right) = \frac{\left|z_S\left(\varphi\right)\right|}{x_S\left(\varphi\right)^{1 + \mu}}\exp\left[-\frac{\nu}{2}\left(x_S\left(\varphi\right)^2 + z_S\left(\varphi\right)^2\right)\right].\tag{D.231}\label{eq:great_deriv_5} \end{align} \]

Now one differentiates Equation (D.231) partially according to $R$:

\[ \begin{align} \frac{2}{z}\frac{\partial z}{\partial R} &= -\frac{\partial x}{\partial R}2x\frac{1}{R^2 - x^2} + 2R\frac{1}{R^2 - x^2} - 2\nu R\frac{1}{1 + \mu + \nu R^2}\nonumber\\ \Leftrightarrow \left(1 + \mu + \nu R^2\right)z\frac{\partial z}{\partial R} &= -\frac{\partial x}{\partial R}\left(1 + \mu + \nu R^2\right)\frac{xz^2}{R^2 - x^2} + \left(1 + \mu + \nu R^2\right)\frac{Rz^2}{R^2 - x^2} - \nu Rz^2. \end{align} \]

If you put Eq. (D.227), you get

\[ \begin{align} \left(1 + \mu + \nu R^2\right)z\frac{\partial z}{\partial R} + x\frac{\partial x}{\partial R} &= R\left(1 - \nu z^2\right).\tag{D.233}\label{eq:great_deriv_6} \end{align} \]

If you differentiate Eq. (D.230) partial to $R$, we get:

\[ \begin{align} \frac{1}{z}\frac{\partial z}{\partial R} &= \left(1 + \mu\right)\frac{1}{x}\frac{\partial x}{\partial R} + \nu\left(x\frac{\partial x}{\partial R} + z\frac{\partial z}{\partial R}\right)\nonumber\\ \Leftrightarrow \left(1 + \mu\right)\frac{1}{x}\frac{\partial x}{\partial R} - \frac{1}{z}\frac{\partial z}{\partial R} &= -\nu\left(x\frac{\partial x}{\partial R} + z\frac{\partial z}{\partial R}\right)\nonumber\\ \Leftrightarrow\left(1 + \mu + \nu R^2\right)\frac{1}{x}\frac{\partial x}{\partial R} - \frac{1}{z}\frac{\partial z}{\partial R} &= -\nu\left(x\frac{\partial x}{\partial R} + z\frac{\partial z}{\partial R}\right) + \nu R^2\frac{1}{x}\frac{\partial x}{\partial R}\nonumber\\ \Leftrightarrow\left(1 + \mu + \nu R^2\right)\frac{1}{x}\frac{\partial x}{\partial R} - \frac{1}{z}\frac{\partial z}{\partial R} &= \frac{\nu}{x}\frac{\partial x}{\partial R}\left(R^2 - x^2\right) - \nu z\frac{\partial z}{\partial R}. \end{align} \]

If you put Eq. (D.227), you get

\[ \begin{align} \Leftrightarrow\left(1 + \mu + \nu R^2\right)\frac{1}{x}\frac{\partial x}{\partial R} - \frac{1}{z}\frac{\partial z}{\partial R} &= \frac{\nu}{x}\frac{\partial x}{\partial R}\left(1 + \mu + \nu R^2\right)z^2 - \nu z\frac{\partial z}{\partial R}\nonumber\\ \Leftrightarrow\left(1 + \mu + \nu R^2\right)\frac{1}{x}\frac{\partial x}{\partial R} - \frac{1}{z}\frac{\partial z}{\partial R} &= \nu z^2\left[\frac{1}{x}\frac{\partial x}{\partial R}\left(1 + \mu + \nu R^2\right)z^2 - \frac{1}{z}\frac{\partial z}{\partial R}\right]\nonumber\\ \Leftrightarrow\left(1 + \mu + \nu R^2\right)\frac{1}{x}\frac{\partial x}{\partial R} - \frac{1}{z}\frac{\partial z}{\partial R} &= 0. \end{align} \]

If you convert this to $\frac{\partial z}{\partial R}$, you get

\[ \begin{align} \frac{\partial z}{\partial R} = \frac{z}{x}\left(1 + \mu + \nu R^2\right)\frac{\partial x}{\partial R}.\tag{D.236}\label{eq:great_deriv_7} \end{align} \]

Putting this into Eq. (D.233), you get

\[ \begin{align} & \left(1 + \mu + \nu R^2\right)^2\frac{z^2}{x}\frac{\partial x}{\partial R} + x\frac{\partial x}{\partial R} = R\left(1 - \nu z^2\right)\nonumber\\ &\Leftrightarrow x\frac{\partial x}{\partial R}\left[1 + \frac{z^2}{x^2}\left(1 + \mu + \nu R^2\right)^2\right] = R\left(1 - \nu z^2\right)\nonumber\\ &\Leftrightarrow\frac{\partial x}{\partial R} = \frac{R\left(1 - \nu z^2\right)x^2}{x\left[x^2 + z^2\left(1 + \mu + \nu R^2\right)^2\right]}\nonumber\\ &\Leftrightarrow\frac{\partial x}{\partial R} = \frac{R\left(1 - \nu z^2\right)\left(R^2 - \left(1 + \mu + \nu R^2\right)z^2\right)}{x\left[R^2 - \left(1 + \mu + \nu R^2\right)z^2 + z^2\left(1 + \mu + \nu R^2\right)^2\right]}\nonumber\\ &\Leftrightarrow\frac{\partial x}{\partial R} = \frac{R\left(1 - \nu z^2\right)\left(R^2 - \left(1 + \mu + \nu R^2\right)z^2\right)}{x\left[R^2 + z^2\left(1 + \mu + \nu R^2\right)\left(\mu + \nu R^2\right)\right]}.\tag{D.237}\label{eq:great_deriv_8} \end{align} \]

Equation was again used. (D.227) is used. If you put Eq. (D.237) in Eq. (D.236) and uses Eq. (D.229), follows

\[ \begin{align} \frac{\partial z}{\partial R} &= \frac{zR\left(1 - \nu z^2\right)\left(1 + \mu + \nu R^2\right)}{R^2 + z^2\left(1 + \mu + \nu R^2\right)\left(\mu + \nu R^2\right)}.\tag{D.238}\label{eq:great_deriv_9} \end{align} \]

So with Eq. (D.227)

\[ \begin{align} h_R &= \sqrt{\left(\frac{\partial x}{\partial R}\right)^2 + \left(\frac{\partial z}{\partial R}\right)^2} = \frac{R\left(1 - \nu z^2\right)\sqrt{z^2\left(1 + \mu + \nu R^2\right)^2 + x^2}}{R^2 + z^2\left(1 + \mu + \nu R^2\right)\left(\mu + \nu R^2\right)}\nonumber \end{align} \]

Now one differentiates Eq. (D.229) partial to $\varphi$:

\[ \begin{align} \frac{2}{z}\frac{\partial z}{\partial\varphi} &= -2x\frac{\frac{\partial x}{\partial\varphi}}{R^2 - x^2}\Leftrightarrow\frac{R^2 - x^2}{z}\frac{\partial z}{\partial\varphi} = -x\frac{\partial x}{\partial\varphi}. \end{align} \]

If you put Eq. (D.227), you get

\[ \begin{align} \left(1 + \mu + \nu R^2\right)z\frac{\partial z}{\partial\varphi} + x\frac{\partial x}{\partial\varphi} &= 0.\tag{D.241}\label{eq:great_deriv_10} \end{align} \]

Proceed analogously with Eq. (D.230), you get

\[ \begin{align} & \frac{1}{z}\frac{\partial z}{\partial\varphi} = \frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi} + \frac{1 + \mu}{x}\frac{\partial x}{\partial\varphi} + \nu x\frac{\partial x}{\partial\varphi} + \nu z\frac{\partial z}{\partial\varphi}\nonumber\\ &\Leftrightarrow\left[1 + \mu + \nu x^2\right]\frac{1}{x}\frac{\partial x}{\partial\varphi} - \left[1 - \nu z^2\right]\frac{1}{z}\frac{\partial z}{\partial\varphi} = -\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}\nonumber\\ &\Leftrightarrow\frac{1 + \mu + \nu x^2}{1 - \nu z^2}\frac{1}{x}\frac{\partial x}{\partial\varphi} - \frac{1}{z}\frac{\partial z}{\partial\varphi} = -\frac{1}{1 - \nu z^2}\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}\nonumber\\ &\Leftrightarrow\left(1 + \mu + \nu R^2\right)\frac{1}{x}\frac{\partial x}{\partial\varphi} - \frac{1}{z}\frac{\partial z}{\partial\varphi} = -\frac{1}{1 - \nu z^2}\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}. \end{align} \]

In the last step, the version of Eq. rescaled using the equations (D.222) - (D.222) was used. (D.220) used. If you convert this to $\frac{\partial z}{\partial\varphi}$, you get

\[ \begin{align} \frac{\partial z}{\partial\varphi} &= z\left(1 + \mu + \nu R^2\right)\frac{1}{x}\frac{\partial x}{\partial\varphi} + \frac{z}{1 - \nu z^2}\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}.\tag{D.243}\label{eq:great_deriv_11} \end{align} \]

Putting this into Eq. (D.241), follows

\[ \begin{align} & \left(1 + \mu + \nu R^2\right)^2z^2\frac{1}{x}\frac{\partial x}{\partial\varphi} + z^2\frac{1 + \mu + \nu R^2}{1 - \nu z^2}\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi} + x\frac{\partial x}{\partial\varphi} = 0\nonumber\\ &\Leftrightarrow\frac{\partial x}{\partial\varphi}\left(x + \frac{z^2}{x}\left(1 + \mu + \nu R^2\right)^2\right) = -\frac{z^2\left(1 + \mu + \nu R^2\right)}{1 - \nu z^2}\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}\nonumber\\ &\Leftrightarrow\frac{\partial x}{\partial\varphi}\left(x^2 + z^2\left(1 + \mu + \nu R^2\right)^2\right) = -\frac{xz^2\left(1 + \mu + \nu R^2\right)}{1 - \nu z^2}\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}\nonumber\\ &\Leftrightarrow\frac{\partial x}{\partial\varphi}\left(x^2 + z^2\left(1 + \mu + \nu R^2\right)\left(\mu + \nu R^2\right) + z^2\left(1 + \mu + \nu R^2\right)\right)\nonumber\\ &= -\frac{xz^2\left(1 + \mu + \nu R^2\right)}{1 - \nu z^2}\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}. \end{align} \]

If you put Eq. (D.227), you get

\[ \begin{align} \frac{\partial x}{\partial\varphi} = -\frac{xz^2\left(1 + \mu + \nu R^2\right)}{\left(1 - \nu z^2\right)\left(R^2 + \left(\mu + \nu R^2\right)\left(1 + \mu + \nu R^2\right)z^2\right)}\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}. \end{align} \]

Putting this into Eq. (D.243) and uses Eq. (D.227), so follows

\[ \begin{align} \frac{\partial z}{\partial\varphi} &= - \frac{z^3\left(1 + \mu + \nu R^2\right)^2}{\left(1 - \nu z^2\right)\left(R^2 + \left(\mu + \nu R^2\right)\left(1 + \mu + \nu R^2\right)z^2\right)}\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi} + \frac{z}{1 - \nu z^2}\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}\nonumber\\ &= \frac{zR^2 + z^3\left(1 + \mu + \nu R^2\right)\left(\mu + \nu R^2\right) - z^3\left(1 + \mu + \nu R^2\right)^2}{\left(1 - \nu z^2\right)\left(R^2 + \left(\mu + \nu R^2\right)\left(1 + \mu + \nu R^2\right)z^2\right)}\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}\nonumber\\ \Leftrightarrow\frac{\partial z}{\partial\varphi} &= \frac{zR^2 - z^3\left(1 + \mu + \nu R^2\right)}{\left(1 - \nu z^2\right)\left(R^2 + \left(\mu + \nu R^2\right)\left(1 + \mu + \nu R^2\right)z^2\right)}\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}\nonumber\\ \Leftrightarrow\frac{\partial z}{\partial\varphi} &= \frac{zx^2}{\left(1 - \nu z^2\right)\left(R^2 + \left(\mu + \nu R^2\right)\left(1 + \mu + \nu R^2\right)z^2\right)}\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}. \end{align} \]

This means that again with Eq. (D.227)

\[ \begin{align} & h_\varphi = \sqrt{\left(\frac{\partial x}{\partial\varphi}\right)^2 + \left(\frac{\partial\varphi}{\partial z\varphi}\right)^2} = \frac{xz\sqrt{x^2 + z^2\left(1 + \mu + \nu R^2\right)^2}}{\left(1 - \nu z^2\right)\left(R^2 + \left(\mu + \nu R^2\right)\left(1 + \mu + \nu R^2\right)z^2\right)}\left|\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}\right|\nonumber\\ &\Leftrightarrow h_\varphi = \frac{xz\sqrt{R^2 + z^2\left(\mu + \nu R^2\right)\left(1 + \mu + \nu R^2\right)}}{\left(1 - \nu z^2\right)\left(R^2 + \left(\mu + \nu R^2\right)\left(1 + \mu + \nu R^2\right)z^2\right)}\left|\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}\right|\nonumber \end{align} \]

If you logarithm both sides of Eq. (D.229), you get

\[ \begin{align} \ln\left(D\right) = \ln\left(\left|z_S\right|\right) - \left(1 + \mu\right)\ln\left(x_S\right) - \frac{\nu}{2}\left(x_S^2 + z_S^2\right). \end{align} \]

If you differentiate this, it follows

\[ \begin{align} & \frac{d\ln\left(D\right)}{d\varphi} = \frac{d\ln\left(\left|z_S\right|\right)}{d\varphi} - \left(1 + \mu\right)\frac{d\ln\left(x_S\right)}{d\varphi} - \nu x_S\frac{dx_S}{d\varphi} - \nu z_S\frac{d z_S}{d\varphi}\nonumber\\ &= \frac{d\ln\left(\left|z_S\right|\right)}{d\varphi}\left(1 - \nu z_S^2\right) - \frac{d\ln\left(x_S\right)}{d\varphi}\left(1 + \mu + \nu x_S^2\right).\tag{D.249}\label{eq:great_deriv_12} \end{align} \]

This equation can be used to express $\frac{d\ln\left(D\right)}{d\varphi}$ in terms of $\frac{d x_S}{d\varphi}$ and $\frac{d z_S}{d\varphi}$:

\[ \begin{align} h_\varphi = \frac{xz}{\left(1 - \nu z^2\right)\left(R^2 + \left(\mu + \nu R^2\right)\left(1 + \mu + \nu R^2\right)z^2\right)^{1/2}}\left|\frac{d\ln\left(\left|z_S\right|\right)}{d\varphi}\left(1 - \nu z_S^2\right) - \frac{d\ln\left(x_S\right)}{d\varphi}\left(1 + \mu + \nu x_S^2\right)\right| \end{align} \]

They apply

\[ \begin{align} x = \sqrt{x^2 + y^2}\cos\lambda, & {} & y = \sqrt{x^2 + y^2}\sin\lambda. \end{align} \]

Thus follows

For the determinant $g$ of the metric tensor of the Great coordinates holds

\[ \begin{align} \sqrt{g} &= h_Rh_\lambda h_\varphi = h_Rr_\perp h_\varphi = h_Rr_\perp\frac{r_\perp z}{\left(1 - \nu z^2\right)\left(R^2 + \left(\mu + \nu R^2\right)\left(1 + \mu + \nu R^2\right)z^2\right)^{1/2}}\left|\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}\right|\nonumber\\ &= h_R\frac{r_\perp^2z}{\left(1 - \nu z^2\right)\left(R^2 + \left(\mu + \nu R^2\right)\left(1 + \mu + \nu R^2\right)z^2\right)^{1/2}}\left|\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}\right|\nonumber\\ &= \frac{R\left(1 - \nu z^2\right)}{\left[R^2 + z^2\left(1 + \mu + \nu R^2\right)\left(\mu + \nu R^2\right)\right]^{1/2}}\frac{r_\perp^2z}{\left(1 - \nu z^2\right)\left(R^2 + \left(\mu + \nu R^2\right)\left(1 + \mu + \nu R^2\right)z^2\right)^{1/2}}\left|\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}\right|\nonumber\\ &= \frac{R}{\left[R^2 + z^2\left(1 + \mu + \nu R^2\right)\left(\mu + \nu R^2\right)\right]^{1/2}}\frac{r_\perp^2z}{\left(R^2 + \left(\mu + \nu R^2\right)\left(1 + \mu + \nu R^2\right)z^2\right)^{1/2}}\left|\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}\right|\nonumber\\ &= \frac{Rr_\perp^2z}{R^2 + \left(\mu + \nu R^2\right)\left(1 + \mu + \nu R^2\right)z^2}\left|\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}\right|\nonumber \end{align} \]

D.4.5 Transformation von Great-Koordinaten in globale Koordinaten

Let a point $\mathbf{r}$ with the Great coordinates $\left(R, \varphi, \lambda\right)$ be given. One wants to determine the global coordinates $\left(x, y, z\right)^T$ of this point. First assume $\lambda = 0$. If you square Eq. (D.228), you get

\[ \begin{align} z^2 = D\left(\varphi\right)^2\left(x^2\right)^{1 + \mu}\exp\left[\nu\left(x^2 + z^2\right)\right]. \end{align} \]

If you put Eq. (D.230), follows

\[ \begin{align} z^2 = D\left(\varphi\right)^2\left[R^2 - z^2\left(1 + \mu + \nu R^2\right)\right]^{1 + \mu}\exp\left[\nu\left(R^2 - z^2\left(\mu + \nu R^2\right)\right)\right].\tag{D.255}\label{eq:trafo_great_zu_glo_deriv_1} \end{align} \]

This equation is implicit and therefore has to be solved iteratively. $D\left(\varphi\right)$ and $R$ are constants, $z^2$ is the only remaining quantity to be determined. About Eq. (D.230) can then be determined $x^2$. The sign of $z$ is the sign of $\varphi$. In the case $\lambda\not = 0$

\[ \begin{align} x^2\to r_\perp^2 \end{align} \]

to replace so that one

\[ \begin{align} x &= r_\perp\cos\left(\lambda\right),\\ y &= r_\perp\sin\left(\lambda\right) \end{align} \]

receives.

D.4.6 Transformation von globalen Koordinaten in Great-Koordinaten

Given a point $\mathbf{r}$ with global coordinates $\left(x, y, z\right)^T$, one wants to determine the Great coordinates $\left(R, \varphi, \lambda\right)$ of $\mathbf{r}$. If you write down Eq. (D.212) with the replacement equations (D.222) - (D.222), one obtains

\[ \begin{align} R^2 = \frac{x^2 + y^2 + \left(1 + \mu\right)z^2}{1 - \nu z^2}. \end{align} \]

$x^2\to x^2 + y^2$ was replaced because the general case $\lambda\not = 0$ is assumed. $\lambda$ is obtained via

\[ \begin{align} \lambda = \sign\left(y\right)\arccos\left(\frac{x}{\sqrt{x^2 + y^2}}\right). \end{align} \]

To determine $\varphi$, one first sets Eq. (D.228) to $D$ and replaces $x^2\to x^2 + y^2$ again:

\[ \begin{align} D = \frac{z}{\sqrt{x^2 + y^2}^{1 + \mu}}\exp\left(-\frac{\nu}{2}\left(x^2 + y^2 + z^2\right)\right)\tag{D.261}\label{eq:trafo_glo_zu_great_deriv_5} \end{align} \]

Now you put in Eq. (D.255) $R = 1$, which corresponds to the surface of the WGS 84, and solves numerically for $z_S^2$, the sign of $z_S$ is identical to the sign of $z$. Now you can use Eq. Use (D.227) with $R = 1$, $z = z_S$ and $x^2\to x_S^2 + y_S^2$ to determine $x_S^2 + y_S^2$. About the equation

\[ \begin{align} \varphi = \arctan\left[\frac{a^2}{c^2}\frac{z_S}{\sqrt{x^2 + y^2}}\right]\tag{D.262}\label{eq:eq:trafo_glo_zu_great_deriv_1} \end{align} \]

one can finally determine $\varphi$. Eq. (D.262) should now be derived. Please consider o. B. d. A. the plane $\lambda = 0$, therefore $y = 0$. The intersection of the WGS 84 with this plane is given by (note the scaling)

\[ \begin{align} x^2 + \frac{z^2}{c^2} = 1.\tag{D.263}\label{eq:def_ell_skal} \end{align} \]

Define

\[ \begin{align} f\left(x, z\right) \coloneqq x^2 + \frac{z^2}{c^2}, \end{align} \]

then applies

\[ \begin{align} \nabla f = \left(\begin{array}{c} 2x\\ \frac{2z}{c^2} \end{array}\right). \end{align} \]

The WGS 84 is given by the contour set $f = 1$. Define the line $\mathbf{g}$ in the xz plane of the global coordinates by

\[ \begin{align} \mathbf{g} = \left(\begin{array}{c} x_S\\ z_S \end{array}\right) + \gamma\left(\begin{array}{c} x_S\\ \frac{z_S}{c^2} \end{array}\right) \end{align} \]

with $\gamma\in \mathbb{R}$. The intersection of $\mathbf{g}$ with the x-axis is given by in the case $z_S\not = 0$

\[ \begin{align} 0 = z_S + \gamma \frac{z_S}{c^2} \Rightarrow \gamma = -c^2. \end{align} \]

For the x coordinate of this point we get

\[ \begin{align} x = x_S\left(1 - c^2\right). \end{align} \]

So it applies

\[ \begin{align} \tan\left(\varphi\right) = \frac{z_S}{x_Sc^2}.\tag{D.269}\label{eq:trafo_glo_zu_great_deriv_2} \end{align} \]

With Eq. (D.263) follows

\[ \begin{align} z_S^2 = c^2\left(1 - x_S^2\right) \end{align} \]

and therefore applies

\[ \begin{align} \tan\left(\varphi\right) &= \frac{c\sqrt{1 - x_S^2}}{x_Sc^2} = \frac{\sqrt{1 - x_S^2}}{x_Sc}\Leftrightarrow \tan^2\left(\varphi\right) = \frac{1 - x_S^2}{x_S^2c^2}\Leftrightarrow x_S^2\left(c^2\tan^2\left(\varphi\right) + 1\right) = 1\nonumber\\ \Leftrightarrow x_S^2 &= \frac{\cos^2\left(\varphi\right)}{c^2\sin^2\left(\varphi\right) + \cos^2\left(\varphi\right)}. \end{align} \]

With $x_S>0$ follows

\[ \begin{align} x_S = \frac{\cos\left(\varphi\right)}{\sqrt{c^2\sin^2\left(\varphi\right) + \cos^2\left(\varphi\right)}}.\tag{D.272}\label{eq:trafo_glo_zu_great_deriv_3} \end{align} \]

If you put Eq. (D.263) to $x_S$ is obtained

\[ \begin{align} x_S = \sqrt{1 - \frac{z^2}{c^2}}. \end{align} \]

Putting this into Eq. (D.269), you get

\[ \begin{align} \tan\left(\varphi\right) &= \frac{z_S}{c^2\sqrt{1 - \frac{z_S^2}{c^2}}}\Leftrightarrow c^4\tan^2\left(\varphi\right)\left(1 - \frac{z_S^2}{c^2}\right) = z_S^2\nonumber\\ \Leftrightarrow z_S^2\left(1 + c^2\tan^2\left(\varphi\right)\right) &= c^4\tan^2\left(\varphi\right)\Leftrightarrow z_S^2 = \frac{c^4\sin^2\left(\varphi\right)}{\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)}\nonumber\\ \Leftrightarrow z_S &= \sign\left(\varphi\right)\frac{c^2\sin\left(\varphi\right)}{\sqrt{\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)}}.\tag{D.274}\label{eq:trafo_glo_zu_great_deriv_4} \end{align} \]

It therefore applies

\[ \begin{align} \frac{z_S}{x_S} = \frac{c^2}{a^2}\tan\left(\varphi\right), \end{align} \]

with the generalization $x_S\to\sqrt{x_S^2 + y_S^2}$ follows Eq. (D.262).

Now we need to find an explicit expression for $D\left(\varphi\right)$. To do this, put the equations (D.272) and (D.274) with $y = 0$ in Eq. (D.261):

\[ \begin{align} & D\left(\varphi\right) = \frac{\left|z_S\right|}{x_S^{1 + \mu}}\exp\left(-\frac{\nu}{2}\left(x_S^2 + z_S^2\right)\right)\nonumber\\ & = \frac{c^2\sin\left(\varphi\right)\left(\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)\right)^{\mu/2}}{\left(\cos\left(\varphi\right)\right)^{1 + \mu}}\exp\left[-\frac{\nu}{2}\left(\frac{\cos^2\left(\varphi\right) + c^4\sin^2\left(\varphi\right)}{\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)}\right)\right]. \end{align} \]

If you take this logarithm, it follows

\[ \begin{align} \ln\left[D\left(\varphi\right)\right] &= 2\ln\left(c\right) + \ln\left(\sin\left(\varphi\right)\right) + \frac{\mu}{2}\ln\left(\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)\right) - \left(1 + \mu\right)\ln\left(\cos\left(\varphi\right)\right)\nonumber\\ & - \frac{\nu}{2}\left(\frac{\cos^2\left(\varphi\right) + c^4\sin^2\left(\varphi\right)}{\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)}\right). \end{align} \]

Thus you get

\[ \begin{align} & \frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi} = \frac{\cos\left(\varphi\right)}{\sin\left(\varphi\right)} + \frac{\mu}{2}\frac{ - 2\sin\left(\varphi\right)\cos\left(\varphi\right) + c^22\cos\left(\varphi\right)\sin\left(\varphi\right)}{\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)} + \left(1 + \mu\right)\frac{\sin\left(\varphi\right)}{\cos\left(\varphi\right)}\nonumber\\ & - \frac{\nu}{2}\frac{\left(-2\cos\left(\varphi\right)\sin\left(\varphi\right) + 2c^4\cos\left(\varphi\right)\sin\left(\varphi\right)\right)\left(\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)\right)}{\left(\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)\right)^2}\nonumber\\ & + \frac{\nu}{2}\frac{\left(\cos^2\left(\varphi\right) + c^4\sin^2\left(\varphi\right)\right)\left(-2\sin\left(\varphi\right)\cos\left(\varphi\right) + c^22\cos\left(\varphi\right)\sin\left(\varphi\right)\right)}{\left(\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)\right)^2}\nonumber\\ &= \frac{1}{\sin\left(\varphi\right)}\left(\cos\left(\varphi\right) + \frac{\sin^2\left(\varphi\right)}{\cos\left(\varphi\right)}\right) + \mu\sin\left(\varphi\right)\cos\left(\varphi\right)\frac{c^2 - 1 + 1 + c^2\tan^2\left(\varphi\right)}{\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)}\nonumber\\ & + \nu\frac{\left(\cos\left(\varphi\right)\sin\left(\varphi\right) - c^4\cos\left(\varphi\right)\sin\left(\varphi\right)\right)\left(\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)\right)}{\left(\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)\right)^2}\nonumber\\ & + \nu\frac{\left(\cos^2\left(\varphi\right) + c^4\sin^2\left(\varphi\right)\right)\left(-\sin\left(\varphi\right)\cos\left(\varphi\right) + c^2\cos\left(\varphi\right)\sin\left(\varphi\right)\right)}{\left(\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)\right)^2}\nonumber\\ &= \frac{1}{\sin\left(\varphi\right)\cos\left(\varphi\right)} + \mu\tan\left(\varphi\right)\frac{c^2}{\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)}\nonumber\\ & + \nu\frac{\sin\left(\varphi\right)\cos\left(\varphi\right)\left(c^2\left(\sin^2\left(\varphi\right) + \cos^2\left(\varphi\right)\right) + c^4\left(-\sin^2\left(\varphi\right) - \cos^2\left(\varphi\right)\right)\right)}{\left(\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)\right)^2}\nonumber \end{align} \]

D.5 Ellipsoidal geometry

D.5.1 Ellipsoidal geodesics

D.5.2 Ellipsoidal triangles

An area $A$ on an ellipsoid is usually calculated using an area-conformal projection. The basic procedure is to transform each vertex $\left(\phi_i,\lambda_i\right)$ into the projection plane, i.e. to determine its coordinates $\left(x_i,y_i\right)$, and then calculate the area $A'$ in plane geometry. Area-conformal projections are exactly those for which $A = A'$ holds.

The Lambert azimuthal projection is used here as the surface-conformal projection.

D.5.3 Ellipsoidal Voronoi lattices

[title = Literatur]