Cycloid finite motion of points on <plate> relative to <plate>

by <your name>

Code by Vince Cronin, revised 12 April 2017

Parts of this code that you will need to change are highlighted in yellow.  You will also need to add appropriate comments and explanations, as indicated by the text enclosed in angle brackets <like this> throughout the notebook, which are not highlighted in yellow.  Delete this cell from your final version of this code that you submit for credit.

1  Introduction

<Explain the purpose of this notebook, what plate boundary section you are working on, which two plates are associated with this boundary, and why you chose to work on the particular plate boundary.  Why is it interesting to you, and why might it be geologically/tectonically interesting?>

<text>

2  Input data

2.1  Angular velocity data for plates relative to NNR reference frame

2.1.1  Instantaneous plate motion in an NNR reference frame from Argus et al. (2011)

<Describe the following data from Argus and others (2011) indicating instantaneous plate motion in an NNR reference frame>

nnr=(({{Antarctica, an, 65.42, -118.11, 0.250}, {Australia, au, 33.86, 37.94, 0.632}, {Caribbean, ca, 35.20, -92.62, 0.286}, {Cocos, co, 26.93, -124.31, 1.198}, {Eurasia, eu, 48.85, -106.50, 0.223}, {NorthAmerica, na, -4.85, -80.64, 0.209}, {Nubia, nb, 47.68, -68.44, 0.292}, {Nazca, nz, 46.23, -101.06, 0.696}, {Pacific, pa, -63.58, 114.70, 0.651}, {SouthAmerica, sa, -22.62, -112.83, 0.109}, {Somalia, sm, 49.95, -84.52, 0.339}}));

2.1.2  Plates and plate boundary segment being studied

<text>

Figure 1.  Boundary segment between <insert plate name> and <insert plate name> from GeoMapApp (http://www.geomapapp.org), showing the approximate current location of the boundary (black polyline).  <insert the rest of the caption here>

In the code below, one plate is designated as plate A and the other is plate B.  (This allows the basic code to be used for different plate pairs without extensive revision.)  Plate A is the frame of reference from which the motion of plate B is observed.

Please identify your choice of plate A:  <insert name here -- example: Australia>.

Please identify your choice of plate B:  <insert name here -- example: Somalia>.

2.1.3  Define the full angular velocity vector eΩplate for each plate relative to an external reference frame (NNR)

<Explain the following...>

eΩA=nnr[[2,5]]*convert2Cart[nnr[[2,3]],nnr[[2,4]]];

<Explain the following...>

eΩB=nnr[[11,5]]*convert2Cart[nnr[[11,3]],nnr[[11,4]]];

2.2  Location data along the boundary

<Reference points along the current plate boundary (Bird, 2003;  GeoMapApp).  In any pair of coordinates, the first one is the latitude and the second is the longitude -- north latitude and east longitude are positive; south latitude and west longitude are negative.  Some example data for the ridge between the Australian and Somalian Plates are provided.>

inputBoundary={{-14.741,66.375},{-13.949,65.883},{-13.343,66.718},{-12.769,66.313},{-12.151,65.610},{-11.618,66.357}};

2.3  Model times

This code is written to produce graphical and numeric answers for evenly spaced model times, beginning at 0 Myr (today) and continuing in the future by some integral number of time steps.  So, for example, if you want to model a boundary in 5 Myr steps, you could begin at 0 Myr and go to +20 Myr (four steps in the future).

timeStep=5;

futureTimeSteps=3;

2.4  Choices that affect the graphic output

2.4.1  Define the length of the north arrow

<Explain the following...>

northArrowScale=0.1;

2.4.2  Choose the distance from Earth’s center from which the boundary is viewed

Variable viewDistance is the distance to the point from which the boundary is viewed looking toward the center of Earth, measured in earth radii.  The value assigned to this variable must be greater than 1.

viewDistance=1.1;

3  User-defined functions

<Supply a sentence or so of introductory text that answers the following question:  What is this section for? >

The user-defined function convert2Cart takes the latitude (lat) and longitude (long) of a point on Earth’s surface and computes the unit location vector to that point in a reference frame whose origin is Earth’s center, the intersection of the equator and the prime meridian (0° latitude, 0° longitude) is along the positive X axis, the positive Y axis passes through 0° latitude and 90°E longitude, and the north pole is along the positive Z axis.

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

<Supply text that answers the following question:  What does the following user-defined function do?>

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

<Supply text that answers the following question:  What does the following user-defined function do?>

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

The user-defined function northArrow takes the average of all the unit location vectors of points along the plate boundary (viewPoint) and a user-defined scale (nArrowScale), generally ranging from about 0.05 to 0.1, and finds a suitable endpoint for a vector pointing toward north in the output graphics.

northArrow[viewPoint_,nArrowScale_]:=Module[{north,center,normal,northward,answer},north={0,0,1};center={0,0,0};normal=Cross[north,viewPoint];northward=Cross[viewPoint,normal]*nArrowScale;answer={center,northward};answer];

<Supply text that answers the following question:  What does the following user-defined function do?>

cycloid[(*list of global variables*)(*obsΩ is the total velocity vector of the observer's plate in the external reference frame*)obsΩ_,(*movΩ is the total velocity vector of the moving plate in the external reference frame*)movΩ_,(*refPt is the unit location vector to the reference point whose finite motion is to be modeled*)refPt_,(*t is time in Myr*)t_]:=Module[{(*list of local variables*)a,b,c,d,e,f,g,h0,h1,h,i0,i1,i,j0,j1,j,k,l,m1, m2, m3, m4, m5,n,p,answer},(*geographic X axis -- equator at prime meridian*)a={1,0,0};(*geographic Y axis*)b={0,1,0};(*geographic Z axis -- north pole*)c={0,0,1};(*eωobs*)d=Norm[obsΩ];(*pole eΩobs -- observer coordinate system (CS) Z axis*)e={(obsΩ[[1]]/d),(obsΩ[[2]]/d),(obsΩ[[3]]/d)};(*eωmov*)f=Norm[movΩ];(*pole eΩmov -- other plate CS Z axis*)g={(movΩ[[1]]/f),(movΩ[[2]]/f),(movΩ[[3]]/f)};(*vector coincident with observer CS Y axis*)h0=e\[Cross]g;(*length of h0*)h1=Norm[h0];
(*unit vector along observer CS Y axis*)h={(h0[[1]]/h1),(h0[[2]]/h1),(h0[[3]]/h1)};(*vector coincident with observer CS X axis*)i0=h\[Cross]e;(*length of i0*)i1=Norm[i0];
(*observer CS X axis*)i={(i0[[1]]/i1),(i0[[2]]/i1),(i0[[3]]/i1)};(*vector coincident with other plate CS X axis*)j0=h\[Cross]g;(*length of j0*)j1=Norm[j0];
(*other plate CS X axis*)j={(j0[[1]]/j1),(j0[[2]]/j1),(j0[[3]]/j1)};(*other plate CS Y axis -- same as observer CS Y axis*)k=h;(*m1 is the 3x3 transformation matrix from the geographic CS to the other plate CS*)m1={{j[[1]],j[[2]],j[[3]]},{k[[1]],k[[2]],k[[3]]},{g[[1]],g[[2]],g[[3]]}};(*m2 is a 3x3 Z-rotation matrix for a rotation of (-eωmov*t) degrees around pole eΩmov*)m2={{Cos[(-f Degree) t],Sin[(-f Degree) t],0},{-Sin[(-f Degree) t],Cos[(-f Degree) t],0},{0,0,1}};(*k is the angle between d and j*)l=VectorAngle[e,j];(*n is a variable that helps to differentiate between supplimentary angles*)n=If[l<(π/2),1,-1];(*p is the angle between f and d that has the proper sign to be used in matrix m3*)p=(n)(VectorAngle[g,e]);(*m3 is a 3x3 Y-rotation matrix for a rotation of p radians around the common Y axis of the observer CS and the other plate CS*)m3={{Cos[p],0,-Sin[p]},{0,1,0},{Sin[p],0,Cos[p]}};(*m4 is a 3x3 Z-rotation matrix for a rotation of (eωobs*t) degrees around pole eΩobs*)m4={{Cos[(d Degree) t],Sin[(d Degree) t],0},{-Sin[(d Degree) t],Cos[(d Degree) t],0},{0,0,1}};(*m5 is the 3x3 transformation matrix from the observer CS to the geographic CS*)m5={{i[[1]],h[[1]],e[[1]]},{i[[2]],h[[2]],e[[2]]},{i[[3]],h[[3]],e[[3]]}};(*answer is the product of multiplying the input location vector to a point on the other plate (refPt) by matrices m1, m2, m3, m4, and m5, resulting in the location vector to the point in "t" Myr in a version of the geographic CS that is fixed to the observer's plate*)answer=m5.m4.m3.m2.m1.refPt;answer];

4  Computation

4.1  Convert geographic location data along plate boundary to unit location vectors

Variable refPts is the set of all points along the boundary of interest, expressed as unit location vectors.  Each point is converted from the input geographic coordinates to Cartesian coordinates using the function convert2Cart.

refPts=Table[convert2Cart[inputBoundary[[i,1]],inputBoundary[[i,2]]],{i,1,Length[inputBoundary]}];

4.3  Define the full angular velocity vector observing.plateΩmoving.plate for one plate relative to the other plate

<Explain the following...>

plateAΩplateB=eΩB-eΩA;

<Explain the following...>

plateAωplateB=Norm[plateAΩplateB];

<Explain the following...>

cartRelPole={plateAΩplateB[[1]]/plateAωplateB,plateAΩplateB[[2]]/plateAωplateB,plateAΩplateB[[3]]/plateAωplateB};

<Explain the following...>

geogRelPole=findGeogCoord[cartRelPole];

4.4  Model time

<Explain the following...>

futureTime=futureTimeSteps*timeStep;

<Explain the following...>

modelTime=Table[i,{i,0,futureTime,timeStep}];

4.5  Compute boundary locations

<Explain the following...>

cycOut = Table[cycloid[eΩA, eΩB, refPts[[i]], modelTime[[j]]], {j,1,Length[modelTime]},{i, 1, Length[refPts]}];

4.7  Prepare the graphics files

4.7.1  Define the view vector

Variable averageVV is the point from which the boundary is viewed looking toward the center of Earth, as plotted using the built-in Mathematica function Graphics3D

averageVV=viewDistance*(Sum[refPts[[i]],{i,Length[refPts]}]/Norm[Sum[refPts[[i]],{i,Length[refPts]}]]);

4.7.2  Compile the graphics file

Note:  Some useful (easily visible) built-in colors for use in the graphics include the following:  Red, Green, Blue, Black, Gray, Cyan, Magenta, Brown, Orange, Pink, and Purple.

<Explain the following...>

boundaryOut = Graphics3D[{Arrow[northArrow[averageVV,northArrowScale]],Black,Line[Table[cycOut[[1,i]],{i,1,Length[refPts]}]],Red, Line[Table[cycOut[[j,i]],{j,2,Length[cycOut]},{i,1,Length[refPts]}]]}, AspectRatio->Automatic, Boxed->False,ViewVector->{averageVV,{0,0,0}}];

trajectoryOut = Graphics3D[{Arrow[northArrow[averageVV,northArrowScale]],Black,Line[Table[cycOut[[1,i]],{i,1,Length[refPts]}]],Orange, Line[Table[cycOut[[j,i]],{i,1,Length[refPts]},{j,1,Length[cycOut]}]],Blue,Dashed,Line[Table[cycOut[[Length[cycOut],i]],{i,1,Length[refPts]}]]}, AspectRatio->Automatic, Boxed->False,ViewVector->{averageVV,{0,0,0}}];

5  Output

<Explain...>

Show[boundaryOut]

Figure 2.  Present boundary (solid black polyline) and modeled location of those same points on the <name of plate B> plate as observed from the <name of plate A> plate at <timeStep> Myr intervals into the future (solid red polylines).    Black arrow indicates direction toward north.  <insert the rest of the caption here>

<Explain...>

Show[trajectoryOut]

Figure 3.  Present boundary (solid black polyline) and modeled trajectories of those same points on the <name of plate B> plate as observed from the <name of plate A> plate at <timeStep> Myr intervals into the future (solid orange polylines) to their modeled positions at <finalTime> Myr (dashed blue polyline).    Black arrow indicates direction toward north.  <insert the rest of the caption here>

Check the relative pole and speed between the two plates.

Compare the computed latitude, longitude, and angular speed of these plates relative to each other today (given below) with those published in DeMets et al (2010).

latRelPole=geogRelPole[[1]]

longRelPole=geogRelPole[[2]]

angularSpeed=plateAωplateB

6  Discussion, Conclusions

<text>

7  References and data sources

Argus, D.F., Gordon, R.G., and DeMets, C., 2011, Geologically current motion of 56 plates relative to the no-net-rotation reference frame:  Geochemistry Geophysics Geosystems, v. 12, Q 1101, 13 p., doi:  10.1029/2011GC003751.

Bird, P., 2003, An updated digital model of plate boundaries:  Geochemistry Geophysics Geosystems, v. 4, no. 3, p. 1027, doi:10.1029/2001GC000252.  Available online via http://peterbird.name/publications/2003_PB2002/2003_PB2002.htm

Cronin, V.S., 2017, A primer on the kinematics of the lithosphere and continental crust [draft]:  accessible on the web via http://CroninProjects.org/Vince/PlateKinematics/KinematicsPrimer/index.htm

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, Ph.D. dissertation, 118 p.

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

DeMets, D. C., Gordon, R.G., and Argus, D.F., 2010, Geologically current plate motions: Geophysical Journal International, v. 181, p. 1–80, doi:10.1111/j.1365-246X.2010.04491.x.

Some web resources to investigate:

http://geoscience.wisc.edu/~chuck/MORVEL/

http://www.virtualocean.org/

http://www.geomapapp.org/

http://worldwind.arc.nasa.gov/java/