MODULE GRL_MAGALI_UM_UTZ_ANISOTROPY
IMPLICIT NONE
CONTAINS

SUBROUTINE voigt_averaging(size3,eul,odf,eul_2,odf_2,perc,Sij0,Sij0_2,Sij)
!!!Voigt averaging routine for a two-phase aggregate based on Mainprice (1990). 
implicit none
integer,intent(in)::size3
double precision,dimension(3,size3),intent(in)::eul,eul_2
double precision,dimension(size3),intent(in)::odf,odf_2
double precision,dimension(6,6),intent(in)::Sij0,Sij0_2
double precision,dimension(6,6),intent(out)::Sij
double precision,dimension(3,3,3,3)::S0ijkl,Sijkl,Sijkl2
double precision,dimension(:),allocatable::phi1,phi2,theta
double precision,dimension(:,:,:),allocatable::g
double precision, intent(in)::perc
integer :: i,j,k,l,m,n,o,p,nu
!!tensor of indices
INTEGER,dimension(3,3)::ijkl
integer,dimension(6)::l1,l2
allocate(phi1(size3),phi2(size3),theta(size3))
allocate(g(3,3,size3))
phi1=0d0;phi2=0d0;theta=0d0
g=0d0
sijkl2=0d0;sijkl=0d0;s0ijkl=0d0;sij=0d0
ijkl=0;l1=0;l2=0
!!! settensor of indices
ijkl(1,1)=1;ijkl(1,2)=6;ijkl(1,3)=5
ijkl(2,1)=6;ijkl(2,2)=2;ijkl(2,3)=4
ijkl(3,1)=5;ijkl(3,2)=4;ijkl(3,3)=3

l1(1)=1;l1(2)=2;l1(3)=3
l1(4)=2;l1(5)=3;l1(6)=1
l2(1)=1;l2(2)=2;l2(3)=3
l2(4)=3;l2(5)=1;l2(6)=2

do i=1,size3
 phi1(i)=(eul(1,i)*4d0*atan(1d0))/180.d0
 theta(i)=(eul(2,i)*4d0*atan(1d0))/180.d0
 phi2(i)=(eul(3,i)*4d0*atan(1d0))/180.d0
end do

!! calculate direction cosines

do i=1,size3
g(1,1,i)=cos(phi2(i))*cos(phi1(i))-cos(theta(i))*sin(phi1(i))*sin(phi2(i))
g(1,2,i)=cos(phi2(i))*sin(phi1(i))+cos(theta(i))*cos(phi1(i))*sin(phi2(i))
g(1,3,i)=sin(phi2(i))*sin(theta(i))

g(2,1,i)=-sin(phi2(i))*cos(phi1(i))-cos(theta(i))*sin(phi1(i))*cos(phi2(i))
g(2,2,i)=-sin(phi2(i))*sin(phi1(i))+cos(theta(i))*cos(phi1(i))*cos(phi2(i))
g(2,3,i)=cos(phi2(i))*sin(theta(i))

g(3,1,i)=sin(theta(i))*sin(phi1(i))
g(3,2,i)=-sin(theta(i))*cos(phi1(i))
g(3,3,i)=cos(theta(i))
end do

!!! Single crystal elastic tensor
do i=1,3;do j=1,3; do k=1,3; do l=1,3
S0ijkl(i,j,k,l)=Sij0(ijkl(i,j),ijkl(k,l))
end do; end do; end do; end do

!! Voigt average the shit out of these tiny crystals
do nu=1,size3

Sijkl2=0d0
 do i=1,3;do j=1,3; do k=1,3; do l=1,3
  do m=1,3;do n=1,3; do o=1,3; do p=1,3
  Sijkl2(i,j,k,l)=Sijkl2(i,j,k,l)+&
  g(m,i,nu)*g(n,j,nu)*g(o,k,nu)*g(p,l,nu)*S0ijkl(m,n,o,p)
  end do; end do; end do; end do
  Sijkl(i,j,k,l)=Sijkl(i,j,k,l)+Sijkl2(i,j,k,l)*odf(nu)*perc
  
end do; end do; end do; end do

end do


!!! Voigt average of the second phase


do i=1,size3
 phi1(i)=(eul_2(1,i)*4d0*atan(1d0))/180.d0
 theta(i)=(eul_2(2,i)*4d0*atan(1d0))/180.d0
 phi2(i)=(eul_2(3,i)*4d0*atan(1d0))/180.d0
end do

!! calculate direction cosines

do i=1,size3
g(1,1,i)=cos(phi2(i))*cos(phi1(i))-cos(theta(i))*sin(phi1(i))*sin(phi2(i))
g(1,2,i)=cos(phi2(i))*sin(phi1(i))+cos(theta(i))*cos(phi1(i))*sin(phi2(i))
g(1,3,i)=sin(phi2(i))*sin(theta(i))

g(2,1,i)=-sin(phi2(i))*cos(phi1(i))-cos(theta(i))*sin(phi1(i))*cos(phi2(i))
g(2,2,i)=-sin(phi2(i))*sin(phi1(i))+cos(theta(i))*cos(phi1(i))*cos(phi2(i))
g(2,3,i)=cos(phi2(i))*sin(theta(i))

g(3,1,i)=sin(theta(i))*sin(phi1(i))
g(3,2,i)=-sin(theta(i))*cos(phi1(i))
g(3,3,i)=cos(theta(i))
end do

!!! Single crystal elastic tensor
do i=1,3;do j=1,3; do k=1,3; do l=1,3
S0ijkl(i,j,k,l)=Sij0_2(ijkl(i,j),ijkl(k,l))
end do; end do; end do; end do

!! Voigt average the shit out of these tiny crystals
do nu=1,size3

Sijkl2=0d0
 do i=1,3;do j=1,3; do k=1,3; do l=1,3
  do m=1,3;do n=1,3; do o=1,3; do p=1,3
  Sijkl2(i,j,k,l)=Sijkl2(i,j,k,l)+&
  g(m,i,nu)*g(n,j,nu)*g(o,k,nu)*g(p,l,nu)*S0ijkl(m,n,o,p)
  end do; end do; end do; end do
  Sijkl(i,j,k,l)=Sijkl(i,j,k,l)+Sijkl2(i,j,k,l)*odf_2(nu)*(1.d0-perc)
  
end do; end do; end do; end do

end do


!!! Return 6x6 stiffnes matrix

do i=1,6;do j=1,6
 Sij(i,j)=Sijkl(l1(i),l2(i),l1(j),l2(j))
end do; end do

deallocate(phi1,phi2,theta,g)
return
end subroutine voigt_averaging

SUBROUTINE ol_single_xtal_calc(P,T,sij0)
IMPLICIT NONE
DOUBLE PRECISION :: p,T
DOUBLE PRECISION,DIMENSION(6,6) :: sij0
sij0=0d0

sij0(1,1) = 320.2 + 7.5*p - 0.5*0.24*p**2 - 0.033*(T-300.0)
sij0(2,2) = 196.5 + 5.5*p -0.5*0.089*p**2 - 0.028*(T-300.0)
sij0(3,3) = 232.3 + 5.4*p -0.5*0.08*p**2 - 0.032*(T-300.0)

sij0(1,2) = 71.0 + 3.8*p +0.5*0.077*p**2- 0.015*(T-300.0)
sij0(1,3) = 71.0 + 3.8*p + 0.5*0.01*p**2- 0.009*(T-300.0)
sij0(2,3) = 76.0 + 3.5*p -0.5*0.059*p**2 - 0.013*(T-300.0)

sij0(4,4) = 63.2 + 2.0*p -0.5*0.099*p**2 - 0.010*(T-300.0)
sij0(5,5) = 76.5 + 2.0*p -0.5*0.14*p**2 - 0.010*(T-300.0)
sij0(6,6) = 79.7 + 2.4*p -0.5*0.122*p**2 - 0.019*(T-300.0)

sij0(2,1)=sij0(1,2)
sij0(3,1)=sij0(1,3)
sij0(3,2)=sij0(2,3)

RETURN
END SUBROUTINE ol_single_xtal_calc


SUBROUTINE wd_single_xtal_calc_wet(P,T,sij0)
IMPLICIT NONE
DOUBLE PRECISION :: p, T
DOUBLE PRECISION,DIMENSION(6,6) :: sij0
sij0=0d0

sij0(1,1) = 341.0 + 7.48*p - 0.5*0.18*p**2 - 0.029*(T-300.0)
sij0(2,2) = 357.0 + 5.66*p - 0.03*(T-300.0)
sij0(3,3) = 244.0 + 8.8*p - 0.5*0.19*p**2 - 0.037*(T-300.0)

sij0(1,2) = 75.72 + 3.68*p  - 0.0046*(T-300.0)
sij0(1,3) = 99.76 + 2.36*p - 0.0083*(T-300.0)
sij0(2,3) = 101.88+ 3.7*p - 0.5*0.0*p**2 - 0.012*(T-300.0)

sij0(4,4) = 106.0 + 0.91*p - 0.0169*(T-300.0)
sij0(5,5) = 109.0 + 0.86*p - 0.011*(T-300.0)
sij0(6,6) = 90.5  + 2.37*p - 0.5*0.034*p**2 - 0.0135*(T-300.0)

sij0(2,1)=sij0(1,2)
sij0(3,1)=sij0(1,3)
sij0(3,2)=sij0(2,3)
RETURN
END SUBROUTINE wd_single_xtal_calc_wet


SUBROUTINE wd_single_xtal_calc_dry(P,T,sij0)
IMPLICIT NONE
DOUBLE PRECISION :: p, T
DOUBLE PRECISION,DIMENSION(6,6) :: sij0
sij0=0d0

sij0(1,1)=335.59+7.12*(p)-0.5*0.09*(p)**2-0.036*(T-300.0)
sij0(2,2)=289.25+7.6*(p)-0.5*0.1*p**2-0.041*(T-300.0)
sij0(3,3)=310.71+7.2*(p)-0.5*0.11*p**2-0.034*(T-300.0)

sij0(1,2)=87.81+3.26*(p)-0.0085*(T-300.0)
sij0(1,3)=91.95+3.25*(p)-0.0087*(T-300.0)
sij0(2,3)=96.47+3.7*(p)-0.012*(T-300.0)

sij0(4,4)=104.91+1.45*(p)-0.5*0.035*p**2-0.014*(T-300.0)
sij0(5,5)=114.96+1.51*(p)-0.5*0.028*p**2-0.013*(T-300.0)
sij0(6,6)=107.51+1.3*(p)-0.5*0.029*(p)**2-0.0123*(T-300.0)

sij0(2,1)=sij0(1,2)
sij0(3,1)=sij0(1,3)
sij0(3,2)=sij0(2,3)

RETURN
END SUBROUTINE wd_single_xtal_calc_dry

subroutine rw_single_xtal_calc(P,T,sij0)
implicit none
double precision :: p, T
double precision,dimension(6,6)::sij0
sij0=0d0

sij0(1,1) = 329.0 + 6.20*p - 0.049*(T-300.0)
sij0(2,2) = sij0(1,1)
sij0(3,3) = sij0(1,1)

sij0(1,2) = 118.0 + 0.8*p - 0.007*(T-300.0)
sij0(1,3) = sij0(1,2)
sij0(2,3) = sij0(1,2)

sij0(4,4) = 130.0 + 2.8*p - 0.013*(T-300.0)
sij0(5,5) = sij0(4,4)
sij0(6,6) = sij0(4,4)

sij0(2,1)=sij0(1,2)
sij0(3,1)=sij0(1,3)
sij0(3,2)=sij0(2,3)

RETURN
END SUBROUTINE rw_single_xtal_calc


SUBROUTINE gt_single_xtal_calc(P,T,sij0)
IMPLICIT NONE
DOUBLE PRECISION :: p, T
DOUBLE PRECISION,DIMENSION(6,6)::sij0
sij0=0d0

sij0(1,1) = 299.10 + 6.54*p - 0.0305*(T-300.0)
sij0(2,2) = sij0(1,1)
sij0(3,3) = sij0(1,1)

sij0(1,2) = 106.70 + 2.87*p - 0.0058*(T-300.0)
sij0(1,3) = sij0(1,2)
sij0(2,3) = sij0(1,2)

sij0(4,4) = 93.7 + 1.72*p - 0.0071*(T-300.0)
sij0(5,5) = sij0(4,4)
sij0(6,6) = sij0(4,4)

sij0(2,1)=sij0(1,2)
sij0(3,1)=sij0(1,3)
sij0(3,2)=sij0(2,3)

RETURN
END SUBROUTINE gt_single_xtal_calc

SUBROUTINE calc_crss_ol(p,tau)
IMPLICIT NONE
DOUBLE PRECISION,DIMENSION(7) :: tau
DOUBLE PRECISION :: p

!(010)[100]
tau(1)=0.0026*p**2-0.0177*p+0.1882

!(010)[001]
tau(2)=0.0006*p**2-0.0218*p+0.3397

!(100)[001]
tau(3)=0.0017*p**2+0.001*p+0.1518

!(001)[100]
tau(4)=0.002*p**2+0.0014*p+0.1837

!(021)[100]
tau(5)=0.0026*p**2-0.0177*p+0.1882

!(110)[001]
tau(6)=0.0017*p**2-0.0655*p+1.0186

!fake slip system {111}<110> + {110}<111>
tau(7)=1.0 !!!

RETURN
END SUBROUTINE calc_crss_ol

END MODULE GRL_MAGALI_UM_UTZ_ANISOTROPY
