Skip to content

Commit 4639be4

Browse files
proway2proway2
authored
Near pole lats (#23)
Near pole values are calculated now.
1 parent ef80f14 commit 4639be4

File tree

5 files changed

+333
-51
lines changed

5 files changed

+333
-51
lines changed

‎README.md‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,21 @@ The output is of type `type IGRFresults struct`. Fields are:
3232
|TotalIntensity|(F)|nT|53814.3||
3333
|TotalSV|(F)|nT/yr|71.8||
3434

35+
### Annual changes (SV values)
36+
37+
The SV values are computed by subtracting the values for the desired input date from corresponding values one year later.
38+
39+
### Values near geographic poles
40+
41+
Unlike in `C` implementation, this software calculates values near geographic poles, e.g. for latitudes higher than 89.999 and -89.999. This is the same approach the reference `FORTRAN` implementaion does have. Be adviced that near pole values are much less accurate.
42+
43+
## Overall accuracy
44+
45+
There are almost 900 unittests and results are compared against those generated from FORTRAN.
46+
- Max allowed absolute error: 0.005
47+
- Max allowed relative error: 0.2 %
48+
Since FORTRAN rounds almost all values, except `D` and `I`, actual results are of much higher accuracy.
49+
3550
## References
3651

3752
- [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

‎igrf/igrf.go‎

Lines changed: 19 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -52,12 +52,12 @@ func IGRF(lat, lon, alt, date float64) (IGRFresults, error) {
5252

5353
// deal with geographic and magnetic poles
5454

55-
// at magnetic poles
56-
if h < 100.0 {
57-
d = math.NaN()
58-
ddot = math.NaN()
59-
/* while rest is ok */
60-
}
55+
// // at magnetic poles
56+
// if h < 100.0 {
57+
// d = math.NaN()
58+
// ddot = math.NaN()
59+
// /* while rest is ok */
60+
// }
6161
// warn_H := 0
6262
// warn_H_val := 99999.0
6363
// var warn_H_strong int
@@ -77,19 +77,19 @@ func IGRF(lat, lon, alt, date float64) (IGRFresults, error) {
7777
// }
7878
// }
7979

80-
// at geographic poles
81-
if 90.0-math.Abs(lat) <= 0.001 {
82-
x = math.NaN()
83-
y = math.NaN()
84-
d = math.NaN()
85-
xdot = math.NaN()
86-
ydot = math.NaN()
87-
ddot = math.NaN()
88-
// warn_P = 1
89-
// warn_H = 0
90-
// warn_H_strong = 0
91-
/* while rest is ok */
92-
}
80+
// // at geographic poles
81+
// if 90.0-math.Abs(lat) <= 0.001 {
82+
// x = math.NaN()
83+
// y = math.NaN()
84+
// d = math.NaN()
85+
// xdot = math.NaN()
86+
// ydot = math.NaN()
87+
// ddot = math.NaN()
88+
// // warn_P = 1
89+
// // warn_H = 0
90+
// // warn_H_strong = 0
91+
// /* while rest is ok */
92+
// }
9393

9494
res := IGRFresults{
9595
Declination: float32(d),

‎igrf/igrf_test.go‎

Lines changed: 45 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -27,8 +27,17 @@ type testsData struct {
2727
}
2828

2929
const dir_path string = "../testdata"
30-
const max_allowed_error = 0.2 // %
31-
// SV fields are just inregers in FORTRAN, so there might be situations where:
30+
31+
// allowed errors
32+
const (
33+
max_allowed_error = 0.2 // %
34+
near_pole_max_allower_error = 3.25 // %
35+
max_absolute_error = 0.005
36+
)
37+
38+
const near_pole_tolerance = 0.001
39+
40+
// SV fields are just integers in FORTRAN, so there might be situations where:
3241
// calculated value 16.47, reference 17
3342
// calculated value 2.52, reference 3
3443
// ...
@@ -43,97 +52,105 @@ func TestIGRFDataCases(t *testing.T) {
4352
t.Errorf("IGRF() error = %v, wantErr %v", err, tt.wantErr)
4453
return
4554
}
55+
allowed_error := getMaxAllowedRelativeError(tt.args.lat)
4656
// Declination
47-
// there are just two digits after the dot in FORTRAN results, so truncating it
48-
dn := convertToPrecision(float64(got.Declination), 2)
49-
compare := compareFloats(dn, float64(tt.want.Declination), max_allowed_error)
57+
compare := compareFloats(float64(got.Declination), float64(tt.want.Declination), allowed_error, max_absolute_error)
5058
if !compare {
5159
t.Errorf("IGRF() Declination = %v, want %v", got.Declination, tt.want.Declination)
5260
}
5361
// Declination SV
54-
compare = compareFloats(math.Round(float64(got.DeclinationSV)), float64(tt.want.DeclinationSV), max_sv_error)
62+
compare = compareFloats(math.Round(float64(got.DeclinationSV)), float64(tt.want.DeclinationSV), max_sv_error, max_absolute_error)
5563
if !compare {
5664
t.Errorf("IGRF() DeclinationSV = %v, want %v", got.DeclinationSV, tt.want.DeclinationSV)
5765
}
5866
// Inclination
59-
in := convertToPrecision(float64(got.Inclination), 2)
60-
compare = compareFloats(in, float64(tt.want.Inclination), max_allowed_error)
67+
compare = compareFloats(float64(got.Inclination), float64(tt.want.Inclination), allowed_error, max_absolute_error)
6168
if !compare {
6269
t.Errorf("IGRF() Inclination = %v, want %v", got.Inclination, tt.want.Inclination)
6370
}
6471
// Inclination SV
65-
compare = compareFloats(math.Round(float64(got.InclinationSV)), float64(tt.want.InclinationSV), max_sv_error)
72+
compare = compareFloats(math.Round(float64(got.InclinationSV)), float64(tt.want.InclinationSV), max_sv_error, max_absolute_error)
6673
if !compare {
6774
t.Errorf("IGRF() InclinationSV = %v, want %v", got.InclinationSV, tt.want.InclinationSV)
6875
}
6976
// Horizontal intensity
70-
compare = compareFloats(float64(got.HorizontalIntensity), float64(tt.want.HorizontalIntensity), max_allowed_error)
77+
compare = compareFloats(float64(got.HorizontalIntensity), float64(tt.want.HorizontalIntensity), allowed_error, max_absolute_error)
7178
if !compare {
7279
t.Errorf("IGRF() Horizontal intensity = %v, want %v", got.HorizontalIntensity, tt.want.HorizontalIntensity)
7380
}
7481
// Horizontal SV
75-
compare = compareFloats(math.Round(float64(got.HorizontalSV)), float64(tt.want.HorizontalSV), max_sv_error)
82+
compare = compareFloats(math.Round(float64(got.HorizontalSV)), float64(tt.want.HorizontalSV), max_sv_error, max_absolute_error)
7683
if !compare {
7784
t.Errorf("IGRF() HorizontalSV = %v, want %v", got.HorizontalSV, tt.want.HorizontalSV)
7885
}
7986
// North component
80-
compare = compareFloats(float64(got.NorthComponent), float64(tt.want.NorthComponent), max_allowed_error)
87+
compare = compareFloats(float64(got.NorthComponent), float64(tt.want.NorthComponent), allowed_error, max_absolute_error)
8188
if !compare {
8289
t.Errorf("IGRF() NorthComponent = %v, want %v", got.NorthComponent, tt.want.NorthComponent)
8390
}
8491
// North SV
85-
compare = compareFloats(math.Round(float64(got.NorthSV)), float64(tt.want.NorthSV), max_sv_error)
92+
compare = compareFloats(math.Round(float64(got.NorthSV)), float64(tt.want.NorthSV), max_sv_error, max_absolute_error)
8693
if !compare {
8794
t.Errorf("IGRF() NorthSV = %v, want %v", got.NorthSV, tt.want.NorthSV)
8895
}
8996
// East Component - FORTRAN results are rounded!!!
9097
new_east_got := math.Round(float64(got.EastComponent))
9198
new_east_want := math.Round(float64(tt.want.EastComponent))
92-
compare = compareFloats(new_east_got, new_east_want, max_allowed_error)
99+
// compare = compareFloats(float64(got.EastComponent), float64(tt.want.EastComponent), allowed_error, max_absolute_error)
100+
compare = compareFloats(new_east_got, new_east_want, allowed_error, max_absolute_error)
93101
if !compare {
94102
t.Errorf("IGRF() EastComponent = %v, want %v", got.EastComponent, tt.want.EastComponent)
95103
}
96104
// EastSV
97-
compare = compareFloats(math.Round(float64(got.EastSV)), float64(tt.want.EastSV), max_sv_error)
105+
compare = compareFloats(math.Round(float64(got.EastSV)), float64(tt.want.EastSV), max_sv_error, max_absolute_error)
98106
if !compare {
99107
t.Errorf("IGRF() EastSV = %v, want %v", got.EastSV, tt.want.EastSV)
100108
}
101109
// Vertical component
102-
compare = compareFloats(float64(got.VerticalComponent), float64(tt.want.VerticalComponent), max_allowed_error)
110+
compare = compareFloats(float64(got.VerticalComponent), float64(tt.want.VerticalComponent), allowed_error, max_absolute_error)
103111
if !compare {
104112
t.Errorf("IGRF() VerticalComponent = %v, want %v", got.VerticalComponent, tt.want.VerticalComponent)
105113
}
106114
// VerticalSV
107-
compare = compareFloats(math.Round(float64(got.VerticalSV)), float64(tt.want.VerticalSV), max_sv_error)
115+
compare = compareFloats(math.Round(float64(got.VerticalSV)), float64(tt.want.VerticalSV), max_sv_error, max_absolute_error)
108116
if !compare {
109117
t.Errorf("IGRF() VerticalSV = %v, want %v", got.VerticalSV, tt.want.VerticalSV)
110118
}
111119
// Total intensity
112-
compare = compareFloats(float64(got.TotalIntensity), float64(tt.want.TotalIntensity), max_allowed_error)
120+
compare = compareFloats(float64(got.TotalIntensity), float64(tt.want.TotalIntensity), allowed_error, max_absolute_error)
113121
if !compare {
114122
t.Errorf("IGRF() Total = %v, want %v", got.TotalIntensity, tt.want.TotalIntensity)
115123
}
116124
// TotalSV
117-
compare = compareFloats(math.Round(float64(got.TotalSV)), float64(tt.want.TotalSV), max_sv_error)
125+
compare = compareFloats(math.Round(float64(got.TotalSV)), float64(tt.want.TotalSV), max_sv_error, max_absolute_error)
118126
if !compare {
119127
t.Errorf("IGRF() TotalSV = %v, want %v", got.TotalSV, tt.want.TotalSV)
120128
}
121129
})
122130
}
123131
}
124132

125-
func convertToPrecision(value float64, precision int) float64 {
126-
format_verbs := fmt.Sprintf("%%.%vf", precision)
127-
vs := fmt.Sprintf(format_verbs, value)
128-
return toFloat64(vs)
133+
func isNearPole(lat float64) bool {
134+
return 90.0-math.Abs(lat) <= near_pole_tolerance
135+
}
136+
137+
func getMaxAllowedRelativeError(lat float64) float64 {
138+
if isNearPole(lat) {
139+
return near_pole_max_allower_error
140+
}
141+
return max_allowed_error
129142
}
130143

131-
func compareFloats(check, base, allowable_error float64) bool {
144+
func compareFloats(check, base, allowable_relative_error, allowable_absolute_error float64) bool {
145+
if math.IsNaN(check) && !math.IsNaN(base) {
146+
return false
147+
}
132148
value1 := math.Abs(check)
133149
value2 := math.Abs(base)
134-
calc_err := math.Abs(100 * ((value1 - value2) / value2))
150+
absolute_error := math.Abs(value1 - value2)
151+
calc_err := 100 * absolute_error / math.Abs(value2)
135152
// allowable error percent
136-
if calc_err > allowable_error {
153+
if calc_err > allowable_relative_error && absolute_error > allowable_absolute_error {
137154
return false
138155
}
139156
return true
@@ -153,17 +170,13 @@ func discoverTestData() []string {
153170

154171
func getTestData() []testsData {
155172
test_data_files := discoverTestData()
156-
tests := make([]testsData, 625)
157-
var index int
173+
tests := make([]testsData, 0)
158174
for _, file := range test_data_files {
159175
f, err := os.Open(file)
160176
defer f.Close()
161177
check(err)
162178
current_file_tests := produceTestsDataFromFile(f)
163-
for _, test := range current_file_tests {
164-
tests[index] = test
165-
index++
166-
}
179+
tests = append(tests, current_file_tests...)
167180
}
168181
return tests
169182
}

0 commit comments

Comments
 (0)