C 
C Program TIAXIS to determine the reference symmetry axis of a generally
C anisotropic medium which is approximately transversely isotropic
C
C Version: 7.20
C Date: 2015, May 4
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/klimes.htm
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Data specifying input file:
C     ANIMOD='string'... Name of the input file containing the
C             density-reduced elastic moduli.
C             Description of file ANIMOD.
C             No default.  ANIMOD must be specified and cannot be blank.
C Data specifying output file:
C     TIAXIS='string'... Name of the output file containing
C             the components of the reference symmetry vector,
C             and the ratio of the norm of the derivative of the
C             stiffness tensor with respect to the angle of rotation
C             and the norm of the stiffness tensor.
C             Default: TIAXIS='tiaxis.out'
C
C                                                  
C Input file ANIMOD with the density-reduced elastic moduli:
C (1) 21 elastic moduli read at once in order
C     A1111,A1122,A1133,A1123,A1113,A1112,
C           A2222,A2233,A2223,A2213,A2212,
C                 A3333,A3323,A3313,A3312,
C                       A2323,A2313,A2312,
C                             A1313,A1312,
C                                   A1212
C
C=======================================================================
C
      EXTERNAL ERROR,RSEP1,RSEP3T,EIGEN
C
C     Filenames and parameters:
      CHARACTER*80 FILE
      INTEGER LU
      PARAMETER (LU=1)
C
C     Stiffness tensor
      REAL A11,A12,A13,A14,A15,A16
      REAL A21,A22,A23,A24,A25,A26
      REAL A31,A32,A33,A34,A35,A36
      REAL A41,A42,A43,A44,A45,A46
      REAL A51,A52,A53,A54,A55,A56
      REAL A61,A62,A63,A64,A65,A66
      REAL A(3,3,3,3)
C
C     Other variables:
      INTEGER I,J,K,L,M,N,MN
      REAL EPS(3,3,3),D(3,3,3,3,3),B(6),E(9),ANORM,RHO1,RHO2,RHO3
C
C     Debugging:
*     INTEGER IJ,KL
C
C-----------------------------------------------------------------------
C
C     Reading main input data:
      WRITE(*,'(A)') '+TIAXIS: Enter input filename: '
      FILE=' '
      READ(*,*) FILE
      IF(FILE.EQ.' ') THEN
C       TIAXIS-01
        CALL ERROR('TIAXIS-01: No input SEP file specified')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      END IF
      CALL RSEP1(LU,FILE)
      WRITE(*,'(A)') '+TIAXIS: Working...            '
C
C     Reading stiffness tensor or density reduced elastic moduli:
      CALL RSEP3T('ANIMOD',FILE,' ')
      IF(FILE.EQ.' ') THEN
C       TIAXIS-02
        CALL ERROR('TIAXIS-02: No input stiffness tensor specified')
C       Name ANIMOD of the input file with the stiffness tensor
C       must be specified and cannot be blank.
      END IF
      OPEN(LU,FILE=FILE,STATUS='OLD')
      READ(LU,*) A11,A12,A13,A14,A15,A16
     *              ,A22,A23,A24,A25,A26
     *                  ,A33,A34,A35,A36
     *                      ,A44,A45,A46
     *                          ,A55,A56
     *                              ,A66
      CLOSE(LU)
C
C     Completing the stiffness tensor A(I,J,K,L) of elastic moduli
      A21=A12
      A31=A13
      A41=A14
      A51=A15
      A61=A16
      A32=A23
      A42=A24
      A52=A25
      A62=A26
      A43=A34
      A53=A35
      A63=A36
      A54=A45
      A64=A46
      A65=A56
      A(1,1,1,1)=A11
      A(2,2,1,1)=A21
      A(3,3,1,1)=A31
      A(2,3,1,1)=A41
      A(3,2,1,1)=A41
      A(1,3,1,1)=A51
      A(3,1,1,1)=A51
      A(1,2,1,1)=A61
      A(2,1,1,1)=A61
      A(1,1,2,2)=A12
      A(2,2,2,2)=A22
      A(3,3,2,2)=A32
      A(2,3,2,2)=A42
      A(3,2,2,2)=A42
      A(1,3,2,2)=A52
      A(3,1,2,2)=A52
      A(1,2,2,2)=A62
      A(2,1,2,2)=A62
      A(1,1,3,3)=A13
      A(2,2,3,3)=A23
      A(3,3,3,3)=A33
      A(2,3,3,3)=A43
      A(3,2,3,3)=A43
      A(1,3,3,3)=A53
      A(3,1,3,3)=A53
      A(1,2,3,3)=A63
      A(2,1,3,3)=A63
      A(1,1,2,3)=A14
      A(2,2,2,3)=A24
      A(3,3,2,3)=A34
      A(2,3,2,3)=A44
      A(3,2,2,3)=A44
      A(1,3,2,3)=A54
      A(3,1,2,3)=A54
      A(1,2,2,3)=A64
      A(2,1,2,3)=A64
      A(1,1,3,2)=A14
      A(2,2,3,2)=A24
      A(3,3,3,2)=A34
      A(2,3,3,2)=A44
      A(3,2,3,2)=A44
      A(1,3,3,2)=A54
      A(3,1,3,2)=A54
      A(1,2,3,2)=A64
      A(2,1,3,2)=A64
      A(1,1,1,3)=A15
      A(2,2,1,3)=A25
      A(3,3,1,3)=A35
      A(2,3,1,3)=A45
      A(3,2,1,3)=A45
      A(1,3,1,3)=A55
      A(3,1,1,3)=A55
      A(1,2,1,3)=A65
      A(2,1,1,3)=A65
      A(1,1,3,1)=A15
      A(2,2,3,1)=A25
      A(3,3,3,1)=A35
      A(2,3,3,1)=A45
      A(3,2,3,1)=A45
      A(1,3,3,1)=A55
      A(3,1,3,1)=A55
      A(1,2,3,1)=A65
      A(2,1,3,1)=A65
      A(1,1,1,2)=A16
      A(2,2,1,2)=A26
      A(3,3,1,2)=A36
      A(2,3,1,2)=A46
      A(3,2,1,2)=A46
      A(1,3,1,2)=A56
      A(3,1,1,2)=A56
      A(1,2,1,2)=A66
      A(2,1,1,2)=A66
      A(1,1,2,1)=A16
      A(2,2,2,1)=A26
      A(3,3,2,1)=A36
      A(2,3,2,1)=A46
      A(3,2,2,1)=A46
      A(1,3,2,1)=A56
      A(3,1,2,1)=A56
      A(1,2,2,1)=A66
      A(2,1,2,1)=A66
C
C     Levi-Civita symbols EPS(I,J,K)
      DO 23 K=1,3
        DO 22 J=1,3
          DO 21 I=1,3
            EPS(I,J,K)=0.
   21     CONTINUE
   22   CONTINUE
   23 CONTINUE
      EPS(1,2,3)=1.
      EPS(3,1,2)=1.
      EPS(2,3,1)=1.
      EPS(3,2,1)=-1.
      EPS(1,3,2)=-1.
      EPS(2,1,3)=-1.
C
C     Derivative D(I,J,K,L,M) of elastic moduli with respect to rotation
      ANORM=0.
      DO 36 I=1,3
        DO 35 J=1,3
          DO 34 K=1,3
            DO 33 L=1,3
              ANORM=ANORM+A(I,J,K,L)*A(I,J,K,L)
              DO 32 M=1,3
                D(I,J,K,L,M)=0.
                DO 31 N=1,3
                  D(I,J,K,L,M)=D(I,J,K,L,M)+EPS(M,I,N)*A(N,J,K,L)
     *                                     +EPS(M,J,N)*A(I,N,K,L)
     *                                     +EPS(M,K,N)*A(I,J,N,L)
     *                                     +EPS(M,L,N)*A(I,J,K,N)
   31           CONTINUE
   32         CONTINUE
   33       CONTINUE
   34     CONTINUE
   35   CONTINUE
   36 CONTINUE
C
C     Quadratic B(MN) form of derivatives D(I,J,K,L,M)
      MN=0
      DO 46 N=1,3
        DO 45 M=1,N
          MN=MN+1
          B(MN)=0.
          DO 44 L=1,3
            DO 43 K=1,3
              DO 42 J=1,3
                DO 41 I=1,3
                  B(MN)=B(MN)+D(I,J,K,L,M)*D(I,J,K,L,N)
   41           CONTINUE
   42         CONTINUE
   43       CONTINUE
   44     CONTINUE
   45   CONTINUE
   46 CONTINUE
C
C     Eigenvalues and eigenvectors of quadratic form B(MN)
      CALL EIGEN(B,E,3,0)
      RHO1=SQRT(B(1)/ANORM)
      RHO2=SQRT(B(3)/ANORM)
      RHO3=SQRT(B(6)/ANORM)
C
C     Output file with the reference symmetry vector:
      CALL RSEP3T('TIAXIS',FILE,'tiaxis.out')
      OPEN(LU,FILE=FILE)
      WRITE(LU,'(6(A,I1,A,F8.6))')  'TI',1,'=',E(7),
     *                             ' TI',2,'=',E(8),
     *                             ' TI',3,'=',E(9),
     *                             ' RHO',3,'=',RHO3,
     *                             ' RHO',2,'=',RHO2,
     *                             ' RHO',1,'=',RHO1
C
C     Derivative of the stiffness tensor for debugging
*     DO 52 KL=1,6
*       IF(KL.LE.3) THEN
*         K=KL
*         L=KL
*       ELSE
*         K=(9-KL-1)/2
*         L=(9-KL+2)/2
*       END IF
*       DO 51 IJ=1,KL
*         IF(IJ.LE.3) THEN
*           I=IJ
*           J=IJ
*         ELSE
*           I=(9-IJ-1)/2
*           J=(9-IJ+2)/2
*         END IF
*         WRITE(LU,'(A,4I1,A,2E14.6)') ' a''',I,J,K,L,'=',
*    *      D(I,J,K,L,1)*E(7)+D(I,J,K,L,2)*E(8)+D(I,J,K,L,3)*E(9)
*  51   CONTINUE
*  52 CONTINUE
C
      CLOSE(LU)
      WRITE(*,'(A)') '+TIAXIS: Done.                 '
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'eigen.for'
C     eigen.for
C
C=======================================================================
C