1  #! /usr/bin/env python
  2  # -*- coding: Latin-1 -*-
  3
  4  __author__ = "Jean-Philippe Durand (jpdurand@free.fr)"
  5  __version__ = "$Revision: 1.0 $"
  6  __date__ = "$Date: jeu mai  4 14:01:02 CEST 2006 $"
  7  __copyright__ = "Copyright (c) 2006 J.-Ph. Durand"
  8  __license__ = "Python"
  9
 10  import os
 11  import sys
 12  #import math
 13  from math import cos
 14  from math import sin
 15  from math import acos
 16  from math import asin
 17  from math import atan2
 18  from math import degrees
 19  from math import radians
 20  from math import sqrt
 21
 22  rt = 6371  # rayon terrestre moyen en km
 23
 24  class Point:
 25    def __init__(self,lat,lon):
 26      self.lat=radians(lat)
 27      self.lon=radians(lon)
 28    def afficher(self):
 29      print "lat=",degrees(self.lat)
 30      print "long=",degrees(self.lon)
 31
 32  # 
 33  #  Calcul de la distance (en km) entre 2 points specifies par leurs
 34  #  latitude/longitude avec la formule de Haversine
 35  # 
 36  #   de : Haversine formula - R. W. Sinnott, "Virtues of the Haversine",
 37  #        Sky and Telescope, vol 68, no 2, 1984
 38  #        http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1
 39  #  
 40
 41
 42
 43  def distHaversine(p1, p2):
 44    dLat  = p2.lat - p1.lat
 45    dLong = p2.lon - p1.lon
 46    a = sin(dLat/2) * sin(dLat/2) + cos(p1.lat) * cos(p2.lat) * sin(dLong/2) * sin(dLong/2)
 47    c = 2 * atan2(sqrt(a), sqrt(1-a))
 48    d = rt * c
 49    return  d
 50
 51  # 
 52  # Calcul de la distance (en km) entre 2 points specifies par leurs
 53  # latitude/longitude en utilisant les fonctions trigonometriques
 54  #
 55  def distCosineLaw(p1, p2):
 56    d = acos(sin(p1.lat)*sin(p2.lat) + cos(p1.lat)*cos(p2.lat)*cos(p2.lon-p1.lon)) * rt
 57    return d
 58
 59
 60  # 
 61  #  calcul du cap (en degres ) entre 2 points
 62  # 
 63  #  de : Ed Williams' Aviation Formulary, http://williams.best.vwh.net/avform.htm#Crs
 64  #  
 65  def cap(p1, p2):
 66     y = sin(p2.lon-p1.lon) * cos(p2.lat)
 67     x = cos(p1.lat)*sin(p2.lat) - sin(p1.lat)*cos(p2.lat)*cos(p2.lon-p1.lon)
 68     return degrees(atan2(y, x))
 69
 70  #
 71  # calcul le point de destination a parit du point de depart, du cap et de la distance
 72  #  voir http://williams.best.vwh.net/avform.htm#LL
 73  #
 74  def destPoint(p1, capinit, dist):
 75     if (dist !=0):
 76       p2=Point(0,0)
 77       d = float(dist)/rt  # d = distance angulaire parcourue sur la surface terrestre
 78       capinit = radians(capinit)
 79       p2.lat = asin( sin(p1.lat)*cos(d) + cos(p1.lat)*sin(d)*cos(capinit) )
 80       p2.lon = p1.lon + atan2(sin(capinit)*sin(d)*cos(p1.lat), cos(d)-sin(p1.lat)*sin(p2.lat))
 81     else:
 82             p2=p1
 83     return p2
 84
 85
 86
 87  p1=Point(30,60)
 88  p2=Point(38,61)
 89  p1.afficher()
 90  p2.afficher()
 91  dist = distHaversine(p1, p2)
 92  print dist
 93  dist = distCosineLaw(p1, p2)
 94  print dist
 95  capinit=cap(p1,p2)
 96  print capinit
 97  p3=destPoint(p1, 35.2643896828, 6671.69559867)
 98  p3.afficher()
 99  print degrees(p3.lat)
100  print degrees(p3.lon)