1 from scipy import *
   2 from numpy import *
   3 from pylab import *
   4 import Image
   5 
   6 
   7 def mapxy(x, y, sgn):
   8     cdr  = 57.29577951
   9     slat = 70.0
  10     re   = 6378.273
  11     ec2   = .006693883
  12     ec    =  sqrt(ec2)
  13     if sgn == 1:
  14         delta = 45.0
  15     else:
  16         delta = 0.0
  17     sl = slat * pi/180.0
  18     rho = sqrt(x**2 + y**2)
  19     cm = cos(sl) / sqrt(1.0 - ec2 * (sin(sl)**2))
  20     t = tan((pi / 4.0) - (sl / 2.0)) / ((1.0 - ec * sin(sl)) / (1.0 + ec * sin(sl)))**(ec / 2.0)
  21     if  (absolute(slat-90.0) < 1.e-5):
  22         t = rho * sqrt((1. + ec)**(1. + ec) * (1. - ec)**(1. - ec)) / 2. / re
  23     else:
  24         t = rho * t / (re * cm)
  25     chi = (pi / 2.0) - 2.0 * arctan(t)
  26     alat = chi + ((ec2 / 2.0) + (5.0 * ec2**2.0 / 24.0) + (ec2**3.0 / 12.0)) * sin(2 * chi) + ((7.0 * ec2**2.0 / 48.0) + (29.0 * ec2**3 / 240.0)) *sin(4.0 * chi)+ (7.0 * ec2**3.0 / 120.0) * sin(6.0 * chi)
  27     alat = sgn * alat
  28     along = arctan2(sgn * x, -sgn * y)
  29     along = sgn * along
  30 
  31     along = along * 180. / pi
  32     alat  = alat * 180. / pi
  33     along = along - delta
  34     return [alat,along]
  35 
  36 def mapll(lat, lon,sgn):
  37     cdr  = 57.29577951
  38     slat = 70.
  39     re   = 6378.273
  40     ec2   = .006693883
  41     ec    =  sqrt(ec2)
  42     if sgn == 1:
  43         delta = 45.0
  44     else:
  45         delta = 0.0
  46     latitude  = absolute(lat) * pi/180.
  47     longitude = (lon + delta) * pi/180.
  48     t =   tan(pi/4-latitude/2)/((1-ec*sin(latitude))/(1+ec*sin(latitude)))**(ec/2)
  49     if ((90-slat) < 1.e-5):
  50         rho = 2.*re*t/sqrt((1.+ec)**(1.+ec)*(1.-ec)**(1.-ec))
  51     else:
  52         sl  = slat*pi/180.
  53         tc  = tan(pi/4.-sl/2.)/((1.-ec*sin(sl))/(1.+ec*sin(sl)))**(ec/2.)
  54         mc  = cos(sl)/sqrt(1.0-ec2*(sin(sl)**2))
  55         rho = re*mc*t/tc
  56 
  57     y=-rho*sgn*cos(sgn*longitude)
  58     x= rho*sgn*sin(sgn*longitude)
  59     return [x,y]
  60 
  61 
  62 lat1,lon1=70,10 
  63 lat2,lon2=75,15
  64 xpc1,ypc1=mapll(lat1,lon1,1) 
  65 xpc2,ypc2=mapll(lat2,lon2,1)
  66 
  67 #print xpc1,ypc1,xpc2,ypc2
  68 
  69 shiftx=(xpc1-xpc2) #in [km]
  70 shifty=(ypc1-ypc2) #in [km]
  71 
  72 print shiftx,shifty