GeoCotrans.jl

Transformations between geocentric reference frames, IGRF-14 geomagnetic field, and field line tracing.

Conventions

• Positions: Cartesian (x, y, z) in Earth radii (R🜨 = 6371.2 km); spherical (r [Re], θ colatitude [rad], φ east longitude [rad]); GDZ (lat [deg], lon [deg], alt [km]).
• Field models take positions in units of their reference radius (Re for IGRF). With Unitful loaded, positions may also be lengths (r, or Cartesian x, y, z).
• Batched input: a vector of 3-vectors, or a matrix of 3-vectors stacked along dim (default: the container's time dimension, else 2, i.e. 3×N) with one t for all samples or a length-N vector of times; a time-keyed matrix (e.g. DimArray) may omit t.

Reference frames

• GEO: Geocentric Geographic (GEO) reference frame.
• GSM: Geocentric Solar Magnetospheric (GSM) reference frame.
• GSE: Geocentric Solar Ecliptic (GSE) reference frame.
• GEI: Geocentric Equatorial Inertial (GEI) reference frame.
• MAG: Geomagnetic (MAG) reference frame.
• SM: Solar Magnetic (SM) reference frame.

GEO(x, y, z) etc. build a CoordinateVector tagged with its frame. Representations Cartesian3, Spherical, Geodetic follow the conventions above; getcsys gives (frame, representation). GDZ (lat, lon, alt) is the geodetic representation of GEO. gdz2sph converts it to spherical.

Coordinate transformations

• transform(from => to, x, t): generic frame transform.
• rotation(from, to, t): the underlying 3×3 rotation matrix.
• a2b(x, t) shortcuts for every pair of frames, e.g. gse2gei, geo2gsm.
• Representation: car2sph, sph2car (radians); car2sphd, sphd2car (degrees).

Geomagnetic field

• IGRF(): IGRF-14 model, called as model(𝐫, t; from, to).
• get_mlt(xGEO, t): magnetic local time [h].

Field line tracing

Requires an OrdinaryDiffEq solver package to be loaded (e.g. using OrdinaryDiffEqTsit5).

• trace(pos, t, solver; model, dir, r0, rlim, maxs): integrate a field line.
• find_magequator(pos, t, solver): |B| minimum along the field line.

References

• Coordinate transformations between geocentric systems
• IGRF - NOAA/NCEI, IGRF-14 evaluation

Transformations between geocentric reference frames, IGRF-14 geomagnetic field, and field line tracing.

Conventions

• Positions: Cartesian (x, y, z) in Earth radii (R🜨 = 6371.2 km); spherical (r [Re], θ colatitude [rad], φ east longitude [rad]); GDZ (lat [deg], lon [deg], alt [km]).
• Field models take positions in units of their reference radius (Re for IGRF). With Unitful loaded, positions may also be lengths (r, or Cartesian x, y, z).
• Batched input: a vector of 3-vectors, or a matrix of 3-vectors stacked along dim (default: the container's time dimension, else 2, i.e. 3×N) with one t for all samples or a length-N vector of times; a time-keyed matrix (e.g. DimArray) may omit t.

Reference frames

• GEO: Geocentric Geographic (GEO) reference frame.
• GSM: Geocentric Solar Magnetospheric (GSM) reference frame.
• GSE: Geocentric Solar Ecliptic (GSE) reference frame.
• GEI: Geocentric Equatorial Inertial (GEI) reference frame.
• MAG: Geomagnetic (MAG) reference frame.
• SM: Solar Magnetic (SM) reference frame.

GEO(x, y, z) etc. build a CoordinateVector tagged with its frame. Representations Cartesian3, Spherical, Geodetic follow the conventions above; getcsys gives (frame, representation). GDZ (lat, lon, alt) is the geodetic representation of GEO. gdz2sph converts it to spherical.

Coordinate transformations

• transform(from => to, x, t): generic frame transform.
• rotation(from, to, t): the underlying 3×3 rotation matrix.
• a2b(x, t) shortcuts for every pair of frames, e.g. gse2gei, geo2gsm.
• Representation: car2sph, sph2car (radians); car2sphd, sphd2car (degrees).

Geomagnetic field

• IGRF(): IGRF-14 model, called as model(𝐫, t; from, to).
• get_mlt(xGEO, t): magnetic local time [h].

Field line tracing

Requires an OrdinaryDiffEq solver package to be loaded (e.g. using OrdinaryDiffEqTsit5).

• trace(pos, t, solver; model, dir, r0, rlim, maxs): integrate a field line.
• find_magequator(pos, t, solver): |B| minimum along the field line.

References

• Coordinate transformations between geocentric systems
• IGRF - NOAA/NCEI, IGRF-14 evaluation

CoordinateVector{F, R}(x, y, z) <: FieldVector{3}

FieldVector tagged with reference frame F and representation R (Cartesian3 by default).

getcsys(v) returns (F(), R()).

Geocentric Equatorial Inertial (GEI) reference frame.

Z-axis is parallel to Earth's rotation axis (positive northward). X-axis points toward the vernal equinox (first point of Aries) of date. As in SPEDAS, mean and true of date (differing by nutation, < 20″, under the 0.01° accuracy) are not distinguished; J2000 maps to true of date.

Geocentric Geographic (GEO) reference frame.

Z-axis is parallel to Earth's rotation axis (positive northward). X-axis passes through the intersection of the equator and the Greenwich meridian. Rotates with Earth.

Geocentric Solar Ecliptic (GSE) reference frame.

X points sunward from Earth's center. Z-axis is perpendicular to the ecliptic plane (positive north).

Geocentric Solar Magnetospheric (GSM) reference frame.

X points sunward from Earth's center. The X-Z plane is defined to contain Earth's dipole axis (positive North).

IGRF()

International Geomagnetic Reference Field (IGRF-14): Earth's main field with coefficients at 5-year epochs from 1965 to 2030, linearly interpolated in time. Valid for years 1965–2030 (errors outside).

(m::IGRF)(𝐫, t; from, to) -> B
(m::IGRF)(r, θ, φ, t; kw...)
(m::IGRF)(R::AbstractMatrix, [times]; dim, kw...)

Magnetic field [nT] at position(s) 𝐫 and time(s) t.

Keyword Arguments

• from: input coordinate system, a frame (GSM/GSM()) or a (frame, representation) tuple. Inferred from 𝐫 for CoordinateVector; defaults to (GEO(), Spherical()).
• to: output coordinate system, defaults to from. Geodetic input yields (Be, Bn, Bu) in East-North-Up.

Examples

m = IGRF()
t = DateTime(2021, 3, 28)
m(1.0, deg2rad(45), deg2rad(45), t)          # (Br, Bθ, Bφ) in spherical GEO
m(GDZ(60.39, 5.32, 0), t)                     # (Be, Bn, Bu), same as m([60.39, 5.32, 0], t; from = GDZ())
m(GSM(3.0, 0.0, 1.0), t)                      # Cartesian GSM in, Cartesian GSM out
m(rand(3, 6), Date.(2015:2020); from = GSM())   # 3×6 matrix, one time per column

J2000 Geocentric Equatorial (J2000, EME2000) reference frame.

Like GEI, but fixed at the mean equator and equinox of the J2000.0 epoch (2000-01-01T12:00 TT). Inertial.

Geomagnetic (MAG) reference frame.

Z-axis is parallel to Earth's magnetic dipole axis (positive northward). Y-axis is perpendicular to the plane containing the dipole and Earth's rotation axis.

Solar Magnetic (SM) reference frame.

Z-axis is parallel to Earth's magnetic dipole axis (positive northward). The X-Z plane contains the Earth-Sun direction.

GDZ(ϕ, λ, h = 0)
GDZ(; lat, lon, alt = 0)

Geodetic coordinate system with:

• ϕ: latitude (north/south)
• λ: longitude (east/west)
• h: ellipsoidal height [km]

https://www.wikipedia.org/wiki/Geodetic_coordinates

Geodetic (GDZ) coordinate system (latitude [deg], longitude [deg], altitude [km]).

Defined using a reference ellipsoid. Both the altitude and latitude depend on the ellipsoid used. GeoCotrans uses the WGS84 reference ellipsoid.

car2sph(x, y, z) -> [r, θ, φ]
car2sph(Bx, By, Bz, x, y, z) -> [Br, Bθ, Bφ]

Cartesian position, or vector B at that position, to spherical, with θ ∈ [0, π] and φ ∈ [0, 2π) in radians (φ = 0 on the polar axis). Also car2sph(xyz), car2sph(B, xyz); likewise for sph2car, car2sphd, sphd2car.

car2sphd(x, y, z) -> [r, θ, φ]

car2sph with angles in degrees.

gdz2sph(lat, lon, alt) -> [r, θ, φ]

Geodetic (lat [deg], lon [deg], alt [km]) to spherical GEO (r [Re], θ [rad], φ [rad]).

gei2geo(x, t)

Transforms coordinate(s) x from Geocentric Equatorial Inertial (GEI) to Geographic (GEO) reference frame at time(s) t.

Equivalent to transform(GEI => GEO, x, t).

gei2gse(x, t)

Transforms coordinate(s) x from Geocentric Equatorial Inertial (GEI) to Geocentric Solar Ecliptic (GSE) reference frame at time(s) t.

Equivalent to transform(GEI => GSE, x, t).

gei2gsm(x, t)

Transforms coordinate(s) x from Geocentric Equatorial Inertial (GEI) to Geocentric Solar Magnetospheric (GSM) reference frame at time(s) t.

Equivalent to transform(GEI => GSM, x, t).

gei2j2000(x, t)

Transforms coordinate(s) x from Geocentric Equatorial Inertial (GEI) to J2000 Geocentric Equatorial (J2000) reference frame at time(s) t.

Equivalent to transform(GEI => J2000, x, t).

gei2mag(x, t)

Transforms coordinate(s) x from Geocentric Equatorial Inertial (GEI) to Geomagnetic (MAG) reference frame at time(s) t.

Equivalent to transform(GEI => MAG, x, t).

gei2sm(x, t)

Transforms coordinate(s) x from Geocentric Equatorial Inertial (GEI) to Solar Magnetic (SM) reference frame at time(s) t.

Equivalent to transform(GEI => SM, x, t).

geo2gei(x, t)

Transforms coordinate(s) x from Geographic (GEO) to Geocentric Equatorial Inertial (GEI) reference frame at time(s) t.

Equivalent to transform(GEO => GEI, x, t).

geo2gse(x, t)

Transforms coordinate(s) x from Geographic (GEO) to Geocentric Solar Ecliptic (GSE) reference frame at time(s) t.

Equivalent to transform(GEO => GSE, x, t).

geo2gsm(x, t)

Transforms coordinate(s) x from Geographic (GEO) to Geocentric Solar Magnetospheric (GSM) reference frame at time(s) t.

Equivalent to transform(GEO => GSM, x, t).

geo2j2000(x, t)

Transforms coordinate(s) x from Geographic (GEO) to J2000 Geocentric Equatorial (J2000) reference frame at time(s) t.

Equivalent to transform(GEO => J2000, x, t).

geo2mag(x, t)

Transforms coordinate(s) x from Geographic (GEO) to Geomagnetic (MAG) reference frame at time(s) t.

Equivalent to transform(GEO => MAG, x, t).

geo2sm(x, t)

Transforms coordinate(s) x from Geographic (GEO) to Solar Magnetic (SM) reference frame at time(s) t.

Equivalent to transform(GEO => SM, x, t).

get_igrf_coeffs(time)

Get IGRF-14 coefficients for a given time.

Similar to IRBEM implementation, but with higher precision (IRBEM uses year as the time unit).

get_mlt(xGEO, time)
get_mlt(xGEO::AbstractMatrix, [times]; dim)

Magnetic local time in hours, [0, 24), of Cartesian GEO position(s) xGEO.

getcsys(x)

Coordinate system (frame, representation) of a position or field model where its evalmodel works in; nothing if unknown.

gse2gei(x, t)

Transforms coordinate(s) x from Geocentric Solar Ecliptic (GSE) to Geocentric Equatorial Inertial (GEI) reference frame at time(s) t.

Equivalent to transform(GSE => GEI, x, t).

gse2geo(x, t)

Transforms coordinate(s) x from Geocentric Solar Ecliptic (GSE) to Geographic (GEO) reference frame at time(s) t.

Equivalent to transform(GSE => GEO, x, t).

gse2gsm(x, t)

Transforms coordinate(s) x from Geocentric Solar Ecliptic (GSE) to Geocentric Solar Magnetospheric (GSM) reference frame at time(s) t.

Equivalent to transform(GSE => GSM, x, t).

gse2j2000(x, t)

Transforms coordinate(s) x from Geocentric Solar Ecliptic (GSE) to J2000 Geocentric Equatorial (J2000) reference frame at time(s) t.

Equivalent to transform(GSE => J2000, x, t).

gse2mag(x, t)

Transforms coordinate(s) x from Geocentric Solar Ecliptic (GSE) to Geomagnetic (MAG) reference frame at time(s) t.

Equivalent to transform(GSE => MAG, x, t).

gse2sm(x, t)

Transforms coordinate(s) x from Geocentric Solar Ecliptic (GSE) to Solar Magnetic (SM) reference frame at time(s) t.

Equivalent to transform(GSE => SM, x, t).

gsm2gei(x, t)

Transforms coordinate(s) x from Geocentric Solar Magnetospheric (GSM) to Geocentric Equatorial Inertial (GEI) reference frame at time(s) t.

Equivalent to transform(GSM => GEI, x, t).

gsm2geo(x, t)

Transforms coordinate(s) x from Geocentric Solar Magnetospheric (GSM) to Geographic (GEO) reference frame at time(s) t.

Equivalent to transform(GSM => GEO, x, t).

gsm2gse(x, t)

Transforms coordinate(s) x from Geocentric Solar Magnetospheric (GSM) to Geocentric Solar Ecliptic (GSE) reference frame at time(s) t.

Equivalent to transform(GSM => GSE, x, t).

gsm2j2000(x, t)

Transforms coordinate(s) x from Geocentric Solar Magnetospheric (GSM) to J2000 Geocentric Equatorial (J2000) reference frame at time(s) t.

Equivalent to transform(GSM => J2000, x, t).

gsm2mag(x, t)

Transforms coordinate(s) x from Geocentric Solar Magnetospheric (GSM) to Geomagnetic (MAG) reference frame at time(s) t.

Equivalent to transform(GSM => MAG, x, t).

gsm2sm(x, t)

Transforms coordinate(s) x from Geocentric Solar Magnetospheric (GSM) to Solar Magnetic (SM) reference frame at time(s) t.

Equivalent to transform(GSM => SM, x, t).

Alias for IGRF()(args...; kw...); see IGRF.

igrf_B(r, θ, φ, t; max_degree=IGRF_degree) -> [Br, Bθ, Bφ]

Calculate IGRF model components in geocentric coordinates (r [Re], θ [rad], φ [rad]) at time t, where r is the radius, θ the colatitude, and φ the longitude.

j20002gei(x, t)

Transforms coordinate(s) x from J2000 Geocentric Equatorial (J2000) to Geocentric Equatorial Inertial (GEI) reference frame at time(s) t.

Equivalent to transform(J2000 => GEI, x, t).

j20002geo(x, t)

Transforms coordinate(s) x from J2000 Geocentric Equatorial (J2000) to Geographic (GEO) reference frame at time(s) t.

Equivalent to transform(J2000 => GEO, x, t).

j20002gse(x, t)

Transforms coordinate(s) x from J2000 Geocentric Equatorial (J2000) to Geocentric Solar Ecliptic (GSE) reference frame at time(s) t.

Equivalent to transform(J2000 => GSE, x, t).

j20002gsm(x, t)

Transforms coordinate(s) x from J2000 Geocentric Equatorial (J2000) to Geocentric Solar Magnetospheric (GSM) reference frame at time(s) t.

Equivalent to transform(J2000 => GSM, x, t).

j20002mag(x, t)

Transforms coordinate(s) x from J2000 Geocentric Equatorial (J2000) to Geomagnetic (MAG) reference frame at time(s) t.

Equivalent to transform(J2000 => MAG, x, t).

j20002sm(x, t)

Transforms coordinate(s) x from J2000 Geocentric Equatorial (J2000) to Solar Magnetic (SM) reference frame at time(s) t.

Equivalent to transform(J2000 => SM, x, t).

mag2gei(x, t)

Transforms coordinate(s) x from Geomagnetic (MAG) to Geocentric Equatorial Inertial (GEI) reference frame at time(s) t.

Equivalent to transform(MAG => GEI, x, t).

mag2geo(x, t)

Transforms coordinate(s) x from Geomagnetic (MAG) to Geographic (GEO) reference frame at time(s) t.

Equivalent to transform(MAG => GEO, x, t).

mag2gse(x, t)

Transforms coordinate(s) x from Geomagnetic (MAG) to Geocentric Solar Ecliptic (GSE) reference frame at time(s) t.

Equivalent to transform(MAG => GSE, x, t).

mag2gsm(x, t)

Transforms coordinate(s) x from Geomagnetic (MAG) to Geocentric Solar Magnetospheric (GSM) reference frame at time(s) t.

Equivalent to transform(MAG => GSM, x, t).

mag2j2000(x, t)

Transforms coordinate(s) x from Geomagnetic (MAG) to J2000 Geocentric Equatorial (J2000) reference frame at time(s) t.

Equivalent to transform(MAG => J2000, x, t).

mag2sm(x, t)

Transforms coordinate(s) x from Geomagnetic (MAG) to Solar Magnetic (SM) reference frame at time(s) t.

Equivalent to transform(MAG => SM, x, t).

rotation(from => to, t)

Rotation matrix taking vectors from reference frame from to to at time t, e.g. rotation(GEI() => GSM(), t).

sm2gei(x, t)

Transforms coordinate(s) x from Solar Magnetic (SM) to Geocentric Equatorial Inertial (GEI) reference frame at time(s) t.

Equivalent to transform(SM => GEI, x, t).

sm2geo(x, t)

Transforms coordinate(s) x from Solar Magnetic (SM) to Geographic (GEO) reference frame at time(s) t.

Equivalent to transform(SM => GEO, x, t).

sm2gse(x, t)

Transforms coordinate(s) x from Solar Magnetic (SM) to Geocentric Solar Ecliptic (GSE) reference frame at time(s) t.

Equivalent to transform(SM => GSE, x, t).

sm2gsm(x, t)

Transforms coordinate(s) x from Solar Magnetic (SM) to Geocentric Solar Magnetospheric (GSM) reference frame at time(s) t.

Equivalent to transform(SM => GSM, x, t).

sm2j2000(x, t)

Transforms coordinate(s) x from Solar Magnetic (SM) to J2000 Geocentric Equatorial (J2000) reference frame at time(s) t.

Equivalent to transform(SM => J2000, x, t).

sm2mag(x, t)

Transforms coordinate(s) x from Solar Magnetic (SM) to Geomagnetic (MAG) reference frame at time(s) t.

Equivalent to transform(SM => MAG, x, t).

sph2car(r, θ, φ) -> [x, y, z]
sph2car(Br, Bθ, Bφ, r, θ, φ) -> [Bx, By, Bz]

Inverse of car2sph.

sphd2car(r, θ, φ) -> [x, y, z]

sph2car with angles in degrees.

transform(from => to, x, t; dim)
transform(from => to, A; dim)
transform(to, x::CoordinateVector, t)

Transform x from reference frame from to to at time(s) t, e.g. transform(GEO() => GSM(), x, t). from can be omitted for x::CoordinateVector (e.g. GEO(x, y, z)). from/to may carry a representation, e.g. GDZ() => GSM or GSM => (GSE(), Spherical()); otherwise components are Cartesian.

x may be a 3-vector or a batched matrix (see GeoCotrans conventions).

Trace magnetic field lines using ODE solvers.

Requires loading SciML ODE solver package (e.g., OrdinaryDiffEqTsit5) before use.

API

• trace, find_magequator

find_magequator(pos, t, solver; kwargs...) -> (; pos, Bmin, s) or nothing

Find the magnetic equator, the |B| minimum along the field line through pos. Traces downhill in |B| and stops at the root of d|B|/ds (a ContinuousCallback), so the result is the nearest local minimum; s is the signed arc length from pos (positive along B). Returns nothing if a boundary (r0, rlim, maxs) is reached first (open field line).

Keyword arguments are those of trace; raise rlim for nightside lines in external models, whose equator can lie beyond 10 Re.

Example

using GeoCotrans, OrdinaryDiffEqTsit5, Dates
eq = find_magequator(GEO(3.0, 0.0, 0.5), DateTime(2020, 1, 1), Tsit5())
eq.pos, eq.Bmin

trace(pos, t, solver; kwargs...) :: ODESolution

Trace a magnetic field line through pos (in from) at time t with a SciML ODE solver, parameterized by arc length.

Keyword Arguments

• model = IGRF(): Magnetic field model to use
• dir = 1: Tracing direction (+1 for parallel to B, -1 for anti-parallel)
• r0 = 1.0: Inner radial boundary (Earth radii)
• rlim = 10.0: Outer radial boundary (Earth radii)
• maxs = 100.0: Maximum arc length for integration
• from = getcsys(pos): Input coordinate system (Reference frame and coordinate representation)
• callback = nothing: extra SciML callback(s), combined with the boundary callbacks; u is in the model's native frame
• Additional keyword arguments are passed to solve()

Example

using GeoCotrans, OrdinaryDiffEqTsit5, Dates
sol = trace(GEO(3.0, 0.0, 0.0), DateTime(2020, 1, 1), Tsit5())
sol.u  # Cartesian positions [Re] in the frame of `from` along the line; sol.t is arc length

GaussCoefficients(g, h)

Schmidt semi-normalized Gauss coefficients; c[n, m] returns (; g, h).

SphericalHarmonicModel(name, coeffs; radius, frame, degree, order) <: InternalFieldModel

Internal field model in frame from coefficients coeffs; positions are in units of radius, and degree/order truncate the expansion.

TimeVaryingGaussCoefficients(epochs, pieces)
TimeVaryingGaussCoefficients(epochs, coefficients)

Gauss coefficients piecewise polynomial in time between sorted epochs and throws outside them. On [epochs[i], epochs[i + 1]], c(t) = Σⱼ pieces[i][j] αʲ⁻¹ with α = (t - epochs[i]) / (epochs[i + 1] - epochs[i]). Given coefficients at the epochs, interpolates linearly.

calc_dipole_angle(g10, g11, h11)

Calculate dipole angle (θ, φ) and dipole strength (b0) from spherical harmonic coefficients g10, g11, h11.

θ: dipole tilt angle (radians) φ: dipole longitude/phase (radians) b0: dipole strength (nT)

calc_dipole_geo(time)

Compute dipole direction in GEO coordinates.

car2gdz(x, y, z; scale = R🜨) -> [lat, lon, alt]

Inverse of gdz2car by Bowring's (1976) formula, as in CoordRefSystems.jl.

dipole_tilt(time)

Compute the dipole tilt angle μ (in radians).

The dipole tilt is the angle between the dipole axis and the GSM Z-axis, positive when the north magnetic pole is tilted toward the Sun.

evalsph(coeffs, r, θ, φ, max_degree, max_order = max_degree) -> [B_r, B_θ, B_φ]

Evaluate the magnetic field at a point in spherical coordinates with Schmidt semi-normalized Gauss coefficients coeffs.g, coeffs.h; r is in units of the model's reference radius.

Mathematical Formulation

The magnetic field is derived from the scalar potential:

\[B_r = -\frac{∂V}{∂r}, \quad B_θ = -\frac{1}{r}\frac{∂V}{∂θ}, \quad B_φ = -\frac{1}{r\sin θ}\frac{∂V}{∂φ}\]

gdz2car(lat, lon, alt; scale = R🜨) -> [x, y, z]

Geodetic to Cartesian GEO in units of scale [km].

get_dipole_terms(g, h)

Compute dipole parameters (θ, φ, x0, y0, z0, b0) from IGRF coefficients.

Returns a named tuple: (θ, φ, x0, y0, z0, b0)

Length unit in km or as a Unitful length. Models defining it accept Unitful positions.