Modeling-Circular-Finite-Motion.nb

Vince Cronin

Compiled on 28-29 March 2019, based on earlier work by the author.

Introduction

The purpose of this code is to rotate a point or a set of points on the surface of a unit sphere (representing Earth) around an arbitrary Euler pole.

Input Data

This code utilizes an Excel spreadsheet to provide some of the input data for this computation.  The initial or current location of the point or points that are to be rotated is contained in a spreadsheet.  Data that are input directly through input lines in this code include the angular speed in degrees per million years and the geographic coordinates of the Euler pole (pole of rotation) around which the points move in a positive (right-handed or anticlockwise) direction.

Establish the location data for the Euler pole and corresponding angular speed

The variables poleLat and poleLong are the latitude and longitude of the Euler pole, and omegaPole is the angular speed around that pole.  A positive value indicates a right-hand or anticlockwise rotation around the pole, as observed from the pole toward Earth’s center.

poleLat=60.5;

poleLong=-108.3;

omegaPole=0.750;

Establish the time interval over which this analysis will run and the associated time steps

Note that positive times indicate time in the future.  All times are expressed in millions of years (Myr).

beginTime=0;

endTime=20;

stepTime=5;

Import the Excel Data File

To modify this code so that it can be used for your dataset containing the latitude and longitude of one or more points, do the following.

1.  In the same folder or directory in which this notebook is located, create a blank Excel spreadsheet and call it “geog_coord_datafile.xlsx.”  In the first column of the new spreadsheet, enter the decimal latitudes of the points you would like to use in this test.  In the second column, enter the corresponding decimal longitudes.  For example, the top of the input-data spreadsheet might look like this:

;

2.  Enter the following as an input line (under the black box, below):

testDataIn = Import[];

3. Put your cursor between the square brackets in the Import[] statement and click to establish the insertion point we will need in the next step.

4. Go to the Insert menu, select File Path, navigate to the correct input data file and choose it, and the correct file path will be inserted at the cursor.  In this example, the Excel file is located on the desktop of Vince’s computer, and so Mathematica will insert the path "/Users/vince/Desktop/geog_coord_datafile.xlsx" between the brackets.  The path will be different for every different file, and on every different computer.

testDataIn = Import[“/Users/vince/Desktop/geog_coord_datafile.xlsx”];

You can make sure it is an input line by clicking on that line’s bracket on the far right edge of this window, going to the Format menu, selecting the Style submenu, and then choosing Input

testDataIn=Import["/Users/vince/Desktop/geog_coord_datafile.xlsx"];

This won' t work unless you follow directions and add the appropriate input line above this line.  Your input line should look something like the example:
testDataIn = Import["/Users/vince/Desktop/geog_coord_datafile.xlsx"];

The imported dataset will have to be modified so that it has the correct dimensions for this analysis, which is accomplished with the following code that defines the list eqData:

inData=Flatten[testDataIn,1];

dims=Dimensions[inData];

noRecords=dims[[1]];

User-Defined Functions

The function unitVect3D operates on a 3D vector (vect_) and outputs the Cartesian coordinates of the corresponding unit vector.

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

zRotation[w_,dT_]:={{Cos[(w dT)Degree],-Sin[(w dT)Degree],0},{Sin[(w dT)Degree],Cos[(w dT)Degree],0},{0,0,1}};

The function geog2Cart operates on two geographic coordinates (lat_ and long_) and outputs the Cartesian coordinates of the unit location vector to that point.

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

The function cart2Geog operates on a 3D unit location vector (vect_) and outputs the geographic coordinates of that point in latitude and longitude.

cart2Geog[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 function distAzFromVectors operates on a 3D unit location vector to a reference point (refPtVect_) and a 3D unit location vector to another point (otherPtVect_), and outputs the distance in km and azimuth in degrees from the reference point to the other point.

distAzFromVectors[refPtVect_,otherPtVect_]:=Module[{radEarth,circEarth,northPole,θ1,θ2,θ3,θ4,normVect1,normVect2,normVect3,distKm,bearing,answer},radEarth=6371.02;circEarth=2*π*radEarth;northPole={0,0,1};normVect1 = Cross[northPole,refPtVect];
normVect2 = Cross[otherPtVect,refPtVect];θ1=VectorAngle[refPtVect,otherPtVect];
distKm=(θ1/(2*π))*circEarth;θ2=(180/π)*(VectorAngle[normVect1,normVect2]);θ3=(180/π)*(VectorAngle[normVect1,otherPtVect]);normVect3=Cross[refPtVect,normVect1];θ4=(180/π)*(VectorAngle[normVect3,otherPtVect]);bearing=If[((θ3==90)∨(θ3==180)),If[(θ4<=90),0,180],If[(θ3>90),(360-θ2),θ2]];answer={distKm,bearing};answer];

The function circMotion operates on a 3D unit location vector (x_), a coordinate transformation matrix (m1_), an angular speed (w_) and a time interval (dT_) in Myr, and outputs the unit location vector to that point after it rotates around its axis.

circMotion[x_,m1_,w_,dT_]:=Module[{m2,m3,answer},m2=zRotation[w,dT];m3=Inverse[m1];answer=m3.m2.m1.x;answer];

The function spherePlot1 operates on a set of 3D unit location vectors of computed data (inData_) and two sets of 3D unit location vectors of colored markers (markers1_ and markers2_) to help establish the coordinate system(s), and outputs a graphic of a unit sphere on which the input points are plotted.

spherePlot1[inData_,markers1_,markers2_]:=Module[{g1,g2,g3,g4,answer},g1=Graphics3D[Sphere[{0,0,0},1],AspectRatio->1,BoxRatios->{1,1,1},PlotRange->All,PlotRangePadding->0.1,ColorOutput->GrayLevel,Lighting->"Neutral"];g2=ListPointPlot3D[inData,AspectRatio->1,BoxRatios->{1,1,1},PlotStyle->Blue,PlotRange->All,PlotRangePadding->0.1];g3=ListPointPlot3D[markers1,AspectRatio->1,BoxRatios->{1,1,1},PlotStyle->Red,PlotRange->All,PlotRangePadding->0.1];g4=ListPointPlot3D[markers2,AspectRatio->1,BoxRatios->{1,1,1},PlotStyle->Green,PlotRange->All,PlotRangePadding->0.1];answer=Show[g1,g2,g3,g4];answer];

Computation

Define the axes of the geographic coordinate system

xGeog={1,0,0};

yGeog={0,1,0};

zGeog={0,0,1};

Define the axes of the Euler pole coordinate system

poleEuler=geog2Cart[poleLat,poleLong];

zPole=poleEuler;

xPole=Cross[zPole,zGeog];

yPole=Cross[zPole,xPole];

Define the angular speed “ω”

ω=omegaPole;

Convert the input data from the Excel spreadsheet from geographic coordinates (latitude and longitude) into Cartesian coordinates in the geographic coordinate system

inDataCart=Table[geog2Cart[inData[[i,1]],inData[[i,2]]],{i,1,noRecords}];

Define the transformation matrix between the geographic coordinate system and the Euler pole coordinate system

j1=(({{Dot[xGeog,xPole], Dot[yGeog,xPole], Dot[zGeog,xPole]}, {Dot[xGeog,yPole], Dot[yGeog,yPole], Dot[zGeog,yPole]}, {Dot[xGeog,zPole], Dot[yGeog,zPole], Dot[zGeog,zPole]}}));

Rotate the points around the Euler pole by an angle equal to the angular speed multiplied by the time interval.

prelimResults=Table[circMotion[inDataCart[[i]],j1,ω,deltaT],{i,1,noRecords},{deltaT,beginTime,endTime,stepTime}];

Flatten the prior result(s) by one dimension.

results1=Flatten[prelimResults,1];

Check the dimensions of the results1 matrix.

dims2=Dimensions[results1];

Establish the number of unit location vectors in the results1 matrix.

noSets=dims2[[1]];

Convert the Cartesian unit location vectors the results1 matrix into geographic coordinates (latitude and longitude).

results2=Table[cart2Geog[results1[[i]]],{i,1,noSets}];

Create some unit location vectors (markers1 and ) to help orient us in the graphic output.  The red markers correspond to the poles associated with the X, Y, and Z axes of the geographic coordinate system (where the Z axis corresponds to the north pole along Earth’s spin axis), and the green marker corresponds to the Euler pole.

markers1={{1,0,0},{0,1,0},{0,0,1}};

markers2={{poleEuler[[1]],poleEuler[[2]],poleEuler[[3]]}};

Output

The following matrix (results1) contains the unit location vectors, expressed in the geographic Cartesian coordinate system, for all of the points whose location was computed using this code.  Each row of the matrix corresponds to a single point location, with the first column holding the X coordinate, the second has the Y coordinate, and the third has the Z coordinate.

MatrixForm[results1]

The following matrix (results2) contains the unit location vectors, expressed in geographic coordinates of latitude and longitude, for all of the points whose location was computed using this code. Positive latitudes are north latitudes, and positive longitudes are east longitudes.  Each row of the matrix corresponds to a single point location, with the first column holding the latitude and the second has the longitude.

MatrixForm[results2]

Graphic output

The following graphic (graphicResult) shows a unit sphere with blue dots corresponding to the rotated points, red dots showing the locations of the poles associated with the axes of the geographic Cartesian coordinate system, and a green dot showing the location of the Euler pole.  Within this Mathematica notebook (but not in any PDF file created from this notebook), the graphic can be rotated or its scale changed to improve the visibility of the points plotted .

graphicResult=spherePlot1[results1,markers1,markers2]

Export data

1.  In the same folder or directory in which this notebook is located, create a blank Excel spreadsheet and call it “Rutabega.xlsx”

2.  Enter the following, formatted as a Mathematica input line, or modify an existing input line (see the blue boxes below):

Export[,outFile1];

3. Put your cursor between the left square bracket and the comma in the Export[,outFile1] statement and click to establish the insertion point we will need in the next step.

4. Go to the Insert menu, select File Path, navigate to the correct input data file and choose it, and the correct file path will be inserted at the cursor.  In this example, the Excel file is located on the desktop of Vince’s computer, and so Mathematica will insert the path "/Users/vince/Desktop/Rutabega.xlsx" between the brackets.  The path will be different for every different file, and on every different computer.

Export[“/Users/vince/Desktop/Rutabega.xlsx”];

You can make sure it is an input line by clicking on that line’s bracket on the far right edge of this window, going to the Format menu, selecting the Style submenu, and then choosing Input

The files created using an Export function will be placed in the same directory on your computer as this notebook file.

IMPORTANT NOTE:  The user must supply a name for the output file that is different from existing file names, or else the existing files may be over-written.

Export["/Users/vince/Desktop/outCartCoords.xlsx", results1];

Export["/Users/vince/Desktop/outGeogCoords.xlsx", results2];

Export["/Users/vince/Desktop/graphic1.jpg", graphicResult];

References

Ahrens, T.J., editor, 1995, Global Earth physics-- a handbook of physical constants : Washington, D.C., American Geophysical Union, AGU Reference Shelf 1, 376 p., ISBN 0 - 87590 - 851 - 9.

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.

Cronin, V.S., 2017, A brief primer on plate kinematics, with Mathematica:  unpublished draft, http://croninprojects.org/Vince/PlateKinematics/KinematicsPrimer/index.htm

Davis, H.F., and Snider, A.D., 1987, Introduction to Vector Analysis [5 th edition] :  Boston, Massachusets, Allyn and Bacon, 365 p., ISBN 0 - 205 - 10263 - 8.

Herring, T., editor, 2009, Geodesy : Amsterdam, Elsevier, Treatise on Geophysics, v.3, 446 p., ISBN 978 - 0 - 444 - 53460 - 6.

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.