Logo Tesseroids 1.0: User Manual and API Documentation

src/c/grav_tess.h File Reference

Functions that calculate the gravitational potential and its first and second derivatives for the tesseroid. More...

#include "utils.h"
#include "glq.h"

Go to the source code of this file.

Functions

double calc_tess_model (TESSEROID *model, int size, double lonp, double latp, double rp, GLQ *glq_lon, GLQ *glq_lat, GLQ *glq_r, double(*field)(TESSEROID, double, double, double, GLQ, GLQ, GLQ))
 Calculates the field of a tesseroid model at a given point.
double calc_tess_model_adapt (TESSEROID *model, int size, double lonp, double latp, double rp, GLQ *glq_lon, GLQ *glq_lat, GLQ *glq_r, double(*field)(TESSEROID, double, double, double, GLQ, GLQ, GLQ))
 Adaptatively calculate the field of a tesseroid model at a given point by splitting the tesseroids if necessary to maintain GLQ stability.
double tess_pot (TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r)
 Calculates potential caused by a tesseroid.
double tess_gx (TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r)
 Calculates gx caused by a tesseroid (Grombein et al., 2010).
double tess_gy (TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r)
 Calculates gy caused by a tesseroid (Grombein et al., 2010).
double tess_gz (TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r)
 Calculates gz caused by a tesseroid (Grombein et al., 2010).
double tess_gxx (TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r)
 Calculates gxx caused by a tesseroid (Grombein et al., 2010).
double tess_gxy (TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r)
 Calculates gxy caused by a tesseroid (Grombein et al., 2010).
double tess_gxz (TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r)
 Calculates gxz caused by a tesseroid (Grombein et al., 2010).
double tess_gyy (TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r)
 Calculates gyy caused by a tesseroid (Grombein et al., 2010).
double tess_gyz (TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r)
 Calculates gyz caused by a tesseroid (Grombein et al., 2010).
double tess_gzz (TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r)
 Calculates gzz caused by a tesseroid (Grombein et al., 2010).

Detailed Description

Functions that calculate the gravitational potential and its first and second derivatives for the tesseroid.

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

\[ g_{\alpha\beta}(r_p,\phi_p,\lambda_p) = G \rho \displaystyle\int_{\lambda_1}^{\lambda_2} \displaystyle\int_{\phi_1}^{\phi_2} \displaystyle\int_{r_1}^{r_2} I_{\alpha\beta}\ d r' d \phi' d \lambda' \ \ \alpha,\beta \in \{1,2,3\} \]

\[ I_{\alpha\beta} = \left(\frac{3\Delta x_i \Delta x_j}{\ell^5} - \frac{\delta_{ij}}{\ell^3} \right)\kappa\ \]

and solved using the Gauss-Legendre Quadrature rule (Asgharzadeh et al., 2007):

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

where $ \rho $ is density, the subscripts 1, 2, and 3 should be interpreted as the x, y, and z axis and

\begin{eqnarray*} \Delta x_1 &=& r' K_{\phi} \\ \Delta x_2 &=& r' \cos \phi' \sin(\lambda' - \lambda_p) \\ \Delta x_3 &=& r' \cos \psi - r_p\\ \ell &=& \sqrt{r'^2 + r_p^2 - 2 r' r_p \cos \psi} \\ \cos\psi &=& \sin\phi_p\sin\phi' + \cos\phi_p\cos\phi' \cos(\lambda' - \lambda_p) \\ K_{\phi} &=& \cos\phi_p\sin\phi' - \sin\phi_p\cos\phi' \cos(\lambda' - \lambda_p)\\ \kappa &=& {r'}^2 \cos \phi' \end{eqnarray*}

$ \phi $ is latitude, $ \lambda $ is longitude, $ r $ is radius. The subscript $ p $ is for the computation point.

The derivatives of the potential are made with respect to the local coordinate system x->North, y->East, z->Up (away from center of the Earth).

To maintain the standard convention, only for component gz the z axis is inverted, so a positive density results in positive gz.

Example:

To calculate the gzz component due to a tesseroid on a regular grid.

#include <stdio.h>
#include "glq.h"r
#include "constants.h"
#include "grav_tess.h"

int main()
{
    TESSEROID tess = {1000, 44, 46, -1, 1, MEAN_EARTH_RADIUS - 100000,
                      MEAN_EARTH_RADIUS};
    GLQ *glqlon, *glqlat, *glqr;
    double lon, lat, r = MEAN_EARTH_RADIUS + 1500000, res;
    int order = 8;

    glqlon = glq_new(order, tess.w, tess.e);
    glqlat = glq_new(order, tess.s, tess.n);
    glqr = glq_new(order, tess.r1, tess.r2);

    for(lat = 20; lat <= 70; lat += 0.5)
    {
        for(lon = -25; lon <= 25; lon += 0.5)
        {
            res = tess_gzz(tess, lon, lat, r, *glqlon, *glqlat, *glqr);
            printf("%g %g %g\n", lon, lat, res);
        }
    }

    glq_free(glqlon);
    glq_free(glqlat);
    glq_free(glqr);

    return 0;
}

References

Todo:

Possible speed up: use pointers for weights and nodes

Allow for tesseroids with depth varying density

Author:
Leonardo Uieda
Date:
27 Jan 2011

Function Documentation

double calc_tess_model ( TESSEROID model,
int  size,
double  lonp,
double  latp,
double  rp,
GLQ glq_lon,
GLQ glq_lat,
GLQ glq_r,
double(*)(TESSEROID, double, double, double, GLQ, GLQ, GLQ field 
)

Calculates the field of a tesseroid model at a given point.

Uses a function pointer to call one of the apropriate field calculating functions:

To pass a function pointer to a function use something like:

calc_tess_model(my_model, 10, 0, 10, 1, glqlon, glqlat, glqr, &tess_gx);

This would calculate the gx effect of the model my_model with 10 tesseroids at lon=0 lat=10 r=1.

Will re-use the same GLQ structures, and therefore the same order, for all the tesseroids.

Parameters:
model TESSEROID array defining the model
size number of tesseroids in the model
lonp longitude of the computation point P
latp latitude of the computation point P
rp radial coordinate of the computation point P
glq_lon pointer to GLQ structure used for the longitudinal integration
glq_lat pointer to GLQ structure used for the latitudinal integration
glq_r pointer to GLQ structure used for the radial integration
field pointer to one of the field calculating functions
Returns:
the sum of the fields of all the tesseroids in the model
double calc_tess_model_adapt ( TESSEROID model,
int  size,
double  lonp,
double  latp,
double  rp,
GLQ glq_lon,
GLQ glq_lat,
GLQ glq_r,
double(*)(TESSEROID, double, double, double, GLQ, GLQ, GLQ field 
)

Adaptatively calculate the field of a tesseroid model at a given point by splitting the tesseroids if necessary to maintain GLQ stability.

See calc_tess_model() for more details.

Will re-use the same GLQ structures, and therefore the same order, for all the tesseroids.

Parameters:
model TESSEROID array defining the model
size number of tesseroids in the model
lonp longitude of the computation point P
latp latitude of the computation point P
rp radial coordinate of the computation point P
glq_lon pointer to GLQ structure used for the longitudinal integration
glq_lat pointer to GLQ structure used for the latitudinal integration
glq_r pointer to GLQ structure used for the radial integration
field pointer to one of the field calculating functions
Returns:
the sum of the fields of all the tesseroids in the model
double tess_gx ( TESSEROID  tess,
double  lonp,
double  latp,
double  rp,
GLQ  glq_lon,
GLQ  glq_lat,
GLQ  glq_r 
)

Calculates gx caused by a tesseroid (Grombein et al., 2010).

\[ g_x(r_p,\phi_p,\lambda_p) = G \rho \displaystyle\int_{\lambda_1}^{\lambda_2} \displaystyle\int_{\phi_1}^{\phi_2} \displaystyle\int_{r_1}^{r_2} \frac{r'K_{\phi}}{\ell^3}\kappa \ d r' d \phi' d \lambda' \]

The derivatives of the potential are made with respect to the local coordinate system x->North, y->East, z->out

Input values in SI units and degrees and returns values in mGal!

Use function glq_new() to create the GLQ parameters required. The integration limits should be set to:

  • glq_lon: lower = tess.w and upper = tess.e (in degrees)
  • glq_lat: lower = tess.s and upper = tess.n (in degrees)
  • glq_r: lower = tess.r1 and upper = tess.r2
Parameters:
tess data structure describing the tesseroid
lonp longitude of the computation point P
latp latitude of the computation point P
rp radial coordinate of the computation point P
glq_lon GLQ structure with the nodes, weights and integration limits set for the longitudinal integration
glq_lat GLQ structure with the nodes, weights and integration limits set for the latitudinal integration
glq_r GLQ structure with the nodes, weights and integration limits set for the radial integration
Returns:
field calculated at P
double tess_gxx ( TESSEROID  tess,
double  lonp,
double  latp,
double  rp,
GLQ  glq_lon,
GLQ  glq_lat,
GLQ  glq_r 
)

Calculates gxx caused by a tesseroid (Grombein et al., 2010).

\[ g_{xx}(r_p,\phi_p,\lambda_p) = G \rho \displaystyle\int_{\lambda_1}^{\lambda_2} \displaystyle\int_{\phi_1}^{\phi_2} \displaystyle\int_{r_1}^{r_2} \frac{3(r' K_{\phi})^2 - \ell^2}{\ell^5}\kappa \ d r' d \phi' d \lambda' \]

The derivatives of the potential are made with respect to the local coordinate system x->North, y->East, z->out

Input values in SI units and degrees and returns values in Eotvos!

Use function glq_new() to create the GLQ parameters required. The integration limits should be set to:

  • glq_lon: lower = tess.w and upper = tess.e (in degrees)
  • glq_lat: lower = tess.s and upper = tess.n (in degrees)
  • glq_r: lower = tess.r1 and upper = tess.r2
Parameters:
tess data structure describing the tesseroid
lonp longitude of the computation point P
latp latitude of the computation point P
rp radial coordinate of the computation point P
glq_lon GLQ structure with the nodes, weights and integration limits set for the longitudinal integration
glq_lat GLQ structure with the nodes, weights and integration limits set for the latitudinal integration
glq_r GLQ structure with the nodes, weights and integration limits set for the radial integration
Returns:
field calculated at P
double tess_gxy ( TESSEROID  tess,
double  lonp,
double  latp,
double  rp,
GLQ  glq_lon,
GLQ  glq_lat,
GLQ  glq_r 
)

Calculates gxy caused by a tesseroid (Grombein et al., 2010).

\[ g_{xy}(r_p,\phi_p,\lambda_p) = G \rho \displaystyle\int_{\lambda_1}^{\lambda_2} \displaystyle\int_{\phi_1}^{\phi_2} \displaystyle\int_{r_1}^{r_2} \frac{3{r'}^2 K_{\phi}\cos\phi'\sin(\lambda' - \lambda_p)}{\ell^5} \kappa \ d r' d \phi' d \lambda' \]

The derivatives of the potential are made with respect to the local coordinate system x->North, y->East, z->out

Input values in SI units and degrees and returns values in Eotvos!

Use function glq_new() to create the GLQ parameters required. The integration limits should be set to:

  • glq_lon: lower = tess.w and upper = tess.e (in degrees)
  • glq_lat: lower = tess.s and upper = tess.n (in degrees)
  • glq_r: lower = tess.r1 and upper = tess.r2
Parameters:
tess data structure describing the tesseroid
lonp longitude of the computation point P
latp latitude of the computation point P
rp radial coordinate of the computation point P
glq_lon GLQ structure with the nodes, weights and integration limits set for the longitudinal integration
glq_lat GLQ structure with the nodes, weights and integration limits set for the latitudinal integration
glq_r GLQ structure with the nodes, weights and integration limits set for the radial integration
Returns:
field calculated at P
double tess_gxz ( TESSEROID  tess,
double  lonp,
double  latp,
double  rp,
GLQ  glq_lon,
GLQ  glq_lat,
GLQ  glq_r 
)

Calculates gxz caused by a tesseroid (Grombein et al., 2010).

\[ g_{xz}(r_p,\phi_p,\lambda_p) = G \rho \displaystyle\int_{\lambda_1}^{\lambda_2} \displaystyle\int_{\phi_1}^{\phi_2} \displaystyle\int_{r_1}^{r_2} \frac{3 r' K_{\phi}(r' \cos\psi - r_p)}{\ell^5}\kappa \ d r' d \phi' d \lambda' \]

The derivatives of the potential are made with respect to the local coordinate system x->North, y->East, z->out

Input values in SI units and degrees and returns values in Eotvos!

Use function glq_new() to create the GLQ parameters required. The integration limits should be set to:

  • glq_lon: lower = tess.w and upper = tess.e (in degrees)
  • glq_lat: lower = tess.s and upper = tess.n (in degrees)
  • glq_r: lower = tess.r1 and upper = tess.r2
Parameters:
tess data structure describing the tesseroid
lonp longitude of the computation point P
latp latitude of the computation point P
rp radial coordinate of the computation point P
glq_lon GLQ structure with the nodes, weights and integration limits set for the longitudinal integration
glq_lat GLQ structure with the nodes, weights and integration limits set for the latitudinal integration
glq_r GLQ structure with the nodes, weights and integration limits set for the radial integration
Returns:
field calculated at P
double tess_gy ( TESSEROID  tess,
double  lonp,
double  latp,
double  rp,
GLQ  glq_lon,
GLQ  glq_lat,
GLQ  glq_r 
)

Calculates gy caused by a tesseroid (Grombein et al., 2010).

\[ g_y(r_p,\phi_p,\lambda_p) = G \rho \displaystyle\int_{\lambda_1}^{\lambda_2} \displaystyle\int_{\phi_1}^{\phi_2} \displaystyle\int_{r_1}^{r_2} \frac{r'\cos\phi'\sin(\lambda'-\lambda)}{\ell^3}\kappa \ d r' d \phi' d \lambda' \]

The derivatives of the potential are made with respect to the local coordinate system x->North, y->East, z->out

Input values in SI units and degrees and returns values in mGal!

Use function glq_new() to create the GLQ parameters required. The integration limits should be set to:

  • glq_lon: lower = tess.w and upper = tess.e (in degrees)
  • glq_lat: lower = tess.s and upper = tess.n (in degrees)
  • glq_r: lower = tess.r1 and upper = tess.r2
Parameters:
tess data structure describing the tesseroid
lonp longitude of the computation point P
latp latitude of the computation point P
rp radial coordinate of the computation point P
glq_lon GLQ structure with the nodes, weights and integration limits set for the longitudinal integration
glq_lat GLQ structure with the nodes, weights and integration limits set for the latitudinal integration
glq_r GLQ structure with the nodes, weights and integration limits set for the radial integration
Returns:
field calculated at P
double tess_gyy ( TESSEROID  tess,
double  lonp,
double  latp,
double  rp,
GLQ  glq_lon,
GLQ  glq_lat,
GLQ  glq_r 
)

Calculates gyy caused by a tesseroid (Grombein et al., 2010).

\[ g_{yy}(r_p,\phi_p,\lambda_p) = G \rho \displaystyle\int_{\lambda_1}^{\lambda_2} \displaystyle\int_{\phi_1}^{\phi_2} \displaystyle\int_{r_1}^{r_2} \frac{3(r'\cos\phi'\sin(\lambda' - \lambda_p))^2 - \ell^2}{\ell^5} \kappa \ d r' d \phi' d \lambda' \]

The derivatives of the potential are made with respect to the local coordinate system x->North, y->East, z->out

Input values in SI units and degrees and returns values in Eotvos!

Use function glq_new() to create the GLQ parameters required. The integration limits should be set to:

  • glq_lon: lower = tess.w and upper = tess.e (in degrees)
  • glq_lat: lower = tess.s and upper = tess.n (in degrees)
  • glq_r: lower = tess.r1 and upper = tess.r2
Parameters:
tess data structure describing the tesseroid
lonp longitude of the computation point P
latp latitude of the computation point P
rp radial coordinate of the computation point P
glq_lon GLQ structure with the nodes, weights and integration limits set for the longitudinal integration
glq_lat GLQ structure with the nodes, weights and integration limits set for the latitudinal integration
glq_r GLQ structure with the nodes, weights and integration limits set for the radial integration
Returns:
field calculated at P
double tess_gyz ( TESSEROID  tess,
double  lonp,
double  latp,
double  rp,
GLQ  glq_lon,
GLQ  glq_lat,
GLQ  glq_r 
)

Calculates gyz caused by a tesseroid (Grombein et al., 2010).

\[ g_{yz}(r_p,\phi_p,\lambda_p) = G \rho \displaystyle\int_{\lambda_1}^{\lambda_2} \displaystyle\int_{\phi_1}^{\phi_2} \displaystyle\int_{r_1}^{r_2} \frac{3 r' \cos\phi' \sin(\lambda' - \lambda_p)(r'\cos\psi - r_p)}{\ell^5} \kappa \ d r' d \phi' d \lambda' \]

The derivatives of the potential are made with respect to the local coordinate system x->North, y->East, z->out

Input values in SI units and degrees and returns values in Eotvos!

Use function glq_new() to create the GLQ parameters required. The integration limits should be set to:

  • glq_lon: lower = tess.w and upper = tess.e (in degrees)
  • glq_lat: lower = tess.s and upper = tess.n (in degrees)
  • glq_r: lower = tess.r1 and upper = tess.r2
Parameters:
tess data structure describing the tesseroid
lonp longitude of the computation point P
latp latitude of the computation point P
rp radial coordinate of the computation point P
glq_lon GLQ structure with the nodes, weights and integration limits set for the longitudinal integration
glq_lat GLQ structure with the nodes, weights and integration limits set for the latitudinal integration
glq_r GLQ structure with the nodes, weights and integration limits set for the radial integration
Returns:
field calculated at P
double tess_gz ( TESSEROID  tess,
double  lonp,
double  latp,
double  rp,
GLQ  glq_lon,
GLQ  glq_lat,
GLQ  glq_r 
)

Calculates gz caused by a tesseroid (Grombein et al., 2010).

\[ g_z(r_p,\phi_p,\lambda_p) = G \rho \displaystyle\int_{\lambda_1}^{\lambda_2} \displaystyle\int_{\phi_1}^{\phi_2} \displaystyle\int_{r_1}^{r_2} \frac{r'\cos\psi - r_p}{\ell^3}\kappa \ d r' d \phi' d \lambda' \]

The derivatives of the potential are made with respect to the local coordinate system x->North, y->East, z->out

Input values in SI units and degrees and returns values in mGal!

Use function glq_new() to create the GLQ parameters required. The integration limits should be set to:

  • glq_lon: lower = tess.w and upper = tess.e (in degrees)
  • glq_lat: lower = tess.s and upper = tess.n (in degrees)
  • glq_r: lower = tess.r1 and upper = tess.r2
Parameters:
tess data structure describing the tesseroid
lonp longitude of the computation point P
latp latitude of the computation point P
rp radial coordinate of the computation point P
glq_lon GLQ structure with the nodes, weights and integration limits set for the longitudinal integration
glq_lat GLQ structure with the nodes, weights and integration limits set for the latitudinal integration
glq_r GLQ structure with the nodes, weights and integration limits set for the radial integration
Returns:
field calculated at P
double tess_gzz ( TESSEROID  tess,
double  lonp,
double  latp,
double  rp,
GLQ  glq_lon,
GLQ  glq_lat,
GLQ  glq_r 
)

Calculates gzz caused by a tesseroid (Grombein et al., 2010).

\[ g_{zz}(r_p,\phi_p,\lambda_p) = G \rho \displaystyle\int_{\lambda_1}^{\lambda_2} \displaystyle\int_{\phi_1}^{\phi_2} \displaystyle\int_{r_1}^{r_2} \frac{3(r'\cos\psi-r_p)^2 - \ell^2}{\ell^5}\kappa \ d r' d \phi' d \lambda' \]

The derivatives of the potential are made with respect to the local coordinate system x->North, y->East, z->out

Input values in SI units and degrees and returns values in Eotvos!

Use function glq_new() to create the GLQ parameters required. The integration limits should be set to:

  • glq_lon: lower = tess.w and upper = tess.e (in degrees)
  • glq_lat: lower = tess.s and upper = tess.n (in degrees)
  • glq_r: lower = tess.r1 and upper = tess.r2
Parameters:
tess data structure describing the tesseroid
lonp longitude of the computation point P
latp latitude of the computation point P
rp radial coordinate of the computation point P
glq_lon GLQ structure with the nodes, weights and integration limits set for the longitudinal integration
glq_lat GLQ structure with the nodes, weights and integration limits set for the latitudinal integration
glq_r GLQ structure with the nodes, weights and integration limits set for the radial integration
Returns:
field calculated at P
double tess_pot ( TESSEROID  tess,
double  lonp,
double  latp,
double  rp,
GLQ  glq_lon,
GLQ  glq_lat,
GLQ  glq_r 
)

Calculates potential caused by a tesseroid.

\[ V(r_p,\phi_p,\lambda_p) = 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' \]

Input and output values in SI units and degrees!

Use function glq_new() to create the GLQ parameters required. The integration limits should be set to:

  • glq_lon: lower = tess.w and upper = tess.e (in degrees)
  • glq_lat: lower = tess.s and upper = tess.n (in degrees)
  • glq_r: lower = tess.r1 and upper = tess.r2
Parameters:
tess data structure describing the tesseroid
lonp longitude of the computation point P
latp latitude of the computation point P
rp radial coordinate of the computation point P
glq_lon GLQ structure with the nodes, weights and integration limits set for the longitudinal integration
glq_lat GLQ structure with the nodes, weights and integration limits set for the latitudinal integration
glq_r GLQ structure with the nodes, weights and integration limits set for the radial integration
Returns:
field calculated at P
Generated on Tue Apr 26 12:17:07 2011 for Tesseroids 1.0: User manual and API documentation by doxygen 1.6.3