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: