Skip to content

Commit 190c245

Browse files
proway2Michael S
proway2
and
Michael S
authored
X, Y, Z are now working (#16)
* Added `nmax` value to tests. * Maximum spherical harmonic degree is returned as well now. * Maximal spherical harmonic degree is handled now. * Fixed typo. * Switching to float64 inputs for `igrf` package. * Fixed automatic tests creation. * shval3 initial revision. * WIP: some calculations are started. * Fixed error with incorrect `date`. * Tests are now working for X, Y, Z * Some clean up * Making linter happy. Co-authored-by: Michael S <[email protected]>
1 parent 0ded6a4 commit 190c245

File tree

6 files changed

+303
-85
lines changed

6 files changed

+303
-85
lines changed

‎calc/shval3.go‎

Lines changed: 151 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,151 @@
1+
package calc
2+
3+
import "math"
4+
5+
// Calculates field components from spherical harmonic (sh) models.
6+
// x - northward component
7+
// y - eastward component
8+
// z - vertically-downward component
9+
func Shval3(flat, flon, elev float64, nmax int, gha, ghb *[]float64) (float64, float64, float64, float64, float64, float64) {
10+
// similar to shval3 from C implementation
11+
var earths_radius float64 = 6371.2
12+
var dtr float64 = 0.01745329
13+
/*
14+
a2,b2 - squares of semi-major and semi-minor axes of
15+
the reference spheroid used for transforming
16+
between geodetic and geocentric coordinates or components
17+
*/
18+
var a2 float64 = 40680631.59 /* WGS84 */
19+
var b2 float64 = 40408299.98 /* WGS84 */
20+
var x, y, z, xtemp, ytemp, ztemp, aa, aa_temp, argument, clat, slat, sd, bb, cc, dd, r, ratio, power, rr, fn, fm float64
21+
var l, n, m, npq int
22+
var sl, cl [14]float64
23+
var p, q [119]float64
24+
argument = flat * dtr
25+
slat = math.Sin(argument)
26+
if (90.0 - flat) < 0.001 {
27+
// 300 ft. from North pole
28+
aa = 89.999
29+
} else {
30+
if (90.0 + flat) < 0.001 {
31+
// 300 ft. from South pole
32+
aa = -89.999
33+
} else {
34+
aa = flat
35+
}
36+
}
37+
argument = aa * dtr
38+
clat = math.Cos(argument)
39+
argument = flon * dtr
40+
// TODO: Why this start from 1?
41+
sl[1] = math.Sin(argument)
42+
cl[1] = math.Cos(argument)
43+
l = 0 // in C index starts from 1
44+
n = 0
45+
m = 1
46+
npq = (nmax * (nmax + 3)) / 2
47+
48+
// this block is for geodetic coordinate system ->
49+
aa = a2 * clat * clat
50+
bb = b2 * slat * slat
51+
cc = aa + bb
52+
argument = cc
53+
dd = math.Sqrt(argument)
54+
argument = elev*(elev+2.0*dd) + (a2*aa+b2*bb)/cc
55+
r = math.Sqrt(argument)
56+
cd := (elev + dd) / r
57+
sd = (a2 - b2) / dd * slat * clat / r
58+
aa = slat
59+
slat = slat*cd - clat*sd
60+
clat = clat*cd + aa*sd
61+
// <- this block is for geodetic coordinate system
62+
ratio = earths_radius / r
63+
argument = 3.0
64+
aa = math.Sqrt(argument)
65+
// TODO: Why all these starts with 1?
66+
p[1] = 2.0 * slat
67+
p[2] = 2.0 * clat
68+
p[3] = 4.5*slat*slat - 1.5
69+
p[4] = 3.0 * aa * clat * slat
70+
q[1] = -clat
71+
q[2] = slat
72+
q[3] = -3.0 * clat * slat
73+
q[4] = aa * (slat*slat - clat*clat)
74+
75+
for k := 1; k <= npq; k++ {
76+
if n < m {
77+
m = 0
78+
n = n + 1
79+
argument = ratio
80+
power = float64(n + 2)
81+
rr = math.Pow(argument, power)
82+
fn = float64(n)
83+
}
84+
fm = float64(m)
85+
if k >= 5 {
86+
if m == n {
87+
argument = (1.0 - 0.5/fm)
88+
aa = math.Sqrt(argument)
89+
j := k - n - 1
90+
p[k] = (1.0 + 1.0/fm) * aa * clat * p[j]
91+
q[k] = aa * (clat*q[j] + slat/fm*p[j])
92+
sl[m] = sl[m-1]*cl[1] + cl[m-1]*sl[1]
93+
cl[m] = cl[m-1]*cl[1] - sl[m-1]*sl[1]
94+
} else {
95+
argument = fn*fn - fm*fm
96+
aa = math.Sqrt(argument)
97+
argument = ((fn - 1.0) * (fn - 1.0)) - (fm * fm)
98+
bb = math.Sqrt(argument) / aa
99+
cc = (2.0*fn - 1.0) / aa
100+
ii := k - n
101+
j := k - 2*n + 1
102+
p[k] = (fn + 1.0) * (cc*slat/fn*p[ii] - bb/(fn-1.0)*p[j])
103+
q[k] = cc*(slat*q[ii]-clat/fn*p[ii]) - bb*q[j]
104+
}
105+
}
106+
aa = rr * (*gha)[l]
107+
aa_temp = rr * (*ghb)[l]
108+
if m == 0 {
109+
x = x + aa*q[k]
110+
z = z - aa*p[k]
111+
xtemp = xtemp + aa_temp*q[k]
112+
ztemp = ztemp - aa_temp*p[k]
113+
l++
114+
} else {
115+
// ->
116+
bb = rr * (*gha)[l+1]
117+
cc = aa*cl[m] + bb*sl[m]
118+
x = x + cc*q[k]
119+
z = z - cc*p[k]
120+
if clat > 0 {
121+
y = y + (aa*sl[m]-bb*cl[m])*
122+
fm*p[k]/((fn+1.0)*clat)
123+
} else {
124+
y = y + (aa*sl[m]-bb*cl[m])*q[k]*slat
125+
}
126+
// <-
127+
128+
// ->
129+
bb_temp := rr * (*ghb)[l+1]
130+
cc_temp := aa_temp*cl[m] + bb_temp*sl[m]
131+
xtemp = xtemp + cc_temp*q[k]
132+
ztemp = ztemp - cc_temp*p[k]
133+
if clat > 0 {
134+
ytemp = ytemp + (aa_temp*sl[m]-bb_temp*cl[m])*fm*p[k]/((fn+1.0)*clat)
135+
} else {
136+
ytemp = ytemp + (aa_temp*sl[m]-bb_temp*cl[m])*q[k]*slat
137+
}
138+
// <-
139+
l = l + 2
140+
}
141+
m++
142+
}
143+
144+
aa = x
145+
x = x*cd + z*sd
146+
z = z*cd - aa*sd
147+
aa = xtemp
148+
xtemp = xtemp*cd + ztemp*sd
149+
ztemp = ztemp*cd - aa*sd
150+
return x, y, z, xtemp, ytemp, ztemp
151+
}

‎coeffs/read.go‎

Lines changed: 13 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -42,41 +42,43 @@ func NewCoeffsData() (*IGRFcoeffs, error) {
4242
return &igrf, nil
4343
}
4444

45-
// Coeffs - returns two sets of SHC coeffs for the given `date`, as well as for `date` plus one year.
46-
func (igrf *IGRFcoeffs) Coeffs(date float64) (*[]float64, *[]float64, error) {
45+
// Coeffs - returns two sets of SHC coeffs for the given `date`, as well as for `date` plus one year. Also returns the maximal spherical harmonic degree.
46+
func (igrf *IGRFcoeffs) Coeffs(date float64) (*[]float64, *[]float64, int, error) {
4747
max_column := len(*igrf.epochs)
4848
min_epoch := (*igrf.epochs)[0]
4949
max_epoch := (*igrf.epochs)[max_column-1]
5050
if date < min_epoch || date > max_epoch {
51-
return nil, nil, errors.New(fmt.Sprintf("Date %v is out of range (%v, %v).", date, min_epoch, max_epoch))
51+
return nil, nil, 0, errors.New(fmt.Sprintf("Date %v is out of range (%v, %v).", date, min_epoch, max_epoch))
5252
}
5353
// calculate coeffs for the requested date
5454
start, end := igrf.findEpochs(date)
55-
coeffs_start := igrf.interpolateCoeffs(start, end, date)
55+
var nmax int
56+
coeffs_start, nmax := igrf.interpolateCoeffs(start, end, date)
5657
// in order to calculate yaerly SV add 1 year to the date
5758
date = date + 1
5859
var coeffs_end *[]float64
5960
if date < max_epoch {
60-
coeffs_end = igrf.interpolateCoeffs(start, end, date)
61+
coeffs_end, nmax = igrf.interpolateCoeffs(start, end, date)
6162
} else {
6263
coeffs_end = igrf.extrapolateCoeffs(start, end, date)
6364
}
64-
return coeffs_start, coeffs_end, nil
65+
return coeffs_start, coeffs_end, nmax, nil
6566
}
6667

67-
func (igrf *IGRFcoeffs) interpolateCoeffs(start_epoch, end_epoch string, date float64) *[]float64 {
68+
func (igrf *IGRFcoeffs) interpolateCoeffs(start_epoch, end_epoch string, date float64) (*[]float64, int) {
6869
factor := findDateFactor(start_epoch, end_epoch, date)
6970
coeffs_start := (*igrf.data)[start_epoch].coeffs
7071
coeffs_end := (*igrf.data)[end_epoch].coeffs
7172
values := make([]float64, len(*coeffs_start))
7273
nmax1 := (*igrf.data)[start_epoch].nmax
7374
nmax2 := (*igrf.data)[end_epoch].nmax
74-
var k, l int
75+
var k, l, nmax int
7576
var interp func(float64, float64, float64) float64
7677
if nmax1 == nmax2 {
7778
// before 2000.0
7879
k = nmax1 * (nmax1 + 2)
7980
l = -100
81+
nmax = nmax1
8082
} else {
8183
if nmax1 > nmax2 {
8284
// the last column has degree of 8
@@ -86,13 +88,15 @@ func (igrf *IGRFcoeffs) interpolateCoeffs(start_epoch, end_epoch string, date fl
8688
interp = func(start, end, f float64) float64 {
8789
return start
8890
}
91+
nmax = nmax1
8992
} else {
9093
// between 1995.0 and 2000.0
9194
k = nmax1 * (nmax1 + 2)
9295
l = nmax2 * (nmax2 + 2)
9396
interp = func(start, end, f float64) float64 {
9497
return f * end
9598
}
99+
nmax = nmax2
96100
}
97101
}
98102
for i := 0; i < coeffs_lines; i++ {
@@ -106,7 +110,7 @@ func (igrf *IGRFcoeffs) interpolateCoeffs(start_epoch, end_epoch string, date fl
106110
}
107111
values[i] = value
108112
}
109-
return &values
113+
return &values, nmax
110114
}
111115

112116
func (igrf *IGRFcoeffs) extrapolateCoeffs(start_epoch, end_epoch string, date float64) *[]float64 {

0 commit comments

Comments
 (0)