*This package is part of the SatelliteToolbox.jl ecosystem.*

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

- The International Geomagnetic Reference Field (IGRF) v13; and
- The simplified dipole model.

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

We have a native Julia implementation of the International Geomagnetic Reference Field
(IGRF) v13 based on **[1]**. This mode can
be accessed by two functions: `irgf`

and `irgfd`

.

`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:

- Distance from the Earth center
`r`

[m]; - Geocentric latitude
`λ`

∈ (-90°, +90°); and - 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)`

.

WarningWe must have`1900 <= date <= 2030`

. A warning message is printed for dates greater than 2025 since the output is not reliable anymore. This message can be suppressed by setting the keyword`show_warnings`

to`false`

.

NoteThe output vector will be represented in the same reference system selected by the parameter`R`

(geocentric or geodetic). The Y-axis of the output reference system always points East. In case ofgeocentric coordinates, the Z-axis points toward the center of Earth and the X-axis completes a right-handed coordinate system. In case ofgeodetic coordinates, the X-axis is tangent to the ellipsoid at the selected location and points toward North, whereas the Z-axis completes a right-hand 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 of lower than 0, it will be set to 1. (**Default**= 13)`show_warnings::Bool`

: Show warnings about the data (**Default**=`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. 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. 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 same input reference (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:

- Distance from the Earth center
`r`

[m]; - Geocentric latitude
`λ`

∈ (-π / 2, +π / 2) [rad]; and - Geocentric longitude
`Ω`

∈ (-π, +π) [rad].

If `R`

is `Val(:geodetic)`

, the input must be **geodetic** coordinates:

- Altitude above the reference ellipsoid
`h`

(WGS-84) \[m]; - Geodetic latitude
`λ`

∈ (-π/2, +π/2) [rad]; and - Geodetic longitude
`Ω`

∈ (-π, +π) [rad].

If `R`

is omitted, it defaults to `Val(:geocentric)`

.

WarningWe must have`1900 <= date <= 2030`

. A warning message is printed for dates greater than 2025 since the output is not reliable anymore. This message can be suppressed by setting the keyword`show_warnings`

to`false`

.

NoteThe output vector will be represented in the same reference system selected by the parameter`R`

(geocentric or geodetic). The Y-axis of the output reference system always points East. In case ofgeocentric coordinates, the Z-axis points toward the center of Earth and the X-axis completes a right-handed coordinate system. In case ofgeodetic coordinates, the X-axis is tangent to the ellipsoid at the selected location and points toward North, whereas the Z-axis completes a right-hand 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 of lower than 0, it will be set to 1. (**Default**= 13)`show_warnings::Bool`

: Show warnings about the data (**Default**=`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. 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. 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 same input reference (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 SVector{3, Float64} with indices SOneTo(3):
15365.787505205588
1274.9958640696996
34201.2182033379
julia> igrfd(2017.12313, 640e3, 50, 25, Val(:geodetic))
3-element SVector{3, Float64} with indices SOneTo(3):
15365.787505205588
1274.9958640696996
34201.2182033379
julia> igrf(2017.12313, 6371e3 + 640e3, 50 * pi / 180, 25 * pi / 180, Val(:geocentric))
3-element SVector{3, Float64} with indices SOneTo(3):
15165.486702524944
1269.7264334427584
34243.04928373083
julia> igrfd(2017.12313, 6371e3 + 640e3, 50, 25, Val(:geocentric))
3-element SVector{3, Float64} with indices SOneTo(3):
15165.486702524944
1269.7264334427584
34243.04928373083
```

```
julia> igrf(2026, 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 2025.
└ @ SatelliteToolboxGeomagneticField ~/tmp/SatelliteToolboxGeomagneticField/src/igrf/igrf.jl:301
3-element SVector{3, Float64} with indices SOneTo(3):
15118.591511098812
1588.129544718569
34668.84185460438
julia> igrfd(2026, 6371e3 + 640e3, 50, 25)
┌ Warning: The magnetic field computed with this IGRF version may be of reduced accuracy for years greater than 2025.
└ @ SatelliteToolboxGeomagneticField ~/tmp/SatelliteToolboxGeomagneticField/src/igrf/igrf.jl:301
3-element SVector{3, Float64} with indices SOneTo(3):
15118.591511098812
1588.129544718569
34668.84185460438
julia> igrf(2026, 6371e3 + 640e3, 50 * pi / 180, 25 * pi / 180; show_warnings = false)
3-element SVector{3, Float64} with indices SOneTo(3):
15118.591511098812
1588.129544718569
34668.84185460438
julia> igrfd(2026, 6371e3+640e3, 50, 25; show_warnings = false)
3-element SVector{3, Float64} with indices SOneTo(3):
15118.591511098812
1588.129544718569
34668.84185460438
```

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:

- The output vector will be represented in the ECEF reference frame.
- The returned vector type is obtained by converting
`T`

to a float. - The south geomagnetic pole position and dipole moment is obtained by interpolating the
values provided in
**[2]**.

```
julia> r_e = [0; 0; R0 + 200e3];
julia> geomag_dipole(r_e)
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
1286.0242861717802
-4232.804339060699
-53444.68086319672
julia> geomag_dipole(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> geomag_dipole(r_e)
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
-2572.0485723435604
-4232.804339060699
26722.34043159836
julia> geomag_dipole(r_e, 1986)
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
-3430.5312142107055
-4964.59806084178
27123.152403574793
```

**[1]****Alken, P.; Thébault, E.; Beggan, C. D. et al. (2021)**.*International Geomagnetic Reference Field: The Thirteenth Generation*.**Earth Planets Space**, v. 73, n. 49.**[2]**http://wdc.kugi.kyoto-u.ac.jp/poles/polesexp.html