Update World Magnetic Model from WMM2020 to WMM2025

README.md

@@ -3,7 +3,7 @@ geomag is an implementation in Go of the NOAA World Magnetic Model.
 
 The World Magnetic Model home is at https://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml.
 
-The coefficients for 2020-2024 can be downloaded at https://www.ngdc.noaa.gov/geomag/WMM/data/WMM2020/WMM2020COF.zip
+The coefficients for 2025-2029 can be downloaded at https://www.ncei.noaa.gov/products/world-magnetic-model
 
 ## Commands
 geomag provides two command line programs, modeled after the command line programs in the official NOAA software.
@@ -64,7 +64,7 @@ h, err := egm96.NewLocationGeodetic(-12.25, 82.75, 1000).HeightAboveMSL()
 ```
 
 ### wmm
-Package wmm provides a representation of the 2020 World Magnetic Model (WMM),
+Package wmm provides a representation of the 2025 World Magnetic Model (WMM),
 a mathematical model of the magnetic field produced by the Earth's core and
 its variation over time.
 
@@ -98,7 +98,7 @@ Please submit an issue on github if you notice any other issues.
 
 ## Updating Coefficients
 Use go-bindata in the root directory to update the coefficients stored in binary form.
-Coefficients are currently updated through the WMM2020 model.
+Coefficients are currently updated through the WMM2025 model.
 First, unzip the new WMM zip file in the assets/wmm directory, then
 `go-bindata -o ../../../pkg/wmm/bindata.go WMM.COF`
 Inside the `bindata.go` file, change the package from `main` to `wmm`.

cmd/wmm_point/main.go

@@ -1,12 +1,12 @@
 // wmm_point estimates the strength and direction of Earth's main Magnetic field for a given point/area.
 //
 // Usage is
-//  wmm_point --cof_file=WMM2020.COF --spherical [latitude] [longitude] [altitude] [date]
+//  wmm_point --cof_file=WMM2025.COF --spherical [latitude] [longitude] [altitude] [date]
 //
-// The World Magnetic Model (WMM) for 2020
+// The World Magnetic Model (WMM) for 2025
 // is a model of Earth's main Magnetic field.  The WMM
 // is recomputed every five (5) years, in years divisible by
-// five (i.e. 2010, 2015, 2020).
+// five (i.e. 2015, 2020, 2025).
 //
 // Information on the model is available at https://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml
 //
@@ -84,7 +84,7 @@ import (
 )
 
 const (
-	usage = "wmm_point --cof_file=WMM2020.COF --spherical [latitude] [longitude] [altitude] [date]"
+	usage = "wmm_point --cof_file=WMM2025.COF --spherical [latitude] [longitude] [altitude] [date]"
 	cofUsage = "COF coefficients file to use, empty for the built-in one"
 	sphericalUsage = "Output spherical values instead of ellipsoidal"
 	lngErr = "Error: Degree input is outside legal range. The legal range is from -180 to 360."

go.mod

@@ -0,0 +1,3 @@
+module github.com/westphae/geomag
+
+go 1.26.1

pkg/wmm/coefficients.go

@@ -58,8 +58,8 @@ func GetWMMCoefficients(n, m int, t time.Time) (gnm, hnm, dgnm, dhnm float64, er
 // Epoch, COFName, and ValidDate.
 // If the passed filename is "", it loads the default (current) coefficients file.
 //
-// The default coefficients file is currently WMM2020.COF, valid from
-// 12/10/2019 until 12/31/2024.
+// The default coefficients file is currently WMM2025.COF, valid from
+// 11/13/2024 until 12/31/2029.
 func LoadWMMCOF(fn string) (err error) {
 	var (
 		data []byte

pkg/wmm/coefficients_test.go

@@ -57,3 +57,28 @@ func TestGetWMM2020Coefficients(t *testing.T) {
 		}
 	}
 }
+
+func TestGetWMM2025Coefficients(t *testing.T) {
+	_ = LoadWMMCOF("testdata/WMM2025.COF")
+	nms := [][]int{{1, 0}, {2, 2}, {5, 1}, {5, 4}, {12, 0}, {12, 6}, {12, 11}}
+	gs := []float64{-29351.8, 1649.3, 368.9, -142.0, -2.0, 0.6, -1.3}
+	hs := []float64{0.0, -815.1, 45.4, 43.0, 0.0, 0.6, 0.1}
+	dgs := []float64{12.0, -8.0, 1.4, 2.2, 0.0, 0.1, 0.0}
+	dhs := []float64{0.0, -12.1, -0.5, 1.7, 0.0, 0.0, 0.0}
+	ts := []time.Time{
+		time.Date(2025, 1, 1, 0, 0, 0, 0, time.UTC),
+		time.Date(2026, 1, 1, 0, 0, 0, 0, time.UTC),
+	}
+
+	for j, tt := range ts {
+		for i, nm := range nms {
+			n := nm[0]
+			m := nm[1]
+			g, h, dg, dh, _ := GetWMMCoefficients(n, m, tt)
+			testDiff(fmt.Sprintf("G(%d,%d)", n, m), g, gs[i]+float64(j)*dgs[i], eps, t)
+			testDiff(fmt.Sprintf("H(%d,%d)", n, m), h, hs[i]+float64(j)*dhs[i], eps, t)
+			testDiff(fmt.Sprintf("DG(%d,%d)", n, m), dg, dgs[i], eps, t)
+			testDiff(fmt.Sprintf("DH(%d,%d)", n, m), dh, dhs[i], eps, t)
+		}
+	}
+}

pkg/wmm/magnetic_field.go

@@ -23,14 +23,14 @@ import (
 
 const (
 	AGeo  = 6371200 // Geomagnetic Reference Radius
-	errX  = 131     // WMM global average X error, nT
-	errY  = 94      // WMM global average Y error, nT
-	errZ  = 157     // WMM global average Z error, nT
-	errH  = 128     // WMM global average H error, nT
-	errF  = 148     // WMM global average F error, nT
-	errI  = 0.21    // WMM global average I error, º
+	errX  = 137     // WMM global average X error, nT
+	errY  = 89      // WMM global average Y error, nT
+	errZ  = 141     // WMM global average Z error, nT
+	errH  = 133     // WMM global average H error, nT
+	errF  = 138     // WMM global average F error, nT
+	errI  = 0.20    // WMM global average I error, º
 	errDA = 0.26    // WMM rough global average D error away from poles, º
-	errDB = 5625    // WMM average H uncertainty scale near the poles, nT
+	errDB = 5417    // WMM average H uncertainty scale near the poles, nT
 )
 
 // MagneticField represents a geomagnetic field and its rate of change.

pkg/wmm/magnetic_field_test.go

@@ -207,6 +207,139 @@ func TestAll2015v2TestValuesFromPaper(t *testing.T) {
 
 }
 
+func TestAll2025TestValuesFromPaper(t *testing.T) {
+	var (
+		date                   DecimalYear
+		height                 float64
+		lat, lon               float64
+		x, y, z                float64
+		h, f, i, d             float64
+		xdot, ydot, zdot       float64
+		hdot, fdot, idot, ddot float64
+		data                   []byte
+		dat                    []string
+		err                    error
+	)
+
+	_ = LoadWMMCOF("testdata/WMM2025.COF")
+
+	data, err = ioutil.ReadFile("testdata/WMM2025_TEST_VALUES.txt")
+	scanner := bufio.NewScanner(bytes.NewReader(data))
+	// Read and parse header
+	if !scanner.Scan() {
+		panic(err)
+	}
+	for scanner.Scan() {
+		// Skip the header lines
+		if scanner.Text()[0]=='#' {
+			continue
+		}
+
+		dat = strings.Fields(scanner.Text())
+
+		dd, err := strconv.ParseFloat(dat[0], 64)
+		if err != nil {
+			panic(err)
+		}
+		date = DecimalYear(dd)
+
+		if dd, err = strconv.ParseFloat(dat[1], 64); err != nil {
+			panic(err)
+		}
+		height = dd*1000
+
+		if dd, err = strconv.ParseFloat(dat[2], 64); err != nil {
+			panic(err)
+		}
+		lat = dd
+
+		if dd, err = strconv.ParseFloat(dat[3], 64); err != nil {
+			panic(err)
+		}
+		lon = dd
+
+		loc := egm96.NewLocationGeodetic(lat,lon,height)
+
+		mag, _ := CalculateWMMMagneticField(loc, date.ToTime())
+		xE, yE, zE, dxE, dyE, dzE := mag.Ellipsoidal()
+
+		if d, err = strconv.ParseFloat(dat[4], 64); err != nil {
+			panic(err)
+		}
+		testDiff("D", mag.D(), d, 0.005, t)
+
+		if i, err = strconv.ParseFloat(dat[5], 64); err != nil {
+			panic(err)
+		}
+		testDiff("I", mag.I(), i, 0.005, t)
+
+		if h, err = strconv.ParseFloat(dat[6], 64); err != nil {
+			panic(err)
+		}
+		testDiff("H", mag.H(), h, 0.05, t)
+
+		if x, err = strconv.ParseFloat(dat[7], 64); err != nil {
+			panic(err)
+		}
+		testDiff("X", xE, x, 0.05, t)
+
+		if y, err = strconv.ParseFloat(dat[8], 64); err != nil {
+			panic(err)
+		}
+		testDiff("Y", yE, y, 0.05, t)
+
+		if z, err = strconv.ParseFloat(dat[9], 64); err != nil {
+			panic(err)
+		}
+		testDiff("Z", zE, z, 0.05, t)
+
+		if f, err = strconv.ParseFloat(dat[10], 64); err != nil {
+			panic(err)
+		}
+		testDiff("F", mag.F(), f, 0.05, t)
+
+		if ddot, err = strconv.ParseFloat(dat[11], 64); err != nil {
+			panic(err)
+		}
+		testDiff("Ddot", MagneticField(mag).DD(), ddot, 0.05, t)
+
+		if idot, err = strconv.ParseFloat(dat[12], 64); err != nil {
+			panic(err)
+		}
+		testDiff("Idot", mag.DI(), idot, 0.05, t)
+
+		if hdot, err = strconv.ParseFloat(dat[13], 64); err != nil {
+			panic(err)
+		}
+		testDiff("Hdot", mag.DH(), hdot, 0.05, t)
+
+		if xdot, err = strconv.ParseFloat(dat[14], 64); err != nil {
+			panic(err)
+		}
+		testDiff("Xdot", dxE, xdot, 0.05, t)
+
+		if ydot, err = strconv.ParseFloat(dat[15], 64); err != nil {
+			panic(err)
+		}
+		testDiff("Ydot", dyE, ydot, 0.05, t)
+
+		if zdot, err = strconv.ParseFloat(dat[16], 64); err != nil {
+			panic(err)
+		}
+		testDiff("Zdot", dzE, zdot, 0.05, t)
+
+		if fdot, err = strconv.ParseFloat(dat[17], 64); err != nil {
+			panic(err)
+		}
+		testDiff("Fdot", mag.DF(), fdot, 0.05, t)
+	}
+
+	if err := scanner.Err(); err != nil {
+		panic(err)
+	}
+
+}
+
 func TestAll2020TestValuesFromPaper(t *testing.T) {
 	var (
 		date                   DecimalYear