Mathematica code to find the unspecified nodal plane, given the other nodal plane from a focal mechanism solution

By Vince Cronin, August 30, 2010

Input Data

The input data is an Excel spreadsheet (.xls file) in which each record/row has three values/columns:  the dip azimuth, the dip angle, and the rake of the slip vector.  The rake is measured in the plane that contains the dip vector and the strike vector, and is measured from the strike vector that is 90° anti-clockwise from the dip direction (i.e., the reference strike).  A positive rake indicates a positive, anti-clockwise, upward rotation of the slip vector from the reference strike;  a negative rake indicates a clockwise (downward) rotation from the reference strike.
The user will need to ensure that the path to the input data file is correctly specified in the first input line,  For example, a correct specification for the file "rawData.xls" located on the desktop of Vince Cronin's office computer would look like this:

mydata = Import["/Users/vince/Desktop/rawData.xls"];

In the code immediately following this text, be certain that the path to the input data file is correctly specified.

mydata=Import["/Users/vince/Desktop/inputData.xls"];

matrixSize=Dimensions[mydata];

Functions

The function vectorNorm computes the length of the specified vector.

vectorNorm[x_]:=Sqrt[Dot[x,x]];

The function vectorAngle computes the angle between the two specified vectors,with output expressed in radian measure.

vectorAngle[a_,b_]:=ArcCos[Dot[a,b]/(vectorNorm[a]*vectorNorm[b])];

The function unitVector takes a vector of arbitrary length and returns the corresponding unit-length vector (i.e., vector whose length is one unit).

unitVector[x_]:={x[[1]]/vectorNorm[x],x[[2]]/vectorNorm[x],x[[3]]/vectorNorm[x]};

Explanation of local variables in the module "findAuxPlane."  Variable "a" is the reference strike azimuth.  Variable "b" is the unit vector toward the reference strike azimuth in a coordinate system in which the +Y axis is north, the +X axis is east, and the +Z axis is up.  Variable "c" is the unit dip vector.  Variable "d" is the vector normal to the input nodal plane, and is the slip vector for the auxilary plane.  Vector "d" points up if the input rake is a positive number, and down if the input rake is a negative number.  Variable "f" is the basis vector for the X axis of the geographic coordinate system, which points east.  Variable "g" is the basis vector for the Y axis of the geographic coordinate system, which points north.  Variable "h" is the basis vector for the Z axis of the geographic coordinate system, which points up.  Variable "i" is the basis vector for the X axis of the nodal plane coordinate system.  If the rake is a negative number, the strike vector is the X axis and the dip vector is the Y axis.  If the rake is a positive vector, the dip vector is the X axis and the strike vector is the Y axis.  Variable "j" is the basis vec  tor for the Y axis of the nodal plane coordinate system.  Variable "d" is the basis vector of the Z axis of the nodal plane coordinate system.  Variable "k" is the unit vector in the nodal plane coordinate system that is coincident with the slip vector on the (input) nodal plane.  Variable "m" is the transformation matrix from the nodal plane coordinate system to the geographic coordinate system.  Variable "rakeVector" is the unit location vector along the slip direction on the nodal plane, expressed in the geographic coordinate system.  Variable "n" is the dip vector of the auxilary plane, expressed in radian measure.  Variable "s" is the unit vector projection of the rake vector onto the X - Y (horizontal) plane of the geographic coordinate system.
Explanation of the result.  The first element is the trend of the dip vector of the auxillary plane, in degrees.  The second element is the dip angle of the auxillary plane, in degrees.  The third element is the trend of the slip vector of the auxillary plane, in degrees.

findAuxPlane[dipAz_, dipAng_, rake_] :=Module[{a,b,c,d,f,g,h,i,j,k,m,rakeVector,n,s,rakeTrend,auxStrike,auxStrikeRad,auxStrikeVector,auxDipAz,auxDipAngle,auxRake},a=If[((dipAz-90)>0),(dipAz-90),(270+dipAz)];b={Sin[a Degree],Cos[a Degree],0};c={Sin[dipAz Degree]Cos[dipAng Degree],Cos[dipAz Degree]Cos[dipAng Degree],-Sin[dipAng Degree]};d=If[(rake<0),(Cross[b,c]),(Cross[c,b])];f={1,0,0};g={0,1,0};h={0,0,1};i=If[(rake<0),b,c];j=If[(rake<0),c,b];k=If[(rake<0),{Cos[rake Degree],-Sin[rake Degree], 0},{-Sin[rake Degree],Cos[rake Degree], 0}];m=(({{Dot[f,i], Dot[f,j], Dot[f,d]}, {Dot[g,i], Dot[g,j], Dot[g,d]}, {Dot[h,i], Dot[h,j], Dot[h,d]}}));rakeVector=m.k;n=ArcCos[Abs[rakeVector[[3]]]];s=unitVector[{rakeVector[[1]],rakeVector[[2]],0}];rakeTrend=If[(rakeVector[[1]]<0),(360-(vectorAngle[s,g]*(180/π))),(vectorAngle[s,g]*(180/π))];auxStrike=If[(rake<0),(If[(rakeTrend>270),(rakeTrend-270),(rakeTrend+90)]),(If[(rakeTrend<90),(rakeTrend+270),(rakeTrend-90)])];auxStrikeRad=auxStrike*(π/180);auxStrikeVector={Sin[auxStrikeRad],Cos[auxStrikeRad],0};auxDipAz=N[IntegerPart[If[(auxStrike>270),(auxStrike-270),(auxStrike+90)]+0.5]];auxDipAngle=N[IntegerPart[((180 n)/π)+0.5]];auxRake=N[IntegerPart[If[(rake<0),(-1*vectorAngle[auxStrikeVector,d]*(180/π)),(vectorAngle[auxStrikeVector,d]*(180/π))]+0.5]];{auxDipAz,auxDipAngle,auxRake}];

Computation

The input data for each earthquake (i.e., each row of the input data file) is passed through the function findAuxPlane to compute the dip azimuth of the corresponding auxilary plane, the dip angle of the auxilary plane, and the rake of the slip vector on the auxilary plane.  The output is a table of data called "preliminary".

result=Table[findAuxPlane[mydata[[1,i,1]],mydata[[1,i,2]],mydata[[1,i,3]]],{i,matrixSize[[2]]}];

Output

The output data is an Excel (.xls) spreadsheet in which each record/row has three values/columns :  the dip azimuth of the auxilary plane, the dip angle of the auxilary plane, and the rake of the slip vector on the auxilary plane.  

The user will need to ensure that the path to the output data file is correctly specified in the input line that follows this text.  For example, a correct specification for the file "outputFile.xls" located on the desktop of Vince Cronin' s office computer would look like this :
  
  Export["/Users/vince/Desktop/outputFile.xls", result];

In the code immediately following this text, be certain that the path to the output data file is correctly specified.

Export["/Users/vince/Desktop/outputData.xls", result];