C
C Program ANIROT to rotate the tensor of elastic moduli
C of an anisotropic elastic medium
C
C Version: 7.30
C Date: 2015, December 2
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 Names of the input and output files:
C     ANIIN='string'... Name of the input file containing the elastic
C             moduli before rotation.
C             Description of file ANIIN.
C             Default: ANIIN='anirot.dat'
C     ANIOUT='string'... Name of the output file containing the elastic
C             moduli after rotation.
C             Description of file ANIOUT.
C             Default: ANIOUT='anirot.out'
C Data describing the rotation:
C     ANIDEG=real... Rotation angle in degrees (not in radians).
C     ANIE1=real, ANIE2=real, ANIE3=real... Three components of the
C             axial vector around which the anisotropic elastic medium
C             is rotated.  The length of the axial vector does not
C             matter.  The medium is rotated according to the right-hand
C             rule in a right-handed coordinate system.  The medium is
C             thus rotated counter-clockwise when looking against the
C             axial vector.
C             Transformation matrix E_ij:
C               The i-th component of the j-th rotated local Cartesian
C               basis vector, expressed with respect to the local
C               Cartesian basis vectors before rotation, is
C                 E_ij = cos(ANIDEG) delta_ij
C                      + [1-cos(ANIDEG)] E_i/E E_j/E
C                      - sin(ANIDEG) epsilon_ijk E_k/E
C               where (E_1,E_2,E_3)=(ANIE1,ANIE2,ANIE3),
C                 E = sqrt(E_k E_k)
C               is the Cartesian length of the rotation vector,
C               delta_ij is the Kronecker delta, and
C               epsilon_ijk is completely skew Levi-Civita's symbol,
C               with epsilon_123=epsilon_231=epsilon_312=1.
C             Transformation of the elastic moduli:
C                 B_ijkl = E_ia E_ib E_ic E_id A_abcd
C               where the Einstein summation is applied,
C               A_ijkl are the elastic moduli before rotation, and
C               B_ijkl are the elastic moduli after rotation.
C             Default:
C               ANIE1=0., ANIE2=0., ANIE3=0. (no rotation)
C Optional parameters specifying the form of the quantities
C written in the output formatted file:
C     MINDIG,MAXDIG=positive integers ... See the description in file
C             forms.for.
C             Recommended value of MAXDIG is 7 or 8 for this program.
C
C                                                   
C Input file ANIIN with the elastic moduli before rotation:
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 Output file ANIOUT with the elastic moduli after rotation:
C (1) Elastic moduli
C     B1111,B1122,B1133,B1123,B1113,B1112
C (2) Elastic moduli
C           B2222,B2233,B2223,B2213,B2212
C (3) Elastic moduli
C                 B3333,B3323,B3313,B3312
C (4) Elastic moduli
C                       B2323,B2313,B2312
C (5) Elastic moduli
C                             B1313,B1312
C (6) Elastic modulus
C                                   B1212
C
C=======================================================================
C
C     Filenames and parameters:
      CHARACTER*80 FSEP,FIN,FOUT
      INTEGER LU1
      PARAMETER (LU1=1)
C
C     Other variables:
      REAL DEG,E1,E2,E3
      REAL RAD,C,S,BMIN,BMAX,AUX1,AUX2,AUX3,AUX4
      REAL AV(6,6),A(3,3,3,3),B(3,3,3,3),E(3,3)
      INTEGER I12,I34,I1,I2,I3,I4,J1,J2,J3,J4
      CHARACTER*17 FORMAT
C
C.......................................................................
C
C     Reading main input data:
      WRITE(*,'(A)') '+ANIROT: Enter input filename: '
      FSEP=' '
      READ (*,*) FSEP
      IF(FSEP.EQ.' ') THEN
C       ANIROT-01
        CALL ERROR('ANIROT-01: No input 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(LU1,FSEP)
      WRITE(*,'(A)') '+ANIROT: Working...            '
C
C     Reading input SEP parameters:
      CALL RSEP3T('ANIIN' ,FIN ,'anirot.dat')
      CALL RSEP3T('ANIOUT',FOUT,'anirot.out')
      CALL RSEP3R('ANIDEG',DEG,0.)
      CALL RSEP3R('ANIE1' ,E1 ,0.)
      CALL RSEP3R('ANIE2' ,E2 ,0.)
      CALL RSEP3R('ANIE3' ,E3 ,0.)
C
C     Transformation matrix E_ij:
      RAD=3.141592654*DEG/180.
      AUX1=SQRT(E1**2+E2**2+E3**2)
      IF(RAD.EQ.0..OR.AUX1.EQ.0.) THEN
        E(1,1)=1.
        E(2,1)=0.
        E(3,1)=0.
        E(1,2)=0.
        E(2,2)=1.
        E(3,2)=0.
        E(1,3)=0.
        E(2,3)=0.
        E(3,3)=1.
      ELSE
        E1=E1/AUX1
        E2=E2/AUX1
        E3=E3/AUX1
        C=COS(RAD)
        S=SIN(RAD)
        AUX1=1.-C
        E(1,1)=AUX1*E1*E1+C
        E(2,1)=AUX1*E2*E1+S*E3
        E(3,1)=AUX1*E3*E1-S*E2
        E(1,2)=AUX1*E1*E2-S*E3
        E(2,2)=AUX1*E2*E2+C
        E(3,2)=AUX1*E3*E2+S*E1
        E(1,3)=AUX1*E1*E3+S*E2
        E(2,3)=AUX1*E2*E3-S*E1
        E(3,3)=AUX1*E3*E3+C
      END IF
C
C     Reading the elastic moduli in Voight notation:
      OPEN(LU1,FILE=FIN,STATUS='OLD')
      READ(LU1,*) AV(1,1),AV(1,2),AV(1,3),AV(1,4),AV(1,5),AV(1,6)
     *                   ,AV(2,2),AV(2,3),AV(2,4),AV(2,5),AV(2,6)
     *                           ,AV(3,3),AV(3,4),AV(3,5),AV(3,6)
     *                                   ,AV(4,4),AV(4,5),AV(4,6)
     *                                           ,AV(5,5),AV(5,6)
     *                                                   ,AV(6,6)
      CLOSE(LU1)
C
C     Transforming the Voight notation to the tensor:
      DO 14 I4=1,3
        DO 13 I3=1,3
          IF(I3.EQ.I4) THEN
            I34=I3
          ELSE
            I34=9-I3-I4
          END IF
          DO 12 I2=1,3
            DO 11 I1=1,3
              IF(I1.EQ.I2) THEN
                I12=I1
              ELSE
                I12=9-I1-I2
              END IF
              IF(I12.LE.I34) THEN
                A(I1,I2,I3,I4)=AV(I12,I34)
              ELSE
                A(I1,I2,I3,I4)=AV(I34,I12)
              END IF
   11        CONTINUE
   12     CONTINUE
   13   CONTINUE
   14 CONTINUE
C
C     Transformation of the elastic moduli:
      BMIN= 1.0E20
      BMAX=-1.0E20
      DO 34 J4=1,3
        DO 33 J3=1,3
          DO 32 J2=1,3
            DO 31 J1=1,3
              AUX1=0.
              DO 24 I4=1,3
                AUX4=E(J4,I4)
                DO 23 I3=1,3
                  AUX3=AUX4*E(J3,I3)
                  DO 22 I2=1,3
                    AUX2=AUX3*E(J2,I2)
                    DO 21 I1=1,3
                      AUX1=AUX1+AUX2*E(J1,I1)*A(I1,I2,I3,I4)
   21               CONTINUE
   22             CONTINUE
   23           CONTINUE
   24         CONTINUE
              B(J1,J2,J3,J4)=AUX1
              BMIN=AMIN1(BMIN,AUX1)
              BMAX=AMAX1(BMAX,AUX1)
   31       CONTINUE
   32     CONTINUE
   33   CONTINUE
   34 CONTINUE
C
C     Setting the output format:
      FORMAT='(    6(F00.0,1X))'
      CALL FORM1(BMIN,BMAX,FORMAT(8:15))
      FORMAT(14:17)= '1X))'
      I1=ICHAR(FORMAT(9:9))-ICHAR('0')
      I1=10*I1+ICHAR(FORMAT(10:10))-ICHAR('0')+1
C
C     Writing the elastic moduli in Voight notation:
      OPEN(LU1,FILE=FOUT)
      WRITE(LU1,FORMAT)
     * B(1,1,1,1),B(1,1,2,2),B(1,1,3,3),B(1,1,2,3),B(1,1,1,3),B(1,1,1,2)
      I2=I1
      I3=I2/10
      I4=I2-I3*10
      FORMAT(2:2)=CHAR(ICHAR('0')+I3)
      FORMAT(3:3)=CHAR(ICHAR('0')+I4)
      FORMAT(4:5)='X,'
      WRITE(LU1,FORMAT)
     *           ,B(2,2,2,2),B(2,2,3,3),B(2,2,2,3),B(2,2,1,3),B(2,2,1,2)
      I2=I2+I1
      I3=I2/10
      I4=I2-I3*10
      FORMAT(2:2)=CHAR(ICHAR('0')+I3)
      FORMAT(3:3)=CHAR(ICHAR('0')+I4)
      WRITE(LU1,FORMAT)
     *                      ,B(3,3,3,3),B(3,3,2,3),B(3,3,1,3),B(3,3,1,2)
      I2=I2+I1
      I3=I2/10
      I4=I2-I3*10
      FORMAT(2:2)=CHAR(ICHAR('0')+I3)
      FORMAT(3:3)=CHAR(ICHAR('0')+I4)
      WRITE(LU1,FORMAT)
     *                                 ,B(2,3,2,3),B(2,3,1,3),B(2,3,1,2)
      I2=I2+I1
      I3=I2/10
      I4=I2-I3*10
      FORMAT(2:2)=CHAR(ICHAR('0')+I3)
      FORMAT(3:3)=CHAR(ICHAR('0')+I4)
      WRITE(LU1,FORMAT)
     *                                            ,B(1,3,1,3),B(1,3,1,2)
      I2=I2+I1
      I3=I2/10
      I4=I2-I3*10
      FORMAT(2:2)=CHAR(ICHAR('0')+I3)
      FORMAT(3:3)=CHAR(ICHAR('0')+I4)
      WRITE(LU1,FORMAT)
     *                                                       ,B(1,2,1,2)
      CLOSE(LU1)
      WRITE(*,'(A)') '+ANIROT: 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 'forms.for'
C     forms.for
C
C=======================================================================
C