Cycloids.nb

An implementation of the cycloid model for finite relative plate motion

Vince Cronin, begun 30 March 2021, last revised 1 April 2021

Input

Import point location data

The location data for the points whose relative motion will be modeled is in an external file.  Each record has two numbers separated by a comma:  decimal latitude in °N and decimal longitude in °E.  The number of records is the same as the number of points on Earth’s surface that are referenced in the imported dataset.  These particular data are from Bird (2003).

The variable points is defined as the list of latitude-longitude data contained in the imported file.

points=Import["/Users/vince/Desktop/IcelandBoundaryBird2003.csv"];

For more information about importing and exporting of data, review the notebook “Exporting-Importing-and-Referencing-Data.nb,” accessible at https://CroninProjects/Vince/PlateKinematics/KinematicsCourse/Exporting-Importing-and-Referencing-Data.nb

Data from GSRM-NNR-2.1 (Kreemer et al., 2014)

The variable nnr is defined as a matrix of data, originally published by Kreemer, Blewitt, and Klein (2014) in their supplementary table 2, with the numerical data from the three columns under “pole NNR (deg/Myr)” column.  In order from left to right, the columns of the nnr matrix contain the plate name, plate abbreviation, pole latitude (°N), pole longitude (°E), and angular speed (°/Myr).

nnr=(({{Africa, af, 49.66, -78.08, 0.285}, {Antarctica, an, 60.08, -120.14, 0.234}, {Caribbean, ca, 37.84, -96.49, 0.290}, {Cocos, co, 27.21, -124.02, 1.169}, {Eurasia, eu, 55.38, -95.41, 0.271}, {JuandeFuca, jf, -37.71, 59.44, 0.977}, {NorthAmerica, na, 2.19, -83.75, 0.219}, {Okhotsk, ok, 28.80, -90.91, 0.209}, {Pacific, pa, -63.09, 109.63, 0.663}, {Rivera, ri, 20.27, -107.10, 4.510}, {SouthAmerica, sa, -14.10, -117.86, 0.123}}));

The element in the third row, fourth column can be accessed using the following reference:  nnr[[3,4]], as in the following:

nnr[[3,4]]

User - defined functions and modules

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

findGeogCoord[vect_]:=Module[{lat,long,a,b,θ},a=ArcSin[vect[[3]]];θ=ArcCos[vect[[1]]/(Sqrt[(((vect[[1]])^2)+((vect[[2]])^2))])];b=If[(vect[[2]]<0),(-θ),(θ)];lat=a(180/π);long=If[((Abs[vect[[1]]]<(1×(10^(-14))))&&(Abs[vect[[2]]]<(1×(10^(-14))))),0,(b(180/π))];{lat,long}];

The following user-defined module cycloid[] acts on the the following input arguments:  the angular velocity vector of the reference plate, the angular velocity vector of the observed plate, the location vector of a point on the observed plate, and the time relative to present where negative values are past time and positive values are future time.  The output will be the unit location vector of the rotated point, which can then be converted back to geographic coordinates using the function findGeogCoord[].

cycloid[refΩ_,obsΩ_,obsPt_,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=Normalize[refΩ];e=Norm[refΩ];f=Normalize[obsΩ];g=Norm[obsΩ];h=Normalize[d\[Cross]f];i=h;j=Normalize[h\[Cross]d];k=Normalize[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.obsPt;answer];

Computation

Angular Velocity Vectors

The variables defined below are angular velocity vectors in an Earth-centered Cartesian geographic coordinate system whose magnitudes are expressed in units inherited from the plate-motion model that is used -- typically °/Myr.  The variable names have a leading letter “e” that indicates that the reference frame is effectively external to the lithosphere (e.g., hotspots, no-net-rotation).  The Greek letter “ω” is a typical symbol used for angular velocity in this application.  The last two letters reflect the plate whose angular velocity is being expressed relative to a plate-external reference frame.

eωaf=nnr[[1,5]]*makeLocVector[nnr[[1,3]],nnr[[1,4]]];

eωan=nnr[[2,5]]*makeLocVector[nnr[[2,3]],nnr[[2,4]]];

eωca=nnr[[3,5]]*makeLocVector[nnr[[3,3]],nnr[[3,4]]];

eωco=nnr[[4,5]]*makeLocVector[nnr[[4,3]],nnr[[4,4]]];

eωeu=nnr[[5,5]]*makeLocVector[nnr[[5,3]],nnr[[5,4]]];

eωjf=nnr[[6,5]]*makeLocVector[nnr[[6,3]],nnr[[6,4]]];

eωna=nnr[[7,5]]*makeLocVector[nnr[[7,3]],nnr[[7,4]]];

eωok=nnr[[8,5]]*makeLocVector[nnr[[8,3]],nnr[[8,4]]];

eωpa=nnr[[9,5]]*makeLocVector[nnr[[9,3]],nnr[[9,4]]];

eωri=nnr[[10,5]]*makeLocVector[nnr[[10,3]],nnr[[10,4]]];

eωsa=nnr[[11,5]]*makeLocVector[nnr[[11,3]],nnr[[11,4]]];

Convert points from geographic coordinates (lat, long) to Cartesian geographic coordinates

The following code line is for user convenience, displaying the imported data represented by the variable points.  The values in the first column are latitudes (°N) while the second column are longitudes (°E).

MatrixForm[N[points]]

The variable ptsCart is the matrix of Cartesian geographic coordinates of the imported points

ptsCart=Table[makeLocVector[points[[i,1]],points[[i,2]]],{i,1,Length[points],1}];

The variable displacedPointsCart is the matrix of Cartesian geographic coordinates of the imported points after they have been displaced by the time steps indicated by the local dummy variable j, which ranges from 0 (present time) to 5 Myr in steps of 1 Myr.  The other local dummy variable, i, indicates the row of the matrix ptsCart that contributes data to the cycloid module at any given step of the process.  Dummy variable i ranges from 1 (the first record in the dataset) to the last record in the dataset (determined using the built-in Mathematica function Length[] acting on the matrix ptsCart) in steps of 1.

displacedPointsCart=Table[cycloid[eωna,eωeu,ptsCart[[i]],j],{j,0,5,1},{i,1,Length[ptsCart],1}];

The variable dimDisplPtsCart contains the dimensions of the matrix displacedPointsCart.  The first value is the number of time steps, the second value is the number of rows in the imported matrix (the number of points whose locations are given in the imported matrix), and the third value is the number of elements in each of the individual results -- the dimensions of each of the resulting location vectors.

dimDisplPtsCart=Dimensions[displacedPointsCart]

The variable resultsCart is used only for display purposes, so that the results for all points at a given time step are shown in the same column.  Hence, the number of columns is the number of time steps.

resultsCart=Transpose[displacedPointsCart];

The variable displacedPointsGeog is the matrix of geographic coordinates of the imported points after they have been displaced by the time steps defined in the input line that defines the variable displacedPointsCart.  The local dummy variable i indicates the row in displacedPointsCart where the target vector is located, and local dummy variable j indicates the column where the target vector is located.  Dummy variable i ranges from 1 (the first row in the dataset) to the last row in the dataset (indicated by the first value of dimDisplPtsCart) in steps of 1, where each row is associated with a specific time step. Dummy variable j ranges from 1 (the first column in the dataset) to the last column in the dataset (indicated by the second value of dimDisplPtsCart) in steps of 1, where each column is associated with a specific point in the imported set of points.

displacedPointsGeog=Table[findGeogCoord[displacedPointsCart[[i,j]]],{i,1,dimDisplPtsCart[[1]],1},{j,1,dimDisplPtsCart[[2]],1}];

The variable resultsGeog is used only for display purposes, so that the results for all points at a given time step are shown in the same column.  Hence, the number of columns is the number of time steps.

resultsGeog=Transpose[displacedPointsGeog];

Output Results

Results Expressed as Unit Location Vectors in Cartesian Geographic Coordinates

The result of using the module cycloid[], acting on input data expressed in geographic coordinates, is given below.  The first number in the set is the decimal latitude (°N) and the second is the decimal longitude (°E).  Each row is a set of 3-element lists (unit location vectors in the Cartesian geographic coordinate system) corresponding to different positions at different times (the time series) for a specific point in the original imported dataset.  Each column is a set of vectors that define the positions of all the points in the imported dataset at a specific time in the time series, with the earliest time in the first column and the latest time in the last column.  So the number of columns of vectors is the number of time steps, and the number of rows of vectors is the number of points we are tracking through time.

MatrixForm[N[resultsCart]]

Results Expressed in Geographic Coordinates (Latitude & Longitude)

The structure of the matrix below is the same as the structure of the preceding matrix, except that the location data are presented in geographic coordinates in which the first number is the latitude and the second number is the longitude.

MatrixForm[N[resultsGeog]]

Export["/Users/vince/Desktop/IcelandResultsCart.csv",resultsCart];

Export["/Users/vince/Desktop/IcelandResultsGeog.csv",resultsGeog];

References

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

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.