Here one chooses a basis $\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3$ situated at the center of the Earth that does not co-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 written as
\[ \begin{align} \mathbf{v} = u_1\mathbf{e}_1 + u_2\mathbf{e}_2 + u_3\mathbf{e}_3 \end{align} \]
The wind field is thereby not equal to the particle motion measured in this system, since the rotation of the Earth is additionally superimposed.
Particle motion in resting coordinates = wind field + Earth's rotation
Here one chooses a basis $\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z$ situated at the center of the Earth just like the basis of the resting coordinates, but co-rotating with it. $\mathbf{e}_x$ points towards the intersection of the prime meridian with the equator, $\mathbf{e}_y$ points towards the intersection of the ninetieth meridian with the equator, and $\mathbf{e}_z$ points towards 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} \]
Alternatively, one may use geographical coordinates $\left(r, \varphi, \lambda\right)$, where $r\geq 0$ is the distance from Earth's center, $-\pi/2\leq \varphi\leq\pi/2$ is the angle between the position vector and the equatorial plane (geographical latitude for a spherical Earth), and $0\leq\lambda<2\pi$ is the angle between the projection onto the equatorial plane and the x-axis (geographical longitude). These coordinates can be transformed to standard spherical coordinates $\left(r, \theta, \phi\right)$ via $\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, spherical coordinates are inconvenient for vector notation. Therefore, one introduces a local right-handed Cartesian system tangent to the sphere: the x-axis points east, the y-axis north, and the z-axis upward. Its basis vectors are denoted by $\mathbf{i}, \mathbf{j}, \mathbf{k}$. The wind field $\mathbf{v}$ is then written as
\[ \begin{align} \mathbf{v} = u\mathbf{i} + v\mathbf{j} + w\mathbf{k}. \end{align} \]
We now determine the shortest surface connection (geodesic) between two points $\left(\varphi_i, \lambda_i\right)$ on a sphere, with $i = 1, 2$. The length $L$ of an arbitrary connecting curve is
\[ \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]$ parametrizes the curve and $\mathbf{r}\left(\tau\right)$ specifies 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)$.
To determine $\varphi\left(\tau\right)$ and $\lambda\left(\tau\right)$, one first works on the unit sphere and ignores the radius. Then
\[ \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
\[ \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
\[ \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,
\[ \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{ in the case }\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)$, one obtains
\[ \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$, this reduces to
\[ \begin{align} d = r\sqrt{2 - 2\left(\cos\left(\varphi_a\right)\cos\left(\varphi_b\right)\cos(\lambda_a - \lambda_b) + \sin\left(\varphi_a\right)\sin\left(\varphi_b\right)\right)}. \end{align} \]
The surface distance between two points on a sphere follows directly. For a sphere of radius $r$, the arc length $\Delta$ between two surface points is
\[ \begin{align} \Delta = r\theta, \end{align} \]
where $\theta$ is the angle between the two position vectors, and
\[ \begin{align} \sin\left(\frac{\theta}{2}\right) = \frac{d}{2r}. \end{align} \]
Since $\frac{\theta}{2}\in\left[0, \frac{\pi}{2}\right]$, the distance $\Delta$ is
Now $\tau$ must be transformed to $\tau'$. We have
\[ \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} \]
Hence,
\[ \begin{align} \tau' &= \frac{1}{2} + \sqrt{\frac{r^2}{d^2} - \frac{1}{4}}\tan\left[\theta\left(\tau - \frac{1}{2}\right)\right].\\ \end{align} \]
For vertical surfaces $A_v$ with base length $L$ and the same notation for the radii, one obtains
\[ \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} \]
Volumes $V$ with base area $A$, inner radius $r_0$, and outer radius $r_1$ are given by
\[ \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} \]
An ellipse in the xz plane is described by
\[ \begin{align} \frac{x^2}{a^2} + \frac{z^2}{c^2} = 1. \end{align} \]
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 is defined as
\[ \begin{align} \epsilon \coloneqq \frac{a - c}{a} \end{align} \]
This gives
\[ \begin{align} c = a\left(1 - \epsilon\right). \end{align} \]
In polar coordinates, with
\[ \begin{align} x = r\cos\left(\varphi\right), & {} & z = r\sin\left(\varphi\right) \end{align} \]
the following equations follow:
\[ \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} \]
this can be rewritten as
\[ \begin{align} r\left(\varphi\right) &= a\left[1 + \frac{2\epsilon - \epsilon^2}{\left(1 - \epsilon\right)^2}\sin^2\left(\varphi\right)\right]^{-1/2}. \end{align} \]
From Eq. (D.29) it further follows
\[ \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} \]
An ellipsoid of rotation, or simply an ellipsoid, is obtained by rotating the ellipse just examined about the z-axis. Its volume $V$ is computed 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} \]
The ellipsoidal coordinates will be oriented on the Earth's geopotential $\phi_g$, for which one has
\[ \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. For the centrifugal potential, one has
\[ \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} \]
since taking the gradient yields
\[ \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 relation
\[ \begin{align} \mathbf{a} = -\nabla\phi_g \end{align} \]
holds 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$ the mass distribution, since for a point mass $M$ at the origin this yields, by 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. It must be solved to determine the gravitational potential of the Earth. In general, one can set up 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 expanded in terms of certain functions. As an ansatz, with $\mu \coloneqq\cos\left(\theta\right)$, a separation ansatz
\[ \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} \]
is made. Substituting this into $\Delta\phi_g = 0$, one obtains, using Eq. (B.91) and the observation
\[ \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$. Thus, for $Q$ one has
\[ \begin{align} Q'' = -m^2Q, \end{align} \]
which 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 one divides Eq. (D.44) by $1 - \mu^2$ and sets the right-hand side equal to $m^2$; one then obtains
\[ \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$. One obtains two differential equations, one for $P$ and one for $U$; the one for $P$ reads
\[ \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 one obtains
\[ \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} \]
Substituting $\lambda = n\left(n + 1\right)$ into Eq. (D.49), one obtains
\[ \begin{align} r^2U'' - n\left(n + 1\right)U = 0. \end{align} \]
$U$ depends on $n$; write $U_n$. One finds 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} \]
are solutions, and therefore the general solution of the ordinary second-order differential equation Eq. (D.52) is
\[ \begin{align} U_n\left(r\right) = a_nr^{n + 1} + \frac{b_n}{r^n}. \end{align} \]
Now solutions 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} \]
have been found; since the Laplace equation (D.40) is linear in $\phi_g$, arbitrary linear combinations of these are again solutions, and one obtains as the 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 particular solution for the case of the rotating ellipsoid is determined by the boundary conditions. It reads:
\[ \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$ thus reads
\[ \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$. Because of the usual boundary condition $\lim\limits_{r\to\infty}\phi_0\left(r\right) = 0$, however, one has
\[ \begin{align} a_0 = a_1 = a_2 = 0, \end{align} \]
furthermore, at infinite distance the solution should tend to $-G\frac{M}{r}$, since the ellipsoidal mass distribution looks like a point mass from afar; one therefore expects
\[ \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 them, one calculates the gravitational potential on the z-axis ($z = r> c$). There, with $\chi = \pi/2$, one has
\[ \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} \]
Here it was assumed that the $b_i$ depend on the eccentricity $\epsilon$. One thus 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 against a simple limiting case, namely the case of a sphere $a = c$. In this case, one obtains
\[ \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 is thus, outside the sphere, equal to that of a point mass at the location of the sphere's center of mass. According to [11] Eq. (3.3.37), one has
\[ \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,
\[ \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 conditions 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&<a\sqrt{\frac{a^2}{c^2} - 1 + \frac{r^2}{c^2}}.\tag{D.76}\label{eq:ungl_schwere} \end{align} \]
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} \]
the derivative of the right-hand expression with respect to $r$ is
\[ \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 equals one for $r = c$ and equals $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
\[ \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 (D.76) is satisfied for $r> c$, 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. One thus has
\[ \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} \]
Substituting $A, B$, and $C$ here, one obtains, using
\[ \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 thus follows, with
\[ \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 calculation
\[ \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 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} \]
can be incorporated:
\[ \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 Eq. (D.67). To achieve this, it is expanded in the small parameter $\epsilon\ll 1$. To this end, one defines 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 directly possible; instead, one must make 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} \]
One can therefore apply L'Hospital's rule to determine $\lim\limits_{\epsilon\to 0}\frac{g\left(\epsilon\right)}{h\left(\epsilon\right)}$. One has
\[ \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} \]
From this follow
\[ \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} \]
From this follow
\[ \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} \]
One thus has
\[ \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} \]
Furthermore,
\[ \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} \]
One thus has
\[ \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 one must again apply L'Hospital's rule; one has
\[ \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} \]
It thus 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 is thus recovered. For the derivative of $F$, one has
\[ \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} \]
By the arithmetic rules for limits of functions, $\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} \]
Using L'Hospital's rule, one has
\[ \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} \]
It thus follows
\[ \begin{align} \lim_{\epsilon\to 0} F'\left(\epsilon\right) = \frac{4}{15}\frac{a^3}{r^3}. \end{align} \]
The higher-order coefficients in $\epsilon$ can be determined analogously. One thus has
\[ \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. Thus, in terms of Eq. (D.67), one has
\[ \begin{align} b_1 &= 0,\\ b_2\left(\epsilon\right) &= \epsilon\frac{2GMa^2}{5}. \end{align} \]
Thus, for the gravitational component of the Earth's gravity potential, one has
\[ \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} \]
One now requires, as an additional condition, that the surface of the ellipsoid of revolution be 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} \]
Neglecting all terms $\mathcal{O}\left(\epsilon^2\right)$ and also $\mathcal{O}\left(\epsilon m\right)$ here, one obtains
\[ \begin{align} - \epsilon + \frac{\epsilon}{5} + \frac{m}{2} &= -\epsilon\frac{2}{5}\Leftrightarrow \frac{m}{2} = \epsilon\frac{2}{5}. \end{align} \]
Thus, in the equilibrium case, one has
\[ \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} \]
One uses
\[ \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:
\[ \begin{align} \phi_g\left(r, \chi\right) &= \frac{GM}{r}\left[\left\lbrace\left(\frac{2\epsilon - m}{2}\right)\frac{a^2}{r^2} + \frac{m}{2}\frac{r^3}{a^3}\right\rbrace\sin^2\left(\chi\right) - 1\right] - \frac{GM}{a}\left\lbrace\left(\frac{2\epsilon - m}{6}\right)\frac{a^3}{r^3} + \frac{m}{2}\frac{r^2}{a^2}\right\rbrace\tag{D.148}\label{eq:schwerepotential} \end{align} \]
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$.
For the construction of a meteorological coordinate system, it is advantageous if surfaces of constant vertical coordinate are equipotential surfaces of the gravity field (see Sect. 13.2). It is also desirable that these surfaces be ellipsoidal. In addition, one must be able to express trajectories orthogonal to all equipotential surfaces in terms of the underlying basis functions. The geopotential approximation should be as physically consistent as possible while still allowing analytical treatment. In particular, it should reproduce the correct Clairaut ratio and satisfy $c = a\left(1 - \epsilon\right)$.
The requirements for a gravity-field approximation are therefore:
possibility of analytically specifying vertical trajectories
ellipsoidal contour surfaces
correct Clairaut ratio
correct semi-axis ratio
The first question is whether all contour surfaces of Eq. (D.148) are ellipsoids. For Earth's surface (without orography), this is already guaranteed to 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 a rotational ellipsoid with semi-major axis $A$ and semi-minor axis $C$ as a function of the geocentric angle $\chi$:
\[ \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: keeping the left-hand side constant and multiplying by $r^3$ yields a fifth-degree polynomial in $r$, which in general can only be solved numerically.
Therefore, one must construct an ellipsoidal approximation to Eq. (D.148). The following discussion is based on [34]. First, one specifies how to label a particular geopotential surface. A convenient choice is its equatorial radius $R$. For this radius,
\[ \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 one now wants to characterize this surface more precisely, namely by its radial distance $r = r\left(\chi\right)$ from Earth's center at geocentric latitude $\phi$, one writes
\[ \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 assumes
\[ \begin{align} \frac{r}{R} = 1 + \frac{r - R}{R} = 1 + \mathcal{O}\left(\epsilon\right), \end{align} \]
Indeed,
\[ \begin{align} \frac{r - R}{R} \lesssim \frac{c - a}{a} = -\epsilon, \end{align} \]
at least within the meteorologically relevant atmosphere. One can then first write
\[ \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} \]
Neglecting terms of order $\mathcal{O}\left(\epsilon^2, m^2, \epsilon m\right)$, this first becomes
\[ \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 further simplifies 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} \]
Its contour surfaces are already ellipsoids. To show this, determine the geocentric distance $r$ as a function of latitude for given 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), one obtains exactly this result. Unfortunately, trajectories orthogonal to the contour surfaces of this potential cannot be determined analytically, so Eq. (D.155) must be approximated further. To this end, the terms $\frac{a^2}{R^2}$ and $\frac{R^3}{a^3}$ are expanded about $R^2 = a^2$, which is justified because $\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} \]
Substituting this into Eq. (D.155) yields
\[ \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} \]
which, using Eq. (D.149), becomes
\[ \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 issue arises here. If one sets $R = a$ in Eq. (D.162), one gets
\[ \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 Earth's surface without orography. One obtains
\[ \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} \]
it follows that
\[ \begin{align} r\left(\pi/2\right) = a - a\epsilon + \mathcal{O}\left(\epsilon^2\right)\approx a\left(1 - \epsilon\right), \end{align} \]
which is 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 one has
\[ \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} \]
Hence,
\[ \begin{align} \newtilde{\epsilon} = \epsilon + \mathcal{O}\left(\epsilon^2\right). \end{align} \]
Thus, to first order, one may replace $\epsilon$ by $\newtilde{\epsilon}$ everywhere. This yields
\[ \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, one likewise obtains 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).
Clairaut's theorem relates the eccentricity $\epsilon$, angular velocity $\omega$, and the ratio of gravitational acceleration $g_a, g_p$ at equator and pole. At the pole, the geopotential is $\phi_g = \phi_g^{(p)}$ (from 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} \]
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,
\[ \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} \]
Evaluating at $r = a$ and $r = c = a\left(1 - \epsilon\right)$ gives
\[ \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} \]
it follows that
\[ \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 one inserts $\epsilon\approx\frac{5m}{4}$ here, one obtains
\[ \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 one has
\[ \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} \]
Hence,
\[ \begin{align} \frac{g_p - g_a}{g_a} \approx 5\frac{m}{4}.\tag{D.188}\label{eq:clairaut_deriv_1} \end{align} \]
Therefore,
This is Clairaut's theorem. We now verify it for Eq. (D.172). It is sufficient to show that Eq. (D.188) still holds. At the equator,
\[ \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,
\[ \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$, Eq. (D.173) gives
\[ \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 $\chi = \pi/2$ and $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} \]
Thus,
\[ \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} \]
It thus follows
\[ \begin{align} \frac{g_p - g_a}{g_a} \approx \frac{5m}{2} - \newtilde{\epsilon}\approx \frac{5m}{2} - \epsilon, \end{align} \]
so Eq. (D.172) also satisfies Clairaut's theorem. Note that Eq. (D.172) is first-order accurate in $\epsilon$ and $m$, and is therefore neither better nor worse than Eq. (D.148) at that order.
Since horizontal derivatives in a suitable geophysical coordinate system are taken along equipotential surfaces of the gravity field, vertical derivatives (for an orthogonal coordinate system) are perpendicular to them and therefore aligned with gravity. It is therefore essential to determine this direction. In this section, trajectories perpendicular to the equipotential surfaces are derived.
Eq. (D.173) describes the equipotential surfaces of the gravity field, more precisely their distance $r$ from Earth's center for given equatorial radius $R$ and geocentric latitude $\chi$:
\[ \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} \]
where
\[ \begin{align} \mu \coloneqq \frac{8\newtilde{\epsilon} - 5m}{2}, \nu \coloneqq \frac{5m - 4\newtilde{\epsilon}}{2} \end{align} \]
are defined. This can be rewritten as
\[ \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} \]
Ellipsoids satisfy
\[ \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$. Substituting Eq. (D.206) gives
\[ \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 satisfied for $c = C$. Therefore, $C$ is the semi-minor axis of the equipotential surfaces.
Since an ellipsoid is rotationally symmetric, it is sufficient to first determine the vertical trajectories in the meridional plane $\lambda = 0$. There, one only needs to consider ellipses, which satisfy
\[ \begin{align} \frac{x^2}{R^2} + \frac{z^2}{C^2} = 1.\tag{D.210}\label{eq:deriv_senkr_ellipse_2d} \end{align} \]
Inserting the definition of $C$ gives
\[ \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} \]
Solving for $R^2$ gives
\[ \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,
\[ \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} \]
For fixed $R$, the coordinate $z$ of an ellipse in one half-space can be treated as a function of $x$. This function $z\left(x\right)$ satisfies Eq. (D.210). Differentiating with respect to $x$ gives
\[ \begin{align} \frac{dz}{dx} = -\frac{C^2}{R^2}\frac{x}{z}. \end{align} \]
This is a directional derivative along the equipotential surface. Geometrically, the local normal to the equipotential surface satisfies
\[ \begin{align} \frac{dz}{dx} = \frac{R^2}{C^2}\frac{z}{x}. \end{align} \]
Using the definition of $C$, one has
\[ \begin{align} \frac{R^2}{C^2} = 1 + \mu + \nu\frac{R^2}{a^2}, \end{align} \]
therefore
\[ \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 by separation of variables as presented in Sect. A.7. Since $z\left(x\right)$ is bijective, the required conditions are satisfied. First,
\[ \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} \]
Thus,
\[ \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 integration constant $D \coloneqq \exp\left(D'\right)>0$, where $D'\in \mathbb{R}$. For given $x$, solving this equation yields the corresponding $z$ on the desired orthogonal trajectory. The specific trajectory is determined by $D$, so $D$ is a function of latitude $\varphi$, i.e. $D = D\left(\varphi\right)$. The exact form of $D\left(\varphi\right)$ is a convention. Here it is fixed by the intersection of the trajectory with WGS 84. Let $\left(x_S, z_S\right)^T$ denote this intersection. Then
\[ \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} \]
Hence, the sought vertical trajectories are determined by
The Great coordinates are described in [34]. They are suitable for expressing differential operators in ellipsoidal geometries, and they are 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
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.
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.
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$.
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$, it is sufficient to work in the meridional plane $\lambda = 0$. Applying the rescaling relations (D.222) - (D.222) to Eqs. (D.211) and (D.219) yields
\[ \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} \]
Taking logarithms gives
\[ \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)$ follows from Eq. (D.220), which in rescaled form reads
\[ \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 differentiate Eq. (D.231) partially with respect 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} \]
Using Eq. (D.227), one obtains
\[ \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} \]
Differentiating Eq. (D.230) partially with respect to $R$ gives
\[ \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} \]
Again using Eq. (D.227), one obtains
\[ \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} \]
Solving for $\frac{\partial z}{\partial R}$ gives
\[ \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} \]
Substituting this into Eq. (D.233) gives
\[ \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} \]
Here, Eq. (D.227) was used again. Substituting Eq. (D.237) into Eq. (D.236) and using Eq. (D.229) yields
\[ \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} \]
Therefore, 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} \]
\[ \begin{align} \Leftrightarrow h_R &= \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}}. \end{align} \]
Now differentiate Eq. (D.229) partially with respect 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} \]
Using Eq. (D.227), one obtains
\[ \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} \]
Proceeding analogously with Eq. (D.230) gives
\[ \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 rescaled version of Eq. (D.220) (using (D.222) - (D.222)) was used. Solving for $\frac{\partial z}{\partial\varphi}$ gives
\[ \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} \]
Substituting this into Eq. (D.241) yields
\[ \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} \]
Using Eq. (D.227), one obtains
\[ \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} \]
Substituting this into Eq. (D.243) and using Eq. (D.227) yields
\[ \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} \]
Hence, 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} \]
\[ \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[D\left(\varphi\right)\right]}{d\varphi}\right|. \end{align} \]
Taking logarithms on both sides of Eq. (D.229) gives
\[ \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} \]
Differentiating this expression yields
\[ \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} \]
We have
\[ \begin{align} x = \sqrt{x^2 + y^2}\cos\lambda, & {} & y = \sqrt{x^2 + y^2}\sin\lambda. \end{align} \]
Thus,
\[ \begin{align} h_\lambda = \sqrt{\left(\frac{\partial x}{\partial \lambda}\right)^2 + \left(\frac{\partial y}{\partial \lambda}\right)^2} = \sqrt{x^2 + y^2} = r_\perp. \end{align} \]
For the determinant $g$ of the metric tensor in Great coordinates, one has
\[ \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} \]
\[ \begin{align} \Leftrightarrow\sqrt{g} &= \frac{Rr_\perp^2z}{R^2 + \left(\mu + \nu R^2\right)\left(1 + \mu + \nu R^2\right)z^2}\cdot\nonumber\\ & \cdot\left|\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)} + \nu\frac{\sin\left(\varphi\right)\cos\left(\varphi\right)c^2\left(1 - c^2\right)}{\left(\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)\right)^2}\right|. \end{align} \]
Let a point $\mathbf{r}$ with Great coordinates $\left(R, \varphi, \lambda\right)$ be given. We seek its global coordinates $\left(x, y, z\right)^T$. First assume $\lambda = 0$. Squaring Eq. (D.228) yields
\[ \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} \]
Substituting Eq. (D.230) gives
\[ \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 must therefore be solved iteratively. Here, $D\left(\varphi\right)$ and $R$ are constants, and $z^2$ is the only unknown. Then $x^2$ can be obtained from Eq. (D.230). 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} \]
must be replaced, so that one obtains
\[ \begin{align} x &= r_\perp\cos\left(\lambda\right),\\ y &= r_\perp\sin\left(\lambda\right). \end{align} \]
Given a point $\mathbf{r}$ with global coordinates $\left(x, y, z\right)^T$, one seeks its Great coordinates $\left(R, \varphi, \lambda\right)$. Writing Eq. (D.212) with the rescalings (D.222) - (D.222) gives
\[ \begin{align} R^2 = \frac{x^2 + y^2 + \left(1 + \mu\right)z^2}{1 - \nu z^2}. \end{align} \]
The replacement $x^2\to x^2 + y^2$ is used because the general case $\lambda\not = 0$ is assumed. Then $\lambda$ is obtained from
\[ \begin{align} \lambda = \sign\left(y\right)\arccos\left(\frac{x}{\sqrt{x^2 + y^2}}\right). \end{align} \]
To determine $\varphi$, one first solves Eq. (D.228) for $D$ and again replaces $x^2\to x^2 + y^2$:
\[ \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 one sets $R = 1$ in Eq. (D.255), which corresponds to the WGS 84 surface, and solves numerically for $z_S^2$; the sign of $z_S$ is identical to the sign of $z$. One can then use Eq. (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$. Finally,
\[ \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 obtains $\varphi$. We now derive Eq. (D.262). Without loss of generality, consider the plane $\lambda = 0$, hence $y = 0$. The intersection of 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
\[ \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}$. For $z_S\not = 0$, the intersection of $\mathbf{g}$ with the x-axis is given by
\[ \begin{align} 0 = z_S + \gamma \frac{z_S}{c^2} \Rightarrow \gamma = -c^2. \end{align} \]
For the x-coordinate of this point, one obtains
\[ \begin{align} x = x_S\left(1 - c^2\right). \end{align} \]
Therefore,
\[ \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} \]
From Eq. (D.263), it follows that
\[ \begin{align} z_S^2 = c^2\left(1 - x_S^2\right) \end{align} \]
and therefore
\[ \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$, it follows that
\[ \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} \]
Solving Eq. (D.263) for $x_S$ gives
\[ \begin{align} x_S = \sqrt{1 - \frac{z^2}{c^2}}. \end{align} \]
Substituting this into Eq. (D.269) gives
\[ \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} \]
Therefore,
\[ \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}$, Eq. (D.262) follows.
Now one needs 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 one takes the logarithm of this, it follows that
\[ \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 one obtains
\[ \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} \]
\[ \begin{align} \Leftrightarrow\frac{d\ln\left[D\left(\varphi\right)\right]}{d\varphi}&= \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)} + \nu\frac{\sin\left(\varphi\right)\cos\left(\varphi\right)c^2\left(1 - c^2\right)}{\left(\cos^2\left(\varphi\right) + c^2\sin^2\left(\varphi\right)\right)^2}\nonumber\\ & \end{align} \]
An area $A$ on an ellipsoid is usually computed using an equal-area projection. The basic procedure is to transform each vertex $\left(\phi_i,\lambda_i\right)$ to the projection plane, i.e. determine coordinates $\left(x_i,y_i\right)$, and then compute the planar area $A'$. Equal-area projections are exactly those for which $A = A'$ holds.
Here, the Lambert azimuthal projection is used as the equal-area projection.
[title = References]