/* Example script on how to use the SH() function to calculate the direction of
* maximum horizontal stress using only the
* directions of the principal stress and R = (S1 - S2)/(S1 - S3).
* The methodology is 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.
*
* Copyright (C) 2007 Bjorn Lund
*
* This program is free software; you have the permission to use, copy, and
* distribute it, either verbatim or with modifications, either gratis or for a
* fee. 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., doi: 10.1111/j.1365-246X.2007.03468.x.
*/

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

#define MAX 400

int main(int ac, char **av) {

   FILE    *infp;
   char    buf[MAX];
   double  S1[3], S2[3], S3[3], R, pi;
   int     i, cnt=0;

   extern double SH(double *,double *,double *,double);

   if (ac<2) {
      fprintf(stderr,"\nGive file with two lines per stress state:\n");
      fprintf(stderr,"\tS1_trend S1_plunge S1_N S1_E S1_D\n");
      fprintf(stderr,"\tS3_trend S3_plunge S3_N S3_E S3_D\n\n");
      exit(-1);
   }

   if ((infp = fopen(av[1],"r")) == NULL) {
      fprintf(stderr,"Cannot find file %s.\n",av[1]);
      exit(-1);
   }
   pi = 4.0*atan(1.0);

   // Read the stress states one by one and calculate SH for 20 R values.
   // The input file must have one line for S1 and one for S3, each line
   // containing trend, plunge, north-, east-, down-coordinates.
   while (fgets(buf,MAX,infp) != NULL) {
      // S1
      //printf("%s",buf);
      if (sscanf(buf,"%*f %*f %lf %lf %lf",&S1[0],&S1[1],&S1[2]) != 3) {
         fprintf(stderr,"Error in reading in S1 coordinates.\n");
         exit(-1);
      }
      if (fgets(buf,MAX,infp) != NULL) {
         // S3
         //printf("%s",buf);
         if (sscanf(buf,"%*f %*f %lf %lf %lf",&S3[0],&S3[1],&S3[2]) != 3) {
            fprintf(stderr,"Error in reading in S3 coordinates.\n");
            exit(-1);
	 }
      }
      cnt += 1;

      // Calculate S2 by S2 = S3xS1
      S2[0] = S3[1]*S1[2] - S3[2]*S1[1];
      S2[1] = S3[2]*S1[0] - S3[0]*S1[2];
      S2[2] = S3[0]*S1[1] - S3[1]*S1[0];

      //printf("%f %f %f\n%f %f %f\n%f %f %f\n",S1[0],S1[1],S1[2],
      //                                  S2[0],S2[1],S2[2],S3[0],S3[1],S3[2]);
 
      // Do 20 R values
      fprintf(stdout," Stress state %d\n  R    SH [deg]\n",cnt);
      for (i=0; i<21; i++) {
         R = (double)i/20.0;
	 fprintf(stdout," %.2f   %.2f\n",R,SH(S1,S2,S3,R)*180.0/pi);
      }
      fprintf(stdout,"\n");
   }
   return 0;
}