function alpha(S1,S2,S3,R)
       implicit none

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


       integer UNDEF
       double precision S1(3),S2(3),S3(3),R,alpha
       double precision EPS,pi
       double precision X,Y,dev2


       EPS = 1.0e-8
       UNDEF = -999
       pi=3.1415926535897932384626



c Calculate the direction of maximum horizontal stress using Eq. 11 in
c 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))

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

c The denominator is non-zero
       else
         alpha = atan(Y/X)/2.0
       endif

c Have we found a minimum or maximum? Use 2nd derivative to find out.
c 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.gt.0) then
c    We found a minimum. Add 90 degrees to get the maximum.
         alpha = alpha + pi/2.0
       endif

c The resulting direction of SH is given as [0,180[ degrees.
       if(alpha.lt.0) then
         alpha = alpha + pi
       endif
       if(alpha.gt.pi.or.abs(alpha-pi).lt.EPS) then
         alpha = alpha - pi
       endif

       return
       end