Skip to content

Commit

Permalink
All values are calculated now. (#21)
Browse files Browse the repository at this point in the history 
All IGRF values are calculated now.
  • Loading branch information
@proway2
proway2 authored Apr 27, 2022
1 parent a3bb55f commit ef80f14
Show file tree
Hide file tree
Showing 5 changed files with 209 additions and 77 deletions.
6 changes: 3 additions & 3 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,11 +23,11 @@ The output is of type `type IGRFresults struct`. Fields are:
|InclinationSV|(I)|arcmin/yr|-2.59||
|HzIntensity|Horizontal intensity (H)|nT|14626.6||
|HorizontalSV|(H)|nT/yr|-16.8||
|NorthComponent|(X)|nT|14181.2|Tested|
|NorthComponent|(X)|nT|14181.2||
|NorthSV|(X)|nT/yr|-28.2||
|EastComponent|(Y)|nT|3581.8|Tested|
|EastComponent|(Y)|nT|3581.8||
|EastSV|(Y)|nT/yr|32.0||
|VerticalComponent|(Z)|nT|51788.4|downward - Tested|
|VerticalComponent|(Z)|nT|51788.4|downward|
|VerticalSV|(Z)|nT/yr|80.1||
|TotalIntensity|(F)|nT|53814.3||
|TotalSV|(F)|nT/yr|71.8||
Expand Down
42 changes: 42 additions & 0 deletions calc/shval3.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,3 +149,45 @@ func Shval3(flat, flon, elev float64, nmax int, gha, ghb *[]float64) (float64, f
ztemp = ztemp*cd - aa*sd
return x, y, z, xtemp, ytemp, ztemp
}

// Computes the geomagnetic d, i, h, and f from x, y, and z.
// D - declination
// I - inclination
// H - horizontal intensity
// F - total intensity
func Dihf(x, y, z float64) (float64, float64, float64, float64) {
var d, i, h, f float64
sn := 0.0001
for j := 1; j <= 1; j++ {
h2 := x*x + y*y
argument := h2
// calculate horizontal intensity
h = math.Sqrt(argument)
argument = h2 + z*z
// calculate total intensity
f = math.Sqrt(argument)
if f < sn {
// If d and i cannot be determined
d = math.NaN()
// set equal to NaN
i = math.NaN()
} else {
argument = z
argument2 := h
i = math.Atan2(argument, argument2)
if h < sn {
d = math.NaN()
} else {
hpx := h + x
if hpx < sn {
d = math.Pi
} else {
argument = y
argument2 = hpx
d = 2.0 * math.Atan2(argument, argument2)
}
}
}
}
return d, i, h, f
}
148 changes: 85 additions & 63 deletions igrf/igrf.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ import (
"errors"
"fmt"
"log"
"math"

"github.com/proway2/go-igrf/calc"
"github.com/proway2/go-igrf/coeffs"
Expand All @@ -18,8 +19,6 @@ func IGRF(lat, lon, alt, date float64) (IGRFresults, error) {
if err := checkInitialConditions(lat, lon, alt); err != nil {
return IGRFresults{}, err
}
// colat := float64(90.0 - lat)
// alt64, colat, sd, cd := gg2geo(float64(alt), float64(colat))
shc, err := coeffs.NewCoeffsData()
if err != nil {
log.Fatal(err.Error())
Expand All @@ -28,16 +27,86 @@ func IGRF(lat, lon, alt, date float64) (IGRFresults, error) {
if err != nil {
return IGRFresults{}, err
}
_ = *start_coeffs
_ = *end_coeffs
x, y, z, xtemp, ytemp, ztemp := calc.Shval3(lat, lon, alt, nmax, start_coeffs, end_coeffs)
_ = x
_ = y
_ = z
_ = xtemp
_ = ytemp
_ = ztemp
res := IGRFresults{NorthComponent: float32(x), EastComponent: float32(y), VerticalComponent: float32(z)}
d, i, h, f := calc.Dihf(x, y, z)

dtemp, itemp, htemp, ftemp := calc.Dihf(xtemp, ytemp, ztemp)

ddot := rad2deg(dtemp - d)
if ddot > 180.0 {
ddot -= 360.0
}
if ddot <= -180.0 {
ddot += 360.0
}
ddot *= 60.0

idot := rad2deg(itemp-i) * 60
d = rad2deg(d)
i = rad2deg(i)
hdot := htemp - h
xdot := xtemp - x
ydot := ytemp - y
zdot := ztemp - z
fdot := ftemp - f

// deal with geographic and magnetic poles

// at magnetic poles
if h < 100.0 {
d = math.NaN()
ddot = math.NaN()
/* while rest is ok */
}
// warn_H := 0
// warn_H_val := 99999.0
// var warn_H_strong int
// warn_H_strong_val := 99999.0
// // warn_P := 0
// if h < 1000.0 {
// // warn_H = 0
// warn_H_strong = 1
// if h < warn_H_strong_val {
// warn_H_strong_val = h
// }
// // } else if h < 5000.0 && !warn_H_strong {
// } else if h < 5000.0 && warn_H_strong != 0 {
// // warn_H = 1
// if h < warn_H_val {
// warn_H_val = h
// }
// }

// at geographic poles
if 90.0-math.Abs(lat) <= 0.001 {
x = math.NaN()
y = math.NaN()
d = math.NaN()
xdot = math.NaN()
ydot = math.NaN()
ddot = math.NaN()
// warn_P = 1
// warn_H = 0
// warn_H_strong = 0
/* while rest is ok */
}

res := IGRFresults{
Declination: float32(d),
DeclinationSV: float32(ddot),
Inclination: float32(i),
InclinationSV: float32(idot),
HorizontalIntensity: float32(h),
HorizontalSV: float32(hdot),
NorthComponent: float32(x),
NorthSV: float32(xdot),
EastComponent: float32(y),
EastSV: float32(ydot),
VerticalComponent: float32(z),
VerticalSV: float32(zdot),
TotalIntensity: float32(f),
TotalSV: float32(fdot),
}
return res, nil
}

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

// // gg2geo - computes geocentric colatitude and radius from geodetic colatitude and height. Uses WGS-84 ellipsoid parameters.
// //
// // Inputs:
// // h - altitude in kilometers
// // gdcolat - geodetic colatitude
// //
// // Outputs:
// // radius - Geocentric radius in kilometers.
// // theta - Geocentric colatitude in degrees.
// // sd - rotate B_X to gd_lat
// // cd - rotate B_Z to gd_lat
// //
// // References:
// // Equations (51)-(53) from "The main field" (chapter 4) by Langel, R. A. in: "Geomagnetism", Volume 1, Jacobs, J. A., Academic Press, 1987.
// // Malin, S.R.C. and Barraclough, D.R., 1981. An algorithm for synthesizing the geomagnetic field. Computers & Geosciences, 7(4), pp.401-405.
// func gg2geo(h, gdcolat float64) (radius, theta, sd, cd float64) {
// eqrad := 6378.137 // equatorial radius
// flat := 1 / 298.257223563
// plrad := eqrad * (1 - flat) // polar radius
// ctgd := math.Cos(deg2rad(gdcolat))
// stgd := math.Sin(deg2rad(gdcolat))

// a2 := eqrad * eqrad
// a4 := a2 * a2
// b2 := plrad * plrad
// b4 := b2 * b2
// c2 := ctgd * ctgd
// s2 := 1 - c2
// rho := math.Sqrt(a2*s2 + b2*c2)

// rad := math.Sqrt(h*(h+2*rho) + (a4*s2+b4*c2)/math.Pow(rho, 2))

// cd = (h + rho) / rad
// sd = (a2 - b2) * ctgd * stgd / (rho * rad)

// cthc := ctgd*cd - stgd*sd // Also: sthc = stgd*cd + ctgd*sd
// theta = rad2deg(math.Acos(cthc)) // arccos returns values in [0, pi]

// return rad, theta, sd, cd
// }

// // deg2rad - converts `degrees` into radians.
// func deg2rad(degrees float64) float64 {
// rad := degrees * math.Pi / 180.0
// return rad
// }

// // rad2deg - converts `radians` into degrees.
// func rad2deg(radians float64) float64 {
// deg := radians * 180.0 / math.Pi
// return deg
// }
// rad2deg - converts `radians` into degrees.
func rad2deg(radians float64) float64 {
deg := radians * 180.0 / math.Pi
return deg
}
2 changes: 1 addition & 1 deletion igrf/igrf_edge_test.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ func TestIGRFEdgeCases(t *testing.T) {
},
// {
// name: "Testing",
// args: args{lat: -59.9, lon: -39.9, alt: -0.5, date: 1915.5},
// args: args{lat: 59.9, lon: 39.9, alt: -0.5, date: 2015.5},
// want: IGRFresults{},
// wantErr: false,
// },
Expand Down
88 changes: 78 additions & 10 deletions igrf/igrf_test.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,11 @@ type testsData struct {

const dir_path string = "../testdata"
const max_allowed_error = 0.2 // %
// SV fields are just inregers in FORTRAN, so there might be situations where:
// calculated value 16.47, reference 17
// calculated value 2.52, reference 3
// ...
const max_sv_error = 50 // %

func TestIGRFDataCases(t *testing.T) {
tests := getTestData()
Expand All @@ -38,28 +43,91 @@ func TestIGRFDataCases(t *testing.T) {
t.Errorf("IGRF() error = %v, wantErr %v", err, tt.wantErr)
return
}
compare_x := compareFloats(float64(got.NorthComponent), float64(tt.want.NorthComponent), max_allowed_error)
if !compare_x {
// Declination
// there are just two digits after the dot in FORTRAN results, so truncating it
dn := convertToPrecision(float64(got.Declination), 2)
compare := compareFloats(dn, float64(tt.want.Declination), max_allowed_error)
if !compare {
t.Errorf("IGRF() Declination = %v, want %v", got.Declination, tt.want.Declination)
}
// Declination SV
compare = compareFloats(math.Round(float64(got.DeclinationSV)), float64(tt.want.DeclinationSV), max_sv_error)
if !compare {
t.Errorf("IGRF() DeclinationSV = %v, want %v", got.DeclinationSV, tt.want.DeclinationSV)
}
// Inclination
in := convertToPrecision(float64(got.Inclination), 2)
compare = compareFloats(in, float64(tt.want.Inclination), max_allowed_error)
if !compare {
t.Errorf("IGRF() Inclination = %v, want %v", got.Inclination, tt.want.Inclination)
}
// Inclination SV
compare = compareFloats(math.Round(float64(got.InclinationSV)), float64(tt.want.InclinationSV), max_sv_error)
if !compare {
t.Errorf("IGRF() InclinationSV = %v, want %v", got.InclinationSV, tt.want.InclinationSV)
}
// Horizontal intensity
compare = compareFloats(float64(got.HorizontalIntensity), float64(tt.want.HorizontalIntensity), max_allowed_error)
if !compare {
t.Errorf("IGRF() Horizontal intensity = %v, want %v", got.HorizontalIntensity, tt.want.HorizontalIntensity)
}
// Horizontal SV
compare = compareFloats(math.Round(float64(got.HorizontalSV)), float64(tt.want.HorizontalSV), max_sv_error)
if !compare {
t.Errorf("IGRF() HorizontalSV = %v, want %v", got.HorizontalSV, tt.want.HorizontalSV)
}
// North component
compare = compareFloats(float64(got.NorthComponent), float64(tt.want.NorthComponent), max_allowed_error)
if !compare {
t.Errorf("IGRF() NorthComponent = %v, want %v", got.NorthComponent, tt.want.NorthComponent)
}
// Fortran results are rounded!!!
// North SV
compare = compareFloats(math.Round(float64(got.NorthSV)), float64(tt.want.NorthSV), max_sv_error)
if !compare {
t.Errorf("IGRF() NorthSV = %v, want %v", got.NorthSV, tt.want.NorthSV)
}
// East Component - FORTRAN results are rounded!!!
new_east_got := math.Round(float64(got.EastComponent))
new_east_want := math.Round(float64(tt.want.EastComponent))
compare_y := compareFloats(new_east_got, new_east_want, max_allowed_error)
if !compare_y {
compare = compareFloats(new_east_got, new_east_want, max_allowed_error)
if !compare {
t.Errorf("IGRF() EastComponent = %v, want %v", got.EastComponent, tt.want.EastComponent)
}
compare_z := compareFloats(float64(got.VerticalComponent), float64(tt.want.VerticalComponent), max_allowed_error)
if !compare_z {
// EastSV
compare = compareFloats(math.Round(float64(got.EastSV)), float64(tt.want.EastSV), max_sv_error)
if !compare {
t.Errorf("IGRF() EastSV = %v, want %v", got.EastSV, tt.want.EastSV)
}
// Vertical component
compare = compareFloats(float64(got.VerticalComponent), float64(tt.want.VerticalComponent), max_allowed_error)
if !compare {
t.Errorf("IGRF() VerticalComponent = %v, want %v", got.VerticalComponent, tt.want.VerticalComponent)
}
// if !reflect.DeepEqual(got, tt.want) {
// t.Errorf("IGRF() = %v, want %v", got, tt.want)
// }
// VerticalSV
compare = compareFloats(math.Round(float64(got.VerticalSV)), float64(tt.want.VerticalSV), max_sv_error)
if !compare {
t.Errorf("IGRF() VerticalSV = %v, want %v", got.VerticalSV, tt.want.VerticalSV)
}
// Total intensity
compare = compareFloats(float64(got.TotalIntensity), float64(tt.want.TotalIntensity), max_allowed_error)
if !compare {
t.Errorf("IGRF() Total = %v, want %v", got.TotalIntensity, tt.want.TotalIntensity)
}
// TotalSV
compare = compareFloats(math.Round(float64(got.TotalSV)), float64(tt.want.TotalSV), max_sv_error)
if !compare {
t.Errorf("IGRF() TotalSV = %v, want %v", got.TotalSV, tt.want.TotalSV)
}
})
}
}

func convertToPrecision(value float64, precision int) float64 {
format_verbs := fmt.Sprintf("%%.%vf", precision)
vs := fmt.Sprintf(format_verbs, value)
return toFloat64(vs)
}

func compareFloats(check, base, allowable_error float64) bool {
value1 := math.Abs(check)
value2 := math.Abs(base)
Expand Down

0 comments on commit ef80f14

Please sign in to comment.