Skip to content

Commit ef80f14

Browse files
authored
All values are calculated now. (#21)
All IGRF values are calculated now.
1 parent a3bb55f commit ef80f14

File tree

5 files changed

+209
-77
lines changed

5 files changed

+209
-77
lines changed

README.md

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -23,11 +23,11 @@ The output is of type `type IGRFresults struct`. Fields are:
2323
|InclinationSV|(I)|arcmin/yr|-2.59||
2424
|HzIntensity|Horizontal intensity (H)|nT|14626.6||
2525
|HorizontalSV|(H)|nT/yr|-16.8||
26-
|NorthComponent|(X)|nT|14181.2|Tested|
26+
|NorthComponent|(X)|nT|14181.2||
2727
|NorthSV|(X)|nT/yr|-28.2||
28-
|EastComponent|(Y)|nT|3581.8|Tested|
28+
|EastComponent|(Y)|nT|3581.8||
2929
|EastSV|(Y)|nT/yr|32.0||
30-
|VerticalComponent|(Z)|nT|51788.4|downward - Tested|
30+
|VerticalComponent|(Z)|nT|51788.4|downward|
3131
|VerticalSV|(Z)|nT/yr|80.1||
3232
|TotalIntensity|(F)|nT|53814.3||
3333
|TotalSV|(F)|nT/yr|71.8||

calc/shval3.go

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -149,3 +149,45 @@ func Shval3(flat, flon, elev float64, nmax int, gha, ghb *[]float64) (float64, f
149149
ztemp = ztemp*cd - aa*sd
150150
return x, y, z, xtemp, ytemp, ztemp
151151
}
152+
153+
// Computes the geomagnetic d, i, h, and f from x, y, and z.
154+
// D - declination
155+
// I - inclination
156+
// H - horizontal intensity
157+
// F - total intensity
158+
func Dihf(x, y, z float64) (float64, float64, float64, float64) {
159+
var d, i, h, f float64
160+
sn := 0.0001
161+
for j := 1; j <= 1; j++ {
162+
h2 := x*x + y*y
163+
argument := h2
164+
// calculate horizontal intensity
165+
h = math.Sqrt(argument)
166+
argument = h2 + z*z
167+
// calculate total intensity
168+
f = math.Sqrt(argument)
169+
if f < sn {
170+
// If d and i cannot be determined
171+
d = math.NaN()
172+
// set equal to NaN
173+
i = math.NaN()
174+
} else {
175+
argument = z
176+
argument2 := h
177+
i = math.Atan2(argument, argument2)
178+
if h < sn {
179+
d = math.NaN()
180+
} else {
181+
hpx := h + x
182+
if hpx < sn {
183+
d = math.Pi
184+
} else {
185+
argument = y
186+
argument2 = hpx
187+
d = 2.0 * math.Atan2(argument, argument2)
188+
}
189+
}
190+
}
191+
}
192+
return d, i, h, f
193+
}

igrf/igrf.go

Lines changed: 85 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@ import (
44
"errors"
55
"fmt"
66
"log"
7+
"math"
78

89
"github.com/proway2/go-igrf/calc"
910
"github.com/proway2/go-igrf/coeffs"
@@ -18,8 +19,6 @@ func IGRF(lat, lon, alt, date float64) (IGRFresults, error) {
1819
if err := checkInitialConditions(lat, lon, alt); err != nil {
1920
return IGRFresults{}, err
2021
}
21-
// colat := float64(90.0 - lat)
22-
// alt64, colat, sd, cd := gg2geo(float64(alt), float64(colat))
2322
shc, err := coeffs.NewCoeffsData()
2423
if err != nil {
2524
log.Fatal(err.Error())
@@ -28,16 +27,86 @@ func IGRF(lat, lon, alt, date float64) (IGRFresults, error) {
2827
if err != nil {
2928
return IGRFresults{}, err
3029
}
31-
_ = *start_coeffs
32-
_ = *end_coeffs
3330
x, y, z, xtemp, ytemp, ztemp := calc.Shval3(lat, lon, alt, nmax, start_coeffs, end_coeffs)
34-
_ = x
35-
_ = y
36-
_ = z
37-
_ = xtemp
38-
_ = ytemp
39-
_ = ztemp
40-
res := IGRFresults{NorthComponent: float32(x), EastComponent: float32(y), VerticalComponent: float32(z)}
31+
d, i, h, f := calc.Dihf(x, y, z)
32+
33+
dtemp, itemp, htemp, ftemp := calc.Dihf(xtemp, ytemp, ztemp)
34+
35+
ddot := rad2deg(dtemp - d)
36+
if ddot > 180.0 {
37+
ddot -= 360.0
38+
}
39+
if ddot <= -180.0 {
40+
ddot += 360.0
41+
}
42+
ddot *= 60.0
43+
44+
idot := rad2deg(itemp-i) * 60
45+
d = rad2deg(d)
46+
i = rad2deg(i)
47+
hdot := htemp - h
48+
xdot := xtemp - x
49+
ydot := ytemp - y
50+
zdot := ztemp - z
51+
fdot := ftemp - f
52+
53+
// deal with geographic and magnetic poles
54+
55+
// at magnetic poles
56+
if h < 100.0 {
57+
d = math.NaN()
58+
ddot = math.NaN()
59+
/* while rest is ok */
60+
}
61+
// warn_H := 0
62+
// warn_H_val := 99999.0
63+
// var warn_H_strong int
64+
// warn_H_strong_val := 99999.0
65+
// // warn_P := 0
66+
// if h < 1000.0 {
67+
// // warn_H = 0
68+
// warn_H_strong = 1
69+
// if h < warn_H_strong_val {
70+
// warn_H_strong_val = h
71+
// }
72+
// // } else if h < 5000.0 && !warn_H_strong {
73+
// } else if h < 5000.0 && warn_H_strong != 0 {
74+
// // warn_H = 1
75+
// if h < warn_H_val {
76+
// warn_H_val = h
77+
// }
78+
// }
79+
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+
}
93+
94+
res := IGRFresults{
95+
Declination: float32(d),
96+
DeclinationSV: float32(ddot),
97+
Inclination: float32(i),
98+
InclinationSV: float32(idot),
99+
HorizontalIntensity: float32(h),
100+
HorizontalSV: float32(hdot),
101+
NorthComponent: float32(x),
102+
NorthSV: float32(xdot),
103+
EastComponent: float32(y),
104+
EastSV: float32(ydot),
105+
VerticalComponent: float32(z),
106+
VerticalSV: float32(zdot),
107+
TotalIntensity: float32(f),
108+
TotalSV: float32(fdot),
109+
}
41110
return res, nil
42111
}
43112

@@ -58,55 +127,8 @@ func checkInitialConditions(lat, lon, alt float64) error {
58127
return nil
59128
}
60129

61-
// // gg2geo - computes geocentric colatitude and radius from geodetic colatitude and height. Uses WGS-84 ellipsoid parameters.
62-
// //
63-
// // Inputs:
64-
// // h - altitude in kilometers
65-
// // gdcolat - geodetic colatitude
66-
// //
67-
// // Outputs:
68-
// // radius - Geocentric radius in kilometers.
69-
// // theta - Geocentric colatitude in degrees.
70-
// // sd - rotate B_X to gd_lat
71-
// // cd - rotate B_Z to gd_lat
72-
// //
73-
// // References:
74-
// // Equations (51)-(53) from "The main field" (chapter 4) by Langel, R. A. in: "Geomagnetism", Volume 1, Jacobs, J. A., Academic Press, 1987.
75-
// // Malin, S.R.C. and Barraclough, D.R., 1981. An algorithm for synthesizing the geomagnetic field. Computers & Geosciences, 7(4), pp.401-405.
76-
// func gg2geo(h, gdcolat float64) (radius, theta, sd, cd float64) {
77-
// eqrad := 6378.137 // equatorial radius
78-
// flat := 1 / 298.257223563
79-
// plrad := eqrad * (1 - flat) // polar radius
80-
// ctgd := math.Cos(deg2rad(gdcolat))
81-
// stgd := math.Sin(deg2rad(gdcolat))
82-
83-
// a2 := eqrad * eqrad
84-
// a4 := a2 * a2
85-
// b2 := plrad * plrad
86-
// b4 := b2 * b2
87-
// c2 := ctgd * ctgd
88-
// s2 := 1 - c2
89-
// rho := math.Sqrt(a2*s2 + b2*c2)
90-
91-
// rad := math.Sqrt(h*(h+2*rho) + (a4*s2+b4*c2)/math.Pow(rho, 2))
92-
93-
// cd = (h + rho) / rad
94-
// sd = (a2 - b2) * ctgd * stgd / (rho * rad)
95-
96-
// cthc := ctgd*cd - stgd*sd // Also: sthc = stgd*cd + ctgd*sd
97-
// theta = rad2deg(math.Acos(cthc)) // arccos returns values in [0, pi]
98-
99-
// return rad, theta, sd, cd
100-
// }
101-
102-
// // deg2rad - converts `degrees` into radians.
103-
// func deg2rad(degrees float64) float64 {
104-
// rad := degrees * math.Pi / 180.0
105-
// return rad
106-
// }
107-
108-
// // rad2deg - converts `radians` into degrees.
109-
// func rad2deg(radians float64) float64 {
110-
// deg := radians * 180.0 / math.Pi
111-
// return deg
112-
// }
130+
// rad2deg - converts `radians` into degrees.
131+
func rad2deg(radians float64) float64 {
132+
deg := radians * 180.0 / math.Pi
133+
return deg
134+
}

igrf/igrf_edge_test.go

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -56,7 +56,7 @@ func TestIGRFEdgeCases(t *testing.T) {
5656
},
5757
// {
5858
// name: "Testing",
59-
// args: args{lat: -59.9, lon: -39.9, alt: -0.5, date: 1915.5},
59+
// args: args{lat: 59.9, lon: 39.9, alt: -0.5, date: 2015.5},
6060
// want: IGRFresults{},
6161
// wantErr: false,
6262
// },

igrf/igrf_test.go

Lines changed: 78 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,11 @@ type testsData struct {
2828

2929
const dir_path string = "../testdata"
3030
const max_allowed_error = 0.2 // %
31+
// SV fields are just inregers in FORTRAN, so there might be situations where:
32+
// calculated value 16.47, reference 17
33+
// calculated value 2.52, reference 3
34+
// ...
35+
const max_sv_error = 50 // %
3136

3237
func TestIGRFDataCases(t *testing.T) {
3338
tests := getTestData()
@@ -38,28 +43,91 @@ func TestIGRFDataCases(t *testing.T) {
3843
t.Errorf("IGRF() error = %v, wantErr %v", err, tt.wantErr)
3944
return
4045
}
41-
compare_x := compareFloats(float64(got.NorthComponent), float64(tt.want.NorthComponent), max_allowed_error)
42-
if !compare_x {
46+
// 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)
50+
if !compare {
51+
t.Errorf("IGRF() Declination = %v, want %v", got.Declination, tt.want.Declination)
52+
}
53+
// Declination SV
54+
compare = compareFloats(math.Round(float64(got.DeclinationSV)), float64(tt.want.DeclinationSV), max_sv_error)
55+
if !compare {
56+
t.Errorf("IGRF() DeclinationSV = %v, want %v", got.DeclinationSV, tt.want.DeclinationSV)
57+
}
58+
// Inclination
59+
in := convertToPrecision(float64(got.Inclination), 2)
60+
compare = compareFloats(in, float64(tt.want.Inclination), max_allowed_error)
61+
if !compare {
62+
t.Errorf("IGRF() Inclination = %v, want %v", got.Inclination, tt.want.Inclination)
63+
}
64+
// Inclination SV
65+
compare = compareFloats(math.Round(float64(got.InclinationSV)), float64(tt.want.InclinationSV), max_sv_error)
66+
if !compare {
67+
t.Errorf("IGRF() InclinationSV = %v, want %v", got.InclinationSV, tt.want.InclinationSV)
68+
}
69+
// Horizontal intensity
70+
compare = compareFloats(float64(got.HorizontalIntensity), float64(tt.want.HorizontalIntensity), max_allowed_error)
71+
if !compare {
72+
t.Errorf("IGRF() Horizontal intensity = %v, want %v", got.HorizontalIntensity, tt.want.HorizontalIntensity)
73+
}
74+
// Horizontal SV
75+
compare = compareFloats(math.Round(float64(got.HorizontalSV)), float64(tt.want.HorizontalSV), max_sv_error)
76+
if !compare {
77+
t.Errorf("IGRF() HorizontalSV = %v, want %v", got.HorizontalSV, tt.want.HorizontalSV)
78+
}
79+
// North component
80+
compare = compareFloats(float64(got.NorthComponent), float64(tt.want.NorthComponent), max_allowed_error)
81+
if !compare {
4382
t.Errorf("IGRF() NorthComponent = %v, want %v", got.NorthComponent, tt.want.NorthComponent)
4483
}
45-
// Fortran results are rounded!!!
84+
// North SV
85+
compare = compareFloats(math.Round(float64(got.NorthSV)), float64(tt.want.NorthSV), max_sv_error)
86+
if !compare {
87+
t.Errorf("IGRF() NorthSV = %v, want %v", got.NorthSV, tt.want.NorthSV)
88+
}
89+
// East Component - FORTRAN results are rounded!!!
4690
new_east_got := math.Round(float64(got.EastComponent))
4791
new_east_want := math.Round(float64(tt.want.EastComponent))
48-
compare_y := compareFloats(new_east_got, new_east_want, max_allowed_error)
49-
if !compare_y {
92+
compare = compareFloats(new_east_got, new_east_want, max_allowed_error)
93+
if !compare {
5094
t.Errorf("IGRF() EastComponent = %v, want %v", got.EastComponent, tt.want.EastComponent)
5195
}
52-
compare_z := compareFloats(float64(got.VerticalComponent), float64(tt.want.VerticalComponent), max_allowed_error)
53-
if !compare_z {
96+
// EastSV
97+
compare = compareFloats(math.Round(float64(got.EastSV)), float64(tt.want.EastSV), max_sv_error)
98+
if !compare {
99+
t.Errorf("IGRF() EastSV = %v, want %v", got.EastSV, tt.want.EastSV)
100+
}
101+
// Vertical component
102+
compare = compareFloats(float64(got.VerticalComponent), float64(tt.want.VerticalComponent), max_allowed_error)
103+
if !compare {
54104
t.Errorf("IGRF() VerticalComponent = %v, want %v", got.VerticalComponent, tt.want.VerticalComponent)
55105
}
56-
// if !reflect.DeepEqual(got, tt.want) {
57-
// t.Errorf("IGRF() = %v, want %v", got, tt.want)
58-
// }
106+
// VerticalSV
107+
compare = compareFloats(math.Round(float64(got.VerticalSV)), float64(tt.want.VerticalSV), max_sv_error)
108+
if !compare {
109+
t.Errorf("IGRF() VerticalSV = %v, want %v", got.VerticalSV, tt.want.VerticalSV)
110+
}
111+
// Total intensity
112+
compare = compareFloats(float64(got.TotalIntensity), float64(tt.want.TotalIntensity), max_allowed_error)
113+
if !compare {
114+
t.Errorf("IGRF() Total = %v, want %v", got.TotalIntensity, tt.want.TotalIntensity)
115+
}
116+
// TotalSV
117+
compare = compareFloats(math.Round(float64(got.TotalSV)), float64(tt.want.TotalSV), max_sv_error)
118+
if !compare {
119+
t.Errorf("IGRF() TotalSV = %v, want %v", got.TotalSV, tt.want.TotalSV)
120+
}
59121
})
60122
}
61123
}
62124

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)
129+
}
130+
63131
func compareFloats(check, base, allowable_error float64) bool {
64132
value1 := math.Abs(check)
65133
value2 := math.Abs(base)

0 commit comments

Comments
 (0)