module invariants

implicit none

contains

!this function computes the trace of a tensor
!Tr[\sigma]=\sum_{i} \sigma_{ii}

function trace (exx,eyy,ezz,exy,eyz,ezx)
double precision trace
double precision exx,eyy,ezz,exy,eyz,ezx
trace=exx+eyy+ezz
end function trace


!this subroutine takes the 6 components of a tensor
!computes the trace and outputs the deviatoric tensor
!\sigma^d=\sigma-\frac{1}{3}Tr[\sigma] {\bf 1}

subroutine deviatoric (exx,eyy,ezz,exy,eyz,ezx)
implicit none
double precision exx,eyy,ezz,exy,eyz,ezx
double precision trace
trace=(exx+eyy+ezz)/3.d0
exx=exx-trace
eyy=eyy-trace
ezz=ezz-trace
end subroutine deviatoric

!this function computes the second invariant:
!J_2=\frac{1}{2} \sum_{ij} e_{ij}e_{ji}
!in the case of a symmetric tensor

function second_invariant (exx,eyy,ezz,exy,eyz,ezx)
double precision second_invariant 
double precision exx,eyy,ezz,exy,eyz,ezx
second_invariant=(exx**2+eyy**2+ezz**2)/2.d0+exy**2+eyz**2+ezx**2
end function second_invariant


!this function computes the second invariant:
!J_3=\frac{1}{3} \sum_{ijk} e_{ij}e_{jk}e_{ji}
!in the case of a symmetric tensor

function third_invariant (exx,eyy,ezz,exy,eyz,ezx)
double precision third_invariant
double precision exx,eyy,ezz,exy,eyz,ezx
third_invariant=(exx*(     exx**2 + 3.d0*exy**2 + 3.d0*ezx**2)       &
                +eyy*(3.d0*exy**2 +      eyy**2 + 3.d0*eyz**2)       &
                +ezz*(3.d0*ezx**2 + 3.d0*eyz**2 +      ezz**2))/3.d0 &
                +2.d0*exy*eyz*ezx
end function third_invariant


!this function computes the Lode angle
!\theta=-\frac{1}{3}\sin^{-1} \left( \frac{3\sqrt{3}}{2}  J_3 / J2^{3/2}  \right)

function lode_angle (J2d,J3d)
use constants
double precision lode_angle,J2d,J3d 
lode_angle=asin(-1.5d0*sqrt3*J3d/J2d**(1.5d0))/3.d0
end function lode_angle



end module invariants