Cycloid Relative Motion Model

Incomplete draft -- subject to revision

Written in Mathematica 4.0 by Vince Cronin, April 2000;  revised 28 February 2012

Introduction

The cycloid relative motion model is the simplest finite motion model that addresses the 3-plate problem of Cox (1973).  It also leads to the development of oceanic fracture zones that are sigmoid in shape, as observed in nature.  Unlike the circular relative motion model, the cycloid model is kinematically admissible for finite relative plate motion in a system in which the instantaneous poles of relative motion are not coaxial/coincident (e.g., Demets and others, 2010).

Assumptions

We assume that Earth is a sphere of unit radius.  We also assume, in the first-order model, that the angular velocity vectors that we input are invariant during the time interval that we are considering.

A word about coordinate systems

All motion is described relative to some reference system.  Several coordinate systems are employed in this analysis, some of which are fixed relative to each other, while others are moving.  All coordinate systems that are used in this analysis are right orthogonal Cartesian coordinate systems that share a common origin at the center of the Earth.  In right orthogonal Cartesian coordinate systems, there are three axes (X, Y and Z) of the sort first described by René Descartes (i.e., they are Cartesian axes) oriented perpendicular to each other (i.e., they are orthogonal) that follow the naming convention consistent with the right-hand rule (i.e., they are right orthogonal Cartesian axes).

Basis vectors are unit vectors that are coincident with each of the coordinate axes.  The set of 3 basis vectors for a given right orthogonal Cartesian coordinate system would have coordinates of {1, 0, 0} for the basis vector associated with the X axis, {0, 1, 0} for the Y-axis basis vector, and {0, 0, 1} for the Z-axis basis vector.

A separate Mathematica notebook will consider the equivalent analysis employing polar, rather than Cartesian, coordinates.  The distance from the center of Earth to Earth's surface along each of these axes will be assumed to be the same and equal to 1 Earth radius.  So we assume that Earth is a unit sphere.

Geographic coordinate system

The geographic coordinate system is defined for the present time, with the positive Z axis extending to the north spin axis of Earth, the positive X axis extending to the intersection of the Equator (latitude = 0) and the Prime Meridian (longitude = 0), and the positive Y axis extending to the intersection of the Equator with latitude 90° east.  Given a point whose latitude and longitude are known, the Cartesian coordinates {x, y, z} of that point in the geographic coordinate system are

({{({{(({{x}, {y}, {z}}))}})=(({{cos[latitude] cos[longitude]}, {cos[latitude] sin[longitude]}, {sin[latitude]}}))}})

or, more explicitly,

({{({{({{x=cos[latitude]*cos[longitude] }, {y=cos[latitude]*sin[longitude]}, {z=sin[latitude]}})}})}})

where I use the asterisk (*) to indicate multiplication.   The resulting coordinates {x, y, z} define the unit location vector to the given point.  A unit location vector is a vector with a length of 1 (i.e., it has unit length) that extends from Earth' s center to that point on Earth' s surface.

Reference coordinate system external to the lithosphere

It is safe to say that I have spent more time thinking about reference frames external to the lithosphere that might be useful for plate-kinematic work than just about anything (other than how to throw a good overhand curveball).  And I am not alone.  As discussed elsewhere in this resource, reference frames external to the lithosphere have been described that utilize hotspots, the constellation of GPS satellites, results from very-long baseline interferometry (VLBI), and inversions of relative-motion data to identify a no-net-rotation (no-net-torque) reference frame (e.g., http://sps.unavco.org/crustal_motion/dxdt/model/).  For the purposes of this discussion, I assume that a reference frame external to the lithosphere has been chosen, so that the present-day instantaneous motion of each individual plate is known relative to that reference frame.

One older example of a reference frame that is external to the lithosphere is provided by the HS2-NUVEL1 model of Gripp and Gordon (1990, Geophys. Res. Lett. 17, 1109-1112).  This reference frame is based on the interpreted track of the Hawaiian hot spot on the Pacific plate, and used the NUVEL1 angular velocities of plates moving relative to each other to define the velocity of each plate relative to the hot spot.  The coordinates of the pole of positive (anticlockwise) rotation of the Pacific and North American plates are as follows:
    Plate                     Latitude      Longitude      Angular Velocity (°/Myr)
North America         -67.2               -11.1                       0.28
Pacific                      -60.192             90.031                   0.9753

Negative values correspond to south latitudes and west longitudes.

Observer's geographic coordinate system

The observer's geographic coordinate system (OGCS) is defined at the present time, when it coincides exactly with the geographic coordinate system (GCS).  What makes the OGCS different from the GCS is that the OGCS is fixed to the observer's plate throughout the finite time interval, so as the observer's plate moves relative to Earth's spin axis and the other plates, so too will the OGCS move relative to the GCS.  From the perspective of the observer's plate, the GCS will rotate around the pole to the observer's plate in a negative (clockwise) manner with the same angular speed as the observer's plate (|(Subscript[ω, a])|).

Other coordinate systems associated with the two plates

Two additional right orthogonal Cartesian coordinate systems that are fixed to the reference coordinate system are used in the cycloid model.  One might think of these as coordinate subsystems of the reference coordinate system, since all three sets of coordinate axes are rigidly fixed to one another throughout the time interval that we are modeling.  
The positive Z axis of the coordinate system associated with the observer's plate (with basis vector called obsZ in the following code) passes through the Euler pole around which the observer's plate rotates in a positive (anticlockwise) direction, as observed in the reference coordinate system.   Similarly, the positive Z axis of the coordinate system associated with the plate that is moving relative to the observer's plate (movZ) passes through the moving plate's Euler pole.  
Let's represent the angular speed of the observer's plate with the symbol |(Subscript[ω, obs])| and the angular speed of the plate that is moving relative to the observer with the symbol |(Subscript[ω, mov])|.  Both of these velocities are measured relative to the reference coordinate system.  If x ≤ x, then the positive Y axis in both coordinate subsystems (i.e., obsY and movY) is coincident with the product of the vector cross product obsZ × movZ .

Specification of some 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 *(π/180)" and "long *(π/180)" in the square brackets 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*(π/180)]*Cos[long*(π/180)],Cos[lat*(π/180)]*Sin[long*(π/180)],Sin[lat*(π/180)]};

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

The functions vectorLat and vectorLong convert a Cartesian location vector back to a geographic coordinate system.

vectorLat[x_]:=ArcSin[x[[3]]]*(180/π);

vectorLong[x_]:=(x[[2]]/Abs[x[[2]]])*(ArcCos[x[[1]]/√((((x[[1]])^2))+(((x[[2]])^2)))])*(180/π);

Input Data

Angular velocity vector for the observer' s plate

Input latitude and longitude data in decimal degrees (i.e., not in degrees, minutes, seconds or in degrees ).  The sign convention follows :  north latitude is positive, south latitude is negative, east longitude is positive, west longitude is negative.

Input variables latObs and longObs are the initial latitude and longitude of the pole to the observer's plate.

latObs = -60.192;

longObs = 90.031;

Input variable obsω is the angular speed around that positive pole, in degrees per million years.

obsω = 0.9753;

Angular velocity vector for the plate that is moving relative to the observer (i.e., the moving plate)

Input variables latMov and longMov are the initial latitude and longitude of the pole to the moving plate.

latMov = -67.2;

longMov = -11.1;

Input variable movω is the angular speed around that positive pole, in degrees per million years.

movω = 0.28;

Initial location of the reference point on the moving plate

Input variables latRP and longRP are the initial latitude and longitude of the reference point Q on the moving plate.

latRP= 34.05;

longRP= -118.24;

Time for which the relative position is to be calculated

Input variables (Subscript[t, n]) is the modeled time, in million years.  Positive values are future time, negative are past time.

modelTime = 1;

Computation

pointRP=makeVector[latRP,longRP];

Establish basis vectors

Specify the coordinates of the basis vectors associated with the geographic coordinate system {geogX, geogY, geogZ}.

geogZ={1,0,0};

geogY={0,1,0};

geogX={0,0,1};

Specify the coordinates of the basis vectors for the coordinate system associated with the observer's plate {obsX, obsY, obsZ} and the coordinate system associated with the moving plate {movX, movY, movZ}, relative to the geographic coordinate system.

obsZ=makeVector[latObs,longObs];

movZ =makeVector[latMov,longMov];

obsY=unitVector[Cross[obsZ,movZ]];

movY = obsY;

obsX = unitVector[Cross[obsY,obsZ]];

movX = unitVector[Cross[movY,movZ]];

Establish transformation matrices

Matrix j1 transforms input vector coordinates expressed in the observer's geographic coordinate system to a the coordinate system associated with the moving plate.

j1 = (({{movX[[1]], movX[[2]], movX[[3]]}, {movY[[1]], movY[[2]], movY[[3]]}, {movZ[[1]], movZ[[2]], movZ[[3]]}}));

MatrixForm[j1]

Matrix j5 transforms input vector coordinates expressed in the coordinate system associated with the observer's plate to the observer's geographic coordinate system.

j5 = (({{obsX[[1]], obsY[[1]], obsZ[[1]]}, {obsX[[2]], obsY[[2]], obsZ[[2]]}, {obsX[[3]], obsY[[3]], obsZ[[3]]}}));

Establish rotation matrices

The matrix j2 reflects the rotation of the moving plate around its pole at an angular speed of movω during the modeled time interval between time zero (today) and some other time (modelTime).

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

We determine the value and sign of the angle γ between the poles to the moving and observer' s plates, measured in the direction from the pole of the moving plate to the pole of the observer's plate.  The range of γ is 180° ≥ γ ≥ -180°.  An anti-clockwise rotation around the vector result of the vector cross product (Subscript[ω, slower]) × (Subscript[ω, faster]) (where (Subscript[ω, slower]) ≤  (Subscript[ω, faster])) yields a positive value of γ.  So if our input values are such that where (Subscript[ω, a]) ≤  (Subscript[ω, b]) , then an anti-clockwise rotation from movZ to obsZ around the vector result of (Subscript[ω, a])× (Subscript[ω, b]) yields a positive value of γ.

a=vectorAngle[obsZ,movX];

If[a<(π/2),(s=1),(s=(-1))];

γ = (s) (vectorAngle[movZ,obsZ]);

The matrix j3 reflects a rotation of coordinate axes around a common Y axis from a system in which the pole to the moving plate lies along the Z axis of a system in which the pole to the observer's plate lies along the Z axis.

j3 =(({{Cos[γ], 0, -Sin[γ]}, {0, 1, 0}, {Sin[γ], 0, Cos[γ]}}));

The matrix j4 reflects the rotation of the observer's plate around its pole at an angular speed of (Subscript[ω, obs]) during the modeled time interval (Δt).

j4 = (({{Cos[(obsω Degree) modelTime], Sin[(obsω Degree) modelTime], 0}, {-Sin[(obsω Degree) modelTime], Cos[(obsω Degree) modelTime], 0}, {0, 0, 1}}));

Compute the cycloid displacement

Multiply the matrices in proper sequence to model the new location of the reference point (pointRP2) relative to the observer's plate.

pointRP2=j5.j4.j3.j2.j1.pointRP;

Convert result to latitude and longitude in the observer's geographic coordinates

Convert the result back to the observer's geographic coordinates.  Of course, the observer's plate will have moved relative to the geographic coordinates that are fixed to Earth's spin axis.  We express the new location relative to the geographic coordinates as they were relative to the observer's plate at time zero (today) for convenience in plotting the motion.  If we wanted to account for the motion of the observer's plate relative to the actual geographic coordinates, we would rotate the geographic coordinates clockwise around the pole of the observer's plate at an angular speed of (Subscript[ω, a]).

endLat=N[IntegerPart[(100*vectorLat[pointRP2])+0.5]/100];

endLong= N[IntegerPart[(100*vectorLong[pointRP2])+0.5]/100];

Output:  Compare the original location with the new location.

originalLocation={latRP,longRP}

otherLocation={endLat,endLong}

10.2  User-defined functions

We will use the following user-defined function developed in a previous chapter.

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

unitVect3D[vect_]:={(vect[[1]]/Norm[vect]),(vect[[2]]/Norm[vect]),(vect[[3]]/Norm[vect])};

findGeogCoord[vect_]:=Module[{lat,long,a,b,c,d,e,f},a=ArcSin[vect[[3]]];b={vect[[1]],vect[[2]],0};c=If[((Abs[vect[[1]]]<(1×(10^(-14))))&&(Abs[vect[[2]]]<(1×(10^(-14))))),{1,1,0},{vect[[1]]/Norm[b],vect[[2]]/Norm[b],0}];d={1,0,0};e=VectorAngle[c,d];f=If[(vect[[2]]<0),(-e),(e)];lat=a(180/π);long=If[((Abs[vect[[1]]]<(1×(10^(-14))))&&(Abs[vect[[2]]]<(1×(10^(-14))))),0,(f(180/π))];{lat,long}];

makeGreatCircle[normal_]:=Module[{a,x1,y1,z1,x2,y2,z2,j1,j2,j3,θ,b},a=Table[{Cos[i Degree],Sin[i Degree],0},{i,0,360,5}];x1={1,0,0};y1={0,1,0};z1={0,0,1};z2=normal;x2=unitVect3D[Cross[z2,z1]];y2=unitVect3D[Cross[z2,x2]];j1={{x1.x2,y1.x2,z1.x2},{x1.y2,y1.y2,z1.y2},{x1.z2,y1.z2,z1.z2}};θ=VectorAngle[z1,z2];j2={{1,0,0},{0,Cos[θ],Sin[θ]},{0,-Sin[θ],Cos[θ]}};j3=Inverse[j1];b=Table[j3.j2.j1.a[[i]],{i,1,Length[a]}];b]

cycloid[obsΩ_,movΩ_,refPt_,t_]:=Module[{a,b,c,d,e,f,g,h,i,j,k,l,m1, m2, m3, m4, m5,n,p,answer},a={1,0,0};b={0,1,0};c={0,0,1};d=unitVect3D[obsΩ];e=Norm[obsΩ];f=unitVect3D[movΩ];g=Norm[movΩ];h=unitVect3D[d\[Cross]f];i=h;j=unitVect3D[h\[Cross]d];k=unitVect3D[i\[Cross]f];m1={{k[[1]],k[[2]],k[[3]]},{i[[1]],i[[2]],i[[3]]}, {f[[1]],f[[2]],f[[3]]}};m2={{Cos[(-g Degree) t],Sin[(-g Degree) t],0},{-Sin[(-g Degree) t],Cos[(-g Degree) t],0},{0,0,1}};l=VectorAngle[d,k];If[l<(π/2),(n=1),(n=(-1))];p = (n) (VectorAngle[f,d]);m3={{Cos[p],0,-Sin[p]},{0,1,0},{Sin[p],0,Cos[p]}};m4={{Cos[(e Degree) t],Sin[(e Degree) t],0},{-Sin[(e Degree) t],Cos[(e Degree) t],0},{0,0,1}};m5={{j[[1]],h[[1]],d[[1]]},{j[[2]],h[[2]],d[[2]]},{j[[3]],h[[3]],d[[3]]}};answer=m5.m4.m3.m2.m1.refPt;answer];

observer=obsω*convert2Cart[latObs,longObs];

moving=movω*convert2Cart[latMov,longMov];

refPoint=convert2Cart[latRP,longRP];

time=1;

cycloid[observer,moving,refPoint,time]

findGeogCoord[cycloid[observer,moving,refPoint,time]]

originalLocation={latRP,longRP}

otherLocation={endLat,endLong}

References Cited and Selected Bibliography

Altamimi, Z., P. Sillard, and C. Boucher, 2002, ITRF2000 -- A new release of the International Terrestrial Reference Frame for earth science applications:  Journal of Geophysical Research, v. 107(B10), 2214, doi:10.1029/2001JB000561.

Argus, D.F. and R.G. Gordon, 1991, No-net-rotation model of current plate velocities incorporating plate motion model NUVEL-1:   Geophysical Research Letters, v. 18, p. 2039-2042.

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.

Cox, A., (Compiler), 1973, Plate tectonics and geomagnetic reversals:  San Francisco, California, W.H. Freeman, 702 p.

Cronin, V.S., 1992, Types and kinematic stability of triple junctions:  Tectonophysics, v. 207, p. 287-301.

Cronin, V.S., 1992, A kinematic perspective on finite relative plate motion, provided by the first-order cycloid model, in Boucot, A.J., Chatterjee, S., Dickins, J.M., and Hotton, N., [editors], Smithsonian Workshop on New Concepts in Global Tectonics:  Lubbock, Texas, Texas Tech University Press, p. 13-21.

Cronin, V.S., 1991, The cycloid relative-motion model and the kinematics of transform faulting, in Hilde, T.W.C., and Carlson, R.L., [editors], Proceedings of the 1987 Geodynamics Symposium:  Tectonophysics, v. 187, p. 215-249;  corrections to publisher errors, 1991, Tectonophysics, v. 192, no. 3/4, p. 401.

Cronin, V.S., 1988, Cycloid tectonics -- a kinematic model of finite relative plate motion:  College Station, Texas A&M University, Department of Geology, dissertation, 118 p.

Cronin, V.S., 1987, Cycloid kinematics of relative plate motion:  Geology, v. 15, p. 1006-1009.

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., R.G. Gordon, D.F. Argus, and S. Stein, 1994, Effect of recent revisions to the geomagnetic reversal time scale on estimates of current plate motions:  Geophysical Research Letters, v. 21, p. 2191-2194.

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  http://www.geology.wisc.edu/~chuck/MORVEL/

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.

Drewes, H., 1998, Combination of VLBI, SLR, and GPS determined station velocities for actual plate kinematic and crustal deformation models, in Geodynamics, M. Feissel (Ed.):  IAG Symposia, Springer.

Drewes, H., and D. Angermann, 2001, The Actual Plate Kinematic and Crustal Deformation Model 2000 (APKIM2000) as a Geodetic Reference System:  AIG 2001 Scientific Assembly, Budapest, 2-8 Sept 2001.

Gripp, A.E., and R.G. Gordon, 1990, Current plate velocities relative to the hopspots incorporating the NUVEL-1 global plate motion model:  Geophysical Research Letters, v. 17, p. 1109-1112.

Gripp, A.E & Gordon R.G., 2002, Young tracks of hotspots and current plate velocities:  Geophysical Journal International, v. 150, 321-361. 

Kreemer, C., W.E. Holt, and A.J. Haines, 2003, An integrated global model of present-day plate motions and plate boundary deformation:  Geophysical Journal International, v. 154, 8-34.

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.

Prawirodirdjo, L., and Y. Bock, 2004, Instantaneous global plate motion model from 12 years of continuous GPS observations:  Journal of Geophysical Research, v. 109, B08405, doi:10.1029/2003JB002944.

Sella, G.F., T.H. Dixon, and A. Mao, 2002, REVEL -- A model for recent plate velocities from space geodesy:  Journal of Geophysical Research, v. 107, B4, doi:10.1029/2000JB000033. 

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.