Tuesday, July 8, 2008

World Magnetic Model

In about a dozen or so math courses I remember having to prove (always by induction) that the sum of N numbers (1+2+3+...+N) is equal to n(n+1)/2 and always thinking, this is pointless. I stand corrected.

At work I was tasked to calculate magnetic declination using java instead of the C implementation NOAA used. The program is pretty much a function: MagDec:(double,double,double)->double (elevation,latitude,longitude)->magnetic_dec. A key component used to calculate the declination is a file with spherical harmonics gathered by satellites. It looks something like this:

          G(m,n)   H(m,n)       G_dot(m,n)  H_dot(m,n)
1 0 -29556.8 0.0 8.0 0.0
1 1 -1671.7 5079.8 10.6 -20.9
2 0 -2340.6 0.0 -15.1 0.0
2 1 3046.9 -2594.7 -7.8 -23.2
2 2 1657.0 -516.7 -0.8 -14.6
3 0 1335.4 0.0 0.4 0.0
3 1 -2305.1 -199.9 -2.6 5.0
3 2 1246.7 269.3 -1.2 -7.0
3 3 674.0 -524.2 -6.5 -0.6
4 0 919.8 0.0 -2.5 0.0
4 1 798.1 281.5 2.8 2.2
4 2 211.3 -226.0 -7.0 1.6
....
....
12 10 -0.1 -0.9 0.0 0.0
12 11 -0.3 -0.4 0.0 0.0
12 12 -0.1 0.8 0.0 0.0


My goal was to create S = G(m,n) - G_dot(m,n) and C = H(m,n) - H_dot(m,n) for each value of m and n.
Since Java IO isn't the most convenient thing in the world, instead of sscanf() I was forced to use java.util.Scanner, which IMO is not nearly as convenient for this type of data organization.
I realized that if I have a flat list it looks something like this:

Col. A Col. B
G(1,0) -29556.8
8.0
G(1,1) -1671.7
10.6
G(2,0) -2340.6
-15.1
G(2,1) 3046.9
-7.8
G(2,2) 1657.0
-0.8
Written in a different way:
G(1) G(2) G(3) ... G(12)
2 3 4 ... 13

which means there are 13(13+1)/2 - 1 coefficients for G. If I am presented with a single column of numbers like Col B. and I want to pick G(2,1) all I need to do is is realize that there are 3(3+1)/2 coefficients up to G(2,2) which means G(2,1) is 3(3+1)/2-1. Of course it would've been a lot easier if I could do the straightforward thing easily like:

WMM = open("WMM.COF")
coeffs = []
try:
for line in WMM:
nline = []
for each in line.split():
nline.append(float(each))
coeffs.append(nline)
finally:
WMM.close()
but it's not to be.

No comments: