geomag.py
Go to the documentation of this file.
1 # geomag.py
2 # by Christopher Weiss cmweiss@gmail.com
3 # From https://code.google.com/p/geomag/
4 # License: LGPL
5 
6 # Adapted from the geomagc software and World Magnetic Model of the NOAA
7 # Satellite and Information Service, National Geophysical Data Center
8 # http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml
9 #
10 # Suggestions for improvements are appreciated.
11 
12 # USAGE:
13 #
14 # >>> gm = geomag.GeoMag("WMM.COF")
15 # >>> mag = gm.GeoMag(80,0)
16 # >>> mag.dec
17 # -6.1335150785195536
18 # >>>
19 
20 import math, os, unittest
21 from datetime import date
22 
23 class GeoMag:
24 
25  def calc(self, dlat, dlon, h=0, time=date.today()): # latitude (decimal degrees), longitude (decimal degrees), altitude (feet), date
26  #time = date('Y') + date('z')/365
27  time = time.year+((time - date(time.year,1,1)).days/365.0)
28  alt = h/3280.8399
29 
30  otime = oalt = olat = olon = -1000.0
31 
32  dt = time - self.epoch
33  glat = dlat
34  glon = dlon
35  rlat = math.radians(glat)
36  rlon = math.radians(glon)
37  srlon = math.sin(rlon)
38  srlat = math.sin(rlat)
39  crlon = math.cos(rlon)
40  crlat = math.cos(rlat)
41  srlat2 = srlat*srlat
42  crlat2 = crlat*crlat
43  self.sp[1] = srlon
44  self.cp[1] = crlon
45 
46  #/* CONVERT FROM GEODETIC COORDS. TO SPHERICAL COORDS. */
47  if (alt != oalt or glat != olat):
48  q = math.sqrt(self.a2-self.c2*srlat2)
49  q1 = alt*q
50  q2 = ((q1+self.a2)/(q1+self.b2))*((q1+self.a2)/(q1+self.b2))
51  ct = srlat/math.sqrt(q2*crlat2+srlat2)
52  st = math.sqrt(1.0-(ct*ct))
53  r2 = (alt*alt)+2.0*q1+(self.a4-self.c4*srlat2)/(q*q)
54  r = math.sqrt(r2)
55  d = math.sqrt(self.a2*crlat2+self.b2*srlat2)
56  ca = (alt+d)/r
57  sa = self.c2*crlat*srlat/(r*d)
58 
59  if (glon != olon):
60  for m in range(2,self.maxord+1):
61  self.sp[m] = self.sp[1]*self.cp[m-1]+self.cp[1]*self.sp[m-1]
62  self.cp[m] = self.cp[1]*self.cp[m-1]-self.sp[1]*self.sp[m-1]
63 
64  aor = self.re/r
65  ar = aor*aor
66  br = bt = bp = bpp = 0.0
67  for n in range(1,self.maxord+1):
68  ar = ar*aor
69 
70  #for (m=0,D3=1,D4=(n+m+D3)/D3;D4>0;D4--,m+=D3):
71  m=0
72  D3=1
73  #D4=(n+m+D3)/D3
74  D4=(n+m+1)
75  while D4>0:
76 
77  # /*
78  # COMPUTE UNNORMALIZED ASSOCIATED LEGENDRE POLYNOMIALS
79  # AND DERIVATIVES VIA RECURSION RELATIONS
80  # */
81  if (alt != oalt or glat != olat):
82  if (n == m):
83  self.p[m][n] = st * self.p[m-1][n-1]
84  self.dp[m][n] = st*self.dp[m-1][n-1]+ct*self.p[m-1][n-1]
85 
86  elif (n == 1 and m == 0):
87  self.p[m][n] = ct*self.p[m][n-1]
88  self.dp[m][n] = ct*self.dp[m][n-1]-st*self.p[m][n-1]
89 
90  elif (n > 1 and n != m):
91  if (m > n-2):
92  self.p[m][n-2] = 0
93  if (m > n-2):
94  self.dp[m][n-2] = 0.0
95  self.p[m][n] = ct*self.p[m][n-1]-self.k[m][n]*self.p[m][n-2]
96  self.dp[m][n] = ct*self.dp[m][n-1] - st*self.p[m][n-1]-self.k[m][n]*self.dp[m][n-2]
97 
98  # /*
99  # TIME ADJUST THE GAUSS COEFFICIENTS
100  # */
101  if (time != otime):
102  self.tc[m][n] = self.c[m][n]+dt*self.cd[m][n]
103  if (m != 0):
104  self.tc[n][m-1] = self.c[n][m-1]+dt*self.cd[n][m-1]
105 
106  # /*
107  # ACCUMULATE TERMS OF THE SPHERICAL HARMONIC EXPANSIONS
108  # */
109  par = ar*self.p[m][n]
110 
111  if (m == 0):
112  temp1 = self.tc[m][n]*self.cp[m]
113  temp2 = self.tc[m][n]*self.sp[m]
114  else:
115  temp1 = self.tc[m][n]*self.cp[m]+self.tc[n][m-1]*self.sp[m]
116  temp2 = self.tc[m][n]*self.sp[m]-self.tc[n][m-1]*self.cp[m]
117 
118  bt = bt-ar*temp1*self.dp[m][n]
119  bp = bp + (self.fm[m] * temp2 * par)
120  br = br + (self.fn[n] * temp1 * par)
121  # /*
122  # SPECIAL CASE: NORTH/SOUTH GEOGRAPHIC POLES
123  # */
124  if (st == 0.0 and m == 1):
125  if (n == 1):
126  self.pp[n] = self.pp[n-1]
127  else:
128  self.pp[n] = ct*self.pp[n-1]-self.k[m][n]*self.pp[n-2]
129  parp = ar*self.pp[n]
130  bpp = bpp + (self.fm[m]*temp2*parp)
131 
132  D4=D4-1
133  m=m+1
134 
135  if (st == 0.0):
136  bp = bpp
137  else:
138  bp = bp/st
139  # /*
140  # ROTATE MAGNETIC VECTOR COMPONENTS FROM SPHERICAL TO
141  # GEODETIC COORDINATES
142  # */
143  bx = -bt*ca-br*sa
144  by = bp
145  bz = bt*sa-br*ca
146  # /*
147  # COMPUTE DECLINATION (DEC), INCLINATION (DIP) AND
148  # TOTAL INTENSITY (TI)
149  # */
150  bh = math.sqrt((bx*bx)+(by*by))
151  ti = math.sqrt((bh*bh)+(bz*bz))
152  dec = math.degrees(math.atan2(by,bx))
153  dip = math.degrees(math.atan2(bz,bh))
154  # /*
155  # COMPUTE MAGNETIC GRID VARIATION IF THE CURRENT
156  # GEODETIC POSITION IS IN THE ARCTIC OR ANTARCTIC
157  # (I.E. GLAT > +55 DEGREES OR GLAT < -55 DEGREES)
158 
159  # OTHERWISE, SET MAGNETIC GRID VARIATION TO -999.0
160  # */
161  gv = -999.0
162  if (math.fabs(glat) >= 55.):
163  if (glat > 0.0 and glon >= 0.0):
164  gv = dec-glon
165  if (glat > 0.0 and glon < 0.0):
166  gv = dec+math.fabs(glon);
167  if (glat < 0.0 and glon >= 0.0):
168  gv = dec+glon
169  if (glat < 0.0 and glon < 0.0):
170  gv = dec-math.fabs(glon)
171  if (gv > +180.0):
172  gv = gv - 360.0
173  if (gv < -180.0):
174  gv = gv + 360.0
175 
176  otime = time
177  oalt = alt
178  olat = glat
179  olon = glon
180 
181  class RetObj:
182  pass
183  retobj = RetObj()
184  retobj.dec = dec
185  retobj.dip = dip
186  retobj.ti = ti
187  retobj.bh = bh
188  retobj.bx = bx
189  retobj.by = by
190  retobj.bz = bz
191  retobj.lat = dlat
192  retobj.lon = dlon
193  retobj.alt = h
194  retobj.time = time
195 
196  return retobj
197 
198  def __init__(self, wmm_filename=None):
199  if not wmm_filename:
200  wmm_filename = os.path.join(os.path.dirname(__file__), 'WMM.COF')
201  wmm=[]
202  with open(wmm_filename) as wmm_file:
203  for line in wmm_file:
204  linevals = line.strip().split()
205  if len(linevals) == 3:
206  self.epoch = float(linevals[0])
207  self.model = linevals[1]
208  self.modeldate = linevals[2]
209  elif len(linevals) == 6:
210  linedict = {'n': int(float(linevals[0])),
211  'm': int(float(linevals[1])),
212  'gnm': float(linevals[2]),
213  'hnm': float(linevals[3]),
214  'dgnm': float(linevals[4]),
215  'dhnm': float(linevals[5])}
216  wmm.append(linedict)
217 
218  z = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
219  self.maxord = self.maxdeg = 12
220  self.tc = [z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13]]
221  self.sp = z[0:14]
222  self.cp = z[0:14]
223  self.cp[0] = 1.0
224  self.pp = z[0:13]
225  self.pp[0] = 1.0
226  self.p = [z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14]]
227  self.p[0][0] = 1.0
228  self.dp = [z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13]]
229  self.a = 6378.137
230  self.b = 6356.7523142
231  self.re = 6371.2
232  self.a2 = self.a*self.a
233  self.b2 = self.b*self.b
234  self.c2 = self.a2-self.b2
235  self.a4 = self.a2*self.a2
236  self.b4 = self.b2*self.b2
237  self.c4 = self.a4 - self.b4
238 
239  self.c = [z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14]]
240  self.cd = [z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14],z[0:14]]
241 
242  for wmmnm in wmm:
243  m = wmmnm['m']
244  n = wmmnm['n']
245  gnm = wmmnm['gnm']
246  hnm = wmmnm['hnm']
247  dgnm = wmmnm['dgnm']
248  dhnm = wmmnm['dhnm']
249  if (m <= n):
250  self.c[m][n] = gnm
251  self.cd[m][n] = dgnm
252  if (m != 0):
253  self.c[n][m-1] = hnm
254  self.cd[n][m-1] = dhnm
255 
256  #/* CONVERT SCHMIDT NORMALIZED GAUSS COEFFICIENTS TO UNNORMALIZED */
257  self.snorm = [z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13]]
258  self.snorm[0][0] = 1.0
259  self.k = [z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13],z[0:13]]
260  self.k[1][1] = 0.0
261  self.fn = [0.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0]
262  self.fm = [0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0]
263  for n in range(1,self.maxord+1):
264  self.snorm[0][n] = self.snorm[0][n-1]*(2.0*n-1)/n
265  j=2.0
266  #for (m=0,D1=1,D2=(n-m+D1)/D1;D2>0;D2--,m+=D1):
267  m=0
268  D1=1
269  D2=(n-m+D1)/D1
270  while (D2 > 0):
271  self.k[m][n] = (((n-1)*(n-1))-(m*m))/((2.0*n-1)*(2.0*n-3.0))
272  if (m > 0):
273  flnmj = ((n-m+1.0)*j)/(n+m)
274  self.snorm[m][n] = self.snorm[m-1][n]*math.sqrt(flnmj)
275  j = 1.0
276  self.c[n][m-1] = self.snorm[m][n]*self.c[n][m-1]
277  self.cd[n][m-1] = self.snorm[m][n]*self.cd[n][m-1]
278  self.c[m][n] = self.snorm[m][n]*self.c[m][n]
279  self.cd[m][n] = self.snorm[m][n]*self.cd[m][n]
280  D2=D2-1
281  m=m+D1
282 
283 class GeoMagTest(unittest.TestCase):
284 
285  d1=date(2010,1,1)
286  d2=date(2012,7,1)
287 
288  test_values = (
289  # date, alt, lat, lon, var
290  (d1, 0, 80, 0, -6.13),
291  (d1, 0, 0, 120, 0.97),
292  (d1, 0, -80, 240, 70.21),
293  (d1, 328083.99, 80, 0, -6.57),
294  (d1, 328083.99, 0, 120, 0.94),
295  (d1, 328083.99, -80, 240, 69.62),
296  (d2, 0, 80, 0, -5.21),
297  (d2, 0, 0, 120, 0.88),
298  (d2, 0, -80, 240, 70.04),
299  (d2, 328083.99, 80, 0, -5.63),
300  (d2, 328083.99, 0, 120, 0.86),
301  (d2, 328083.99, -80, 240, 69.45),
302  )
303 
304  def test_declination(self):
305  gm = GeoMag()
306  for values in self.test_values:
307  calcval=gm.GeoMag(values[2], values[3], values[1], values[0])
308  self.assertAlmostEqual(values[4], calcval.dec, 2, 'Expected %s, result %s' % (values[4], calcval.dec))
309 
310 if __name__ == '__main__':
311  unittest.main()
geomag.geomag.GeoMag.__init__
def __init__(self, wmm_filename=None)
Definition: geomag.py:198
geomag.geomag.GeoMagTest.test_declination
def test_declination(self)
Definition: geomag.py:304
geomag.geomag.GeoMag.calc
def calc(self, dlat, dlon, h=0, time=date.today())
Definition: geomag.py:25

declination
Author(s):
autogenerated on Sat May 23 2020 03:41:45