InstRelVel
(Direction and speed of instantaneous relative motion)

Coded in Mathematica v. 7.0.0 by Vince Cronin
Revised March 29, 2010

Introduction

The purpose of this Mathematica notebook is to explain how to compute the direction and magnitude of the vector of instantaneous relative plate motion at a given location on Earth's surface, and to perform the computation using data provided by the user.  That is, we seek to understand where a point on one plate is going relative to another plate.

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 magnitude of the tangential velocity vector (i.e., the tangential speed) at a point.

Instantaneous Relative Velocity Models

Geoscientists have generated several models of the instantaneous relative velocities of plates.  By "instantaneous" we mean "over a short time interval."  The available data allow us to define a velocity vector for a given point (that is, the direction and speed of the point), but not the actual path of motion (or finite trajectory) of that point over an extended time interval.
The current standard model for the instantaneous relative motion of plates is the MORVEL model of DeMets and others (2010).  MORVEL is the latest in a series of models developed since the early 1970s, including the models of Chase (1972, 1978), Minster and others (1974), Minster and Jordan (1978), and the NUVEL and NUVEL-1A models of Charles DeMets and others (1990, 1994).  In the MORVEL model (DeMets and others, 2010, table 2), the North American plate moves with an anti-clockwise rotation around a pole located at latitude 48.9°N and longitude 71.7°W with an angular velocity of 0.750°/Myr relative to the Pacific plate.

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 instantaneous pole of relative motion is to be determined.

refPtLat=39.292;

refPtLong=-120;

The variables rotPoleLat and rotPoleLong are the latitude and longitude of the pole around which the observed plate is moving in a positive (anticlockwise) direction.  These data can be obtained from the MORVEL plate-motion model of DeMets and others (2010).  For the motion of the Pacific Plate relative to the stable interior of the North American Plate, the pole is located at 48.9°S latitude (-48.9) and 108.3°E longitude (108.3).

rotPoleLat=45;

rotPoleLong=-120;

The variable angVelDeg gives the angular velocity in degrees per million years.  For the motion of the Pacific Plate relative to the stable interior of the North American Plate, the angular velocity is 0.750°/Myr (Demets and others, 2010).

angVelDeg=0.750;

Operate on the input data

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]};

The function vectorAngle provides the angle between two vectors, in radians.  Recall that 360° =  (2 π) radians, so to convert from radians to degrees, multiply the radian value by (180/π).  
This function uses the vector dot product (also called an inner product or scalar product).   The notation "|A|" is the length of vector A, so if we note that vector A has components {(Subscript[a, x]),  (Subscript[a, y]), (Subscript[a, z])}, its length is

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

Given two vectors of arbitrary length, A and B, the dot product of A and B is given by

A · B = |A| |B| Cos[θ] = ((Subscript[a, x]) (Subscript[b, x])) + ((Subscript[a, y]) (Subscript[b, y])) + ((Subscript[a, z]) (Subscript[b, z]) ).

Consequently,

θ = ArcCos[(A · B)/(|A| |B|)] = ArcCos[(((Subscript[a, x]) (Subscript[b, x])) + ((Subscript[a, y]) (Subscript[b, y])) + ((Subscript[a, z]) (Subscript[b, z]) ))/(|A| |B|)].

For more information, consult any textbook on vector analysis (e.g., Davis and Snider, 1987).

vectorAngle[a_,b_]:=ArcCos[Dot[a,b]/(Norm[a]*Norm[b])];

Computation

The circumference of a circle, whose radius (r) is known, is given by

circumference = 2 * π * r

so given Earth's mean radius, in kilometers (Yoder, 1995)

meanRadiusEarth=6371;

Earth's circumference, in kilometers, is

circumferenceEarth=2*π*meanRadiusEarth;

Define the unit location vector associated with Earth' s north pole (northPole)

northPole={0,0,1};

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

refPtVect=makeVector[refPtLat,refPtLong];

rotPoleVect=makeVector[rotPoleLat,rotPoleLong];

The tangengial velocity (tanVel) varies with angular distance (θ1) from the pole of instantaneous relative motion (rotPoleVect).

θ1=vectorAngle[refPtVect,rotPoleVect];

When the point at which the velocity is to be computed is 90 ° from the pole, the velocity is equal to the product of the angular velocity (angVelDeg) multiplied by the degrees-to-kilometers conversion factor (circumferenceEarth/360  km/°).  When the point at which the velocity is to be computed is 0 ° from the pole (i.e., at the pole), the velocity is zero.  The velocity varies in accord with the sine function of the angle formed (with apex at Earth' s center) between the reference point and the rotation pole.  Recalling that Sin[0] = 0 and Sin[90 °] = 1, the full expression for the tangential velocity at any angular distance θ is given by

tangential velocity = Sin(θ) * angVelDeg * (circumferenceEarth/360).

tanVelRaw=Sin[θ1]*angVelDeg*(circumferenceEarth/360);

The built-in Mathematica function Round rounds a value to the appropriate integer, so if a value V is given in units of km-per-million years, the code (0.01 * Round[(10*V)]) would return a value in centimeters per year with two digits to the right of the decimal.

tanVel=0.01Round[10*tanVelRaw];

Find the azimuth from the reference point to the instantaneous pole of relative motion

Find the vector that is normal to the plane that contains the given two vectors, using a vector cross product.  The built-in Mathematica function Cross performs this computation.  If you had two vectors A and B, respectively specified by components {(Subscript[a, x]),  (Subscript[a, y]), (Subscript[a, z])} and {(Subscript[b, x]), (Subscript[b, y]), (Subscript[b, z])}, the result of the vector cross product A × B would be a vector C with components {(Subscript[c, x]), (Subscript[c, y]), (Subscript[c, z])} = {(((Subscript[a, y]) (Subscript[b, z]))-((Subscript[a, z]) (Subscript[b, y]))), (((Subscript[a, z]) (Subscript[b, x]))-((Subscript[a, x]) (Subscript[b, z]))), (((Subscript[a, x]) (Subscript[b, y]))-((Subscript[a, y]) (Subscript[b, x])))} that is perpendicular to the the plane containing A and B.  Given A × B = C, vector B is a right-handed rotation around C starting at A.

vectorAlpha = Cross[northPole,refPtVect];

vectorBeta = Cross[rotPoleVect,refPtVect];

Utilize the user-defined function vectorAngle to find the angle θ2  between the two vectors vectorAlpha and vectorBeta, which is also the angular difference between [1] the azimuth from the reference point toward north, and [2] the azimuth from the reference point toward the relative-motion pole.

θ2=(180/π)*(vectorAngle[vectorAlpha,vectorBeta]);

Find the angular distance (θ3), in degrees, between the location vector to the pole of instantaneous relative motion (rotPoleVect) and the vectorAlpha.  This angle helps us define the bearing from the reference point to the pole of instantaneous relative motion.

θ3=(180/π)*(vectorAngle[vectorAlpha,rotPoleVect]);

Imagine that Earth is divided into two hemispheres bounded by the plane that contains Earth's spin axis and the reference point.  We will specify that the north pole is "up" in this frame of reference, so one hemisphere would be to the right of the reference point and the other would be to the left.
If the pole of instantaneous relative motion is located along the same circle of longitude as the reference point, then θ3 will equal 90° or 180°, and the bearing from the reference point to the pole will either be 0° or 180°.
If the pole of instantaneous relative motion is located on the hemisphere to the right of the reference point, then θ3 will be less than 90°.  In that case, the azimuth or bearing from the reference point to the instantaneous pole of relative motion is equal to θ2.  Bearings are expressed in degrees relative to north, where due north is 0°, due east is 90°, due south is 180°, and due west is 270°.
If the pole of instantaneous relative motion is located on the hemisphere to the left of the reference point, then θ3 will be greater than 90°.  In that case, the bearing from the reference point to the instantaneous pole of relative motion is equal to (360°- θ2).

The vector vectorGamma and the angle θ4 are used when the rotational pole is on the same meridian as the reference point, to differentiate between a bearing of 0° (north) or 180° (south) from the reference point to the rotational pole.  If θ4 is less than 90°,  the pole of instantaneous relative motion is north of the reference point.  If θ4 is equal to 90°,  the rotational pole is either coincident with the reference point or 180° from the reference point, in which case the bearing from the reference point to the rotational pole is undefined.  For simplicity, when θ4=90°, we set the bearing to 0°.

vectorGamma=Cross[refPtVect,vectorAlpha];

θ4=(180/π)*(vectorAngle[vectorGamma,rotPoleVect]);

The code fragment "If[((θ3 == 90) ∨ (θ3 == 180))" in the following line  literally means "if θ3 equals 90 or if θ3 equals 180."

bearing=If[((θ3==90)∨(θ3==180)),If[(θ4<=90),0,180],If[(θ3>90),(360-θ2),θ2]];

The direction of relative motion is 90 ° clockwise from the azimuth to the pole.  This is because the pole defines the location of the axis through Earth's center around which the plate being observed is instantaneously rotating in an anti-clockwise direction relative to an observer on the other plate.  By the "right-hand" convention used in science and engineering, clockwise rotations are considered negative rotations, and anti-clockwise (counter-clockwise) rotations are considered positive rotations.
Once the bearing from the reference point to the pole of instantaneous relative motion is defined, the azimuth of instantaneous motion is equal to (bearing + 90°).  If this sum is greater than 360°, subtract 360° to obtain a sensible (and valid) azimuth.  For example, an azimuth computed to be 395° is correctly stated as an azimuth of 35°.

rawAzimuth=If[(bearing+90>360),((bearing+90)-360),(bearing+90)];

Round the raw azimuth  (rawAzimuth) of instantaneous relative motion to the nearest whole degree.

azimuth=0.01 Round[(100*rawAzimuth)];

Output

The direction of tangential relative motion is given in degrees relative to north, where due north is 0°, due east is 90 °, due south is 180°, and due west is 270 °:

azimuth

The speed of tangential relative motion is given in centimeters per year:

tanVel

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.
Chase, C.G., 1972, The n-plate problem of plate tectonics:  Geophysical Journal of the Royal Astronomical Society, v. 29, p. 117–122.
Chase, C.G., 1978, Plate kinematics: the Americas, East Africa, and the rest of the world:  Earth and Planetary Science Letters, v. 37, p. 355–368.
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.
DeMets, C., Gordon, R.G., Argus, D.F. & Stein, S., 1990, Current plate motions:  Geophysical Journal International, v. 101, p. 425–478.
DeMets, C., Gordon, R.G., and Argus, D.F., 2010, Geologically current plate motions:  Geophysical Journal International, v.  121, p. x, doi: 10.1111/j.1365-246X.2009.04491.x
DeMets, C., Gordon, R.G., Argus, D.F., and Stein, S., 1994, Effect of recent revisions to the geomagnetic reversal timescale:  Geophysical Research Letters, v. 21, p. 2191–2194.
Minster, J.B., and Jordan, T.H., 1978, Present-day plate motions:  Journal of Geophysical Research, v. 83, p. 5331–5354.
Minster, J.B., Jordan, T.H., Molnar, P., and Haines, E., 1974, Numerical modeling of instantaneous plate tectonics:  Geophysical Journal of the Royal Astronomical Society, v. 36, p. 541–576.
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.