1  #! /usr/bin/env python
 2  # -*- coding: Latin-1 -*-
 3  """ geomag
 4
 5   Conversion de coordonnees geographiques (latitude, longitude) 
 6   en coordonnees geomagnetiques (latitude, longitude). 
 7   Les altitudes ne changent pas.
 8   Les latitudes et longitudes sont exprimees en degres.
 9
10   Exemple :
11    coordonnees geographiques de Halley (-75.5 °N, -26,6 °E)
12  > -65.8 °N, 24.3 °E
13
14   Source : 
15    http://www.nerc-bas.ac.uk/public/uasd/instrums/magnet/gmrot.html
16
17  """
18  #import os 
19  #import sys
20  #import re
21  from math import sin
22  from math import cos
23  from math import acos
24  from math import atan2
25  from math import degrees
26  from math import radians
27  from math import pi
28
29  def usage():
30      print __doc__
31
32  class Point:
33    def __init__(self,lat,lon):
34      self.lat=radians(lat)
35      self.lon=radians(lon)
36    def afficher(self):
37      print "lat=",degrees(self.lat)
38      print "long=",degrees(self.lon)
39
40  def geomag(p1):
41     alpha=11.5
42     beta=-69.0
43     beta=radians(beta)
44     alpha=radians(alpha)
45     # colatitude geographique
46     thetag=pi/2.0-p1.lat
47     # longitude geographique
48     phig=p1.lon
49     #longitude magnetique
50     phim=atan2((cos(beta)*sin(thetag)*sin(phig)-sin(beta)*sin(thetag)*cos(phig)),(cos(alpha)*cos(beta)*sin(thetag)*cos(phig)+cos(alpha)*sin(beta)*sin(thetag)*sin(phig)-sin(alpha)*cos(thetag)))
51     # colatitude magnetique
52     thetam=acos(sin(alpha)*cos(beta)*sin(thetag)*cos(phig)+sin(alpha)*sin(beta)*sin(thetag)*sin(phig)+cos(alpha)*cos(thetag))
53     # latitude magnetique
54     thetam=degrees(pi/2-thetam)
55     phim= degrees(phim)
56     p2=Point(thetam,phim)
57     return p2
58
59  p1=Point(-75.5, -26.6)
60  print "Coordonnees geographiques : "
61  p1.afficher()
62  print "Coordonnees geomagnetiques : "
63  geomag(p1).afficher()