Some geologic utilities

PLEASE DO NOT REDISTRIBUTE WITHOUT PERMISSION
2014 February 14 by Joshua R. Davis, jdavis@carleton.edu

Miscellany

Clear[arcCos, directionAgreement, nGradient, nHessian, determinantOne, traceZero, frobenius, symmetricNormSquared, argMin, ellipsoidPower];

(* Like ArcCos, but ignores slight numerical errors. *)
arcCos[x_] :=
	If[x ≥ 1, 0,
		If[x ≤ -1, Pi,
			ArcCos[x]]];

(* Cosine squared of the angle between the two vectors. *)
(* Equals 0 if perpendicular, 1 if parallel (or antiparallel). *)
directionAgreement[v_, w_] := 
(Dot[v, w])^2 / (Dot[v, v] * Dot[w, w]);

(* Numerically approximates the gradient of f at x = vals, using step size h. *)
nGradient[f_, vals_, h_] :=
Map[Function[d, (Apply[f, vals + d] - Apply[f, vals - d]) / (2 h)], h IdentityMatrix[Length[vals]]];

(* Numerically approximates the Hessian of f at x = vals, using step size h. *)
nHessian[f_, vals_, h_] :=
	Module[{d, res, i, j},
		d = h IdentityMatrix[Length[vals]];
		res = Table[Table[
			(Apply[f, vals + d[[i]] + d[[j]]] - Apply[f, vals + d[[i]] - d[[j]]] - Apply[f, vals - d[[i]] + d[[j]]] + Apply[f, vals - d[[i]] - d[[j]]]) / (4 h^2),
			{j, 1, i}], {i, 1, Length[vals]}];
		For[i = 1, i < Length[vals], i += 1,
			For[j = i + 1, j ≤ Length[vals], j += 1,
				res[[i]] = Append[res[[i]], res[[j, i]]]]];
		res];

(* Returns the matrix scaled to have determinant 1. *)
determinantOne[mat_] :=
mat / (Det[mat])^(1 / Length[mat]);

(* Subtracts from the matrix's diagonals, to give it trace zero. *)
traceZero[mat_] :=
	mat - (Tr[mat] / 3) IdentityMatrix[Length[mat]];

(* Returns the Frobenius inner product tr(Mt N) = tr(M Nt). *)
frobenius[m_, n_] :=
	Tr[Transpose[m].n];

(* Like Frobenius norm, but assuming the matrix is symmetric and only half should be measured. *)
symmetricNormSquared[e_] :=
	e[[1, 1]]^2 + e[[1, 2]]^2 + e[[1, 3]]^2 + e[[2, 2]]^2 + e[[2, 3]]^2 + e[[3, 3]]^2;

(* Returns a list of indices into the array, where the minimum value in the array occurs. *)
argMin[array_] :=
	Module[{i, j},
		If[VectorQ[array],
			If[Length[array] == 1,
				{1},
				i = argMin[array[[2;;]]];
				If[array[[1]] <= array[[i[[1]] + 1]],
					{1},
					{i[[1]] + 1}]],
			If[Length[array] == 1,
				Join[{1}, argMin[array[[1]]]],
				i = argMin[array[[1]]];
				j = argMin[array[[2;;]]];
				If[Apply[Part, Join[{array, 1}, i]] <= Apply[Part, Join[{array, j[[1]] + 1}, j[[2;;]]]],
					Join[{1}, i],
					Join[{j[[1]] + 1}, j[[2;;]]]]]]];

(* Given positive-semi-definite symmetric matrix E *)
(* and real number x, computes E^x. *)
ellipsoidPower[e_, x_] :=
	Module[{vals, vecs},
		{vals, vecs} = Eigensystem[e];
		vecs = Map[unit, vecs];
		Transpose[vecs].DiagonalMatrix[vals^x].vecs];

Vectors

Clear[sphericalFromCartesian, cartesianFromSpherical, cylindricalFromCartesian, cartesianFromCylindrical, cartesianFromHorizontal, sphericalFromTrendPlungeDeg, trendPlungeDegFromSpherical, cartesianFromTrendPlungeDeg, trendPlungeDegFromCartesian, lowerHemisphereTrendPlungeDeg, sphericalFromStrikeDipDeg, strikeDipDegFromSpherical, cartesianFromStrikeDipDeg, strikeDipDegFromCartesian, strikeLineFromPole, cartesianFromStrikeDipRakeDeg, strikeDipRakeDegFromCartesian, randomSpherePoint, unit, lowerHemisphere, orthogonalVector, mostUpwardOrthogonal];

(* Convert between spherical and Cartesian coordinates. *)
sphericalFromCartesian[{x_, y_, z_}] :=
	Module[{rho, phi, rhoSinPhi, theta},
		rho = Sqrt[x^2 + y^2 + z^2];
		If[rho == 0,
	{0, 0, 0},
	phi = ArcCos[z / rho];
	rhoSinPhi = rho Sin[phi];
	If[rhoSinPhi == 0,
				If[z >= 0,
					{z, 0, 0},
					{-z, Pi, 0}],
				{rho, phi, ArcTan[x, y]}]]];
cartesianFromSpherical[{rho_, phi_, theta_}] := 
{rho * Sin[phi] * Cos[theta], rho * Sin[phi] * Sin[theta], rho * Cos[phi]};

(* Convert between cylindrical and spherical coordinates. *)
cylindricalFromCartesian[{x_, y_, z_}] :=
	{Sqrt[x^2 + y^2], ArcTan[x, y], z};
cartesianFromCylindrical[{r_, theta_, z_}] :=
	{r Cos[theta] Sqrt[1 - z^2], r Sin[theta] Sqrt[1 - z^2], z}
(* Horizontal coordinates are cylindrical coordinates on the unit sphere. *)
(* As a map between sphere and cylinder, they are area-preserving. *)
horizontalFromCartesian[{x_, y_, z_}] :=
	If[x == 0 && y == 0,
		{0, z},
		{ArcTan[x, y], z}];
cartesianFromHorizontal[{theta_, z_}] :=
	{Cos[theta] Sqrt[1 - z^2], Sin[theta] Sqrt[1 - z^2], z};

(* Trend and plunge describe the orientation (direction) of a ray. *)
(* Trend is measured in degrees, clockwise from the positive y axis. *)
(* Plunge is measured in degrees, down from the x-y plane. *)
sphericalFromTrendPlungeDeg[{trendDeg_, plungeDeg_}] :=
{1, (plungeDeg + 90) Degree, Mod[(90 - trendDeg) Degree, 2 Pi]};
trendPlungeDegFromSpherical[{rho_, phi_, theta_}] :=
	{Mod[90 - theta / Degree, 360], phi / Degree - 90};
cartesianFromTrendPlungeDeg[x_] :=
cartesianFromSpherical[sphericalFromTrendPlungeDeg[x]];
trendPlungeDegFromCartesian[x_] :=
	trendPlungeDegFromSpherical[sphericalFromCartesian[x]];
lowerHemisphereTrendPlungeDeg[{trendDeg_, plungeDeg_}] :=
	If[plungeDeg >= 0,
		{trendDeg, plungeDeg},
		{Mod[trendDeg + 180, 360], -plungeDeg}];

(* Strike and dip describe the orientation of a plane. *)
(* Strike is measured in degrees, clockwise from the positive y axis. *)
(* Dip is measured in degrees, down from the x-y plane. *)
(* These functions assume the right-hand rule: *)
(* the plane dips to the right, when walking along strike. *)
(* Spherical and Cartesian coordinates are for the POLE to the plane. *)
sphericalFromStrikeDipDeg[{strikeDeg_, dipDeg_}] :=
{1, dipDeg Degree, (360 - strikeDeg) Degree};
strikeDipDegFromSpherical[{rho_, phi_, theta_}] :=
	If[phi ≤ Pi / 2,
		{Mod[360 - theta / Degree, 360], phi / Degree},
		strikeDipDegFromSpherical[{rho, Pi - phi, Pi + theta}]];
cartesianFromStrikeDipDeg[x_] :=
cartesianFromSpherical[sphericalFromStrikeDipDeg[x]];
strikeDipDegFromCartesian[x_] :=
	strikeDipDegFromSpherical[sphericalFromCartesian[x]];
strikeLineFromPole[pole_] :=
	Module[{theta},
		theta = If[pole[[3]] < 0,
			ArcTan[pole[[2]], -pole[[1]]],
			ArcTan[-pole[[2]], pole[[1]]]];
		{Cos[theta], Sin[theta], 0}];

(* Strike, dip, and rake describe the orientation of a plane and a ray *)
(* in that plane. Strike and dip are discussed above. Rake is measured *)
(* down-dip from the strike line. Returns {rakeVector, poleVector, *)
(* strikeVector}. Pole is orthogonal to strike and rake. Rake and *)
(* strike are not orthogonal, unless rake = 90 degrees or 270 degrees. *)
cartesianFromStrikeDipRakeDeg[sdrDeg_] :=
	Module[{s, d, r, matrix},
		s = -sdrDeg[[1]] Degree;
		s = {{Cos[s], -Sin[s], 0}, {Sin[s], Cos[s], 0}, {0, 0, 1}};
		d = sdrDeg[[2]] Degree;
		d = {{Cos[d], 0, Sin[d]}, {0, 1, 0}, {-Sin[d], 0, Cos[d]}};
		r = -sdrDeg[[3]] Degree + Pi / 2;
		r = {{Cos[r], -Sin[r], 0}, {Sin[r], Cos[r], 0}, {0, 0, 1}};
		matrix = Transpose[s.d.r];
		{matrix[[1]], matrix[[3]], Transpose[s][[2]]}];
strikeDipRakeDegFromCartesian[{rakeVec_, poleVec_, strikeVec_}] :=
	Module[{strikeDeg, dipDeg, rakeDeg, s, d},
		strikeDeg = trendPlungeDegFromCartesian[strikeVec][[1]];
		dipDeg = 90 - trendPlungeDegFromCartesian[-poleVec][[2]];
		s = {{Cos[strikeDeg Degree], -Sin[strikeDeg Degree], 0}, {Sin[strikeDeg Degree], Cos[strikeDeg Degree], 0}, {0, 0, 1}};
		d = {{Cos[-dipDeg Degree], 0, Sin[-dipDeg Degree]}, {0, 1, 0}, {-Sin[-dipDeg Degree], 0, Cos[-dipDeg Degree]}};
		rakeDeg = trendPlungeDegFromCartesian[d.s.rakeVec][[1]];
		{strikeDeg, dipDeg, rakeDeg}];

(* Returns a uniformly chosen point on the unit sphere. *)
randomSpherePoint[] :=
	cartesianFromHorizontal[{RandomReal[{0, 2 Pi}], RandomReal[{-1, 1}]}];

(* Normalizes the vector, assuming it is nonzero. *)
unit[v_] :=
	v / Norm[v];

(* Returns the parallel vector on the lower unit hemisphere. *)
lowerHemisphere[v_] := 
v / If[v[[-1]] <= 0, Norm[v], -Norm[v]];

(* Returns some unit vector orthogonal to the given nonzero vector. *)
(* This function is not continuous, as a map from the unit sphere to itself. *)
orthogonalVector[u_] :=
	If[Abs[u[[1]]] < 0.5,
		Orthogonalize[{u, {1, 0, 0}}][[2]],
		Orthogonalize[{u, {0, 1, 0}}][[2]]];

(* Among all vectors perpendicular to v, returns the one pointing most up. *)
mostUpwardOrthogonal[v_] :=
	Module[{rpt},
		rpt = sphericalFromCartesian[v];
		If[rpt[[2]] ≤ Pi / 2,
			cartesianFromSpherical[{rpt[[1]], Pi / 2 - rpt[[2]], rpt[[3]] + Pi}],
			cartesianFromSpherical[{rpt[[1]], rpt[[2]] - Pi / 2, rpt[[3]]}]]];

Hemispherical plots

The equal angle (i.e. stereographic) projection maps to a disk of radius 1; the equal area projection maps to a disk of radius sqrt(2). Both assume that their arguments lie on the unit sphere.
	The hemisphericalPlot function draws an equal-area or equal-area hemispherical plot. The first argument is a list of points to be drawn as open circles. (All points are 3D Cartesian; this function pre-processes them to the unit lower hemisphere for you.) The second argument is a list of points to be drawn as filled circles. You can insert some Mathematica graphics commands into these point lists; for example, inserting Red switches the drawing color to red. The third argument is a list of curves; each curve is a list of points. You may not insert graphics commands into the curve list. There are also two options; the pointSize option is a positive real number, by default 0.02. The projection option is either equalAreaProjection (the default) or equalAngleProjection.

Clear[equalAngleProjection, equalAreaProjection, fringePrimitives, pointPrimitives, hemisphericalPlot];

(* Returns the 2D stereographic projection of the given 3D point. *)
equalAngleProjection[{x_, y_, z_}] := 
	{x / (1 - z), y / (1 - z)};
(* Returns the 2D Lambert azimuthal equal-area projection of the 3D point. *)
equalAreaProjection[{x_, y_, z_}] := 
	Module[{scaling}, 
		If[x^2 + y^2 == 0, {0, 0},
			scaling = Sqrt[(2 (1 + z)) / (x^2 + y^2)];
			{scaling * x, scaling * y}]];

(* Helper functions for hemisphericalPlot. *)
fringePrimitives[plotRadius_, padding_] :=
	{Circle[{0.0, 0.0}, plotRadius],
	Table[Line[{
			{plotRadius * Cos[angle], plotRadius * Sin[angle]},
{(plotRadius + padding) * Cos[angle], (plotRadius + padding) * Sin[angle]}}],
{angle, 0, 2 Pi, Pi / 18}],
	Text["000°",{0,plotRadius + 2 * padding}],
	Text["090°",{plotRadius + 2 * padding,0}],
	Text["180°",{0,-(plotRadius + 2 * padding)}],
	Text["270°",{-(plotRadius + 2 * padding),0}]};

pointPrimitives[projection_, pointList_, primitive_] :=
	Map[
		Function[item,
			If[ListQ[item],
			primitive[projection[lowerHemisphere[item]]],
			item]],
		pointList];

Options[hemisphericalPlot] = {pointSize -> 0.02, projection -> equalAreaProjection};
(* Returns an equal-area or equal-angle plot with circles, disks, and curves. *)
hemisphericalPlot[circleList_, diskList_, curveList_, OptionsPattern[]] :=
	Graphics[Join[
			{Black},
			fringePrimitives[If[OptionValue[projection] == equalAreaProjection, Sqrt[2], 1], 0.05],
			{Black},
			pointPrimitives[OptionValue[projection], circleList, Function[point, Circle[point, OptionValue[pointSize]]]],
			{Black},
			pointPrimitives[OptionValue[projection], diskList, Function[point, Disk[point, OptionValue[pointSize]]]],
			{Black},
			Map[Function[curve, Line[Map[OptionValue[projection], Map[lowerHemisphere, curve]]]], curveList]](*,
		AspectRatio -> 1*)];

Miscellaneous statistics

The BHEP statistic is a multivariate normality test. See Henze (2002, p. 486.). This implementation has not been extensively tested.
	If you have a distribution of lines in 3D, and want to know the center of the distribution and nothing else, then use binghamCenter. The input is a list of vectors (3-dimensional Cartesian, unit lower hemisphere). The output is a single vector (3-dimensional Cartesian, unit lower hemisphere), the center of the data. In typical applications it is not necessary to project to the lower hemisphere, but it’s done anyway for safety. The binghamPlot function plots a set of lines with their Bingham center in blue, the “second center” in green, and the “third center” in red. If the binghamPlot function plotted the confidence ellipse, it would be centered on blue, pointing toward green.
	The ajne function computes the Ajne statistic on a sphere --- a measure of how anisotropic a set of lines is. Lines are given as unit, 3D, Cartesian vectors, not necessarily lower-hemisphere. 1 ≤ ajne ≤ n. Small values mean that the lines are unaligned (isotropic). See Gine’ (1975), p. 1262. Other accounts, such as Mardia (1972) and Beran (1965), have errors.

Clear[bootstrapSample, bhepHelper, bhep, binghamCenter, binghamPlot, ajne];

(* Given some data, resamples the data with replacement and returns *)
(* a pseudo-sample of the same size. *)
(* to get a sample from the distribution for the stat. *)
bootstrapSample[data_] :=
	Table[data[[RandomInteger[{1, Length[data]}]]], {Length[data]}];

(* Multivariate normality test statistic "BHEP". See Henze (2002, p. 486). *)
(* xs is a list of n d-dimensional vectors. *)
(* beta > 0 is a smoothing parameter. Try beta = 0.5? *)
(* You can pass a list of betas to get a list of BHEP test statistics. *)
bhepHelper[beta_, n_, d_, first_, second_] :=
	Module[{firstTerm, secondTerm},
	firstTerm = (n^-2) Total[E^((-beta^2 / 2) first)];
	secondTerm = (2 (n^-1) (1 + beta^2)^(-d / 2)) Total[E^((-beta^2 / (2 (1 + beta^2))) second)];
	firstTerm - secondTerm + (1 + 2 beta^2)^(-d / 2)];
bhep[xs_, beta_] :=
	Module[{n, d, mean, s, vals, vecs, sNegHalf, ys, squaresFirst, squaresSecond},
		n = Length[xs];
		d = Length[xs[[1]]];
		mean = Mean[xs];
		(* S is the sample covariance, like the ML estimator: n^-1, not (n - 1)^-1. *)
		s = (n - 1) / n * Covariance[xs];
		{vals, vecs} = Eigensystem[s];
		vecs = Map[unit, vecs];
		sNegHalf = Transpose[vecs].DiagonalMatrix[vals^(-1 / 2)].vecs;
		(* The Ys are the scaled residuals of the data. *)
		ys = Table[sNegHalf.(xs[[j]] - mean), {j, 1, n}];
		squaresFirst = Flatten[Table[Table[
			(ys[[j]] - ys[[k]]).(ys[[j]] - ys[[k]]),
			{j, 1, n}], {k, 1, n}]];
		squaresSecond = Table[ys[[j]].ys[[j]], {j, 1, n}];
		(* Compute the BHEP test statistic for one or more betas. *)
		If[NumericQ[beta],
			bhepHelper[beta, n, d, squaresFirst, squaresSecond],
			Table[
				bhepHelper[beta[[i]], n, d, squaresFirst, squaresSecond],
				{i, 1, Length[beta]}]]];

(* Returns the mean of a set of points on the unit sphere. *)
(* Inputs need not be lower-hemisphere, and output is not necessarily. *)
binghamCenter[vectorList_] := 
	Module[{t, center},
		t = Total[Map[
			Function[vector, Outer[Times, vector, vector]],
			vectorList]];
		center = (N[Eigenvectors[t, 1]])[[1]];
		unit[center]];

(* Returns a plot of the three principal directions of the Bingham *)
(* distribution for the given sample of points on the unit sphere. *)
binghamPlot[vectorList_] :=
	Module[{t, eig},
		t = Total[Map[
			Function[vector, Outer[Times, vector, vector]],
			vectorList]];
		eig = Eigenvectors[t];
		hemisphericalPlot[vectorList, {Red, eig[[3]], Green, eig[[2]], Blue, eig[[1]]}, {}]];

(* Returns the Ajne statistic (Gine', 1975) measuring the anisotropy of *)
(* a set of N lines. Each line is given as a unit 3D Cartesian vector. *)
(* 1 ≤ A ≤ N. Point clusters have A ~ N. Isotropy has A ~ 0.37 N? *)
ajne[us_] :=
	Module[{n, accum, i, j},
		n = Length[us];
		accum = 0;
		For[i = 1, i < n, i += 1,
			For[j = i + 1, j <= n, j += 1,
				accum += ArcCos[Abs[Dot[us[[i]], us[[j]]]]]]];
		n - 4 accum / (Pi n)];