IGRF 14th model added. (#41)

* Updated to IGRF 14th model

Author: proway2

.gitignore

@@ -13,3 +13,6 @@
 
 # Dependency directories (remove the comment below to include it)
 # vendor/
+
+# VS Code
+.vscode

README.md

@@ -2,7 +2,7 @@
 [![golangci-lint](https://github.com/proway2/go-igrf/actions/workflows/golangci-lint.yml/badge.svg)](https://github.com/proway2/go-igrf/actions/workflows/golangci-lint.yml)
 [![Go Reference](https://pkg.go.dev/badge/github.com/proway2/go-igrf.svg)](https://pkg.go.dev/github.com/proway2/go-igrf)
 [![Go Report Card](https://goreportcard.com/badge/github.com/proway2/go-igrf)](https://goreportcard.com/report/github.com/proway2/go-igrf)
-[![IGRF badge](https://badgen.net/static/IGRF%20model/13th/green)](https://www.ncei.noaa.gov/products/international-geomagnetic-reference-field)
+[![IGRF badge](https://badgen.net/static/IGRF%20model/14th/green)](https://www.ncei.noaa.gov/products/international-geomagnetic-reference-field)
 
 # go-igrf
 Pure Go IGRF (International Geomagnetic Reference Field). This is based on the existing `C` implementation. This package computes values for the geomagnetic field and secular variation for a given set of coordinates and date.
@@ -20,9 +20,9 @@ The output is of type `type IGRFresults struct`. Fields are:
 
 |Field name|Fortran name|Unit of measurement|Sample value|Notes|
 |-|-|-|-|-|
-|Declination|(D)|decimal degrees(°)|14.175 °|(+ve east)|
+|Declination|(D)|decimal degrees(°)|14.175|degrees east|
 |DeclinationSV|(D)|arcmin/yr|-8.92||
-|Inclination|(I)|decimal degrees(°)|74.229|(+ve down)
+|Inclination|(I)|decimal degrees(°)|74.229|degrees down|
 |InclinationSV|(I)|arcmin/yr|-2.59||
 |HzIntensity|Horizontal intensity (H)|nT|14626.6||
 |HorizontalSV|(H)|nT/yr|-16.8||
@@ -75,6 +75,7 @@ func main() {
 
 ## References
 
+- [IAGA V-MOD WG and main IGRF website](https://www.ncei.noaa.gov/products/international-geomagnetic-reference-field)
 - [Alken, P., Thébault, E., Beggan, C.D. et al. International Geomagnetic Reference Field: the thirteenth generation. Earth Planets Space 73, 49 (2021).](https://rdcu.be/cKqZv) https://doi.org/10.1186/s40623-020-01288-x
-
+- [International Association of Geomagnetism and Aeronomy. (2024). IGRF-14. Zenodo.](https://doi.org/10.5281/zenodo.14012302)
 - Implementation is based on [geomag70.c](https://www.ngdc.noaa.gov/IAGA/vmod/geomag70_linux.tar.gz). [License information](https://www.ngdc.noaa.gov/IAGA/vmod/geomag70_license.html).

calc/shval3.go

@@ -6,11 +6,11 @@ import "math"
 // The calculation is performed for two sets of coeffs for a single location,
 // thus it returns two sets of X, Y, Z.
 //
-// X - northward component
+// # X - northward component
 //
-// Y - eastward component
+// # Y - eastward component
 //
-// Z - vertically-downward component
+// # Z - vertically-downward component
 func Shval3(flat, flon, elev float64, nmax int, gha, ghb *[]float64) (float64, float64, float64, float64, float64, float64) {
 	// similar to shval3 from C implementation
 	var earths_radius float64 = 6371.2
@@ -157,13 +157,13 @@ func Shval3(flat, flon, elev float64, nmax int, gha, ghb *[]float64) (float64, f
 
 // Computes the geomagnetic D, I, H, and F from X, Y, and Z.
 //
-// D  - declination
+// # D  - declination
 //
-// I  - inclination
+// # I  - inclination
 //
-// H  - horizontal intensity
+// # H  - horizontal intensity
 //
-// F  - total intensity
+// # F  - total intensity
 func Dihf(x, y, z float64) (float64, float64, float64, float64) {
 	var d, i, h, f float64
 	sn := 0.0001

coeffs/coeffs.go

@@ -49,7 +49,7 @@ func (igrf *IGRFcoeffs) Coeffs(date float64) (*[]float64, *[]float64, int, error
 	min_epoch := (*igrf.epochs)[0]
 	max_epoch := (*igrf.epochs)[max_column-1]
 	if date < min_epoch || date > max_epoch {
-		return nil, nil, 0, fmt.Errorf("Date %v is out of range (%v, %v).", date, min_epoch, max_epoch)
+		return nil, nil, 0, fmt.Errorf("date %v is out of range (%v, %v)", date, min_epoch, max_epoch)
 	}
 	// calculate coeffs for the requested date
 	start, end := igrf.findEpochs(date)
@@ -90,7 +90,7 @@ func (igrf *IGRFcoeffs) interpolateCoeffs(start_epoch, end_epoch string, date fl
 	} else {
 		if nmax1 > nmax2 {
 			// the last column has degree of 8
-			// now it's anything after 2020.0
+			// now it's anything after 2025.0
 			k = nmax2 * (nmax2 + 2)
 			l = nmax1 * (nmax1 + 2)
 			interp = func(start, end, f float64) float64 {
@@ -101,7 +101,7 @@ func (igrf *IGRFcoeffs) interpolateCoeffs(start_epoch, end_epoch string, date fl
 			// between 1995.0 and 2000.0
 			k = nmax1 * (nmax1 + 2)
 			l = nmax2 * (nmax2 + 2)
-			interp = func(start, end, f float64) float64 {
+			interp = func(_, end, f float64) float64 {
 				return f * end
 			}
 			nmax = nmax2
@@ -212,7 +212,7 @@ func getEpochs(reader <-chan string) (*[]string, *[]float64, error) {
 		names, epochs := parseHeader(line, line2)
 		return &names, &epochs, nil
 	}
-	return nil, nil, errors.New("Unable to get epochs.")
+	return nil, nil, errors.New("unable to get epochs")
 }
 
 // Parses the header of the coeffs. Usually it's the first two non-comment lines.
@@ -257,7 +257,7 @@ func (igrf *IGRFcoeffs) getCoeffsForEpochs(provider <-chan string) error {
 		line_data := space_re.Split(line, -1)
 		line_coeffs, err := parseArrayToFloat(line_data[3:])
 		if err != nil {
-			return errors.New("Unable to parse coeffs.")
+			return errors.New("unable to parse coeffs")
 		}
 		igrf.loadCoeffs(i, line_coeffs)
 		i++
@@ -284,7 +284,7 @@ func nMaxForEpoch(epoch string) (int, error) {
 	}
 	if epoch_f < 2000.0 {
 		nmax = 10
-	} else if epoch_f > 2020.0 {
+	} else if epoch_f > 2025.0 {
 		nmax = 8
 	} else {
 		nmax = 13

coeffs/read.go

@@ -98,7 +98,7 @@ func parseArrayToFloat(raw_data []string) (*[]float64, error) {
 	for index, token := range raw_data {
 		real_data, err := strconv.ParseFloat(token, 32)
 		if err != nil {
-			return nil, errors.New("Unable to parse coeffs.")
+			return nil, errors.New("unable to parse coeffs")
 		}
 		if index == len(raw_data)-1 {
 			// real value calculated for the SV column

coeffs/read_test.go

@@ -34,7 +34,7 @@ func New() *IGRFdata {
 // date - decimal date (1900.00 to 2025).
 func (igd *IGRFdata) IGRF(lat, lon, alt, date float64) (IGRFresults, error) {
 	if igd.shc == nil {
-		return IGRFresults{}, errors.New("IGRFdata structure is not initialized.")
+		return IGRFresults{}, errors.New("IGRFdata structure is not initialized")
 	}
 	if err := checkInitialConditions(lat, lon, alt); err != nil {
 		return IGRFresults{}, err

coeffs/utils.go

@@ -92,7 +92,12 @@ func TestIGRFdata_IGRF(t *testing.T) {
 			}
 			// Horizontal intensity
 			allowed_abs_error := getMaxAllowedAbsoluteTolerance(tt.want.HorizontalIntensity)
-			compare = isClose(got.HorizontalIntensity, tt.want.HorizontalIntensity, allowed_relative_error, allowed_abs_error)
+			is_close := isClose(got.HorizontalIntensity, tt.want.HorizontalIntensity, allowed_relative_error, allowed_abs_error)
+			is_close_no_frac := noFracComparison(got.HorizontalIntensity, tt.want.HorizontalIntensity)
+			compare = is_close || is_close_no_frac
+			// For lat=-64.081, lon=135.866, alt=0, data=2021.5, horizontal intensity is 74.38, whereas reference value is 74.
+			// Fortran rounds/truncates almost all values, therefore it's hard to use either relative or absolute tolerance.
+			// noFracComparison compares a value with it's rounded/truncated value.
 			if !compare {
 				t.Errorf("IGRF() Horizontal intensity = %v, want %v", got.HorizontalIntensity, tt.want.HorizontalIntensity)
 			}
@@ -185,6 +190,17 @@ func getMaxAllowedAbsoluteTolerance(value float64) float64 {
 	return abs_tol
 }
 
+func noFracComparison(val, ref float64) bool {
+	// Since almost all values produces by Fortran are rounded/truncated, it's sometimes needed to compare
+	// calculated value with rounded/truncated value. Rounding doesn't always work, as it's unclear whether Fortran values
+	// are actually rounded or truncated.
+	// Therefore values are truncated and compared to be within certain range.
+	min := ref - 1.0
+	max := ref + 1.0
+	val_trunc := math.Floor((val))
+	return min <= val && val_trunc <= max
+}
+
 func isClose(a, b, rel_tol, abs_tol float64) bool {
 	abs_diff := math.Abs(a - b)
 	a_abs := math.Abs(a)
@@ -222,8 +238,8 @@ func produceTestsDataFromFile(file_descr *os.File) []testsData {
 	scanner := bufio.NewScanner(file_descr)
 	split_regex := regexp.MustCompile(`\s+`)
 	var lat, lon, alt float64
-	tests := make([]testsData, 125)
-	for num, i := 0, 0; scanner.Scan(); num++ {
+	var tests []testsData
+	for num := 0; scanner.Scan(); num++ {
 		line := scanner.Text()
 		if num == 1 {
 			// this is just a column names
@@ -244,8 +260,7 @@ func produceTestsDataFromFile(file_descr *os.File) []testsData {
 			want:    igrf_res,
 			wantErr: false,
 		}
-		tests[i] = current_test
-		i++
+		tests = append(tests, current_test)
 	}
 	return tests
 }

igrf/igrf.go

@@ -1,4 +1,4 @@
-Lat  59.900 geodetic    Long -109.900     1.100 km                     
+Lat  59.900 geodetic    Long -109.900     1.100 km
    DATE       D   SV      I  SV      H    SV       X    SV       Y    SV       Z    SV      F    SV
  1900.5   33.01    -2  82.33   0   8436     4    7074     7    4596    -3   62650   -20  63216   -19
  1901.5   32.97    -2  82.33   0   8440     4    7081     7    4593    -3   62630   -20  63196   -19
@@ -115,13 +115,18 @@ Lat  59.900 geodetic    Long -109.900     1.100 km
  2012.5   14.63   -17  80.21  -4  10020    43    9695    54    2531   -37   58088  -111  58946  -102
  2013.5   14.35   -17  80.15  -4  10063    43    9749    54    2494   -37   57977  -111  58844  -102
  2014.5   14.07   -17  80.09  -4  10106    43    9803    54    2457   -37   57867  -111  58743  -101
- 2015.5   13.87    -7  80.03  -4  10152    48    9855    51    2434    -9   57760  -103  58645   -93
- 2016.5   13.75    -7  79.97  -4  10199    48    9907    51    2425    -9   57657  -103  58553   -93
- 2017.5   13.64    -7  79.90  -4  10247    48    9958    51    2416    -9   57555  -103  58460   -93
- 2018.5   13.52    -7  79.84  -4  10295    48   10009    51    2406    -9   57452  -103  58367   -93
- 2019.5   13.40    -7  79.78  -4  10342    48   10061    51    2397    -9   57350  -103  58275   -93
- 2020.5   13.29    -6  79.71  -4  10391    49   10112    52    2389    -7   57250   -96  58186   -85
- 2021.5   13.19    -6  79.65  -4  10439    49   10164    52    2382    -7   57155   -96  58100   -85
- 2022.5   13.09    -6  79.58  -4  10488    49   10215    52    2376    -7   57059   -96  58015   -85
- 2023.5   12.99    -6  79.52  -4  10537    49   10267    52    2369    -7   56963   -96  57930   -85
- 2024.5   12.90    -6  79.46  -4  10585    49   10318    52    2362    -7   56868   -96  57844   -85
+ 2015.5   13.87    -7  80.03  -4  10152    47    9855    51    2434   -10   57760  -103  58645   -93
+ 2016.5   13.75    -7  79.97  -4  10199    47    9907    51    2425   -10   57657  -103  58553   -93
+ 2017.5   13.63    -7  79.91  -4  10247    48    9958    51    2415   -10   57555  -103  58460   -93
+ 2018.5   13.51    -7  79.84  -4  10294    48   10009    51    2406   -10   57452  -103  58367   -93
+ 2019.5   13.40    -7  79.78  -4  10342    48   10060    51    2396   -10   57349  -103  58274   -93
+ 2020.5   13.28    -7  79.71  -4  10391    50   10113    54    2387    -8   57251   -94  58186   -83
+ 2021.5   13.17    -7  79.65  -4  10441    50   10166    54    2380    -8   57157   -94  58103   -83
+ 2022.5   13.07    -6  79.58  -4  10491    50   10220    54    2372    -8   57063   -94  58020   -83
+ 2023.5   12.96    -6  79.52  -4  10542    50   10273    54    2364    -8   56970   -94  57937   -83
+ 2024.5   12.85    -6  79.45  -4  10592    50   10327    54    2356    -8   56876   -94  57854   -83
+ 2025.5   12.73    -9  79.38  -4  10643    52   10382    57    2345   -15   56778  -102  57767   -91
+ 2026.5   12.58    -9  79.31  -4  10696    52   10439    57    2330   -15   56676  -102  57676   -91
+ 2027.5   12.44    -9  79.24  -4  10748    52   10496    57    2314   -15   56574  -102  57586   -90
+ 2028.5   12.29    -9  79.17  -4  10800    52   10552    57    2299   -15   56472  -102  57495   -90
+ 2029.5   12.15    -8  79.10  -4  10852    52   10609    57    2284   -15   56370  -102  57405   -90