magneticField: Earth magnetic declination, inclination, and intensity

Description

Implements the 12th generation International Geomagnetic Reference Field (IGRF), based on a reworked version of a Fortran program downloaded from a NOAA website [1]. The code (subroutine igrf12syn) seems to have been written by Susan Macmillan of the British Geological Survey. Comments in the code indicate that it employs coefficients agreed to in December 2014 by the IAGA Working Group V-MOD. Comments in the igrf12syn source code also suggest that the valid time interval is from years 1900 to 2020, with only the values from 1945 to 2010 being considered definitive.

Usage

magneticField(longitude, latitude, time)

Arguments

longitude

longitude in degrees east (negative for degrees west). The dimensions must conform to lat.

latitude

latitude in degrees north, a number, vector, or matrix.

time

either a decimal year or a POSIX time corresponding to the longitude and latitude values, or a vector or matrix matching these location values.

Value

A list containing declination, inclination, and intensity.

References

1. The underlying Fortran code is from igrf12.f, downloaded the NOAA website (https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html) on June 7, 2015. (That website suggests that there have been no update to the algorithm as of March 29, 2017.)

Examples

Run this code

# NOT RUN {
library(oce)
# Halifax NS
magneticField(-(63+36/60), 44+39/60, 2013)

# }
# NOT RUN {
## map of North American values
data(coastlineWorld)
mapPlot(coastlineWorld, longitudelim=c(-130,-55), latitudelim=c(35,60),
        projection="+proj=lcc +lat_0=20 +lat_1=60 +lon_0=-100")
lon <- seq(-180, 180, 1)
lat <- seq(-90, 90)
lonm <- rep(lon, each=length(lat))
latm <- rep(lat, times=length(lon))
## Note the counter-intuitive nrow argument
decl <- matrix(magneticField(lonm, latm, 2013)$declination,
               nrow=length(lon), byrow=TRUE)
mapContour(lon, lat, decl, col='red', levels=seq(-90, 90, 5))
incl <- matrix(magneticField(lonm, latm, 2013)$inclination,
               nrow=length(lon), byrow=TRUE)
mapContour(lon, lat, incl, col='blue', levels=seq(-90, 90, 5))
# }
# NOT RUN {
# }