igrf – International Geomagnetic Reference Field

This module computes the compass deviation (or variance) at a given point. It allows mapping from true bearings to magnetic bearings.

This is a migration of the Geomag 7.0 software from Fortran to Python.

It relies on the coefficients table, provided by NOAA. See https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html for more information.

Background

From the web page.

\[V(r,\theta,\phi,t) = a \sum_{n=1}^{N}\sum_{m=0}^{n}\Big(\frac{a}{r}\Big)^{n+1}[g_n^m(t)\cos (m \phi) + h_n^m(t)\sin (m \phi)] P_n^m(\cos \theta)\]

With \(g_n^m\) and \(h_n^m\) are the coefficients to determine phase and amplitude. These are provided by the associated igrf13coeffs.txt file.

\(a\) is the radius of the Earth, 6371.2 km.

\(N\) is the degree of truncation, 13.

The position, \(r,\theta,\phi,t\), is the altitude, latitude, longitude and time. At the surface, \(r = a\).

\(P_n^m\) are associated Legendre functions of the first kind of degree \(n\), order \(m\).

\[P_n^m(z) = \frac{1}{\Gamma(1-m)}\Big[\frac{1+z}{1-z}\Big]^{m / 2} {{}_{2}F_{1}}\left(-n, n+1; 1-m; \frac{1-z}{2}\right)\]

\({}_{2}F_{1}(a,b;c;z)\) is the Hypergeometric function.

\[{}_{2}F_{1}(a,b;c;z)=\sum_{n=0}^{\infty} {\frac {(a)_{n}(b)_{n}}{(c)_{n}}}{\frac {z^{n}}{n!}}=1+{\frac {ab}{c}}{\frac {z}{1!}}+{\frac {a(a+1)b(b+1)}{c(c+1)}}{\frac {z^{2}}{2!}}+\cdots\]

The Gamma function is \(\Gamma (n)=(n-1)!\), or \(\Gamma (z)=\int _{0}^{\infty }x^{z-1}e^{-x}\,dx,\ \qquad \Re (z)>0\)

http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html

Testing

This is tested with IGRF-11, and IGRF-13 coefficients.

Alternative

See https://pypi.org/project/pyIGRF/ for an implementation that can be forked to remove a numpy dependency.

Implementation

The core model.

class navtools.igrf.IGRF(model: str = './igrf11coeffs.txt')

The IGRF model.

The IGRF calculation requires a complex table of coefficients loaded from ./igrf11coeffs.txt

We defer coefficient loading until a call is made. Eager loading when an instance is created leads to odd errors when creating documentation.

prepare() → None

Locate the coefficients file and load them. This is a lazy loader, checked for each computation.

This will look in the following three places:

1. The given name, unmodified.

2. The installed directory for the given name.

3. (Since this may be running from the distribution kit) the parent of the installed directory.

static load_coeffs(file_path: pathlib.Path) → tuple

Returns two dictionaries of g and h coefficients by year and [n,m] index, and the year for which extrapolation must start.

Coefficient loading includes disentangling the unpleasant sparse matrix optimizations. The Fortran program had a number of very clever techniques to ruthlessly minimize the memory footprint.

Years 1900 through 1990 have 120 values. Degree of 10. Years 1995 to 2012 have 195 values. Degree of 13.

The file has a number of comment lines which begin with #. These are simply skipped.

The file has heading lines. The last of these starts with g/h. This line provides useful column titles.

Parameters

file_name – the path to the igrf11coeffs.txt file. Ideally this is a pathlib.Path, but a string will do.

Todo

Use urllib.request

We can then use urllib instead of simple Path.open() so that we could use “file:///path/to/file” or “https://server/path/to/file”.

Even though most use cases involve disconnected computers, and rely on a statically downloaded file. It’s easy to cover all the bases by opening a URL instead of assuming a local file.

The official file: http://www.ngdc.noaa.gov/IAGA/vmod/igrf11coeffs.txt

navtools.igrf.igrfsyn(date: float, nlat: float, elong: float, alt: float = 0.0, coord: str = 'D') → tuple[float, float, float, float]

The IGRF model.

The IGRF calculation requires a complex table of coefficients loaded from ./igrf11coeffs.txt

We defer coefficient loading until a call is made. Eager loading when an instance is created leads to odd errors when creating documentation.

A slightly simplified function that allows a client to get today’s declination at a specific point.

navtools.igrf.declination(nlat: float, elong: float, date: Optional[datetime.date] = None) → float

IGRF model for current declination.

Parameters

• nlat – north latitude as floating-point degrees

• elog – east longitude as floating-point degrees

• date – datetime.date in question, default is today.

Returns

declination as degrees.

A utility to convert simple degrees to (degree, minute) two-tuples.

navtools.igrf.deg2dm(deg: float) → tuple

Converts a simple degree value into a proper (degree,minute) pair with appropriate signs.

Parameters

deg – Degrees to convert

Returns

tuple of (deg, min) with proper signs.