World Magnetic Model
The World Magnetic Model (WMM) is a data model, which describes the Earth’s magnetic field.
The magnetic field lines originate at the magnetic north pole, which is located somewhere north of Canada. So the magnetic north pole is not identical to the geographic north pole. For this reason a compass does not exactly point to the north pole. For example the deviation of the field lines in central Europe is about 2 degrees east.
The WMM describes the above deviation of the magnetic field lines with respect to the geographic north pole. To make use of the WMM we first download the C source code of the WMM Geomagnetism Library from here (we download the file WMM2015_Linux.tar.gz):
https://www.ngdc.noaa.gov/geomag/WMM/soft.shtml
After unpacking the tar ball, we compile the “grid” software in the src folder on the Unix command line like so:
gcc -c wmm_grid.c gcc -c GeomagnetismLibrary.c gcc wmm_grid.o GeomagnetismLibrary.o -o grid
The grid executable needs to reside in the same directory as the WMM.COF file. Here is an example run of the grid utility:
Welcome to the World Magnetic Model (WMM) C-Program
of the US National Geophysical Data Center
--- Grid Calculation Program ----
--- Model Release Date: 15 Dec 2014 ---
--- Software Release Date: 08 Dec 2014 ---
This program may be used to generate a grid/volume of magnetic field values
over latitude, longitude, altitude and time axes. To skip an axis,
keep the start and end values the same and enter zero for the step size.
Enter grid parameters
Please Enter Minimum Latitude (in decimal degrees):
-90
Please Enter Maximum Latitude (in decimal degrees):
90
Please Enter Minimum Longitude (in decimal degrees):
-180
Please Enter Maximum Longitude (in decimal degrees):
180
Please Enter Step Size (in decimal degrees):
5
Select height (default : above MSL)
1. Above Mean Sea Level
2. Above WGS-84 Ellipsoid
1
Please Enter Minimum Height above MSL (in km):
0
Please Enter Maximum Height above MSL (in km):
0
Please Enter height step size (in km):
0
Please Enter the decimal year starting time:
2015.5
Please Enter the decimal year ending time:
2015.5
Please Enter the time step size:
0
Enter a geomagnetic element to print. Your options are:
1. Declination 9. Ddot
2. Inclination 10. Idot
3. F 11. Fdot
4. H 12. Hdot
5. X 13. Xdot
6. Y 14. Ydot
7. Z 15. Zdot
8. GV 16. GVdot
For gradients enter: 17
1
Select output :
1. Print to a file
2. Print to Screen
2
-90.00 -180.00 0.0000 2015.50 149.88 0.40
-90.00 -175.00 0.0000 2015.50 144.88 0.40
-90.00 -170.00 0.0000 2015.50 139.88 0.40
...snip...
90.00 170.00 0.0000 2015.50 164.63 2.85
90.00 175.00 0.0000 2015.50 169.63 2.85
90.00 180.00 0.0000 2015.50 174.63 2.85
Press any key to exit...
According to the above output, the declination of the magnetic field at the equator (with latitude zero) varies from −20 to +10 degrees.
Now we textually format the above declination data into a matrix and put it as an array into a C++ module just like we did with the geoid data:
static const int size_y=37;
static const float data[size_y][size_x] = {
{149.88,144.88,139.88,...
...,164.63,169.63,174.63}
};
And add a linear interpolation function to the above data:
#include <math.h>
#include <limits> // for nan
#ifndef NAN
# define NAN (std::numeric_limits<double>::quiet_NaN())
#endif
// check whether of not a floating point value is not a number
inline int isNAN(const double v) {return(v!=v);}
// bi-linearily sample data matrix at normalized position (x,y)
inline float interpolate(float x, float y,
const float *data, unsigned int width, unsigned int height)
{
unsigned int i,j;
float val1,val2,val3,val4;
const float *floatptr;
if (x>=0.0f && x<=1.0f && y>=0.0f && y<=1.0f)
{
x*=width-1; y*=height-1;
i=(unsigned int)(x+0.5f); j=(unsigned int)(y+0.5f);
x-=i; y-=j;
if (i==width-1 && width>1) {i=width-2; x=1.0f;}
if (j==height-1 && height>1) {j=height-2; y=1.0f;}
floatptr=&data[i+j*width];
val1=floatptr[0];
if (isNAN(val1)) return(NAN);
if (width<2)
return(val1);
else
{
val2=floatptr[1];
if (isNAN(val2)) return(NAN);
if (height<2)
return((1.0f-x)*val1+x*val2);
else
{
val3=floatptr[width];
val4=floatptr[width+1];
if (isNAN(val3) || isNAN(val4)) return(NAN);
return((1.0f-y)*((1.0f-x)*val1+x*val2)+
y*((1.0f-x)*val3+x*val4));
}
}
}
return(NAN);
}
static const int size_x=73;
static const int size_y=37;
static const double LLWGS84_swx=-180;
static const double LLWGS84_swy=-90;
static const double LLWGS84_nwx=-180;
static const double LLWGS84_nwy=90;
static const double LLWGS84_nex=180;
static const double LLWGS84_ney=90;
static const double LLWGS84_sex=180;
static const double LLWGS84_sey=-90;
static const float data[size_y][size_x] = {
{149.88,144.88,139.88,...
...,164.63,169.63,174.63}
};
// sample the magnetic declication at a geographic position
float sample_mag(double latitude, double longitude)
{
// calculate normalized coordinates from geographic position
double nx = (longitude-LLWGS84_swx)/(LLWGS84_sex-LLWGS84_swx);
double ny = (latitude-LLWGS84_swy)/(LLWGS84_nwy-LLWGS84_swy);
// adjust wraparound longitudes
nx = nx - floor(nx);
// clamp out-of-range latitudes
if (ny<0.0) ny = 0.0;
else if (ny>1.0) ny = 1.0;
// interpolate magnetic declication
return(interpolate(nx, ny, (const float *)data, size_x, size_y));
}
Full source code is here:
Here is a usage example:
// example of calculating the wmm declination at Lat/Lon 49/11
double lat = 49.0, lon = 11.0;
double decl = sample_mag(lat, lon);
Note: The above example uses a subsampled coarse grid with 5 degrees spacing to estimate the magnetic declination. The error of the subsampled grid is somewhere in the order of ±1 degree. This is sufficient if you want to correct a magnetic compass of a smartphone, for example, which has a much higher systemic error than the approximation error of the grid. For use cases, where you need a smaller approximation error, the same principle as outlined above applies, but we just use a denser spacing of the samples.
Please do me a favour: If you happen to use the above code, please drop me a note and link to this page!