#!/usr/bin/python
# encoding: latin-1
#
# 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., v.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, 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.

import sys, math, SH, BL

if len(sys.argv) != 5:
   print '\n\t%s S1_az S1_pl S3_az S3_pl\n' % (sys.argv[0])
   print '    Calculate the direction of the maximum horizontal stress'
   print '    for a range of R-values from 0 to 1.'
   print '    Input the stress tensor as the azimuth (0-360) and'
   print '    plunge (0-90) of S1 and S3. The orientation of S2 is'
   print '    computed from S1 and S3.\n'
   sys.exit()

S1az = float(sys.argv[1])
S1pl = float(sys.argv[2])
S3az = float(sys.argv[3])
S3pl = float(sys.argv[4])

deg2rad = math.pi/180.0

# Calculate coordinates in the North, East, Down (NED) system from the
# trend and plunge angles.
S1 = BL.TrendPlunge2Coords([S1az*deg2rad,S1pl*deg2rad])
S3 = BL.TrendPlunge2Coords([S3az*deg2rad,S3pl*deg2rad])

# Make sure the stress vectors are orthogonal
S1,S3 = BL.Orthogonalize(S1,S3)

# Recalculate trend and plunge for the adjusted stress vectors
a = BL.Coords2TrendPlunge(S1)
S1az = a[0]/deg2rad;   S1pl = a[1]/deg2rad
b = BL.Coords2TrendPlunge(S3)
S3az = b[0]/deg2rad;   S3pl = b[1]/deg2rad
print '\nAdjusted S1: Trend %.2f Plunge %.2f' % (S1az,S1pl)
print 'Adjusted S3: Trend %.2f Plunge %.2f' % (S3az,S3pl)

# Calculate S2
S2 = BL.Cross(S3,S1)
tp2 = BL.Coords2TrendPlunge(S2)
S2az = tp2[0]/deg2rad
S2pl = tp2[1]/deg2rad
print 'S2: Trend %.2f Plunge %.2f' % (S2az,S2pl)
print '\nCoordinates in (N,E,D) of:'
print 'S1:',S1
print 'S2:',S2
print 'S3:',S3
print

N1 = 20
np1 = float(N1)
sh = []
# Calculate SH az vs R
for i in range(N1+1):
   tmp = []
   tmp.append(float(i)/np1)
   tmp.append(SH.SH(S1,S2,S3,tmp[0])/deg2rad)
   if tmp[1] < -998:
      tmp[1] = -999
   sh.append(tmp)

# Write out SH
print ' R:  SH azimuth:'
for d in sh:
   print '%.2f   %.2f' % (d[0],d[1])
print