program main
       implicit none

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

       character*100    infil
       double precision alpha, S1(3), S2(3), S3(3), R, pi
       double precision sh, dum1, dum2
       integer          i, cnt, ac, err

       ac = COMMAND_ARGUMENT_COUNT()
       if (ac.LT.1.OR.ac.GT.1) then
          write(*,*) '\nGive file with two lines per stress state:'
          write(*,*) '\tS1_trend S1_plunge S1_N S1_E S1_D'
          write(*,*) '\tS3_trend S3_plunge S3_N S3_E S3_D'
          stop
       else
         CALL GET_COMMAND_ARGUMENT(1,infil)
       endif

       open(10,FILE=infil,IOSTAT=err,STATUS='old')
       if (err.NE.0) then
          write(*,200) 'Error reading file ',infil
 200      format(A,A60)
	  stop
       endif

       pi = 4.0*atan(1.0d0)

C Read the stress states one by one and calculate SH for 20 R values.
C The input file must have one line for S1 and one for S3, each line
C containing trend, plunge, north-, east-, down-coordinates.
       cnt = 0
       do while (err.eq.0)
          read(10,*,iostat=err) dum1,dum2,S1(1),S1(2),S1(3)
	  if (err.NE.0) then
	     stop
	  else
	     read(10,*,iostat=err) dum1,dum2,S3(1),S3(2),S3(3)
	  end if
	  cnt = cnt + 1
C         write(*,*) S1,S3,cnt
C Calculate S2 by S2 = S3xS1
          S2(1) = S3(2)*S1(3) - S3(3)*S1(2)
          S2(2) = S3(3)*S1(1) - S3(1)*S1(3)
          S2(3) = S3(1)*S1(2) - S3(2)*S1(1)
C Do 20 R values
          write(*,205) cnt
 205      format(' Stress state ',i2)
          write(*,*) '  R     SH [deg]'
	  do i = 0,20,1
	     R = i/20.0d0
	     sh = alpha(S1,S2,S3,R)
	     write(*,210) R, sh*180./pi
 210         format(2x,F4.2,2x,F7.2)
	  end do
       end do

       end