declination: geomag.py Source File

Go to the documentation of this file.
00001 # geomag.py
00002 # by Christopher Weiss cmweiss@gmail.com
00003 # From https://code.google.com/p/geomag/
00004 # License: LGPL
00005 
00006 # Adapted from the geomagc software and World Magnetic Model of the NOAA
00007 # Satellite and Information Service, National Geophysical Data Center
00008 # http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml
00009 #
00010 # Suggestions for improvements are appreciated.
00011 
00012 # USAGE:
00013 #
00014 # >>> gm = geomag.GeoMag("WMM.COF")
00015 # >>> mag = gm.GeoMag(80,0)
00016 # >>> mag.dec
00017 # -6.1335150785195536
00018 # >>>
00019 
00020 import math, os, unittest
00021 from datetime import date
00022 
00023 class GeoMag:
00024 
00025     def calc(self, dlat, dlon, h=0, time=date.today()): # latitude (decimal degrees), longitude (decimal degrees), altitude (feet), date
00026         #time = date('Y') + date('z')/365
00027         time = time.year+((time - date(time.year,1,1)).days/365.0)
00028         alt = h/3280.8399
00029 
00030         otime = oalt = olat = olon = -1000.0
00031 
00032         dt = time - self.epoch
00033         glat = dlat
00034         glon = dlon
00035         rlat = math.radians(glat)
00036         rlon = math.radians(glon)
00037         srlon = math.sin(rlon)
00038         srlat = math.sin(rlat)
00039         crlon = math.cos(rlon)
00040         crlat = math.cos(rlat)
00041         srlat2 = srlat*srlat
00042         crlat2 = crlat*crlat
00043         self.sp[1] = srlon
00044         self.cp[1] = crlon
00045 
00046         #/* CONVERT FROM GEODETIC COORDS. TO SPHERICAL COORDS. */
00047         if (alt != oalt or glat != olat):
00048             q = math.sqrt(self.a2-self.c2*srlat2)
00049             q1 = alt*q
00050             q2 = ((q1+self.a2)/(q1+self.b2))*((q1+self.a2)/(q1+self.b2))
00051             ct = srlat/math.sqrt(q2*crlat2+srlat2)
00052             st = math.sqrt(1.0-(ct*ct))
00053             r2 = (alt*alt)+2.0*q1+(self.a4-self.c4*srlat2)/(q*q)
00054             r = math.sqrt(r2)
00055             d = math.sqrt(self.a2*crlat2+self.b2*srlat2)
00056             ca = (alt+d)/r
00057             sa = self.c2*crlat*srlat/(r*d)
00058 
00059         if (glon != olon):
00060             for m in range(2,self.maxord+1):
00061                 self.sp[m] = self.sp[1]*self.cp[m-1]+self.cp[1]*self.sp[m-1]
00062                 self.cp[m] = self.cp[1]*self.cp[m-1]-self.sp[1]*self.sp[m-1]
00063 
00064         aor = self.re/r
00065         ar = aor*aor
00066         br = bt = bp = bpp = 0.0
00067         for n in range(1,self.maxord+1):
00068             ar = ar*aor
00069             
00070             #for (m=0,D3=1,D4=(n+m+D3)/D3;D4>0;D4--,m+=D3):
00071             m=0
00072             D3=1
00073             #D4=(n+m+D3)/D3
00074             D4=(n+m+1)
00075             while D4>0:
00076 
00077         # /*
00078                 # COMPUTE UNNORMALIZED ASSOCIATED LEGENDRE POLYNOMIALS
00079                 # AND DERIVATIVES VIA RECURSION RELATIONS
00080         # */
00081                 if (alt != oalt or glat != olat):
00082                     if (n == m):
00083                         self.p[m][n] = st * self.p[m-1][n-1]
00084                         self.dp[m][n] = st*self.dp[m-1][n-1]+ct*self.p[m-1][n-1]
00085 
00086                     elif (n == 1 and m == 0):
00087                         self.p[m][n] = ct*self.p[m][n-1]
00088                         self.dp[m][n] = ct*self.dp[m][n-1]-st*self.p[m][n-1]
00089 
00090                     elif (n > 1 and n != m):
00091                         if (m > n-2):
00092                             self.p[m][n-2] = 0
00093                         if (m > n-2):
00094                             self.dp[m][n-2] = 0.0
00095                         self.p[m][n] = ct*self.p[m][n-1]-self.k[m][n]*self.p[m][n-2]
00096                         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]
00097 
00098         # /*
00099                 # TIME ADJUST THE GAUSS COEFFICIENTS
00100         # */
00101                 if (time != otime):
00102                     self.tc[m][n] = self.c[m][n]+dt*self.cd[m][n]
00103                     if (m != 0):
00104                         self.tc[n][m-1] = self.c[n][m-1]+dt*self.cd[n][m-1]
00105 
00106         # /*
00107                 # ACCUMULATE TERMS OF THE SPHERICAL HARMONIC EXPANSIONS
00108         # */
00109                 par = ar*self.p[m][n]
00110                 
00111                 if (m == 0):
00112                     temp1 = self.tc[m][n]*self.cp[m]
00113                     temp2 = self.tc[m][n]*self.sp[m]
00114                 else:
00115                     temp1 = self.tc[m][n]*self.cp[m]+self.tc[n][m-1]*self.sp[m]
00116                     temp2 = self.tc[m][n]*self.sp[m]-self.tc[n][m-1]*self.cp[m]
00117 
00118                 bt = bt-ar*temp1*self.dp[m][n]
00119                 bp = bp + (self.fm[m] * temp2 * par)
00120                 br = br + (self.fn[n] * temp1 * par)
00121         # /*
00122                     # SPECIAL CASE:  NORTH/SOUTH GEOGRAPHIC POLES
00123         # */
00124                 if (st == 0.0 and m == 1):
00125                     if (n == 1):
00126                         self.pp[n] = self.pp[n-1]
00127                     else:
00128                         self.pp[n] = ct*self.pp[n-1]-self.k[m][n]*self.pp[n-2]
00129                     parp = ar*self.pp[n]
00130                     bpp = bpp + (self.fm[m]*temp2*parp)
00131                     
00132                 D4=D4-1
00133                 m=m+1
00134 
00135         if (st == 0.0):
00136             bp = bpp
00137         else:
00138             bp = bp/st
00139         # /*
00140             # ROTATE MAGNETIC VECTOR COMPONENTS FROM SPHERICAL TO
00141             # GEODETIC COORDINATES
00142         # */
00143         bx = -bt*ca-br*sa
00144         by = bp
00145         bz = bt*sa-br*ca
00146         # /*
00147             # COMPUTE DECLINATION (DEC), INCLINATION (DIP) AND
00148             # TOTAL INTENSITY (TI)
00149         # */
00150         bh = math.sqrt((bx*bx)+(by*by))
00151         ti = math.sqrt((bh*bh)+(bz*bz))
00152         dec = math.degrees(math.atan2(by,bx))
00153         dip = math.degrees(math.atan2(bz,bh))
00154         # /*
00155             # COMPUTE MAGNETIC GRID VARIATION IF THE CURRENT
00156             # GEODETIC POSITION IS IN THE ARCTIC OR ANTARCTIC
00157             # (I.E. GLAT > +55 DEGREES OR GLAT < -55 DEGREES)
00158 
00159             # OTHERWISE, SET MAGNETIC GRID VARIATION TO -999.0
00160         # */
00161         gv = -999.0
00162         if (math.fabs(glat) >= 55.):
00163             if (glat > 0.0 and glon >= 0.0):
00164                 gv = dec-glon
00165             if (glat > 0.0 and glon < 0.0):
00166                 gv = dec+math.fabs(glon);
00167             if (glat < 0.0 and glon >= 0.0):
00168                 gv = dec+glon
00169             if (glat < 0.0 and glon < 0.0):
00170                 gv = dec-math.fabs(glon)
00171             if (gv > +180.0):
00172                 gv = gv - 360.0
00173             if (gv < -180.0):
00174                 gv = gv + 360.0
00175 
00176         otime = time
00177         oalt = alt
00178         olat = glat
00179         olon = glon
00180 
00181         class RetObj:
00182             pass
00183         retobj = RetObj()
00184         retobj.dec = dec
00185         retobj.dip = dip
00186         retobj.ti = ti
00187         retobj.bh = bh
00188         retobj.bx = bx
00189         retobj.by = by
00190         retobj.bz = bz
00191         retobj.lat = dlat
00192         retobj.lon = dlon
00193         retobj.alt = h
00194         retobj.time = time
00195 
00196         return retobj
00197 
00198     def __init__(self, wmm_filename=None):
00199         if not wmm_filename:
00200             wmm_filename = os.path.join(os.path.dirname(__file__), 'WMM.COF')
00201         wmm=[]
00202         with open(wmm_filename) as wmm_file:
00203             for line in wmm_file:
00204                 linevals = line.strip().split()
00205                 if len(linevals) == 3:
00206                     self.epoch = float(linevals[0])
00207                     self.model = linevals[1]
00208                     self.modeldate = linevals[2]
00209                 elif len(linevals) == 6:
00210                     linedict = {'n': int(float(linevals[0])),
00211                     'm': int(float(linevals[1])),
00212                     'gnm': float(linevals[2]),
00213                     'hnm': float(linevals[3]),
00214                     'dgnm': float(linevals[4]),
00215                     'dhnm': float(linevals[5])}
00216                     wmm.append(linedict)
00217 
00218         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]
00219         self.maxord = self.maxdeg = 12
00220         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]]
00221         self.sp = z[0:14]
00222         self.cp = z[0:14]
00223         self.cp[0] = 1.0
00224         self.pp = z[0:13]
00225         self.pp[0] = 1.0
00226         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]]
00227         self.p[0][0] = 1.0
00228         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]]
00229         self.a = 6378.137
00230         self.b = 6356.7523142
00231         self.re = 6371.2
00232         self.a2 = self.a*self.a
00233         self.b2 = self.b*self.b
00234         self.c2 = self.a2-self.b2
00235         self.a4 = self.a2*self.a2
00236         self.b4 = self.b2*self.b2
00237         self.c4 = self.a4 - self.b4
00238 
00239         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]]
00240         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]]
00241         
00242         for wmmnm in wmm:
00243             m = wmmnm['m']
00244             n = wmmnm['n']
00245             gnm = wmmnm['gnm']
00246             hnm = wmmnm['hnm']
00247             dgnm = wmmnm['dgnm']
00248             dhnm = wmmnm['dhnm']
00249             if (m <= n):
00250                 self.c[m][n] = gnm
00251                 self.cd[m][n] = dgnm
00252                 if (m != 0):
00253                     self.c[n][m-1] = hnm
00254                     self.cd[n][m-1] = dhnm
00255 
00256         #/* CONVERT SCHMIDT NORMALIZED GAUSS COEFFICIENTS TO UNNORMALIZED */
00257         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]]
00258         self.snorm[0][0] = 1.0
00259         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]]
00260         self.k[1][1] = 0.0
00261         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]
00262         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]
00263         for n in range(1,self.maxord+1):
00264             self.snorm[0][n] = self.snorm[0][n-1]*(2.0*n-1)/n
00265             j=2.0
00266             #for (m=0,D1=1,D2=(n-m+D1)/D1;D2>0;D2--,m+=D1):
00267             m=0
00268             D1=1
00269             D2=(n-m+D1)/D1
00270             while (D2 > 0):
00271                 self.k[m][n] = (((n-1)*(n-1))-(m*m))/((2.0*n-1)*(2.0*n-3.0))
00272                 if (m > 0):
00273                     flnmj = ((n-m+1.0)*j)/(n+m)
00274                     self.snorm[m][n] = self.snorm[m-1][n]*math.sqrt(flnmj)
00275                     j = 1.0
00276                     self.c[n][m-1] = self.snorm[m][n]*self.c[n][m-1]
00277                     self.cd[n][m-1] = self.snorm[m][n]*self.cd[n][m-1]
00278                 self.c[m][n] = self.snorm[m][n]*self.c[m][n]
00279                 self.cd[m][n] = self.snorm[m][n]*self.cd[m][n]
00280                 D2=D2-1
00281                 m=m+D1
00282 
00283 class GeoMagTest(unittest.TestCase):
00284 
00285     d1=date(2010,1,1)
00286     d2=date(2012,7,1)
00287     
00288     test_values = (
00289         # date, alt, lat, lon, var
00290         (d1, 0, 80, 0, -6.13),
00291         (d1, 0, 0, 120, 0.97),
00292         (d1, 0, -80, 240, 70.21),
00293         (d1, 328083.99, 80, 0, -6.57),
00294         (d1, 328083.99, 0, 120, 0.94),
00295         (d1, 328083.99, -80, 240, 69.62),
00296         (d2, 0, 80, 0, -5.21),
00297         (d2, 0, 0, 120, 0.88),
00298         (d2, 0, -80, 240, 70.04),
00299         (d2, 328083.99, 80, 0, -5.63),
00300         (d2, 328083.99, 0, 120, 0.86),
00301         (d2, 328083.99, -80, 240, 69.45),
00302     )
00303     
00304     def test_declination(self):
00305         gm = GeoMag()
00306         for values in self.test_values:
00307             calcval=gm.GeoMag(values[2], values[3], values[1], values[0])
00308             self.assertAlmostEqual(values[4], calcval.dec, 2, 'Expected %s, result %s' % (values[4], calcval.dec))
00309 
00310 if __name__ == '__main__':
00311     unittest.main()