# Theorerical background¶

## What is a tesseroid anyway?¶

A tesseroid, or spherical prism, is segment of a sphere. It is delimited by:

1. 2 meridians, $$\lambda_1$$ and $$\lambda_2$$

2. 2 parallels, $$\phi_1$$ and $$\phi_2$$

3. 2 spheres of radii $$r_1$$ and $$r_2$$

Original images (licensed CC-BY) at doi:10.6084/m9.figshare.1495537.

The figure bellow shows a tesseroid, the global coordinate system (X, Y, Z), and the local coordinate system ($$x,\ y,\ z$$) of a point P. View of a tesseroid, the integration point Q, the global coordinate system (X, Y, Z), the computation P and it’s local coordinate system ($$x,\ y,\ z$$). $$r,\ \phi,\ \lambda$$ are the radius, latitude, and longitude, respectively, of point P. Original image (licensed CC-BY) at doi:10.6084/m9.figshare.1495525.

The global system has origin on the center of the Earth and Z axis aligned with the Earth’s mean rotation axis. The X and Y axis are contained on the equatorial parallel with X intercepting the mean Greenwich meridian and Y completing a right-handed system.

The local system has origin on the computation point P. It’s $$z$$ axis is oriented along the radial direction and points away from the center of the Earth. The $$x$$ and $$y$$ axis are contained on a plane normal to the $$z$$ axis. $$x$$ points North and $$y$$ East.

The gravitational attraction and gravity gradient tensor of a tesseroid are calculated with respect to the local coordinate system of the computation point P.

Warning

The $$g_z$$ component is an exception to this. In order to conform with the regular convention of z-axis pointing toward the center of the Earth, this component ONLY is calculated with an inverted z axis. This way, gravity anomalies of tesseroids with positive density are positive, not negative.

## Gravitational fields of a tesseroid¶

The gravitational potential of a tesseroid can be calculated using the formula

$V(r,\phi,\lambda) = G \rho \displaystyle\int_{\lambda_1}^{\lambda_2} \displaystyle\int_{\phi_1}^{\phi_2} \displaystyle\int_{r_1}^{r_2} \frac{1}{\ell} \kappa \ d r' d \phi' d \lambda'$

The gravitational attraction can be calculated using the formula (Grombein et al., 2013):

$g_{\alpha}(r,\phi,\lambda) = G \rho \displaystyle\int_{\lambda_1}^{\lambda_2} \displaystyle\int_{\phi_1}^{\phi_2} \displaystyle\int_{r_1}^{r_2} \frac{\Delta_{\alpha}}{\ell^3} \kappa \ d r' d \phi' d \lambda' \ \ \alpha \in \{x,y,z\}$

The gravity gradients can be calculated using the general formula (Grombein et al., 2013):

$g_{\alpha\beta}(r,\phi,\lambda) = G \rho \displaystyle\int_{\lambda_1}^{\lambda_2} \displaystyle\int_{\phi_1}^{\phi_2} \displaystyle\int_{r_1}^{r_2} I_{\alpha\beta}({r'}, {\phi'}, {\lambda'} ) \ d r' d \phi' d \lambda' \ \ \alpha,\beta \in \{x,y,z\}$
$I_{\alpha\beta}({r'}, {\phi'}, {\lambda'}) = \left( \frac{3\Delta_{\alpha} \Delta_{\beta}}{\ell^5} - \frac{\delta_{\alpha\beta}}{\ell^3} \right) \kappa\ \ \ \alpha,\beta \in \{x,y,z\}$

where $$\rho$$ is density, $$\{x, y, z\}$$ correspond to the local coordinate system of the computation point P (see the tesseroid figure), $$\delta_{\alpha\beta}$$ is the Kronecker delta, and

\begin{eqnarray*} \Delta_x &=& r' K_{\phi} \\ \Delta_y &=& r' \cos \phi' \sin(\lambda' - \lambda) \\ \Delta_z &=& r' \cos \psi - r\\ \ell &=& \sqrt{r'^2 + r^2 - 2 r' r \cos \psi} \\ \cos\psi &=& \sin\phi\sin\phi' + \cos\phi\cos\phi' \cos(\lambda' - \lambda) \\ K_{\phi} &=& \cos\phi\sin\phi' - \sin\phi\cos\phi' \cos(\lambda' - \lambda)\\ \kappa &=& {r'}^2 \cos \phi' \end{eqnarray*}

$$\phi$$ is latitude, $$\lambda$$ is longitude, and $$r$$ is radius.

Note

The gravitational attraction and gravity gradient tensor are calculated with respect to $$(x, y, z)$$, the local coordinate system of the computation point P.

## Numerical integration¶

The above integrals are solved using the Gauss-Legendre Quadrature rule (Asgharzadeh et al., 2007):

$g_{\alpha\beta}(r,\phi,\lambda) \approx G \rho \frac{(\lambda_2 - \lambda_1)(\phi_2 - \phi_1)(r_2 - r_1)}{8} \displaystyle\sum_{k=1}^{N^{\lambda}} \displaystyle\sum_{j=1}^{N^{\phi}} \displaystyle\sum_{i=1}^{N^r} W^r_i W^{\phi}_j W^{\lambda}_k I_{\alpha\beta}({r'}_i, {\phi'}_j, {\lambda'}_k ) \ \alpha,\beta \in \{1,2,3\}$

where $$W_i^r$$, $$W_j^{\phi}$$, and $$W_k^{\lambda}$$ are weighting coefficients and $$N^r$$, $$N^{\phi}$$, and $$N^{\lambda}$$ are the number of quadrature nodes (i.e., the order of the quadrature), for the radius, latitude, and longitude, respectively.

Tesseroids implements a modified version the adaptive discretization algorithm of Li et al (2011). This helps guarantee that the numerical integration will achieve a maximum error of 0.1%.

Warning

The integration error may be larger than this if the computation points are closer than 1km of the tesseroids. This effect is more significant in the gravity gradient components.

## Gravitational fields of a prism in spherical coordinates¶

The gravitational potential and its first and second derivatives for the right rectangular prism can be calculated in Cartesian coordinates using the formula of Nagy et al. (2000).

However, several transformations have to made in order to calculate the fields of a prism in a global coordinate system using spherical coordinates (see this figure). View of a right rectangular prism with its corresponding local coordinate system ($$x^*,\ y^*,\ z^*$$), the global coordinate system (X, Y, Z), the computation P and it’s local coordinate system ($$x,\ y,\ z$$). $$r,\ \phi,\ \lambda$$ are the radius, latitude, and longitude, respectively.

The formula of Nagy et al. (2000) require that the computation point be given in the Cartesian coordinates of the prism ($$x^*,\ y^*,\ z^*$$ in this figure). Therefore, we must first transform the spherical coordinates ($$r,\ \phi,\ \lambda$$) of the computation point P into $$x^*,\ y^*,\ z^*$$. This means that we must convert vector $$\bar{e}$$ (from this other figure) to the coordinate system of the prism. We must first obtain vector $$\bar{e}$$ in the global Cartesian coordinates (X, Y, Z):

$\bar{e}^g = \bar{E} - \bar{E}^*$

where $$\bar{e}^g$$ is the vector $$\bar{e}$$ in the global Cartesian coordinates and

$\begin{split}\bar{E} = \begin{bmatrix} r \cos\phi\cos\lambda \\ r \cos\phi\sin\lambda \\ r \sin\phi \end{bmatrix}\end{split}$
$\begin{split}\bar{E}^* = \begin{bmatrix} r^* \cos\phi^*\cos\lambda^* \\ r^* \cos\phi^*\sin\lambda^* \\ r^* \sin\phi^* \end{bmatrix}\end{split}$ The position vectors involved in the coordinate transformations. $$\bar{E}^*$$ is the position vector of point Q in the global coordinate system, $$\bar{E}$$ is the position vector of point P in the global coordinate system, and $$\bar{e}$$ is the position vector of point P in the local coordinate system of the prism ($$x^*,\ y^*,\ z^*$$).

Next, we transform $$\bar{e}^g$$ to the local Cartesian system of the prism by

$\bar{e} = \underbrace{ \bar{\bar{P}}_y \bar{\bar{R}}_y(90^\circ - \phi^*) \bar{\bar{R}}_z(180^\circ - \lambda^*) }_{ \bar{\bar{W}} } \bar{e}^g$

where $$\bar{\bar{P}}_y$$ is a deflection matrix of the y axis, $$\bar{\bar{R}}_y$$ and $$\bar{\bar{R}}_z$$ are counterclockwise rotation matrices around the y and z axis, respectively (see Wolfram MathWorld).

$\begin{split}\bar{\bar{P}}_y = \begin{bmatrix} 1 & 0 & 0\\ 0 & -1 & 0\\ 0 & 0 & 1\\ \end{bmatrix}\end{split}$
$\begin{split}\bar{\bar{R}}_y(\alpha) = \begin{bmatrix} \cos\alpha & 0 & \sin\alpha\\ 0 & 1 & 0\\ -\sin\alpha & 0 & \cos\alpha\\ \end{bmatrix}\end{split}$
$\begin{split}\bar{\bar{R}}_z(\alpha) = \begin{bmatrix} \cos\alpha & -\sin\alpha & 0\\ \sin\alpha & \cos\alpha & 0\\ 0 & 0 & 1\\ \end{bmatrix}\end{split}$
$\begin{split}\bar{\bar{W}} = \begin{bmatrix} \cos(90^\circ - \phi^*)\cos(180^\circ - \lambda^*) & -\cos(90^\circ - \phi^*)\sin(180^\circ - \lambda^*) & \sin(90^\circ - \phi^*) \\ -\sin(180^\circ - \lambda^*) & -\cos(180^\circ - \lambda^*) & 0 \\ -\sin(90^\circ - \phi^*)\cos(180^\circ - \lambda^*) & \sin(90^\circ - \phi^*)\sin(180^\circ - \lambda^*) & \cos(90^\circ - \phi^*) \end{bmatrix}\end{split}$

Which gives us

$\begin{split}\bar{e} = \begin{bmatrix} x\\y\\z \end{bmatrix}\end{split}$

Note

Nagy et al. (2000) use the z axis pointing down, so we still need to invert the sign of $$z$$.

Vector $$\bar{e}$$ can then be used with the Nagy et al. (2000) formula. These formula give us the gravitational attraction and the gravity gradient tensor calculated with respect to the coordinate system of the prism (i.e., $$x^*,\ y^*,\ z^*$$). However, we need them in the coordinate system of the observation point P, where they are measured by GOCE and calculated for the tesseroids. We perform these transformations via the global Cartesian system (tip: the rotation matrices are orthogonal). $$\bar{g}^*$$ is the gravity vector in the coordinate system of the prism, $$\bar{g}^g$$ is the gravity vector in the global coordinate system, and $$\bar{g}$$ is the gravity vector in the coordinate system of computation point P.

$\bar{g}^g = \bar{\bar{R}}_z(\lambda^* - 180^\circ) \bar{\bar{R}}_y(\phi^* - 90^\circ) \bar{\bar{P}}_y \bar{g}^*$
$\bar{g} = \bar{\bar{P}}_y \bar{\bar{R}}_y(90^\circ - \phi) \bar{\bar{R}}_z(180^\circ - \lambda)\bar{g}^g$
$\bar{g} = \bar{\bar{P}}_y \bar{\bar{R}}_y(90^\circ - \phi) \underbrace{ \bar{\bar{R}}_z(180^\circ - \lambda) \bar{\bar{R}}_z(\lambda^* - 180^\circ) }_{ \bar{\bar{R}}_z(\lambda^* - \lambda) } \bar{\bar{R}}_y(\phi^* - 90^\circ) \bar{\bar{P}}_y \bar{g}^*$
$\bar{g} = \bar{\bar{R}} \bar{g}^*$
$\bar{\bar{R}} = \bar{\bar{P}}_y \bar{\bar{R}}_y(90^\circ - \phi) \bar{\bar{R}}_z(\lambda^* - \lambda) \bar{\bar{R}}_y(\phi^* - 90^\circ) \bar{\bar{P}}_y$
$\begin{split}\bar{\bar{R}} = \begin{bmatrix} \cos\beta\cos\alpha\cos\gamma - \sin\alpha\sin\gamma & \sin\beta\cos\alpha & \cos\beta\cos\alpha\sin\gamma + \sin\alpha\cos\gamma \\ -\sin\beta\cos\gamma & \cos\beta & -\sin\beta\sin\gamma \\ -\cos\beta\sin\alpha\cos\gamma - \cos\alpha\sin\gamma & -\sin\beta\sin\alpha & -\cos\beta\sin\alpha\sin\gamma + \cos\alpha\cos\gamma \end{bmatrix}\end{split}$
$\begin{split}\bar{\bar{R}} = \begin{bmatrix} \cos\beta\sin\phi\sin\phi^* + \cos\phi\cos\phi^* & \sin\beta\sin\phi & -\cos\beta\sin\phi\cos\phi^* + \cos\phi\sin\phi^* \\ -\sin\beta\sin\phi^* & \cos\beta & \sin\beta\cos\phi^* \\ -\cos\beta\cos\phi\sin\phi^* + \sin\phi\cos\phi^* & -\sin\beta\cos\phi & \cos\beta\cos\phi\cos\phi^* + \sin\phi\sin\phi^* \end{bmatrix}\end{split}$

where

\begin{eqnarray*} &\alpha = 90^\circ - \phi \\ &\beta = \lambda^* - \lambda \\ &\gamma = \phi^* - 90^\circ \\ &\cos\alpha = \sin\phi \\ &\sin\alpha = \cos\phi \\ &\cos\gamma = \sin\phi^* \\ &\sin\gamma = -\cos\phi^* \end{eqnarray*}

Likewise, transformation for the gravity gradient tensor $$T$$ is

$\bar{\bar{T}} = \bar{\bar{R}} \bar{\bar{T}}^* \bar{\bar{R}}^T$