C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%   Subroutine "trelastensor" is used to convert the       %%%%
C%%%%   elastic-constants tensor from one reference system     %%%%
C%%%%   with Cartesian coordinates to another one xi with      %%%%
C%%%%   transformed coordinates Xi                             %%%%
C%%%%                        REFERENCE                         %%%%
C%%%%   1\ Golesorkhtabar R, Pavone P, Spitaler J, et al.      %%%%
C%%%%      ElaStic : A tool for calculating second-order       %%%%
C%%%%      elastic constants from first principles[J].         %%%%
C%%%%      Computer Physics Communications, 2013, 184(8):      %%%%
C%%%%      1861-1873.                                          %%%%
C%%%%   2\ Marmier A, Lethbridge Z A D, Walton R I, et al.     %%%%
C%%%%      ElAM: A computer program for the analysis and       %%%%
C%%%%      representation of anisotropic elastic properties[J].%%%%
C%%%%      Computer Physics Communications, 2010, 181(12):     %%%%
C%%%%      2102-2115.                                          %%%%
C%%%%   3\ Gaillac R, Pullumbi P, Coudert F X. ELATE: an       %%%%
C%%%%      open-source online application for analysis and     %%%%
C%%%%      visualization of elastic tensors[J]. Journal of     %%%%
C%%%%      Physics: Condensed Matter, 2016, 28(27): 275201.    %%%%      
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
C            >>>   Fourth-order elastic constants   <<<
C      
C       Tensor notation:  11  22  33  23,32  31,13  12,21
C      Voigt's notation:   1   2   3    4      5      6    
C      
C       1\ Cijmn=Cjimn, Cijmn=Cijnm
C       2\ Cijmn=Cpq if p and q are 1, 2, 3 only
C          Cijmn=1/2*Cpq if either p or q are 4, 5, 6 (and the
C          other is 1, 2, or 3)
C          Cijmn=1/4*Cpq if p and q are 4, 5, 6 only
C
C      
C           C1111    C1122    C1133  2*C1123  2*C1131  2*C1112
C           C2211    C2222    C2233  2*C2223  2*C2231  2*C2212
C           C3311    C3322    C3333  2*C3323  2*C2231  2*C3312
C         2*C2311  2*C2322  2*C2333  4*C2323  4*C2331  4*C2312
C         2*C3111  2*C3122  2*C3133  4*C3123  4*C3131  4*C3112
C         2*C1211  2*C1222  2*C1233  4*C1223  4*C1231  4*C1212
C
C                                ==
C
C                   c11  c12  c13  c14  c15  c16
C                   c12  c22  c23  c24  c25  c26
C                   c13  c23  c33  c34  c35  c36
C                   c14  c24  c34  c44  c45  c46
C                   c15  c25  c35  c45  c55  c56
C                   c16  c26  c36  c46  c56  c66
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        subroutine trelastensor( Celasnew, Celas, cosmat )
        double precision Celasnew, Celas, cosmat
        dimension Celasnew(6,6), Celas(6,6), cosmat(3,3)

        integer i,j,k,l
        integer ii,jj,kk,ll
        integer m,n,XX,YY
        double precision CCold, CCnew
        dimension CCold(3,3,3,3), CCnew(3,3,3,3)

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%      sedond-order to fourth-order elastic constants      %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        do i=1,3
           do j=1,3
              XX=i*j
              select case (XX)
                case(1)
                       m=1
                case(4)
                       m=2
                case(9)
                       m=3
                case(6)
                       m=4
                case(3)
                       m=5
                case(2)
                       m=6
              end select

              do k=1,3
                 do l=1,3
                    YY=k*l
                    select case (YY)
                      case(1)
                             n=1
                      case(4)
                             n=2
                      case(9)
                             n=3
                      case(6)
                             n=4
                      case(3)
                             n=5
                      case(2)
                             n=6
                    end select

                    if ( m .LE. 3 .AND. n .LE. 3) THEN
                            CCold(i,j,k,l)=Celas(m,n)
                    endif
                    if ( m .GE. 4 .OR. n .GE. 4) THEN
                            CCold(i,j,k,l)=1.00/2.00*Celas(m,n)
                    endif
                    if ( m .GE. 4 .AND. n .GE. 4) THEN
                            CCold(i,j,k,l)=1.00/4.00*Celas(m,n)
                    endif

                 enddo
              enddo
           enddo
        enddo

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%   convert the elastic-constants tensor from one          %%%%  
C%%%%   reference system with Cartesian coordinates to another %%%%
C%%%%   one "x" with transformed coordinates "X"               %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C       X(i)=sum( x(j)*cosmat(i,j), j=1,3 ), x=1,3

C       cosmat(i,j): the cosine of the angle between the directions
C       of X(i) and x(j)

        do i=1,3
           do j=1,3
              do k=1,3
                 do l=1,3
                    
                    CCnew(i,j,k,l)=0.0000
                    do ii=1,3
                       do jj=1,3
                          do kk=1,3
                             do ll=1,3
                                CCnew(i,j,k,l)=CCnew(i,j,k,l)
     &                          +cosmat(i,ii)*cosmat(j,jj)
     &                          *cosmat(k,kk)*cosmat(l,ll)
     &                          *CCold(ii,jj,kk,ll)
                             enddo
                          enddo
                       enddo
                    enddo 

                 enddo
              enddo
           enddo
        enddo                 

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%      fourth-order to second-order elastic constants      %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        do m=1,6

           select case (m)
             case(1)
                    i=1, j=1
             case(2)
                    i=2, j=2
             case(3)
                    i=3, j=3
             case(4)
                    i=2, j=3
             case(5)
                    i=1, j=3
             case(6)
                    i=1, j=2
           end select

           do n=1,6
              
              select case (n)
                case(1)
                       k=1, l=1
                case(2)
                       k=2, l=2
                case(3)
                       k=3, l=3
                case(4)
                       k=2, l=3
                case(5)
                       k=1, l=3
                case(6)
                       k=1, l=2
              end select              

              if ( m .LE. 3 .AND. n .LE. 3) THEN
                      Celasnew(m,n)=CCnew(i,j,k,l)
              endif
              if ( m .GE. 4 .OR. n .GE. 4) THEN
                      Celasnew(m,n)=2.00*CCnew(i,j,k,l)
              endif
              if ( m .GE. 4 .AND. n .GE. 4) THEN
                      Celasnew(m,n)=4.00*CCnew(i,j,k,l)
              endif

           enddo
        enddo      

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%                   end of trelastensor.f                  %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
