startTime=AbsoluteTime[];

SLAMCompositeCentroid.nb

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

Begun 5 August 2004;  Revised December 20, 2020.

Introduction

This version of the SLAM code is meant to project a single plane 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 (in the blue 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/vincecronin/Desktop/Strasser_Input_EQ_datafile.xls”];

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/Dec18-BFP-C-mean/Data4BestFitPlaneSLAM.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:
eqDataIn = Import["/Users/vincecronin/Desktop/Strasser_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];

Specify which record in the input earthquake datafile will be acted upon -- which earthquake nodal plane will be used

The value that you input for counter determines which record (which row) of the earthquake dataset is evaluated.

counter=10;

Modify the digital elevation model (DEM) so that the resulting hillshade map has an esthetically pleasing range of tones or colors

The range of tones or colors on the hillshade map generated by this code varies with the difference (relief) between the smallest and largest elevation value in the DEM dataset.  If there are exceptionally small “nodata values” in the elevation dataset, the hillshade map might have a very small range of tones or colors -- it might look flat.  If there are no “nodata values”, a relatively flat ground surface might look like a mountain range on the hillshade map.  The trick is to intercede and re-define the nodata values (if any) as being much closer to the ground-surface elevations, or to define one of the ground-surface elevations as a value that is much lower than the rest but not as low as the original “nodata value.”

In the fragment of the DEM dataset shown above, which has no “nodata” values in the elevation matrix, the very first value in the elevation data (row 7, column 1) was changed to 0.  By trial and error, you can adjust this bogus value to produce the desired range of tones or colors in the hillshade map.

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/Dec18-BFP-C-mean/DogValleyFull27.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.  The zone meridian for work in the Tahoe area is -123

zoneMeridian=-123;

Input the geographic coordinates of the corners of the map area that you are particularly interested in analyzing.  That might be the same as the extent of the input DEM, or it might be a smaller area within that extent.  These data are used to compute the gridConvergence:  the angle between true north and grid north at/near the center of the map area.

nwLat=39.5801; nwLong=-120.2732;
swLat=39.3375; swLong=-120.2727;
seLat=39.3375; seLong=-119.9875;
neLat=39.5801;neLong=-119.9878;

Remember to change the names of the output files

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

The module epicenterEval defines the location of the epicenter on the hillshade map.

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

Computation

Adjustment for input geographic coordinates based on true north to the relevant UTM grid north

midLat=((nwLat-swLat)/2)+swLat;

N[%]

midLongN=((neLong-nwLong)/2)+nwLong;

midLongS=((seLong-swLong)/2)+swLong;

midLong=((midLongS-midLongN)/2)+midLongN;

N[%]

The gridConvergence 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 gridConvergence is positive;  otherwise, it is negative.  
Imagine you have a horizontal vector pointing toward true north (azimuth 0°).  If the longitude of the point-of-action (or the middle of the map area) is greater than the longitude of the central meridian of the UTM zone (i.e., > -123 for Zone 10_), the azimuth of that vector in the UTM grid at that point-of-action will be the true azimuth minus the gridConvergence.

gridConvergence=ArcTan[Tan[(midLong-zoneMeridian)Degree]*Sin[midLat Degree]]*(180/π)

Miscellanious utilitarian parameters that are not usually changed

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=50;

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;

Read the input data

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-gridConvergence)<0),(360+(eh1Azimuth-gridConvergence)),If[((eh1Azimuth-gridConvergence)>360),((eh1Azimuth-gridConvergence)-360),(eh1Azimuth-gridConvergence)]];

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

Make the hillshade image

Select and label the matrix of elevation data

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

Work with the orientation data for a given fault - plane solution

npDipTrend is the dip trend relative to geographic north.

npDipTrend = eqData[[record, 15]]

npDipTr is the dip trend relative to grid north, expressed in degrees azimuth.

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

npDipTrRad is the dip trend relative to grid north, expressed in radians azimuth.

npDipTrRad=npDipTr*(π/180);

Set the lighting angle

lightAzDeg=If[((npDipTrend+180)>360),(npDipTrend-180),(npDipTrend+180)];

lightSourceFrom=If[(npDipTrend≥180),(npDipTrend-180),(npDipTrend+180)]

e is a conversion parameter used to convert from a grid-north CS to a grid-east CS, which is necessary because of the way Mathematica specifies illumination directions (relative to east with positive angles reflecting an anticlockwise or right-handed rotation).

e=90-lightSourceFrom+gridConvergence;

lightDir is the azimuth (in degrees) from which the light shines relative to grid east.

lightDir=If[(e>360),(e-360),If[(e<0),(360+e),e]];

lightDirection is the azimuth (in radians) from which the light shines relative to grid east, in which positive values reflect an anticlockwise (right-handed) rotation such that 90° is grid north, 180° is grid west, 270° is grid south.

lightDirection=lightDir*(π/180);

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

ClearAll[elev3];

Computations specific to the selected plane

Make the trace of the plane

npDipPlunge=eqData[[record,16]];

npDipTrendUncert=eqData[[record,18]];

npDipAngUncert=eqData[[record,19]];

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

Make the mark signifying the epicenter or centroid

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

Output

outFile1=Show[hillshadeMap1];

Show[%]

lightSourceFrom

lightElevation=(π/12)*(180/π);

N[%]

outFile2=Show[traceImgFileMThin,traceImgFileE];

Show[%]

outFile3=Show[outFile1,outFile2];

Show[%]

The black square in each of the maps above 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/Dec18-BFP-C-mean/OutputFiles-CompositeMeanSLAM/10.CompositeMeanV4-map1.jpg", outFile1];

Export["/Users/vince/Desktop/Dec18-BFP-C-mean/OutputFiles-CompositeMeanSLAM/10.CompositeMeanV4-map2.jpg", outFile2];

Export["/Users/vince/Desktop/Dec18-BFP-C-mean/OutputFiles-CompositeMeanSLAM/10.CompositeMeanV4-map3.jpg", outFile3];

ClearAll[outFile1,outFile2,outFile3];

Closing

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

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

N[%]

References

Colley, S.J., 2002, Vector calculus [2nd edition]:  Upper Saddle River, New Jersey, Prentice-Hall, 558 p., (pp. 43-50), ISBN 0-13-041531-6.

Cronin, V.S., 2014, Seismo-Lineament Analysis Method (SLAM), using earthquake focal mechanisms to help recognize seismogenic faults, in Grutzner, C., Choi, J - H., Edwards, P., and Kim, Y - .S., [editors], Proceeding of the 5th International INQUA Meeting on Paleoseismology, Active Tectonics and Archeoseismology, 21-27 September 2014, p.28 - 31, ISBN 9791195344109 93450; available via http://CroninProjects.org/Vince/SLAM/CroninINQUA_PATA14.pdf

Cronin, V.S., and Sverdrup, K.A., 1998, Preliminary assessment of the seismicity of the Malibu Coast Fault Zone, southern California, and related issues of philosophy and practice, in Welby, C.W., and Gowan, M.E. [editors], A Paradox of Power--Voices of Warning and Reason in the Geosciences:  Geological Society of America, Reviews in Engineering Geology, p. 123-155.

Cronin, V.S., Byars, B.W., and Gammill, T., 2003, Developing techniques for regional structural interpretation using GIS and DEM-based terrain analysis:  Geological Society of America, Abstracts with Programs, v. 34, no. 7, p. 261.  http://gsa.confex.com/gsa/2003AM/finalprogram/abstract_61851.htm

Cronin, V.S., Millard, M.A., Seidman, L.E., and Bayliss, B.G., 2008, The Seismo-Lineament Analysis Method [SLAM] -- A reconnaissance tool to help find seismogenic faults: Environmental and Engineering Geoscience, v. 14, no. 3, p. 199-219.

Davis, H.F., and Snider, A.D., 1987, Introduction to vector analysis [5th edition]:  Boston, Allyn and Bacon, 365 p. (pp. 34-37), ISBN 0-205-10263-8.

Gammill, T., Cronin, V.S., and Byars, B.W., 2004, Combining earthquake focal data and digital map analysis in reconnaissance for active faults, central Santa Monica Mountains and northern Santa Monica Bay, California:  Geological Society of America, Abstracts with Programs,v. 36, no. 5, http://gsa.confex.com/gsa/2004AM/finalprogram/abstract_80352.htm

Haneberg, W.C., 2004, Computational geosciences with Mathematica:  Berlin, Springer-Verlag, 381 p., CD, ISBN 3-540-40245-4.

Snyder, J.P., 1982, Map projections used by the U.S. Geological Survey:  U.S. Geological Survey Bulletin 1532, pp. 63-69 and 233-235.

Other references related to this work are available online at http://croninprojects.org/DVFZ/index.htm and http://croninprojects.org/Vince/SLAM/index.htm