Circular Rotation Modules . nb

Vince Cronin, begun March 28, 2021

Explanation of the module

When the pole (the point on Earth’s surface around which our point moves in a positive or anti-clockwise direction as you are looking down on the pole for beyond Earth) is not coincident with the the North Pole.
The axes of the Cartesian geographic coordinate system are
   X = {1,0,0} and X is located at the intersection of the Equator and the Prime Meridian
   Y = {0,1,0} and Y is at the point latitude 0 longitude 90°E
   Z = {0,0,1} and Z is located at the North Pole.
Note that in a right orthogonal Cartesian coordinate system, 
   X is a right-handed rotation from Y to Z,
   Y is a right-handed rotation from Z to X, and
   Z is a right-handed rotation from x to Y.
Or, more usefully employing unit-vector cross products, 
      (OverHat[x]) = (OverHat[y]) × (OverHat[z])
      (OverHat[y]) = (OverHat[z]) × (OverHat[x])
      (OverHat[z]) = (OverHat[x]) × (OverHat[y])

The axes of the Cartesian pole-geographic coordinate system are based on the Cartesian geographic coordinates of the pole of rotation.  The Z axis of the Cartesian pole-geographic coordinate system coincides with the location vector to the pole of rotation

The X axis of the Cartesian pole-geographic coordinate system is the normalized vector result of the cross product of the pole location vector and the location vector to the North Pole.

The Y axis of the Cartesian pole-geographic coordinate system is the vector result of the cross product of the Z axis of the Cartesian pole-geographic coordinate system and the X axis of the Cartesian pole-geographic coordinate system.

The transformation matrix from the Cartesian geographic coordinate system to the Cartesian pole-geographic coordinate system is called j1, where

j1=(({{xPoleGeog[[1]], xPoleGeog[[2]], xPoleGeog[[3]]}, {yPoleGeog[[1]], yPoleGeog[[2]], yPoleGeog[[3]]}, {zPoleGeog[[1]], zPoleGeog[[2]], zPoleGeog[[3]]}}));

This transformation matrix is simplified relative to the general Cartesian-coordinate transformation matrix because the basis vectors of the positive geographic axes all have components of either 0 or 1.

The inverse of j1 is j3, which is the transformation matrix from the Cartesian pole-geographic coordinate system back to the Cartesian geographic coordinate system, and our old friend, the z-rotation matrix, is here called j2:

j2=(({{Cos[θ Degree], -Sin[θ Degree], 0}, {Sin[θ Degree], Cos[θ Degree], 0}, {0, 0, 1}}));

The Cartesian coordinates of a point (ptInitLocVect) rotated by θ degrees around the pole of rotation (poleVect) is ptRotLocVect, where

ptRotLocVect  =  j3 · j2 · j1 · ptInitLocVect

User - defined functions and modules

makeLocationVector[lat_, long_]:={Cos[lat Degree] Cos[long Degree], Cos[lat Degree] Sin[long Degree], Sin[lat Degree]};

findGeogCoord[vect_]:=Module[{lat,long,a,b,θ},a=ArcSin[vect[[3]]];θ=ArcCos[vect[[1]]/(Sqrt[(((vect[[1]])^2)+((vect[[2]])^2))])];b=If[(vect[[2]]<0),(-θ),(θ)];lat=a(180/π);long=If[((Abs[vect[[1]]]<(1×(10^(-14))))&&(Abs[vect[[2]]]<(1×(10^(-14))))),0,(b(180/π))];{lat,long}];

The following user-defined module circularRotationFromGeog[] acts on the the following input arguments:  the rotation angle and the latitude and longitude of the pole and a point.  The output will be the unit location vector of the rotated point, which can then be converted back to geographic coordinates using the function findGeogCoord[].

circularRotationFromGeog[ptLat_,ptLong_,poleLat_,poleLong_,rotAngleDeg_]:=Module[{xGeog,yGeog,zGeog,xPole,yPole,zPole,ptCartGeog,poleCartGeog,j1,j2,j3,resultVect},xGeog={1,0,0};yGeog={0,1,0};zGeog={0,0,1};ptCartGeog={Cos[ptLat Degree] Cos[ptLong Degree], Cos[ptLat Degree] Sin[ptLong Degree], Sin[ptLat Degree]};poleCartGeog={Cos[poleLat Degree] Cos[poleLong Degree], Cos[poleLat Degree] Sin[poleLong Degree], Sin[poleLat Degree]};zPole=poleCartGeog;
xPole=Cross[zPole,zGeog];yPole=Cross[zPole,xPole];j1=(({{xPole[[1]], xPole[[2]], xPole[[3]]}, {yPole[[1]], yPole[[2]], yPole[[3]]}, {zPole[[1]], zPole[[2]], zPole[[3]]}}));
j2=(({{Cos[rotAngleDeg Degree], -Sin[rotAngleDeg Degree], 0}, {Sin[rotAngleDeg Degree], Cos[rotAngleDeg Degree], 0}, {0, 0, 1}}));j3=Inverse[j1];resultVect=j3.j2.j1.ptCartGeog;resultVect];

The following user-defined module circularRotationFromVect[] acts on the the following input arguments:  the rotation angle and the unit location vectors to the pole and a point, respectively.  The output will be the unit location vector of the rotated point, which can then be converted back to geographic coordinates using the function findGeogCoord[].

circularRotationFromVect[ptCartGeog_,poleCartGeog_,rotAngleDeg_]:=Module[{xGeog,yGeog,zGeog,xPole,yPole,zPole,j1,j2,j3,resultVect},xGeog={1,0,0};yGeog={0,1,0};zGeog={0,0,1};zPole=poleCartGeog;
xPole=Cross[zPole,zGeog];yPole=Cross[zPole,xPole];j1=(({{xPole[[1]], xPole[[2]], xPole[[3]]}, {yPole[[1]], yPole[[2]], yPole[[3]]}, {zPole[[1]], zPole[[2]], zPole[[3]]}}));
j2=(({{Cos[rotAngleDeg Degree], -Sin[rotAngleDeg Degree], 0}, {Sin[rotAngleDeg Degree], Cos[rotAngleDeg Degree], 0}, {0, 0, 1}}));j3=Inverse[j1];resultVect=j3.j2.j1.ptCartGeog;resultVect];

Input Test Data

The variables ptLat and ptLong together define the location of the reference point, expressed in geographic coordinates in decimal degrees.

ptLat=35;

ptLong=30;

The variables poleLat and poleLong together define the location of the pole of rotation, expressed in geographic coordinates in decimal degrees.

poleLat=65;

poleLong=70;

The variable θ is the angle of rotation, expressed in degrees.  A positive value indicates a positive (right-handed or anti-clockwise) rotation around the given pole as viewed looking at the pole from outside of Earth toward the Earth’s center.

θ=30;

Computation

The variable poleLocVect is the location vector to the pole of rotation, expressed in Cartesian geographic coordinates.  It is determined using the user-defined function makeLocationVector[], acting on the pole location expressed in geographic coordinates.

poleLocVect=makeLocationVector[poleLat,poleLong];

The variable ptLocVect is the location vector to a given reference point, expressed in Cartesian geographic coordinates.  It is determined using the user-defined function makeLocationVector[], acting on the point location expressed in geographic coordinates.

ptLocVect=makeLocationVector[ptLat,ptLong];

resultVector1=circularRotationFromGeog[ptLat,ptLong,poleLat,poleLong,θ];

resultGeog1=findGeogCoord[resultVector1];

resultVector2=circularRotationFromVect[ptLocVect,poleLocVect,θ];

resultGeog2=findGeogCoord[resultVector2];

Output Results

The result of using the module circularRotationFromGeog[], acting on input data expressed in geographic coordinates, is given below.  The first number in the set is the decimal latitude (°N) and the second is the decimal longitude (°E).

N[resultGeog1]

The result of using the module circularRotationFromVect[], acting on input data expressed in Cartesian geographic coordinates, is given below.  The first number in the set is the decimal latitude (°N) and the second is the decimal longitude (°E).

N[resultGeog2]