startTime=AbsoluteTime[];

SLAM-code-2017.nb

This code is copyright © 2004-2017 by Vincent S. Cronin.

Begun August 5, 2004;  Revised September 18, 2017

Introduction

This version of the SLAM code is meant to project six planes to the ground surface and to map the intersection of those planes and the ground surface.

Prepare the Earthquake Data File

The eqDataIn is an Excel (.xls) spreadsheet that contains the earthquake data necessary for the analysis.  The first record (row) of the spreadsheet contains words that describe the contents of each column, and each other record (row) contains the data for one earthquake focal mechanism.  The data types are as follows, starting with the left-most column (column A): year, month, day, hour, minute, seconds, decimal latitude (°), decimal longitude (°), depth (km), greatest horizontal error (EH1; km), least horizontal error (EH2; km), azimuth of EH1 axis (°), vertical error (km), dip direction of the fault plane (°), dip angle of  the fault plane (°), rake of hanging-wall slip vector on the fault plane (°), strike uncertainty (°), dip uncertainty (°), rake uncertainty (°), and optional other columns containing ID or other information.

Input

Import the Earthquake Data File

To modify this code so that it can be used for your dataset, do the following.

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

eqDataIn = Import[];

2. 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.

3. 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/Input_EQ_datafile.xls" between the brackets.  The path will be different for every different file, and on every different computer.

eqDataIn = Import[“/Users/vince/Desktop/Input_EQ_datafile.xls”];

Add an input line below this line.  Type  eqDataIn = Import[];  and follow the rest of the instructions above.  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

eqDataIn=Import["/Users/vince/Desktop/Input_EQ_datafile.xls"];

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:
eqDataIn = Import["/Users/vince/Desktop/Input_EQ_datafile.xls"];

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:

eqData=Flatten[eqDataIn,1];

Import the digital elevation model (DEM)

The ground surface is represented by a digital elevation model.  The DEM file is an ASCII text file (that uses the generic data  tag “.dat”) that has a 6 line 2 column header followed by a matrix of many rows and columns of elevation data.  The geographic grid used is a UTM grid (meters).  The horizontal position of any given elevation datum is provided by its position in the matrix, relative to the UTM coordinates of the lower left corner and the cellsize -- the horizontal distance between adjacent matrix nodes.  The gridded DEM datafile looks like the following:

In a DEM that has no elements that have a null datapoint equal to the “NODATA_value,” we use the first position in the data grid below the header to condition the dataset so that it displays a hillshade image with sufficient contrast.  The DEM sampled above has a zero in the first grid position.

mydata=Import["/Users/vince/Desktop/thinDVFZCrop.dat"];

Set various operating parameters for the code

The zoneMeridian is the longitude at the center of the UTM zone that contains the epicentral region.  If the epicenter is near a zone boundary, base the UTM grid for this analysis in the zone that is to the west.

zoneMeridian=-123;

The gridNorthAdjustment is the angle between true north and grid north at/near the center of the map area.  If grid north is found by a clockwise rotation from true north, the sign of the gridNorthAdjustment is negative;  otherwise, it is positive.  The gridNorthAdjustment is approximately equal to (zoneMeridian – (longitude at center of map area))*Sin[latitude at center of map area].

gridNorthAdjustment=(zoneMeridian-(-120.15))*Sin[37.47Degree]

The epicenterWidth is the half-width of the square that will be put over the top of the epicenter on the maps generated by this code, in meters

epicenterWidth=250;

The widthFactor is a multiplier that allows for the adjustment of the half-width of the zone between parallel planes that generate the trace of the mean nodal plane on the maps generated by this code, in meters.  A value of 1.5 for the widthFactor sets the distance between the planes at 1.5 times the cell size of the DEM.

widthFactor=1.5;

The topoBase value forces one point of the DEM to have a particular value, which in turn allows the user to adjust the tonal range of the hillshade map.

topoBase=0;

The constant minusOne is used as input for the pointEvaluator module, in the place for the null value.

minusOne=-1;

The meanMultiplier increases the width of the zone between parallel planes around the mean nodal plane, and is used to simplify interpretation of the seismo-lineament.  Care must be taken so that the meanMultiplier value does not become so large that the ground-surface trace of the resulting planes extends beyond the boundaries of the seismo-lineament..

meanMultiplier=1;

Input for the first earthquake in the set

counter=2;

Remember to change the names of the output files

Evaluate the Notebook (Run the Program)

Select the Evaluation drop-down menu from the horizontal menu bar at the top of the window/screen, and choose Evaluate Notebook.

Some user-defined functions

makeVector[plunge_,trend_]:={Cos[plunge Degree] Sin[trend Degree],Cos[plunge Degree] Cos[trend Degree],-Sin[plunge Degree]};

vectorNorm[x_]:=(Sqrt[x.x]);

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

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

The lat/long to UTM conversion is after Snyder (1982), and assumes use of the North American Datum of 1927 (NAD27).  Another excellent resource for converting from lat/long to UTM and back, using a variety of datums (e.g., NAD27, WGS84 and so on) is provided by Steven Dutch at http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm, and by an online calculator at http://www.rcn.montana.edu/resources/tools/coordinates.aspx

convertToUTM[inLat_,inLong_,centMerid_]:=Module[{c1,c2,c3,c4,c5,v1,v2,v3,v4,v5,v6,utmX,utmY},c1=6378206.4;c2=0.00676866;c3=0;c4=centMerid;c5=0.9996;v1=c2/(1-c2);v2=c1/Sqrt[1-(c2*(((Sin[inLat Degree])^2)))];v3=((Tan[inLat Degree])^2);v4=v1*(((Cos[inLat Degree])^2));v5=(Cos[inLat Degree])*((inLong-c4) ((π/180)));v6=(111132.0894*inLat)-(16216.94*Sin[2*(inLat Degree)])+(17.21*Sin[4*(inLat Degree)])-(0.02*Sin[6*(inLat Degree)]);utmX=(c5*v2*(v5+(((1-v3+v4)*(v5^3))/6)+((5-(18*v3)+((v3^2))+(72*v4)-(58*v1))*(v5^5))/120))+500000;utmY=c5*(v6 - 0+(v2*Tan[inLat Degree]*(((v5^2)/2)+(((5-v3+(9*v4)+(4*((v4^2))))*(v5^4))/24)+(((61-(58*v3)+((v3^2))+(600*v4)-(330*v1))*(v5^6))/720))));{utmX,utmY}];

The following module differentiates between points that are within "width" meters from the fault plane and those that are further away.  Given unit vector N that is normal to the fault plane that passes through the origin of the coordinate system, the distance from an arbitrary point (whose position vector is P) to that plane is given by |N·P|.

pointEvaluator[xCoord_,yCoord_,zCoord_,width_,fltUNrml_,nulData_]:=Module[{locVect,distToFlt,result},locVect={xCoord,yCoord,zCoord};distToFlt=Abs[Dot[fltUNrml,locVect]];result=If[(distToFlt≤width),10.,nulData];result];

Text

epicenterEval[xCoord_,yCoord_,zCoord_,width_,nulData_]:=Module[{locVect,distToEpicenter,result},result=If[(Abs[xCoord]<width),If[(Abs[yCoord]<width),zCoord,nulData],nulData];result];

The module findHalfWidth is used to determine the half-width of the uncertainty envelope given the lengths of the semi-axes of the triaxial uncertainty ellipsoid around the hypocenter, the azimuth of the greatest horizontal semi-axis, and the dip azimuth and dip angle of a nodal plane.  The uncertainty envelope is bounded by parallel planes that are tangent to the uncertainty ellipsoid.
Input Parameters
•  dAz is the trend or azimuth of the dip vector of the nodal plane, measured in degrees clockwise from north.  All input azimuths are measured relative to the geographic coordinate system.  Hence, an azimuth of 90° corresponds to due east.
•  dAng is the plunge of the dip vector of the nodal plane, measured in degrees down from horizontal.  Hence, a dip angle of 90° corresponds to a vertical vector pointing down.
•  hEr1 is the length of the a semi-axis of the 95% CI uncertainty ellipsoid surrounding the mean hypocenter, measured in kilometers.  Hence, a value of 1 for hEr1 means that the surface of the 95% CI uncertainty ellipse is 1 km from the mean hypocenter measured horizontally in the direction hEr1Az.
•  hEr1Az is the trend of the e1h uncertainty axis (i.e., of the a semi-axis of the 95% CI uncertainty ellipsoid surrounding the mean hypocenter), measured in degrees clockwise from north.  All input azimuths are measured relative to the geographic coordinate system.  Hence, an azimuth of 90° corresponds to due east.
•  hEr2 is the length of the b semi-axis of the 95% CI uncertainty ellipsoid surrounding the mean hypocenter, measured in kilometers.  Hence, a value of 1 for hEr2 means that the surface of the 95% CI uncertainty ellipse is 1 km from the mean hypocenter measured horizontally 90° from the direction hEr1Az.
•  vEz is the length of the c semi-axis of the 95% CI uncertainty ellipsoid surrounding the mean hypocenter, measured in kilometers.  Hence, a value of 1 for vEz means that the surface of the 95% CI uncertainty ellipse is 1 km from the mean hypocenter measured vertically.
Local Variables
•  var1 is the angle, in degrees, through which the geographic coordinate system (in which north is {0,1,0} and east is {1,0,0}) is transformed by rotation around the common z axis to the coordinate systen that is fixed to the axes of the triaxial uncertainty ellipsoid around the hypocenter (in which the x axis coincides with the a ellipsoidal axis and with hEr1, the y axis coincides with the b ellipsoidal axis and with hEr2, and the z axis is vertical and coincides with the c ellipsoidal axis and with vEr).
•  var2 is the unit vector parallel to the dip vector, in the geographic coordinate system in which north is {0,1,0} and east is {1,0,0}.
•  var3 is the unit vector parallel to the right - hand - rule strike, in the geographic coordinate system in which north is {0,1,0} and east is {1,0,0}.
•  var4 is the 3x3 rotation matrix in which a positive angle results in an anticlockwise rotation around the positive z axis.
•  var5 is the unit vector that is normal to the nodal plane in the geographic coordinate system in which north is {0,1,0} and east is {1,0,0}.
•  var6 is the unit vector that is normal to the nodal plane in the ellipsoid-centered coordinate system that is rotated var1 degrees around the vertical (z) axis.
•  var7-9 are intermediate values that are ultimately used to compute the value of the Lagrangian multiplier  lambda.
•  var10 is the Lagrangian multiplier lambda
•  var11 is the set of components of the location vector to one of the tangent points, expressed in the ellipsoid-centered coordinate system that is rotated var1 degrees around the vertical (z) axis.
•  var12 is the distance, in km, between the a nodal plane through the hypocenter and a parallel nodal plane through one of the tangent points to the triaxial uncertainty ellipsoid around the hypocenter.
Output
•  result is the distance, in km, between the a nodal plane through the hypocenter and a parallel nodal plane through one of the tangent points to the triaxial uncertainty ellipsoid around the hypocenter..

findHalfWidth[dipAz_,dipAng_,hEr1_,hEr1Az_,hEr2_,vEr_]:=Module[{var1,var2,var3,var4,var5,var6,var7,var8,var9,var10,var11,var12,result},var1=hEr1Az-90;var2={Cos[dipAng Degree]Sin[dipAz Degree],Cos[dipAng Degree]Cos[dipAz Degree],-Sin[dipAng Degree]};var3={Sin[(dipAz-90)Degree],Cos[(dipAz-90)Degree],0};var4={{Cos[var1 Degree],-Sin[var1 Degree],0},{Sin[var1 Degree],Cos[var1 Degree],0},{0,0,1}};var5=Cross[var2,var3];var6=var4.var5;var7=(1/(hEr1^2))(((var6[[1]]((hEr1^2)/2)))^2);var8=(1/(hEr2^2))(((var6[[2]]((hEr2^2)/2)))^2);var9=(1/(vEr^2))(((var6[[3]]((vEr^2)/2)))^2);var10=(Sqrt[(1/(var7+var8+var9))]);var11={var10*var6[[1]]*((hEr1^2)/2),var10*var6[[2]]*((hEr2^2)/2),var10*var6[[3]]*((vEr^2)/2)};var12=var6.var11;result=var12];

Computations specific to an individual hypocenter, for which there might be multiple focal mechanism solutions available

Read and interpret the DEM header information

headerData=Table[mydata[[i,j]],{i,6},{j,2}];

ncols=headerData[[1,2]];

nrows=headerData[[2,2]];

xllcorner=headerData[[3,2]];

yllcorner=headerData[[4,2]];

cellsize=headerData[[5,2]];

nodataValue=headerData[[6,2]];

Read and interpret the hypocenter location data

record=counter;

focalLat=eqData[[record,7]];

focalLong=eqData[[record,8]];

focalDepthKm=eqData[[record,9]];

The UTM coordinates that specify the horizontal position of the earthquake focus (i.e., the epicenter location) are given in meters.  The variable focalDepth gives the depth of the earthquake focus, in meters.

focalDepth=focalDepthKm*(-1000);

The reported vertical and horizontal uncertainties in the location of the earthquake focus are given in kilometers.  The variables horizError and vertError convert these to meters.

eh1=eqData[[record,10]]*1000;

eh2=eqData[[record,11]]*1000;

eh1Azimuth=eqData[[record,12]];

eh1Az=If[((eh1Azimuth+gridNorthAdjustment)<0),(360+eh1Azimuth+gridNorthAdjustment),If[((eh1Azimuth+gridNorthAdjustment)>360),((eh1Azimuth+gridNorthAdjustment)-360),(eh1Azimuth+gridNorthAdjustment)]];

ez=eqData[[record,13]]*1000;

Convert the input coordinates of the earthquake focus to the same coordinate system as the DEM data

utmCoordinates=convertToUTM[focalLat,focalLong,zoneMeridian];

The "focus" is the (x, y, z) coordinates in meters of the reported earthquake focus, in a coordinate system in which the origin is the lower-left (southwest) corner of the DEM at sea level.

focusMean={(utmCoordinates[[1]]-xllcorner),(utmCoordinates[[2]]-yllcorner),focalDepth};

demImageFileM = Table[mydata[[i + 6, j]] - focusMean[[3]], {i, nrows}, {j, ncols}];

rawElevDataM = Table[mydata[[i + 6, j]] - focusMean[[3]], {i, nrows}, {j, ncols}];

demImageFileM = Table[rawElevDataM[[i, j]], {i, nrows}, {j, ncols}];

Computations specific to the selected nodal plane

npDipTrend=eqData[[record,15]];

npDipTr=If[((npDipTrend+gridNorthAdjustment)<0),(360+npDipTrend+gridNorthAdjustment),If[((npDipTrend+gridNorthAdjustment)>360),((npDipTrend+gridNorthAdjustment)-360),(npDipTrend+gridNorthAdjustment)]];

npDipPlunge=eqData[[record,16]];

npDipTrendUncert=eqData[[record,18]];

npDipAngUncert=eqData[[record,19]];

Make the hillshade image

rawElevDataInput = Table[mydata[[i + 6, j]] , {i, nrows}, {j, ncols}];

elev3 = Table[demImageFileM[[i, j]], {i, nrows, 1, -1}, {j, ncols}];

npDipTrRad=npDipTr*(π/180);

lightDirection = If[npDipTrRad>(3π/2),(2π)-npDipTrRad,(π/2)-npDipTrRad];

hillshadeMap1 = ReliefPlot[elev3, ColorFunction -> "GrayTones", LightingAngle -> {lightDirection, Pi/12}];

ClearAll[elev3];

Swath 1 of 7:  The middle of the road

Convert the fault dip vector from trend and plunge to a unit location vector

npDipVect=unitVector[makeVector[npDipPlunge,npDipTr]];

Find the strike vector defined using the right-hand rule

npStrike=If[(npDipTr<90),(npDipTr+270),(npDipTr-90)];

npStrikeVect = {Sin[npStrike Degree], Cos[npStrike Degree], 0};

Find the unit vector that is normal to the fault plane and is directed upwards

npNormalVect=unitVector[Cross[npDipVect,npStrikeVect]];

zoneHalfWidthMThin=widthFactor*cellsize*Sin[npDipPlunge Degree];

N[%]

uncertSwath=Table[pointEvaluator[((cellsize*(j-1))-focusMean[[1]]),((cellsize*(nrows-i))-focusMean[[2]]),demImageFileM[[i,j]],zoneHalfWidthMThin,npNormalVect,minusOne],{i,nrows},{j,ncols}];

elevData=Table[uncertSwath[[i,j]],{i,nrows,1,-1},{j,ncols}];

traceImgFileMThin = ListContourPlot[elevData, ContourShading -> False, AspectRatio -> Automatic, Contours -> {topoBase}];

elevDataTemp1=elevData;

ClearAll[elevData,uncertSwath];

Swath 2 of 7:  Gentle dip, smaller strike azimuth (GS)

Subtract the dip-angle uncertainty from the mean dip angle (i.e., from the plunge of the mean dip vector)

gentleDipAng=If[(npDipPlunge-npDipAngUncert)<0,Abs[(npDipPlunge-npDipAngUncert)],(npDipPlunge-npDipAngUncert)];

Subtract the dip-trend uncertainty from the mean dip azimuth (i.e., from the mean dip trend)

gentleDipAzS=If[(npDipPlunge-npDipAngUncert)<0,If[((npDipTr-npDipTrendUncert)>180),((npDipTr-npDipTrendUncert)-180),((npDipTr-npDipTrendUncert)+180)],If[((npDipTr-npDipTrendUncert)<0),360+(npDipTr-npDipTrendUncert),(npDipTr-npDipTrendUncert)]];

Convert the fault dip vector from trend and plunge to a unit location vector

dipVectGS=unitVector[makeVector[gentleDipAng,gentleDipAzS]];

Find the strike vector defined using the right-hand rule

strikeGS=If[(gentleDipAzS<90),(gentleDipAzS+270),(gentleDipAzS-90)];

strikeVectGS={Sin[strikeGS Degree],Cos[strikeGS Degree],0};

Find the unit vector that is normal to the fault plane and is directed upwards

normalVectGS=unitVector[Cross[dipVectGS,strikeVectGS]];

Find the width of the nodal-plane uncertainty zone with the given dip angle, dip azimuth, and triaxial hypocenter uncertainty ellipsoid

zoneHalfWidthGS = findHalfWidth[gentleDipAzS, gentleDipAng, eh1, eh1Az, eh2, ez];

N[zoneHalfWidthGS]

Find the boundaries of the uncertainty swath as it intersects the ground surface

uncertSwath=Table[pointEvaluator[((cellsize*(j-1))-focusMean[[1]]),((cellsize*(nrows-i))-focusMean[[2]]),demImageFileM[[i,j]],zoneHalfWidthGS,normalVectGS,minusOne],{i,nrows},{j,ncols}];

elevData=Table[uncertSwath[[i,j]],{i,nrows,1,-1},{j,ncols}];

elevDataTemp2=elevData+elevDataTemp1;

ClearAll[elevData,elevDataTemp1,uncertSwath];

Swath 3 of 7:  Gentle dip, larger strike azimuth (GL)

Add the dip-trend uncertainty from the mean dip azimuth (i.e., from the mean dip trend)

gentleDipAzL=If[(npDipPlunge-npDipAngUncert)<0,If[((npDipTr+npDipTrendUncert)>180),((npDipTr+npDipTrendUncert)-180),((npDipTr+npDipTrendUncert)+180)],If[((npDipTr+npDipTrendUncert)>360),(npDipTr+npDipTrendUncert)-360,(npDipTr+npDipTrendUncert)]];

Convert the fault dip vector from trend and plunge to a unit location vector

dipVectGL=unitVector[makeVector[gentleDipAng,gentleDipAzL]];

Find the strike vector defined using the right-hand rule

strikeGL=If[(gentleDipAzL<90),(gentleDipAzL+270),(gentleDipAzL-90)];

strikeVectGL={Sin[strikeGL Degree],Cos[strikeGL Degree],0};

Find the unit vector that is normal to the fault plane and is directed upwards

normalVectGL=unitVector[Cross[dipVectGL,strikeVectGL]];

Find the width of the nodal-plane uncertainty zone with the given dip angle, dip azimuth, and triaxial hypocenter uncertainty ellipsoid

zoneHalfWidthGL = findHalfWidth[gentleDipAzL, gentleDipAng, eh1, eh1Az, eh2, ez];

N[zoneHalfWidthGL]

Find the boundaries of the uncertainty swath as it intersects the ground surface

uncertSwath=Table[pointEvaluator[((cellsize*(j-1))-focusMean[[1]]),((cellsize*(nrows-i))-focusMean[[2]]),demImageFileM[[i,j]],zoneHalfWidthGL,normalVectGL,minusOne],{i,nrows},{j,ncols}];

elevData=Table[uncertSwath[[i,j]],{i,nrows,1,-1},{j,ncols}];

elevDataTemp1=elevData+elevDataTemp2;

ClearAll[elevData,elevDataTemp2,uncertSwath];

Dimensions[elevDataTemp1]

Swath 4 of 7:  Steep dip, smaller strike azimuth (SS)

Add the dip-angle uncertainty to the mean dip angle (i.e., from the plunge of the mean dip vector)

steepDipAng=If[(npDipPlunge+npDipAngUncert)>90,(180-(npDipPlunge+npDipAngUncert)),(npDipPlunge+npDipAngUncert)];

Subtract the dip-trend uncertainty from the mean dip azimuth (i.e., from the mean dip trend)

steepDipAzS=If[(npDipPlunge+npDipAngUncert)>90,If[((npDipTr-npDipTrendUncert)>180),((npDipTr-npDipTrendUncert)-180),((npDipTr-npDipTrendUncert)+180)],If[((npDipTr-npDipTrendUncert)<0),360+(npDipTr-npDipTrendUncert),(npDipTr-npDipTrendUncert)]];

Convert the fault dip vector from trend and plunge to a unit location vector

dipVectSS=unitVector[makeVector[steepDipAng,steepDipAzS]];

Find the strike vector defined using the right-hand rule

strikeSS=If[(steepDipAzS<90),(steepDipAzS+270),(steepDipAzS-90)];

strikeVectSS={Sin[strikeSS Degree],Cos[strikeSS Degree],0};

Find the unit vector that is normal to the fault plane and is directed upwards

normalVectSS=unitVector[Cross[dipVectSS,strikeVectSS]];

Find the width of the nodal-plane uncertainty zone with the given dip angle, dip azimuth, and triaxial hypocenter uncertainty ellipsoid

zoneHalfWidthSS = findHalfWidth[steepDipAzS, steepDipAng, eh1, eh1Az, eh2, ez];

N[zoneHalfWidthSS]

Find the boundaries of the uncertainty swath as it intersects the ground surface

uncertSwath=Table[pointEvaluator[((cellsize*(j-1))-focusMean[[1]]),((cellsize*(nrows-i))-focusMean[[2]]),demImageFileM[[i,j]],zoneHalfWidthSS,normalVectSS,minusOne],{i,nrows},{j,ncols}];

elevData=Table[uncertSwath[[i,j]],{i,nrows,1,-1},{j,ncols}];

elevDataTemp2=elevData+elevDataTemp1;

ClearAll[elevData,elevDataTemp1,uncertSwath];

Swath 5 of 7:  Steep dip, larger strike azimuth (SL)

Add the dip-trend uncertainty from the mean dip azimuth (i.e., from the mean dip trend)

steepDipAzL=If[(npDipPlunge+npDipAngUncert)>90,If[((npDipTr+npDipTrendUncert)>180),((npDipTr+npDipTrendUncert)-180),((npDipTr+npDipTrendUncert)+180)],If[((npDipTr+npDipTrendUncert)>360),(npDipTr+npDipTrendUncert)-360,(npDipTr+npDipTrendUncert)]];

Convert the fault dip vector from trend and plunge to a unit location vector

dipVectSL=unitVector[makeVector[steepDipAng,steepDipAzL]];

Find the strike vector defined using the right-hand rule

strikeSL=If[(steepDipAzL<90),(steepDipAzL+270),(steepDipAzL-90)];

strikeVectSL={Sin[strikeSL Degree],Cos[strikeSL Degree],0};

Find the unit vector that is normal to the fault plane and is directed upwards

normalVectSL=unitVector[Cross[dipVectSL,strikeVectSL]];

Find the width of the nodal-plane uncertainty zone with the given dip angle, dip azimuth, and triaxial hypocenter uncertainty ellipsoid

zoneHalfWidthSL = findHalfWidth[steepDipAzL, steepDipAng, eh1, eh1Az, eh2, ez];

N[zoneHalfWidthSL]

Find the boundaries of the uncertainty swath as it intersects the ground surface

uncertSwath=Table[pointEvaluator[((cellsize*(j-1))-focusMean[[1]]),((cellsize*(nrows-i))-focusMean[[2]]),demImageFileM[[i,j]],zoneHalfWidthSL,normalVectSL,minusOne],{i,nrows},{j,ncols}];

elevData=Table[uncertSwath[[i,j]],{i,nrows,1,-1},{j,ncols}];

elevDataTemp1=elevData+elevDataTemp2;

ClearAll[elevData,elevDataTemp2,uncertSwath];

Swath 6 of 7:  Steep dip, mean strike azimuth (SM)

Add the dip-trend uncertainty from the mean dip azimuth (i.e., from the mean dip trend)

steepDipAzM=If[(npDipPlunge+npDipAngUncert)>90,If[(npDipTr>180),(npDipTr-180),(npDipTr+180)],npDipTr];

Convert the fault dip vector from trend and plunge to a unit location vector

dipVectSM=unitVector[makeVector[steepDipAng,steepDipAzM]];

Find the strike vector defined using the right-hand rule

strikeSM=If[(steepDipAzM<90),(steepDipAzM+270),(steepDipAzM-90)];

strikeVectSM={Sin[strikeSM Degree],Cos[strikeSM Degree],0};

Find the unit vector that is normal to the fault plane and is directed upwards

normalVectSM=unitVector[Cross[dipVectSM,strikeVectSM]];

Find the width of the nodal-plane uncertainty zone with the given dip angle, dip azimuth, and triaxial hypocenter uncertainty ellipsoid

zoneHalfWidthSM = findHalfWidth[steepDipAzM, steepDipAng, eh1, eh1Az, eh2, ez];

N[zoneHalfWidthSM]

Find the boundaries of the uncertainty swath as it intersects the ground surface

uncertSwath=Table[pointEvaluator[((cellsize*(j-1))-focusMean[[1]]),((cellsize*(nrows-i))-focusMean[[2]]),demImageFileM[[i,j]],zoneHalfWidthSM,normalVectSM,minusOne],{i,nrows},{j,ncols}];

elevData=Table[uncertSwath[[i,j]],{i,nrows,1,-1},{j,ncols}];

elevDataTemp2=elevData+elevDataTemp1;

ClearAll[elevData,elevDataTemp1,uncertSwath];

Swath 7 of 7:  The middle of the road

Convert the fault dip vector from trend and plunge to a unit location vector

npDipVect=unitVector[makeVector[npDipPlunge,npDipTr]];

Find the strike vector defined using the right-hand rule

npStrike=If[(npDipTr<90),(npDipTr+270),(npDipTr-90)];

npStrikeVect = {Sin[npStrike Degree], Cos[npStrike Degree], 0};

Find the unit vector that is normal to the fault plane and is directed upwards

npNormalVect=unitVector[Cross[npDipVect,npStrikeVect]];

Find the width of the nodal-plane uncertainty zone with the given dip angle, dip azimuth, and triaxial hypocenter uncertainty ellipsoid

zoneHalfWidthMean = findHalfWidth[npDipTr, npDipPlunge, eh1, eh1Az, eh2, ez];

N[zoneHalfWidthMean]

uncertSwath=Table[pointEvaluator[((cellsize*(j-1))-focusMean[[1]]),((cellsize*(nrows-i))-focusMean[[2]]),demImageFileM[[i,j]],zoneHalfWidthMean,npNormalVect,minusOne],{i,nrows},{j,ncols}];

elevData=Table[uncertSwath[[i,j]],{i,nrows,1,-1},{j,ncols}];

elevDataTemp1=elevData+elevDataTemp2;

summary = Table[If[elevDataTemp1[[i, j]] < 0, minusOne, 1], {i, nrows}, {j, ncols}];

traceImgFileSum=ListContourPlot[summary,ContourShading->False,AspectRatio->Automatic,Contours->{topoBase}];

ClearAll[elevData,elevDataTemp1,elevDataTemp2,uncertSwath];

Epicenter

epicenterDot = Table[epicenterEval[((cellsize*(j - 1)) - focusMean[[1]]), ((cellsize*(nrows - i)) - focusMean[[2]]), demImageFileM[[i, j]], epicenterWidth,minusOne], {i, nrows}, {j, ncols}];

elevData = Table[epicenterDot[[i, j]], {i, nrows, 1, -1}, {j, ncols}];

traceImgFileE = ListContourPlot[elevData, ContourShading -> False, AspectRatio -> Automatic, Contours -> {topoBase}];

ClearAll[elevData, epicenterDot];

ClearAll[npDipTr,npDipTrend,npDipPlunge,npDipVect,npStrike,npStrikeVect,npNormalVect,zoneHalfWidthMean,gentleDipAng,gentleDipAzS,dipVectGS,strikeGS,strikeVectGS,normalVectGS,zoneHalfWidthGS,gentleDipAzL,dipVectGL,strikeGL,strikeVectGL,normalVectGL,zoneHalfWidthGL,steepDipAng,steepDipAzS,dipVectSS,strikeSS,strikeVectSS,normalVectSS,zoneHalfWidthSS,steepDipAzL,dipVectSL,strikeSL,strikeVectSL,normalVectSL,zoneHalfWidthSL,steepDipAzM,dipVectSM,strikeSM,strikeVectSM,normalVectSM,zoneHalfWidthSM];

Output

outFile1=Show[hillshadeMap1,traceImgFileMThin,traceImgFileE];

Show[%]

outFile2=Show[traceImgFileMThin,traceImgFileSum,traceImgFileE];

Show[%]

outFile3=Show[hillshadeMap1,outFile2];

Show[%]

The black square in the "outFile2" image marks the epicenter of the earthquake.  The length of each side is approximately "epicenterWidth" in meters.

ClearAll[traceImgFileMean,traceImgFileGS,traceImgFileGL,traceImgFileSS,traceImgFileSL,traceImgFileSM,traceImgFileMThin];

Export data and image files

The file created in the next line is JPEG image that contains the graphic showing the hillshade image developed from the DEM, curves marking the boundaries of the uncertainty regions, and a circle with 1 km diameter centered on the epicenter.

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/test/map1.jpg", outFile1];

Export["/Users/vince/Desktop/test/map2.jpg", outFile2];

Export["/Users/vince/Desktop/test/map3.jpg", outFile3];

ClearAll[outFile1,outFile2,outFile3];

Closing

How long did this program take to run, in minutes?

minutesForProcessing=(AbsoluteTime[]-startTime)/60;

N[%]