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.30
C Date: 2016, March 25
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 files:
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     ANITI='string'... Optional name of the output file containing the
C             density-reduced elastic moduli of the reference
C             transversely isotropic medium.
C             This output file has the same form as input file ANIMOD.
C             Default: ANITI=' ' (no output)
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     Variables related to optional reference TI stiffness tensor:
      REAL T(3),Z(3,3),C(3,3),S(3,3),ATI(3,3,3,3)
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
C-----------------------------------------------------------------------
C
C     Optional reference TI stiffness tensor:
C
C     Reference symmetry vector
      T(1)=E(7)
      T(2)=E(8)
      T(3)=E(9)
C
C     Generating transformation matrices Z, C and S
      DO 62 J=1,3
        DO 61 I=1,3
          Z(I,J)=T(I)*T(J)
          IF(I.EQ.J) THEN
            C(I,J)=1.-Z(I,J)
          ELSE
            C(I,J)=  -Z(I,J)
            K=6-I-J
            IF(I.LT.J) THEN
              S(I,J)=T(K)
            ELSE
              S(I,J)=-T(K)
            END IF
            IF(K.EQ.2) THEN
              S(I,J)=-S(I,J)
            END IF
          END IF
   61   CONTINUE
   62 CONTINUE
C
C     Initializing transversely isotropic stiffness tensor ATI
      DO 74 L=1,3
        DO 73 K=1,3
          DO 72 J=1,3
            DO 71 I=1,3
              ATI(I,J,K,L)=0.
   71       CONTINUE
   72     CONTINUE
   73   CONTINUE
   74 CONTINUE
C
C     Creating transversely isotropic stiffness tensor ATI
      CALL TRANSA(1.000,Z,Z,Z,Z,A,ATI)
C
      CALL TRANSA(0.500,Z,Z,C,C,A,ATI)
      CALL TRANSA(0.500,Z,C,Z,C,A,ATI)
      CALL TRANSA(0.500,Z,C,C,Z,A,ATI)
      CALL TRANSA(0.500,C,Z,Z,C,A,ATI)
      CALL TRANSA(0.500,C,Z,C,Z,A,ATI)
      CALL TRANSA(0.500,C,C,Z,Z,A,ATI)
C
      CALL TRANSA(0.500,Z,Z,S,S,A,ATI)
      CALL TRANSA(0.500,Z,S,Z,S,A,ATI)
      CALL TRANSA(0.500,Z,S,S,Z,A,ATI)
      CALL TRANSA(0.500,S,Z,Z,S,A,ATI)
      CALL TRANSA(0.500,S,Z,S,Z,A,ATI)
      CALL TRANSA(0.500,S,S,Z,Z,A,ATI)
C
      CALL TRANSA(0.125,C,C,S,S,A,ATI)
      CALL TRANSA(0.125,C,S,C,S,A,ATI)
      CALL TRANSA(0.125,C,S,S,C,A,ATI)
      CALL TRANSA(0.125,S,C,C,S,A,ATI)
      CALL TRANSA(0.125,S,C,S,C,A,ATI)
      CALL TRANSA(0.125,S,S,C,C,A,ATI)
C
      CALL TRANSA(0.375,C,C,C,C,A,ATI)
      CALL TRANSA(0.375,S,S,S,S,A,ATI)
C
C     Relative Frobenius norm of the difference
      RNORM=0.
      DO 84 L=1,3
        DO 83 K=1,3
          DO 82 J=1,3
            DO 81 I=1,3
              RNORM=RNORM+(ATI(I,J,K,L)-A(I,J,K,L))**2
   81       CONTINUE
   82     CONTINUE
   83   CONTINUE
   84 CONTINUE
      RNORM=SQRT(RNORM/ANORM)
      WRITE(LU,'(A,F8.6)')  'Relative difference=',RNORM
      CLOSE(LU)
C
C     Writing transversely isotropic stiffness tensor ATI
      CALL RSEP3T('ANITI',FILE,' ')
      IF(FILE.EQ.' ') THEN
        GO TO 90
      END IF
      OPEN(LU,FILE=FILE)
      WRITE(LU,'(    6(F9.6,1X))')ATI(1,1,1,1),ATI(1,1,2,2),ATI(1,1,3,3)
     *                           ,ATI(1,1,2,3),ATI(1,1,1,3),ATI(1,1,1,2)
      WRITE(LU,'(10X,5(F9.6,1X))')             ATI(2,2,2,2),ATI(2,2,3,3)
     *                           ,ATI(2,2,2,3),ATI(2,2,1,3),ATI(2,2,1,2)
      WRITE(LU,'(20X,4(F9.6,1X))')                          ATI(3,3,3,3)
     *                           ,ATI(3,3,2,3),ATI(3,3,1,3),ATI(3,3,1,2)
      WRITE(LU,'(30X,3(F9.6,1X))')ATI(2,3,2,3),ATI(2,3,1,3),ATI(2,3,1,2)
      WRITE(LU,'(40X,2(F9.6,1X))')             ATI(1,3,1,3),ATI(1,3,1,2)
      WRITE(LU,'(50X,1(F9.6,1X))')                          ATI(1,2,1,2)
      CLOSE(LU)
C
   90 CONTINUE
      WRITE(*,'(A)') '+TIAXIS: Done.                 '
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE TRANSA(FACT,B1,B2,B3,B4,A1,A2)
      REAL FACT,B1(3,3),B2(3,3),B3(3,3),B4(3,3),A1(3,3,3,3),A2(3,3,3,3)
C
C-----------------------------------------------------------------------
C
      INTEGER I1,I2,I3,I4,J1,J2,J3,J4
      REAL B34,B234,AUX
C
      DO 80 J4=1,3
        DO 70 J3=1,3
          DO 60 J2=1,3
            DO 50 J1=1,3
              AUX=0.
              DO 40 I4=1,3
                DO 30 I3=1,3
                  B34=B3(J3,I3)*B4(J4,I4)
                  DO 20 I2=1,3
                    B234=B2(J2,I2)*B34
                    DO 10 I1=1,3
                      AUX=AUX+B1(J1,I1)*B234*A1(I1,I2,I3,I4)
   10               CONTINUE
   20             CONTINUE
   30           CONTINUE
   40         CONTINUE
              A2(J1,J2,J3,J4)=A2(J1,J2,J3,J4)+FACT*AUX
   50       CONTINUE
   60     CONTINUE
   70   CONTINUE
   80 CONTINUE
C
      RETURN
      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