Правильны ли мои уравнения?Поверните на сфере от широты / долготы от A до B, где будет точка C? - PullRequest
Новая
       44

Правильны ли мои уравнения?Поверните на сфере от широты / долготы от A до B, где будет точка C?

2 голосов
/ 30 марта 2012

Я написал приведенный ниже скрипт на Python.Идея состоит в том, чтобы вычислить новое местоположение точки C после вращения земного шара из точки A в точку B. Сначала я вычисляю точку P, которая является полюсом вращения.С расчетной точкой P уже что-то не так.Со следующим вводом fe я бы предположил, что точка P будет иметь широту 90 или –90.

Я задавал этот вопрос раньше здесь: Поворот сферы от скоординированной1 до координат2, где будет координата3? Но я решил, что лучше спросить еще раз с включенным скриптом;) 

# GreatCircle can be downloaded from: http://www.koders.com/python/fid0A930D7924AE856342437CA1F5A9A3EC0CAEACE2.aspx?s=coastline
from GreatCircle import *
from math import *

# Points A and B defining the rotation:
LonA = radians(0)
LatA = radians(1)
LonB = radians(45)
LatB = radians(1)

# Point C which will be translated:
LonC = radians(90)
LatC = radians(1)

# The following equation is described here: http://articles.adsabs.harvard.edu//full/1953Metic...1...39L/0000040.000.html
# It calculates the rotation pole at point P of the Great Circle defined by point A and B.
# According to http://www.tutorialspoint.com/python/number_atan2.htm
# atan2(x, y) = atan(y / x)
LonP = atan2(((sin(LonB) * tan(LatA)) - (sin(LonA) * tan(LatB))), ((cos(LonA) * tan(LatB)) - (cos(LonB) * tan(LatA))))
LatP = atan2(-tan(LatA),(cos(LonP - LonA)))
print degrees(LonP), degrees(LatP)

# The equations to calculate the translated point C location were found here: http://www.uwgb.edu/dutchs/mathalgo/sphere0.htm
# The Rotation Angle in radians:
gcAP = GreatCircle(1,1,degrees(LonA),degrees(LatA),degrees(LonP),degrees(LatP))
gcBP = GreatCircle(1,1,degrees(LonB),degrees(LatB),degrees(LonP),degrees(LatP))
RotAngle = abs(gcAP.azimuth12 - gcBP.azimuth12)

# The rotation pole P in Cartesian coordinates:
Px = cos(LatP) * cos(LonP)
Py = cos(LatP) * sin(LonP)
Pz = sin(LatP)

# Point C in Cartesian coordinates:
Cx = cos(radians(LatC)) * cos(radians(LonC))
Cy = cos(radians(LatC)) * sin(radians(LonC))
Cz = sin(radians(LatC))

# The translated point P in Cartesian coordinates:
NewCx = (Cx * cos(RotAngle)) + (1 - cos(RotAngle)) * (Px * Px * Cx + Px * Py * Cy + Px * Pz * Cz) + (Py * Cz - Pz * Cy) * sin(RotAngle)
NewCy = (Cy * cos(RotAngle)) + (1 - cos(RotAngle)) * (Py * Px * Cx + Py * Py * Cy + Py * Pz * Cz) + (Pz * Cx - Px * Cz) * sin(RotAngle)
NewCz = (Cz * cos(RotAngle)) + (1 - cos(RotAngle)) * (Pz * Px * Cx + Pz * Py * Cy + Pz * Pz * Cz) + (Px * Cy - Py * Cx) * sin(RotAngle)

# The following equation I got from http://rbrundritt.wordpress.com/2008/10/14/conversion-between-spherical-and-cartesian-coordinates-systems/
# The translated point P in lat/long:
Cr = sqrt((NewCx*NewCx) + (NewCy*NewCy) + (NewCz*NewCz))
NewCLat = degrees(asin(NewCz/Cr))
NewCLon = degrees(atan2(NewCy, NewCx))

# Output:
print str(NewCLon) + "," + str(NewCLat)
...