I’ve written the below python script. The idea is to calculate the new location of point C after you rotate the globe from point A to point B. I first calculate point P, which is the rotation pole. With calculating point P already something goes wrong. With the following input f.e. I would assume point P to be having latitude 90 or –90.
I asked this question before here: Rotate a sphere from coord1 to coord2, where will coord3 be? But I figured it's better to ask again with the script included ;)
# 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)