#!/usr/bin/env python
"""
ref_ellipsoid.py
Written by Tyler Sutterley (12/2022)
Computes parameters for a reference ellipsoid
CALLING SEQUENCE
wgs84 = ref_ellipsoid('WGS84', UNITS='MKS')
INPUT:
refell - reference ellipsoid name
CLK66 = Clarke 1866
GRS67 = Geodetic Reference System 1967 (IAU ellipsoid)
GRS80 = Geodetic Reference System 1980
WGS72 = World Geodetic System 1972
WGS84 = World Geodetic System 1984
ATS77 = Quasi-earth centred ellipsoid for ATS77
NAD27 = North American Datum 1927
NAD83 = North American Datum 1983
INTER = International
KRASS = Krassovsky (USSR)
MAIRY = Modified Airy (Ireland 1965/1975)
TOPEX = TOPEX/POSEIDON ellipsoid
EGM96 = EGM 1996 gravity model
HGH80 = Hughes 1980 Ellipsoid used in some NSIDC data
OPTIONS:
UNITS: output units
MKS: meters, kilograms, seconds
CGS: centimeters, grams, seconds
OUTPUTS:
a: semimajor semi-axis (m)
b: semiminor semi-axis (m)
f: flattening
rad_e: mean radius of ellipsoid having the same volume
rad_p: Polar radius of curvature
C20: Normalized C20 harmonic
J2: Oblateness coefficient
norm_a: Normal gravity at the equator
norm_b: Normal gravity at the pole
U0: Normal potential at the ellipsoid
dk: ratio between gravity at pole versus gravity at equator
m: m parameter (m)
lin_ecc: Linear eccentricity
ecc1: First eccentricity
ecc2: Second eccentricity
area: Area of the ellipsoid
volume: Volume of the ellipsoid
rho_e: Average density
PYTHON DEPENDENCIES:
numpy: Scientific Computing Tools For Python
https://numpy.org
https://numpy.org/doc/stable/user/numpy-for-matlab-users.html
REFERENCE:
Hofmann-Wellenhof and Moritz (2006)
UPDATE HISTORY:
Updated 12/2022: output average Earth's density
Updated 04/2022: updated docstrings to numpy documentation format
Updated 11/2020: added function docstrings
Updated 09/2017: added parameters for Hughes (1980) ellipsoid
Updated 07/2017: output J2 harmonic (in addition to C20)
Updated 06/2017: added parameters for EGM1996. CGS now True/False
Updated 06/2016: using __future__ print function, formatted output
Updated 08/2015: changed sys.exit to raise ValueError
Updated 04/2015: added command line option. little updates to comments
Updated 06/2014: changed message to sys.exit
Updated 02/2014: minor update to if statements
"""
import numpy as np
[docs]
def ref_ellipsoid(refell, UNITS='MKS'):
"""
Computes parameters for a reference ellipsoid
:cite:p:`HofmannWellenhof:2006hy`
Parameters
----------
refell: str
Reference ellipsoid name
- ``'CLK66'``: Clarke 1866
- ``'GRS67'``: Geodetic Reference System 1967
- ``'GRS80'``: Geodetic Reference System 1980
- ``'HGH80'``: Hughes 1980 Ellipsoid
- ``'WGS72'``: World Geodetic System 1972
- ``'WGS84'``: World Geodetic System 1984
- ``'ATS77'``: Quasi-earth centred ellipsoid for ATS77
- ``'NAD27'``: North American Datum 1927
- ``'NAD83'``: North American Datum 1983
- ``'INTER'``: International
- ``'KRASS'``: Krassovsky (USSR)
- ``'MAIRY'``: Modified Airy (Ireland 1965/1975)
- ``'TOPEX'``: TOPEX/POSEIDON ellipsoid
- ``'EGM96'``: EGM 1996 gravity model
UNITS: str, default 'MKS'
Output units
- ``'MKS'``: meters, kilograms, seconds
- ``'CGS'``: centimeters, grams, seconds
Returns
-------
a: float
semimajor semi-axis (m)
b: float
semiminor semi-axis (m)
f: float
flattening
rad_e: float
Mean radius of ellipsoid having the same volume
rad_p: float
Polar radius of curvature
C20: float
Normalized C20 harmonic
J2: float
Oblateness coefficient
norm_a: float
Normal gravity at the equator
norm_b: float
Normal gravity at the pole
U0: float
Normal potential at the ellipsoid
dk: float
ratio between gravity at pole versus gravity at equator
m: float
m parameter (m)
lin_ecc: float
Linear eccentricity
ecc1: float
First eccentricity
ecc2: float
Second eccentricity
area: float
Area of the ellipsoid
volume: float
Volume of the ellipsoid
rho_e: float
Average density
"""
# validate units
assert UNITS in ('MKS', 'CGS')
# set parameters for ellipsoid
if refell.upper() in ('CLK66', 'NAD27'):
# Clarke 1866
a_axis = 6378206.4 # [m] semimajor axis of the ellipsoid
flat = 1.0 / 294.9786982 # flattening of the ellipsoid
elif refell.upper() in ('GRS80', 'NAD83'):
# Geodetic Reference System 1980
# North American Datum 1983
a_axis = 6378135.0 # [m] semimajor axis of the ellipsoid
flat = 1.0 / 298.26 # flattening of the ellipsoid
GM = 3.986005e14 # [m^3/s^2] Geocentric Gravitational Constant
elif refell.upper() == 'GRS67':
# Geodetic Reference System 1967
# International Astronomical Union (IAU ellipsoid)
a_axis = 6378160.0 # [m] semimajor axis of the ellipsoid
flat = 1.0 / 298.247167427 # flattening of the ellipsoid
GM = 3.98603e14 # [m^3/s^2] Geocentric Gravitational Constant
omega = 7292115.1467e-11 # angular velocity of the Earth [rad/s]
elif refell.upper() == 'WGS72':
# World Geodetic System 1972
a_axis = 6378135.0 # [m] semimajor axis of the ellipsoid
flat = 1.0 / 298.26 # flattening of the ellipsoid
elif refell.upper() == 'WGS84':
# World Geodetic System 1984
a_axis = 6378137.0 # [m] semimajor axis of the ellipsoid
flat = 1.0 / 298.257223563 # flattening of the ellipsoid
elif refell.upper() == 'ATS77':
# Quasi-earth centred ellipsoid for ATS77
a_axis = 6378135.0 # [m] semimajor axis of the ellipsoid
flat = 1.0 / 298.257 # flattening of the ellipsoid
elif refell.upper() == 'KRASS':
# Krassovsky (USSR)
a_axis = 6378245.0 # [m] semimajor axis of the ellipsoid
flat = 1.0 / 298.3 # flattening of the ellipsoid
elif refell.upper() == 'INTER':
# International
a_axis = 6378388.0 # [m] semimajor axis of the ellipsoid
flat = 1 / 297.0 # flattening of the ellipsoid
elif refell.upper() == 'MAIRY':
# Modified Airy (Ireland 1965/1975)
a_axis = 6377340.189 # [m] semimajor axis of the ellipsoid
flat = 1 / 299.3249646 # flattening of the ellipsoid
elif refell.upper() == 'TOPEX':
# TOPEX/POSEIDON ellipsoid
a_axis = 6378136.3 # [m] semimajor axis of the ellipsoid
flat = 1.0 / 298.257 # flattening of the ellipsoid
GM = 3.986004415e14 # [m^3/s^2]
elif refell.upper() == 'EGM96':
# EGM 1996 gravity model
a_axis = 6378136.3 # [m] semimajor axis of the ellipsoid
flat = 1.0 / 298.256415099 # flattening of the ellipsoid
GM = 3.986004415e14 # [m^3/s^2]
elif refell.upper() == 'HGH80':
# Hughes 1980 Ellipsoid used in some NSIDC data
a_axis = 6378273.0 # [m] semimajor axis of the ellipsoid
flat = 1.0 / 298.279411123064 # flattening of the ellipsoid
else:
raise ValueError('Incorrect reference ellipsoid Name')
if refell.upper() not in ('GRS80', 'GRS67', 'NAD83', 'TOPEX', 'EGM96'):
# for ellipsoids not listing the Geocentric Gravitational Constant
GM = 3.986004418e14 # [m^3/s^2]
if refell.upper() not in ('GRS67'):
# for ellipsoids not listing the angular velocity of the Earth
omega = 7292115e-11 # [rad/s]
# universal gravitational constant [N*m^2/kg^2]
G = 6.67430e-11
# convert units to CGS
if UNITS == 'CGS':
a_axis *= 100.0
GM *= 1e6
G *= 1000.0 # [dyn*cm^2/g^2]
# DERIVED PARAMETERS:
# mean radius of the Earth having the same volume
# (4pi/3)R^3 = (4pi/3)(a^2)b = (4pi/3)(a^3)(1 - f)
rad_e = a_axis * (1.0 - flat) ** (1.0 / 3.0)
# semiminor axis of the ellipsoid
b_axis = (1.0 - flat) * a_axis # [m]
# Ratio between ellipsoidal axes
ratio = 1.0 - flat
# Polar radius of curvature
pol_rad = a_axis / (1.0 - flat)
# Linear eccentricity
lin_ecc = np.sqrt((2.0 * flat - flat**2) * a_axis**2)
# first numerical eccentricity
ecc1 = lin_ecc / a_axis
# second numerical eccentricity
ecc2 = lin_ecc / b_axis
# m parameter [omega^2*a^2*b/(GM)]
# p. 70, Eqn.(2-137)
mp = omega**2 * ((1 - flat) * a_axis**3) / GM
# q, q_0
# p. 67, Eqn.(2-113)
q = 0.5 * ((1.0 + 3.0 / (ecc2**2)) * np.arctan(ecc2) - 3.0 / ecc2)
q_0 = (
3 * (1.0 + 1.0 / (ecc2**2)) * (1.0 - 1.0 / ecc2 * np.arctan(ecc2)) - 1.0
)
# J_2 p. 75 Eqn.(2-167), p. 76 Eqn.(2-172)
j_2 = (ecc1**2) * (1.0 - 2.0 * mp * ecc2 / (15.0 * q)) / 3.0
# Normalized C20 terms.
# p. 60, Eqn.(2-80)
C20 = -j_2 / np.sqrt(5.0)
# Normal gravity at the equator.
# p. 71, Eqn.(2-141)
ga = GM / (a_axis * b_axis) * (1.0 - mp - mp * ecc2 * q_0 / (6.0 * q))
# Normal gravity at the pole.
# p. 71, Eqn.(2-142)
gb = GM / (a_axis**2.0) * (1.0 + mp * ecc2 * q_0 / (3.0 * q))
# ratio between gravity at pole versus gravity at equator
dk = b_axis * gb / (a_axis * ga) - 1.0
# Normal potential at the ellipsoid
# p. 68, Eqn.(2-123)
U0 = GM / lin_ecc * np.arctan(ecc2) + (1.0 / 3.0) * omega**2 * a_axis**2
# Surface area of the reference ellipsoid [m^2]
area = (
np.pi
* a_axis**2.0
* (2.0 + ((1.0 - ecc1**2) / ecc1) * np.log((1.0 + ecc1) / (1.0 - ecc1)))
)
# Volume of the reference ellipsoid [m^3]
vol = (4.0 * np.pi / 3.0) * (a_axis**3.0) * (1.0 - ecc1**2.0) ** 0.5
# Average density [kg/m^3]
rho_e = GM / (G * vol)
return {
'a': a_axis,
'b': b_axis,
'f': flat,
'rad_p': pol_rad,
'rad_e': rad_e,
'ratio': ratio,
'GM': GM,
'omega': omega,
'C20': C20,
'J2': j_2,
'U0': U0,
'dk': dk,
'norm_a': ga,
'norm_b': gb,
'mp': mp,
'q': q,
'q0': q_0,
'ecc': lin_ecc,
'ecc1': ecc1,
'ecc2': ecc2,
'area': area,
'volume': vol,
'rho_e': rho_e,
}