NubiaAtAfarTJ.nb
The direction and tangential speed of plate motion at the Afar triple junction, given the Argus et al 2011 NNR data

Coded in Mathematica v. 11.2.0 by Vince Cronin
January 30, 2019

Introduction

The purpose of this Mathematica notebook is to test and describe a module to compute the direction and speed of the tangential velocity at a given point on Earth’s surface, given the assumption that Earth is spherical.

Assumptions

We assume that Earth is a sphere of unit radius in all but one part of this analysis.  Of course, Earth is not really a sphere.  Earth has an irregular surface with deep ocean basins and tall mountains.  The average radius of Earth is around 6,371.01±0.02  km (Yoder, 1995, p. 8).  The deepest part of the ocean basins, called the Challenger Deep in the Mariana Trench, is ~11.0 km below sea level (that is, ~0.2% of Earth's mean radius), and the summit of Mt. Everest is ~8.848 km above sea level (~0.1% of Earth's mean radius).  Earth's polar radius (6356.7530 km) differs from its mean equatorial radius (6378.1363 km) by only ~21 km, so Earth's shape is different from a sphere by about 1 part in 300 or about 0.3 % (e.g., Cazenave, 1995, p. 36).  

Assuming that Earth is a sphere of unit radius makes the mathematics we use in this problem much simpler than if we used either an ellipsoid or the geoid.  For just one example, the radius of a sphere is the same measured at any point on the sphere's surface, and this is not true for either the reference ellipsoid or the geoid.  Although Earth is a little wider at its equator than along its spin axis, we do not generate serious errors in our calculations by assuming Earth's radius is the same in all directions -- that it is a sphere.   

When we say that Earth is assumed to have unit radius, that means that we set Earth's radius to 1 unit of length.  The one part of this analysis in which we must use the actual size of the Earth in kilometers is in the determination of the circumferencial distance between the reference point and the observed pointt.

Specify the input data

Input only decimal degrees (i.e., not degrees, minutes, seconds).  The sign convention follows:  north latitude is positive, south latitude is negative, east longitude is positive, west longitude is negative.  
The variables refPtLat and refPtLong are the latitude and longitude of the reference point from which the distance and azimuth to the observed point is to be determined.

refPtLat=11.317;

refPtLong=41.442;

The variables poleLat and poleLong are the latitude and longitude of the Euler pole, and omegaPole is the angular velocity around that pole.  A positive value indicates a right-hand or anticlockwise rotation around the pole, as observed from the pole toward Earth’s center.

poleLat=47.68;

poleLong=-68.44;

omegaPole=0.292;

The variable movingPlate is the name of the plate that is rotating around the Euler pole in a positive (anticlockwise) direction with respect to the reference point, which is assumed to be either fixed to the “other” plate or to a reference frame external to the lithosphere (e.g., NNR or a hotspot reference frame).  The plate name must be enclosed in quotation marks, as in “Pacific”.

movingPlate="Nubia";

User-defined functions

The function makeVector computes the unit location vector of a given point on an assumed-spherical Earth, given that point's latitude and longitude.  The x coordinate is equal to Cos[latitude]*Cos[longitude].  The y coordinate is equal to Cos[latitude]*Sin[longitude].  The z coordinate is equal to Sin[latitude].  The input latitude and longitude values must be in decimal degrees (i.e., not degrees, minutes, seconds).  The forms "lat Degree" and "long Degree" in the square brackets invokes the built-in Mathematica function Degree to convert the input data from degrees to radians, as required for the built-in Mathematica trigonometric functions (e.g., Sin, Cos, Tan, ArcSin, ArcCos, ArcTan).

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

The function unitVector finds the corresponding vector with a length of 1 (i.e., unit length) given a vector of arbitrary length.  It uses a built-in Mathematica function Norm, which returns the length of the input vector.  Given a vector A with components {(Subscript[a, x]),  (Subscript[a, y]), (Subscript[a, z])}, the length or norm of A is equal to

Norm[A] = (Sqrt[((((Subscript[a, x])))^2)+((((Subscript[a, y])))^2)+((((Subscript[a, z])))^2)]).

The unit-vector that corresponds to A has components {((Subscript[a, x])/Norm[A]), ((Subscript[a, y])/Norm[A]), ((Subscript[a, z])/Norm[A])}.

unitVector[a_]:={a[[1]]/Norm[a],a[[2]]/Norm[a],a[[3]]/Norm[a]};

<insert the descriptive text here>

tangVelocity[refPtVect_,poleVect_,angVel_]:=Module[{radEarth,circEarth,northPole,θ1,θ2,θ3,θ4,normVect1,normVect2,normVect3,distKm,bearing,speed,direction, answer},radEarth=6371.02;circEarth=2*π*radEarth;northPole={0,0,1};normVect1 = Cross[northPole,refPtVect];
normVect2 = Cross[poleVect,refPtVect];
θ1=VectorAngle[refPtVect,poleVect];
distKm=(θ1/(2*π))*circEarth;θ2=(180/π)*(VectorAngle[normVect1,normVect2]);θ3=(180/π)*(VectorAngle[normVect1,poleVect]);normVect3=Cross[refPtVect,normVect1];θ4=(180/π)*(VectorAngle[normVect3,poleVect]);bearing=If[((θ3==90)∨(θ3==180)),If[(θ4<=90),0,180],If[(θ3>90),(360-θ2),θ2]];
speed=(angVel*Sin[θ1])*(circEarth/360);
direction=If[(bearing≤270),bearing+90,(bearing-270)];answer={direction,speed,θ1,bearing,circEarth};answer];

Computation

Utilize the user-defined function makeVector to convert the input latitude and longitude data into unit location vectors

refPtVect=makeVector[refPtLat,refPtLong];

poleVect=makeVector[poleLat,poleLong];

Utilize the user-defined function tangVelocity to find the instantaneous tangential velocity of the plate at the reference point.

answer=tangVelocity[refPtVect,poleVect,omegaPole];

Output

Input Summary

Reference point {lat, long}

{refPtLat,refPtLong}

Euler pole {lat, long}

{poleLat,poleLong}

Instantaneous angular velocity around the Euler pole (degrees per Myr)

omegaPole

Name of the moving plate.

movingPlate

Output Summary

The answer is a list with the following elements:  direction or azimuth of the tangential motion (deg), tangential speed (km/Myr), angular distance between the reference point and the Euler pole θ1 (radian), angular velocity around the Euler pole angVel (deg/Myr), bearing or azimuth from the reference point to the Euler pole (deg), and the circumference of the Earth (km).

answer

direction or azimuth of the tangential motion (deg)

answer[[1]]

tangential speed (km/Myr)

answer[[2]]

References

Cazenave, A., 1995, Geoid, topography and distribution of landforms, in Ahrens, T.J., [editor], Global Earth Physics -- A Handbook of Physical Constants:  American Geophysical Union Reference Shelf 1, p. 32-39, ISBN 0-87590-851-9.
Davis, H.F., and Snider, A.D., 1987, Introduction to Vector Analysis [5th edition]:  Boston, Massachusets, Allyn and Bacon, 365 p., ISBN 0-205-10263-8.
Yoder, C.F., 1995, Astrometric and geodetic properties of Earth and the solar system, in Ahrens, T.J., [editor], Global Earth Physics -- A Handbook of Physical Constants:  American Geophysical Union Reference Shelf 1, p. 1-31, ISBN 0-87590-851-9.