TangentialVelocity.nb

Written in Mathematica 12.2.0.0 Vince Cronin
Begun March 9, 2021, with some material from earlier documents;  revised March 14.

The Question

What is the instantaneous tangential velocity of the Pacific plate, relative to the North American plate, at Wallace Creek along the San Andreas fault?

A Chosen Point

We will look at a point along the San Andreas fault in central California at Wallace Creek:  latitude 35.271641, longitude -119.827653.

refPtLat=35.271641;
refPtLong=-119.827653;

A Model of the Instantaneous Motion of Plates Relative To Each Other

We will use data from Kreemer, Blewitt, and Klein (2014), table 2, for the instantaneous motion of the North American plate relative to the Pacific plate:  pole latitude 49.30°, pole longitude -76.01°, angular velocity 0.791°/Myr.  The antipole of this pole is the pole of instantaneous motion of the Pacific plate relative to North America.  The antipole latitude (-49.30°) is the inverse of the pole latitude (49.30°).  The antipole longitude (103.99°) is 180° from the pole longitude (-76.01°).

paRelNALat=-49.30;
paRelNALong=103.99;

omegaPArelNA=0.791;

User-Defined Functions and Modules

The function makeLocationVector[] converts location data from geographic coordinates (latitude and longitude) to Cartesian geographic coordinates.

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

The module findAzimuthVectInput[] determines the azimuth from a reference point (location vector refPtVect) to a target point (location vector otherPtVect).

findAzimuthVectInput[refPtVect_,otherPtVect_]:=Module[{northPole,vectorA,vectorB,θ1,θ2,vectorC,θ3,azimuth},northPole={0,0,1};vectorA = Cross[northPole,refPtVect];vectorB = Cross[otherPtVect,refPtVect];θ1=VectorAngle[vectorA,vectorB](180/π);θ2=VectorAngle[vectorA,otherPtVect](180/π);vectorC=Cross[refPtVect,vectorA];θ3=VectorAngle[vectorC,otherPtVect](180/π);azimuth=If[((θ2==90)∨(θ2==180)),If[(θ3<=90),0,180],If[(θ2>90),(360-θ1),θ1]];azimuth ];

The module findAzimuthGeogInput[] determines the azimuth from a reference point (location vector refPtLat and refPtLong ) to a target point (location vector otherPtLat and otherPtLong).

findAzimuthGeogInput[refPtLat_,refPtLong_,otherPtLat_,otherPtLong_]:=Module[{refPtVect,otherPtVect,northPole,vectorA,vectorB,θ1,θ2,vectorC,θ3,azimuth},refPtVect={Cos[refPtLat Degree] Cos[refPtLong Degree],Cos[refPtLat Degree] Sin[refPtLong Degree],Sin[refPtLat Degree]};otherPtVect={Cos[otherPtLat Degree] Cos[otherPtLong Degree],Cos[otherPtLat Degree] Sin[otherPtLong Degree],Sin[otherPtLat Degree]};northPole={0,0,1};vectorA = Cross[northPole,refPtVect];vectorB = Cross[otherPtVect,refPtVect];θ1=VectorAngle[vectorA,vectorB](180/π);θ2=VectorAngle[vectorA,otherPtVect](180/π);vectorC=Cross[refPtVect,vectorA];θ3=VectorAngle[vectorC,otherPtVect](180/π);azimuth=If[((θ2==90)∨(θ2==180)),If[(θ3<=90),0,180],If[(θ2>90),(360-θ1),θ1]];azimuth ];

The module findTangVelGeog[] determines the tangential speed (in mm/yr) and azimuth of relative motion at a particular point (geographic coordinates ptLat and ptLong ), given the pole of positive (anticlockwise) motion (geographic coordinates poleLat and poleLong) and the angular velocity (angVel) expressed in degrees per million years.
•  This module calls (and so requires) the user-defined module findAzimuthVectInput[]
•  The output of this module is a list of two numbers.  The first number is the tangential speed, expressed in millimeters per year.  The second number is the azimuth, expressed in degrees.

findTangVelGeog[ptLat_,ptLong_,poleLat_,poleLong_,angVel_]:=Module[{ptCart,poleCart,phi,tangSpeed,azToPole,azimuth,answer},ptCart={Cos[ptLat Degree]*Cos[ptLong Degree],Cos[ptLat Degree]*Sin[ptLong Degree],Sin[ptLat Degree]};poleCart={Cos[poleLat Degree]*Cos[poleLong Degree],Cos[poleLat Degree]*Sin[poleLong Degree],Sin[poleLat Degree]};phi=VectorAngle[ptCart,poleCart];tangSpeed=(Sin[phi](angVel*((2π 6371010000)/360)))/1000000;azToPole=findAzimuthVectInput[ptCart,poleCart];azimuth=If[azToPole>270,azToPole-270,azToPole+90];answer={tangSpeed,azimuth};answer];

The module findTangVelCart[] determines the tangential speed (in mm/yr) and azimuth of relative motion at a particular point (location vector ptCart), given the pole of positive (anticlockwise) motion (location vector poleCart) and the angular velocity (angVel) expressed in degrees per million years.
•  This module calls (and so requires) the user-defined module findAzimuthVectInput[]
•  The output of this module is a list of two numbers.  The first number is the tangential speed, expressed in millimeters per year.  The second number is the azimuth, expressed in degrees.

findTangVelCart[ptCart_,poleCart_,angVel_]:=Module[{phi,tangSpeed,azToPole,azimuth,answer},phi=VectorAngle[ptCart,poleCart];tangSpeed=(Sin[phi](angVel*((2π 6371010000)/360)))/1000000;azToPole=findAzimuthVectInput[ptCart,poleCart];azimuth=If[azToPole>270,azToPole-270,azToPole+90];answer={tangSpeed,azimuth};answer];

Computation

Convert input geographic coordinates to Cartesian geographic coordinates

refPtCart=makeLocationVector[refPtLat,refPtLong];

paRelNACart=makeLocationVector[paRelNALat,paRelNALong];

Solve the problem “longhand”

The variable phi is the angle (in radians) between the unit location vectors to the reference point and the pole of rotation.

phi=VectorAngle[ptCart,poleCart];

The variable tangSpeed is the magnitude of the tangential velocity vector at the reference point

tangSpeed=(Sin[phi](omegaPArelNA*((2π 6371010000)/360)))/1000000;

The variable azToPole is the azimuth (in degrees relative to north, measured clockwise from north) from the reference point to the pole of rotation.

azToPole=findAzimuthVectInput[ptCart,poleCart];

The variable azimuth is the azimuth (in degrees relative to north, measured clockwise from north) of the tangential motion vector at the reference point

azimuth=If[azToPole>270,azToPole-270,azToPole+90];

The variable answer1 is a list of two numbers.  The first number is the tangential speed, expressed in millimeters per year.  The second number is the azimuth, expressed in degrees.

answer1={tangSpeed,azimuth};

Solve the problem using findTangVelGeog[ ]

In this case, we use only the raw input data that was expressed in geographic coordinates.

answer2=findTangVelGeog[refPtLat,refPtLong,paRelNALat,paRelNALong,omegaPArelNA];

Solve the problem using findTangVelCart[ ]

In this case, we use the Cartesian vectors that we computed using the raw input data that was expressed in geographic coordinates.

answer3=findTangVelCart[refPtCart,paRelNACart,omegaPArelNA];

Output

Output from “longhand” solution

Tangential speed, expressed in millimeters per year.

N[answer1[[1]]]

Azimuth of the tangential velocity vector, expressed in degrees.

N[answer1[[2]]]

Output solution using findTangVelGeog[ ]

Tangential speed, expressed in millimeters per year.

N[answer2[[1]]]

Azimuth of the tangential velocity vector, expressed in degrees.

N[answer2[[2]]]

Output solution using findTangVelCart[ ]

Tangential speed, expressed in millimeters per year.

N[answer3[[1]]]

Azimuth of the tangential velocity vector, expressed in degrees.

N[answer3[[2]]]

References

Boyd, S., and Vandenberghe, L., 2018, Introduction to Applied Linear Algebra — Vectors, Matrices, and Least Squares:  Cambridge, UK, Cambridge University Press, 463 p.

Davis, H.F., and Snider, A.D., 1987, Introduction to Vector Analysis [5th edition]:  Boston, Massachusetts, Allyn and Bacon, 365 p.

Goldstein, H., 1980, Classical Mechanics [2nd edition]:  Reading, Massachusetts, Addison-Wesley Publishing Company, 672 p.

Kreemer, C., Blewitt, G., and Klein, E.C., 2014, A geodetic plate motion and Global Strain Rate Model:  Geochemistry, Geophysics, Geosystems, v. 15, p. 3849-3889, doi:10.1002/2014GC005407.
Supplementary table S2 is accessible as an Excel spreadsheet via https://CroninProject.org/Vince/PlateKinematics/KinematicsCourse/KreemerEtAl2014TableS2.xls