GeoCotrans.jl

Coordinate systems, transformations, geomagnetic field models, and field line tracing.

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.

Utility functions for specifying coordinate representaion in GEO reference frames.

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

Coordinate transformations

Direct transformations

• geo2gei, gei2geo: Transform between GEO and GEI reference frames.
• geo2gsm, gsm2geo: Transform between GEO and GSM reference frames.
• gei2gsm, gsm2gei: Transform between GEI and GSM reference frames.
• gse2gsm, gsm2gse: Transform between GSE and GSM reference frames.
• geo2mag, mag2geo: Transform between GEO and MAG reference frames.
• gsm2sm, sm2gsm: Transform between GSM and SM reference frames.

Chain transformations

• gei2sm, sm2gei: Transform between GEI and SM reference frames.

• geo2sm, sm2geo: Transform between GEO and SM reference frames.

• gei2mag, mag2gei: Transform between GEI and MAG reference frames.

• sph2car, car2sph: Transform between spherical and cartesian coordinate representations.

References

• Coordinate transformations between geocentric systems

International Geomagnetic Reference Field (IGRF)

The International Geomagnetic Reference Field (IGRF) is a standard mathematical description of the Earth's main magnetic field. It is used widely in studies of the Earth's deep interior, crust, ionosphere, and magnetosphere.

API

• IGRF / igrf: Compute the geomagnetic field (IGRF-14)

Examples

using GeoCotrans, Dates
r, θ, φ = 1.0, deg2rad(45), deg2rad(45)
t = Date(2015)
Br, Bθ, Bφ = igrf(r, θ, φ, t) # ≈ [-45469, -21942, 2933]

# Input position in geodetic coordinates, output magnetic field in East-North-Up (ENU) coordinates
lat, lon = 60, 5
Be, Bn, Bu = igrf(GDZ(lat, lon), t) # [430, 15184, -48864]

# Input position as a matrix in GSM Cartesian coordinates, output magnetic field as a matrix
pos = rand(3, 6)
times = Date.(2015:2020)
Bgsm = igrf(pos, times; in = GSM())

References

• IAGA - NOAA/NCEI
• IGRF-14 Evaluation

GeoCotrans.jl provides functionality to trace magnetic field lines.

This feature is provided via a package extension and requires an ODE solver package to be loaded in the active environment. For example, load OrdinaryDiffEqTsit5 (or another SciML ODE solver) before calling trace.

API

• trace, FieldLineProblem, FieldLineCallback

Geocentric Equatorial Inertial (GEI) reference frame.

To Construct a CoordinateVector in GEI reference frame with cartesian representation.

GEI(x, y, z, t = nothing)
GEI(𝐫, t = nothing)

Geocentric Geographic (GEO) reference frame.

To Construct a CoordinateVector in GEO reference frame with cartesian representation.

GEO(x, y, z, t = nothing)
GEO(𝐫, t = nothing)

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).

To Construct a CoordinateVector in GSM reference frame with cartesian representation.

GSM(x, y, z, t = nothing)
GSM(𝐫, t = nothing)

Geocentric Solar Ecliptic (GSE) reference frame.

To Construct a CoordinateVector in GSE reference frame with cartesian representation.

GSE(x, y, z, t = nothing)
GSE(𝐫, t = nothing)

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.

To Construct a CoordinateVector in MAG reference frame with cartesian representation.

MAG(x, y, z, t = nothing)
MAG(𝐫, t = nothing)

Solar Magnetic (SM) reference frame.

To Construct a CoordinateVector in SM reference frame with cartesian representation.

SM(x, y, z, t = nothing)
SM(𝐫, t = nothing)

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.

Coordinate systems, transformations, geomagnetic field models, and field line tracing.

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.

Utility functions for specifying coordinate representaion in GEO reference frames.

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

Coordinate transformations

Direct transformations

• geo2gei, gei2geo: Transform between GEO and GEI reference frames.
• geo2gsm, gsm2geo: Transform between GEO and GSM reference frames.
• gei2gsm, gsm2gei: Transform between GEI and GSM reference frames.
• gse2gsm, gsm2gse: Transform between GSE and GSM reference frames.
• geo2mag, mag2geo: Transform between GEO and MAG reference frames.
• gsm2sm, sm2gsm: Transform between GSM and SM reference frames.

Chain transformations

• gei2sm, sm2gei: Transform between GEI and SM reference frames.

• geo2sm, sm2geo: Transform between GEO and SM reference frames.

• gei2mag, mag2gei: Transform between GEI and MAG reference frames.

• sph2car, car2sph: Transform between spherical and cartesian coordinate representations.

References

• Coordinate transformations between geocentric systems

International Geomagnetic Reference Field (IGRF)

The International Geomagnetic Reference Field (IGRF) is a standard mathematical description of the Earth's main magnetic field. It is used widely in studies of the Earth's deep interior, crust, ionosphere, and magnetosphere.

API

• IGRF / igrf: Compute the geomagnetic field (IGRF-14)

Examples

using GeoCotrans, Dates
r, θ, φ = 1.0, deg2rad(45), deg2rad(45)
t = Date(2015)
Br, Bθ, Bφ = igrf(r, θ, φ, t) # ≈ [-45469, -21942, 2933]

# Input position in geodetic coordinates, output magnetic field in East-North-Up (ENU) coordinates
lat, lon = 60, 5
Be, Bn, Bu = igrf(GDZ(lat, lon), t) # [430, 15184, -48864]

# Input position as a matrix in GSM Cartesian coordinates, output magnetic field as a matrix
pos = rand(3, 6)
times = Date.(2015:2020)
Bgsm = igrf(pos, times; in = GSM())

References

• IAGA - NOAA/NCEI
• IGRF-14 Evaluation

CoordinateVector{F, T}

3-element FieldVector in a specific coordinate frame F using representation R.

Geocentric Equatorial Inertial (GEI) reference frame.

To Construct a CoordinateVector in GEI reference frame with cartesian representation.

GEI(x, y, z, t = nothing)
GEI(𝐫, t = nothing)

Geocentric Geographic (GEO) reference frame.

To Construct a CoordinateVector in GEO reference frame with cartesian representation.

GEO(x, y, z, t = nothing)
GEO(𝐫, t = nothing)

Geocentric Solar Ecliptic (GSE) reference frame.

To Construct a CoordinateVector in GSE reference frame with cartesian representation.

GSE(x, y, z, t = nothing)
GSE(𝐫, t = nothing)

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).

To Construct a CoordinateVector in GSM reference frame with cartesian representation.

GSM(x, y, z, t = nothing)
GSM(𝐫, t = nothing)

IGRF(;)

Load the International Geomagnetic Reference Field (IGRF) model.

IGRF is a time-varying model of Earth's main magnetic field with coefficients at 5-year epochs from 1965 to 2030, linearly interpolated between epochs.

By default the input coordinate system in is (GEO(), Spherical()).

Examples

m = IGRF()
r, θ, φ = 1.0, deg2rad(45), deg2rad(45)
t = DateTime("2021-03-28")
m(r, θ, φ, t)

# When input is in different coordinate system, specify `in` or decorate the input
# By default, the output for GDZ input is (Be, Bn, Bu) in East-North-Up (ENU) frame
lat, lon = 60.39299, 5.32415
m2 = IGRF(; in = GDZ())
m2(lat, lon, 0, t)
m(GDZ(lat, lon, 0), t)

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.

To Construct a CoordinateVector in MAG reference frame with cartesian representation.

MAG(x, y, z, t = nothing)
MAG(𝐫, t = nothing)

Solar Magnetic (SM) reference frame.

To Construct a CoordinateVector in SM reference frame with cartesian representation.

SM(x, y, z, t = nothing)
SM(𝐫, t = nothing)

GDZ(ϕ, λ, h)

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, θ, φ)

Convert Cartesian coordinates to spherical coordinates with angles in radians.

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

Convert Cartesian coordinate to spherical coordinate with angles in degrees.

gdz2sph(lat, lon, alt)

Convert (lat [deg], lon [deg], alt [km]) in Geodetic coordinate to Spherical geocentric coordinate (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.

gei2gsm(x, t)

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

gei2mag(x, t)

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

gei2sm(x, t)

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

geo2gei(x, t)

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

geo2gsm(x, t)

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

geo2mag(x, t)

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

geo2sm(x, t)

Transforms coordinate(s) x from Geographic (GEO) to Solar Magnetic (SM) reference frame at time(s) 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)

Compute magnetic local time (MLT) in hours in the range [0, 24).

MLT is computed from the difference between the magnetic longitudes of the position and the subsolar direction in MAG coordinates.

• IRBEM implementation

gse2gsm(x, t)

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

gsm2gei(x, t)

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

gsm2geo(x, t)

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

gsm2gse(x, t)

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

gsm2sm(x, t)

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

igrf(𝐫, t; kw...) = IGRF()(𝐫, t; kw...)

A convenience function for IGRF().

Calculate IGRF model given coordinate(s) 𝐫 at time(s) t.

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.

r is the radius in Earth radii (Re), θ is the colatitude in degrees, and φ is the longitude in degrees.

mag2gei(x, t)

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

mag2geo(x, t)

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

mag2sm(x, t)

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

sm2gei(x, t)

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

sm2geo(x, t)

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

sm2gsm(x, t)

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

sm2mag(x, t)

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

Calculates cartesian field components from spherical ones

theta and phi are spherical angles of the point in radians

sph2car(r, θ, φ) -> (x, y, z)

Convert spherical coordinates to Cartesian coordinates.

Arguments

• r: radial distance
• θ: colatitude in radians [0, π]
• φ: east longitude in radians [0, 2π]

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

Convert spherical coordinates to Cartesian coordinates, with angles in degrees.

GeoCotrans.jl provides functionality to trace magnetic field lines.

This feature is provided via a package extension and requires an ODE solver package to be loaded in the active environment. For example, load OrdinaryDiffEqTsit5 (or another SciML ODE solver) before calling trace.

API

• trace, FieldLineProblem, FieldLineCallback

FieldLineCallback(; r0=1.0, rlim=10.0)

Create a callback for terminating field line integration at boundaries.

Keyword Arguments

• r0 = 1.0: Inner radial boundary (Earth radii)
• rlim = 10.0: Outer radial boundary (Earth radii)

FieldLineProblem(pos, tspan, t; model=IGRF(), dir=1)

Create an ODEProblem for tracing a magnetic field line in model at time t.

• dir::Int = 1: Tracing direction (+1 for parallel to B, -1 for anti-parallel)

Example

using GeoCotrans, OrdinaryDiffEqTsit5, Dates

t = DateTime(2020, 1, 1)
pos = [3.0, 0.0, 0.0]
prob = FieldLineProblem(pos, (0.0, 50.0), t)
sol = solve(prob, Tsit5())

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

Trace a magnetic field line using the specified SciML solver.

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
• in = getcsys(pos): Input coordinate system (Reference frame and coordinate representation)
• Additional keyword arguments are passed to solve()

Example

using GeoCotrans, OrdinaryDiffEqTsit5
sol = trace([3.0, 0.0, 0.0], t, Tsit5())

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.

calc_sun_gei(ra, dec)

Calculate sun direction vector in GEI given right ascension ra and declination dec.

calculate_gst_alt(time)

Alternative implementation of Greenwich sidereal time calculation based on the reference algorithm from pyspedas.cotrans_tools.csundir_vect.

car2gdz(x, y, z; scale=R🜨)

Convert Cartesian GEO coordinates (x, y, z) to Geodetic coordinates (φ [deg], λ [deg], h [km]).

Uses Bowring's formula (1976). The input coordinates are assumed to be in units of scale (default: Earth radii).

Reference:

• https://en.wikipedia.org/wiki/Geodetic_coordinates
• https://github.com/JuliaEarth/CoordRefSystems.jl/blob/main/src/crs/geographic/geodetic.jl#L197

csundir(time) -> (gst, ra, dec, elong, obliq)

Calculate the direction of the sun, returns a tuple of (gst, ra, dec, elong, obliq) in radians.

• gst: Greenwich mean sidereal time
• ra: Right ascension
• dec: Declination of the sun
• elong: ecliptic longitude
• obliq: Inclination of Earth's axis

See also AstroLib.sunpos.

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, θ, φ)

Evaluate the magnetic field at a point in spherical coordinates.

Arguments

• coeffs: The magnetic field coefficients
• r: Radial distance in planetary radii (dimensionless)
• θ: Colatitude in radians [0, π]
• φ: East longitude in radians [0, 2π]

Returns

• SVector{3,Float64}: Magnetic field vector [Br, Bθ, B_φ] in nanoTesla

Coordinate System

• Spherical coordinates (r, θ, φ) where:
  • r: radial distance from planet center
  • θ: colatitude (0 at north pole, π at south pole)
  • φ: east longitude (0 at prime meridian, increases eastward)

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(φ, λ, h; scale=R🜨)

Convert (φ [deg], λ [deg], h [km]) in Geodetic coordinate to Cartesian GEO coordinate (x [Re], y [Re], z [Re]).

Uses the standard conversion formula from Wikipedia:

• X = (N + h) * cos(φ) * cos(λ)
• Y = (N + h) * cos(φ) * sin(λ)
• Z = (N * (b²/a²) + h) * sin(φ)

where N is the prime vertical radius of curvature: N = a² / √(a² cos²(φ) + b² sin²(φ))

gei2gsm_mat(time)

Compute the GEI to GSM transformation matrix.

First axis: sun direction (xS, yS, zS) Second axis: y-axis (cross product of dipole and sun, normalized) Third axis: z-axis (cross product of sun and y-axis, normalized)

geo2mag_mat(time)

Compute the GEO to MAG transformation matrix.

This follows Hapgood (1992) / UCL geo_tran convention:

T5 = <lat-90, Y> * <long, Z>

where lat/long are the dipole pole coordinates from the IGRF dipole terms. With dipole colatitude θ = π/2 - lat and dipole longitude φ = long, this can be written as R_y(θ) * R_z(-φ).

get_dipole_terms(g, h)

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

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

gse2gsm_mat(time)

Compute the GSE to GSM transformation matrix T3 = <- psi,X>, where the rotation angle psi is the GSE-GSM angle.

This transformation is a rotation in the GSE YZ plane from the GSE Z axis to the GSM Z axis.

gsm2sm_mat(time)

Compute the GSM to SM transformation matrix.

The transformation is a simple rotation around the Y-axis by the dipole tilt angle μ. GSM and SM share the same Y-axis (perpendicular to the dipole-sun plane).