function alpha = SH(S1,S2,S3,R)

% SH() Matlab function for the calculation of the direction of maximum
%      horizontal stress.
%
%   Copyright (C) 1998-2007 Bjorn Lund
%
% This program is free software; you have the permission to use, copy, and
% distribute, either verbatim or with modifications, either withot cost or for a
% fee, as long as the copyright notice above is distributed with the software
% and you clearly state modifications made to the original code.
% If you use this algorithm for work that is to be published, please cite
% the paper:
% Lund and Townend, (2007). Calculating horizontal stress orientations
%    with full or partial knowledge of the tectonic stress tensor,
%    Geophys. J. Int., 170, 1328-1335, doi: 10.1111/j.1365-246X.2007.03468.x.
%
%---------------------------------------------------------------------------
%
%   SH()
%      Calculate the direction of maximum horizontal stress using only the
%      directions of the principal stress and R = (S1 - S2)/(S1 - S3).
%      This function uses the methodology described in
%      Lund and Townend, (2007). Calculating horizontal stress orientations
%      with full or partial knowledge of the tectonic stress tensor,
%      Geophys. J. Int., 170, 1328-1335, doi: 10.1111/j.1365-246X.2007.03468.x.
%      In particular, SH() uses Equations 11 and 10 from the paper.
%
%   Input:
%      S1, S2, S3 are the principal stress orientations.
%                 The variables holds the coordinates in the North, East
%                 and Down geographical coordinate system as vectors
%                 of three floats, e.g. S1 = [s1N s1E s1D]
%      R is a single number describing the relative magnitude of S2 with
%        respect to S1 and S3. R = (S1 - S2)/(S1 - S3)
%
%   Returns:
%      The direction of SH from North, angle in radians
%

EPS = 1.0e-8;
UNDEF = -999;

% Calculate the direction of maximum horizontal stress using Eq. 11 in
% Lund and Townend (2007)
Y = 2.0*(S1(1)*S1(2) + (1.0 - R)*(S2(1)*S2(2)));
X = S1(1)*S1(1) - S1(2)*S1(2) + (1.0 - R)*(S2(1)*S2(1) - S2(2)*S2(2));

% Is the denominator (X here) from Eq. 11 in Lund and Townend (2007) zero?
if abs(X) < EPS
   % If so, the first term in Eq. 10 is zero and we are left with the
   % second term, which must equal zero for a stationary point.
   % The second term is zero either if
   % s1Ns1E + (1-R)*s2Ns2E = 0   (A)
   % or if
   % cos(2*alpha) = 0            (B)
   % If (A) holds, the direction of SH is undefined since Eq. 10 is zero
   % irrespective of the value of alpha. We therefore check for (A) first.
   % If (A) holds, R = 1 + s1Ns1E/s2Ns2E unless s2Ns2E = 0, in which case
   % s1Ns1E also has to be zero for (A) to hold.
   %
   if abs(S2(1)*S2(2)) < EPS
      % s2Ns2E = 0
      if abs(S1(1)*S1(2)) < EPS
	 alpha = UNDEF;
	 return
      else
	 alpha = pi/4.0;
      end
   else
      if abs(R - (1.0 + S1(1)*S1(2)/S2(1)*S2(2))) < EPS
	 alpha = UNDEF;
	 return
      else
	 alpha = pi/4.0;
      end
   end

% The denominator is non-zero
else
   alpha = atan(Y/X)/2.0;
end

% Have we found a minimum or maximum? Use 2nd derivative to find out.
% A negative 2nd derivative indicates a maximum, which is what we want.
dev2 = -2.0*X*cos(2.0*alpha) - 2.0*Y*sin(2.0*alpha);
if dev2 > 0
   % We found a minimum. Add 90 degrees to get the maximum.
   alpha = alpha + pi/2.0;
end

% The resulting direction of SH is given as [0,180[ degrees.
if alpha < 0
   alpha = alpha + pi;
end
if alpha > pi or abs(alpha - pi) < EPS
   alpha = alpha - pi;
end

return