/* SH() 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 it, either verbatim or with modifications, either without 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., v.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., 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 three
*                 component varables of 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
*/

#include <stdio.h>
#include <math.h>

double SH(double *S1, double *S2, double *S3, double R) {

   double EPS = 1.0e-8, UNDEF = -999;
   double X, Y, alpha, dev2, pi=4.0*atan(1.0), tmp;

   // Check input, should be unit vectors that are orthogonal
   if ((tmp = sqrt(S1[0]*S1[0] + S1[1]*S1[1] + S1[2]*S1[2])) > EPS) {
      S1[0] /= tmp;
      S1[0] /= tmp;
      S1[0] /= tmp;
   }
   if ((tmp = sqrt(S2[0]*S2[0] + S2[1]*S2[1] + S2[2]*S2[2])) > EPS) {
      S2[0] /= tmp;
      S2[0] /= tmp;
      S2[0] /= tmp;
   }
   if ((tmp = sqrt(S3[0]*S3[0] + S3[1]*S3[1] + S3[2]*S3[2])) > EPS) {
      S3[0] /= tmp;
      S3[0] /= tmp;
      S3[0] /= tmp;
   }
   if ((tmp = S1[0]*S2[0] + S1[1]*S2[1] + S1[2]*S2[2]) > 1.0e-4)
      fprintf(stderr,"SH(): S1 and S2 are not orthogonal! %f\n",tmp);
   if ((tmp = S1[0]*S3[0] + S1[1]*S3[1] + S1[2]*S3[2]) > 1.0e-4)
      fprintf(stderr,"SH(): S1 and S3 are not orthogonal! %f\n",tmp);
   if ((tmp = S2[0]*S3[0] + S2[1]*S3[1] + S2[2]*S3[2]) > 1.0e-4)
      fprintf(stderr,"SH(): S2 and S3 are not orthogonal! %f\n",tmp);
   // Check R
   if (R<0 || R>1) {
      fprintf(stderr,"SH(): R is out of range!\n");
      return UNDEF;
   }

   // Calculate the denominator and numerator in Eq. 11
   X = S1[0]*S1[0] - S1[1]*S1[1] + (1.0 - R)*(S2[0]*S2[0] - S2[1]*S2[1]);
   Y = 2.0*(S1[0]*S1[1] + (1.0 - R)*(S2[0]*S2[1]));

   // Is the denominator (X here) from Eq. 11 in Lund and Townend (2007) zero?
   if (fabs(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 (fabs(S2[0]*S2[1]) < EPS) {
         // s2Ns2E = 0
	 if (fabs(S1[0]*S1[1]) < EPS)
	    return UNDEF;
	 else
	    alpha = pi/4.0;
      }
      else {
         if (fabs(R - (1.0 + S1[0]*S1[1]/S2[0]*S2[1])) < EPS)
	    return UNDEF;
	 else
	    alpha = pi/4.0;
      }
   }

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

   // 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 += pi/2.0;

   // The resulting direction of SH is given as [0,180[ degrees.
   if (alpha < 0)
      alpha += pi;
   if (alpha > pi || fabs(alpha - pi) < EPS)
      alpha -= pi;

   return alpha;
}