This package is part of the SatelliteToolbox.jl ecosystem.

SatelliteToolboxGeomagneticField.jl

This package contains models to compute the geomagnetic field vector. We currently have two models implemented:

1. The International Geomagnetic Reference Field (IGRF) v14; and
2. The simplified dipole model.

Installation

julia> using Pkg
julia> Pkg.add("SatelliteToolboxGeomagneticField")

Usage

IGRF

We have a native Julia implementation of the International Geomagnetic Reference Field (IGRF) v14 based on [1]. This model can be accessed by two functions: igrf and igrfd.

function igrfd(date::Number, <r, h>::Number, λ::Number, Ω::Number[, R]; kwargs...)

Compute the geomagnetic field vector [nT] at the date date [Year A.D.] and position (r or h, λ, Ω).

The position representation is defined by R. If R is Val(:geocentric), the input must be geocentric coordinates:

1. Distance from the Earth center r [m];
2. Geocentric latitude λ ∈ [-90°, +90°]; and
3. Geocentric longitude Ω ∈ [-180°, +180°].

If R is Val(:geodetic), the input must be geodetic coordinates:

1. Altitude above the reference ellipsoid h (WGS-84) [m];
2. Geodetic latitude λ ∈ [-90°, +90°]; and
3. Geodetic longitude Ω ∈ [-180°, +180°].

If R is omitted, it defaults to Val(:geocentric).

Warning We must have 1900 <= date <= 2035. A warning message is printed for dates greater than 2030 since the output is not reliable anymore. This message can be suppressed by setting the keyword show_warnings to Val(false).

Note The output vector will be represented in the North-East-Down (NED) reference system associated with the representation selected by the parameter R (geocentric or geodetic). The Y-axis always points East. In case of geocentric coordinates, the Z-axis points toward the center of Earth and the X-axis completes a right-handed coordinate system. In case of geodetic coordinates, the X-axis is tangent to the ellipsoid at the selected location and points toward North, whereas the Z-axis completes a right-handed coordinate system.

The following keywords are available:

• max_degree::Int: Maximum degree used in the spherical harmonics when computing the geomagnetic field. If it is higher than the available number of coefficients in the IGRF matrices, it will be clamped. If it is equal to or lower than 0, it will be set to 1. (Default = 13)
• show_warnings::Val: If it is Val(true), a warning is printed using @warn for dates greater than 2030, when the accuracy of the model is reduced. If it is Val(false), the code related to the warning is removed at compile time, enabling allocation-free calls (Default = Val(true)).
• P::Union{Nothing, AbstractMatrix}: An optional matrix that must contain at least max_degree + 1 × max_degree + 1 real numbers that will be used to store the Legendre coefficients, reducing the allocations. Its element type should be the output type T to avoid loss of precision. If it is nothing, the matrix will be created when calling the function.
• dP::Union{Nothing, AbstractMatrix}: An optional matrix that must contain at least max_degree + 1 × max_degree + 1 real numbers that will be used to store the Legendre derivative coefficients, reducing the allocations. Its element type should be the output type T to avoid loss of precision. If it is nothing, the matrix will be created when calling the function.

This function returns a SVector{3, T}, which is the geomagnetic field vector [nT] at the desired location represented in the NED reference system associated with the input representation (geocentric or geodetic). Notice that the output type T is obtained by promoting T1, T2, and T3 to a float.

function igrf(date::Number, <r, h>::Number, λ::Number, Ω::Number[, R]; kwargs...)

Compute the geomagnetic field vector [nT] at the date date [Year A.D.] and position (r or h, λ, Ω).

The position representation is defined by R. If R is Val(:geocentric), the input must be geocentric coordinates:

1. Distance from the Earth center r [m];
2. Geocentric latitude λ ∈ [-π / 2, +π / 2] [rad]; and
3. Geocentric longitude Ω ∈ [-π, +π] [rad].

If R is Val(:geodetic), the input must be geodetic coordinates:

1. Altitude above the reference ellipsoid h (WGS-84) [m];
2. Geodetic latitude λ ∈ [-π / 2, +π / 2] [rad]; and
3. Geodetic longitude Ω ∈ [-π, +π] [rad].

If R is omitted, it defaults to Val(:geocentric).

Warning We must have 1900 <= date <= 2035. A warning message is printed for dates greater than 2030 since the output is not reliable anymore. This message can be suppressed by setting the keyword show_warnings to Val(false).

Note The output vector will be represented in the North-East-Down (NED) reference system associated with the representation selected by the parameter R (geocentric or geodetic). The Y-axis always points East. In case of geocentric coordinates, the Z-axis points toward the center of Earth and the X-axis completes a right-handed coordinate system. In case of geodetic coordinates, the X-axis is tangent to the ellipsoid at the selected location and points toward North, whereas the Z-axis completes a right-handed coordinate system.

The following keywords are available:

• max_degree::Int: Maximum degree used in the spherical harmonics when computing the geomagnetic field. If it is higher than the available number of coefficients in the IGRF matrices, it will be clamped. If it is equal to or lower than 0, it will be set to 1. (Default = 13)
• show_warnings::Val: If it is Val(true), a warning is printed using @warn for dates greater than 2030, when the accuracy of the model is reduced. If it is Val(false), the code related to the warning is removed at compile time, enabling allocation-free calls (Default = Val(true)).
• P::Union{Nothing, AbstractMatrix}: An optional matrix that must contain at least max_degree + 1 × max_degree + 1 real numbers that will be used to store the Legendre coefficients, reducing the allocations. Its element type should be the output type T to avoid loss of precision. If it is nothing, the matrix will be created when calling the function.
• dP::Union{Nothing, AbstractMatrix}: An optional matrix that must contain at least max_degree + 1 × max_degree + 1 real numbers that will be used to store the Legendre derivative coefficients, reducing the allocations. Its element type should be the output type T to avoid loss of precision. If it is nothing, the matrix will be created when calling the function.

This function returns a SVector{3, T}, which is the geomagnetic field vector [nT] at the desired location represented in the NED reference system associated with the input representation (geocentric or geodetic). Notice that the output type T is obtained by promoting T1, T2, and T3 to a float.

julia> igrf(2017.12313, 640e3, 50 * pi / 180, 25 * pi / 180, Val(:geodetic))
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 15365.317549816931
  1274.459965624728
 34200.321000410804

julia> igrfd(2017.12313, 640e3, 50, 25, Val(:geodetic))
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 15365.317549816931
  1274.459965624728
 34200.321000410804

julia> igrf(2017.12313, 6371e3 + 640e3, 50 * pi / 180, 25 * pi / 180, Val(:geocentric))
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 15165.026159783198
  1269.1937290447652
 34242.14698863755

julia> igrfd(2017.12313, 6371e3 + 640e3, 50, 25, Val(:geocentric))
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 15165.026159783198
  1269.1937290447652
 34242.14698863755

julia> igrf(2031, 6371e3 + 640e3, 50 * pi/180, 25 * pi/180)
┌ Warning: The magnetic field computed with this IGRF version may be of reduced accuracy for years greater than 2030.
└ @ SatelliteToolboxGeomagneticField ~/.julia/packages/SatelliteToolboxGeomagneticField/src/igrf/igrf.jl:212
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 15121.462463886119
  1653.8244319548928
 34799.38328175518

julia> igrfd(2031, 6371e3 + 640e3, 50, 25)
┌ Warning: The magnetic field computed with this IGRF version may be of reduced accuracy for years greater than 2030.
└ @ SatelliteToolboxGeomagneticField ~/.julia/packages/SatelliteToolboxGeomagneticField/src/igrf/igrf.jl:212
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 15121.462463886119
  1653.8244319548928
 34799.38328175518

julia> igrf(2031, 6371e3 + 640e3, 50 * pi / 180, 25 * pi / 180; show_warnings = Val(false))
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 15121.462463886119
  1653.8244319548928
 34799.38328175518

julia> igrfd(2031, 6371e3+640e3, 50, 25; show_warnings = Val(false))
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 15121.462463886119
  1653.8244319548928
 34799.38328175518

Simplified dipole model

This model assumes that the Earth geomagnetic field is a perfect dipole. The approximation is not good but it can be sufficient for some analysis, such as those carried out at the Pre-Phase A of a space mission when the uncertainties are high.

The following function computes the geomagnetic field using the simplified dipole model:

function geomagnetic_dipole_field(r_e::AbstractVector{T}, year::Number = 2020) where T

Compute the geomagnetic field [nT] using the simplified dipole model at position r_e (ECEF reference frame) [m]. This function uses the year year to obtain the position of the South geomagnetic pole (which lies in the North hemisphere) and the dipole moment. If year is omitted, it defaults to 2020.

Notice that:

1. The output vector will be represented in the ECEF reference frame.
2. The returned vector type is obtained by converting T to a float.
3. The south geomagnetic pole position and dipole moment are obtained by interpolating the values provided in [2].

julia> R0 = 6378.137e3;

julia> r_e = [0; 0; R0 + 200e3];

julia> geomagnetic_dipole_field(r_e)
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
   1301.6719081615724
  -4179.187349236483
 -53460.041620751326

julia> geomagnetic_dipole_field(r_e, 1986)
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
   1715.2656071053527
  -4964.59806084178
 -54246.30480714959

julia> r_e = [R0 + 200e3; 0; 0];

julia> geomagnetic_dipole_field(r_e)
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 -2603.343816323145
 -4179.187349236483
 26730.020810375663

julia> geomagnetic_dipole_field(r_e, 1986)
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 -3430.5312142107055
 -4964.59806084178
 27123.152403574793

References

• [1] Beggan, C. D.; Kloss, C.; Amblard, P. et al. (2026). International Geomagnetic Reference Field: The Fourteenth Generation. Earth, Planets and Space, v. 78, n. 127.
• [2] http://wdc.kugi.kyoto-u.ac.jp/poles/polesexp.html