217 lines
8.0 KiB
Fortran
217 lines
8.0 KiB
Fortran
module math
|
|
|
|
implicit none
|
|
|
|
contains
|
|
|
|
function math_axisAngleToR(axis,omega) result(math_axisAngleToR1)
|
|
|
|
implicit none
|
|
|
|
real, dimension(3), intent(in) :: axis
|
|
real, intent(in) :: omega
|
|
real, dimension(3) :: n
|
|
real :: norm,s,c,c1
|
|
|
|
real, dimension(3,3), parameter :: &
|
|
I3 = real(reshape([&
|
|
1, 0, 0, &
|
|
0, 1, 0, &
|
|
0, 0, 1 &
|
|
],shape(I3))) !< 3x3 Identity
|
|
|
|
real, dimension(3,3) :: math_axisAngleToR1
|
|
|
|
norm = norm2(axis)
|
|
wellDefined: if (norm > 1.0e-8) then
|
|
n = axis/norm ! normalize axis to be sure
|
|
|
|
s = sin(omega)
|
|
c = cos(omega)
|
|
c1 = 1.0 - c
|
|
|
|
math_axisAngleToR1(1,1) = c + c1*n(1)**2.0
|
|
math_axisAngleToR1(1,2) = c1*n(1)*n(2) - s*n(3)
|
|
math_axisAngleToR1(1,3) = c1*n(1)*n(3) + s*n(2)
|
|
|
|
math_axisAngleToR1(2,1) = c1*n(1)*n(2) + s*n(3)
|
|
math_axisAngleToR1(2,2) = c + c1*n(2)**2.0
|
|
math_axisAngleToR1(2,3) = c1*n(2)*n(3) - s*n(1)
|
|
|
|
math_axisAngleToR1(3,1) = c1*n(1)*n(3) - s*n(2)
|
|
math_axisAngleToR1(3,2) = c1*n(2)*n(3) + s*n(1)
|
|
math_axisAngleToR1(3,3) = c + c1*n(3)**2.0
|
|
else wellDefined
|
|
math_axisAngleToR1 = I3
|
|
endif wellDefined
|
|
|
|
end function
|
|
|
|
end module math
|
|
|
|
|
|
program corresponcence_matrix
|
|
use math
|
|
implicit none
|
|
|
|
integer, dimension(4) :: &
|
|
active = [6,6,6,6], & !< number of active twin systems
|
|
potential = [6,6,6,6] !< all the potential twin systems
|
|
|
|
real, dimension(3) :: &
|
|
direction, normal
|
|
|
|
real, dimension(3,24) :: normal_vector, direction_vector
|
|
|
|
real, dimension(3,3,24) :: SchmidMatrix, corresponcenceMatrix
|
|
|
|
real, dimension(24) :: characteristicShear
|
|
|
|
real :: cOverA = 1.6235
|
|
|
|
real :: pi = 3.14159274
|
|
|
|
real, dimension(8,24) :: &
|
|
system = reshape(real([&
|
|
! <-10.1>{10.2} systems, shear = (3-(c/a)^2)/(sqrt(3) c/a)
|
|
! tension in Co, Mg, Zr, Ti, and Be; compression in Cd and Zn
|
|
-1, 0, 1, 1, 1, 0, -1, 2, & !
|
|
0, -1, 1, 1, 0, 1, -1, 2, &
|
|
1, -1, 0, 1, -1, 1, 0, 2, &
|
|
1, 0, -1, 1, -1, 0, 1, 2, &
|
|
0, 1, -1, 1, 0, -1, 1, 2, &
|
|
-1, 1, 0, 1, 1, -1, 0, 2, &
|
|
! <11.6>{-1-1.1} systems, shear = 1/(c/a)
|
|
! tension in Co, Re, and Zr
|
|
-1, -1, 2, 6, 1, 1, -2, 1, &
|
|
1, -2, 1, 6, -1, 2, -1, 1, &
|
|
2, -1, -1, 6, -2, 1, 1, 1, &
|
|
1, 1, -2, 6, -1, -1, 2, 1, &
|
|
-1, 2, -1, 6, 1, -2, 1, 1, &
|
|
-2, 1, 1, 6, 2, -1, -1, 1, &
|
|
! <10.-2>{10.1} systems, shear = (4(c/a)^2-9)/(4 sqrt(3) c/a)
|
|
! compression in Mg
|
|
1, 0, -1, -2, 1, 0, -1, 1, &
|
|
0, 1, -1, -2, 0, 1, -1, 1, &
|
|
-1, 1, 0, -2, -1, 1, 0, 1, &
|
|
-1, 0, 1, -2, -1, 0, 1, 1, &
|
|
0, -1, 1, -2, 0, -1, 1, 1, &
|
|
1, -1, 0, -2, 1, -1, 0, 1, &
|
|
! <11.-3>{11.2} systems, shear = 2((c/a)^2-2)/(3 c/a)
|
|
! compression in Ti and Zr
|
|
1, 1, -2, -3, 1, 1, -2, 2, &
|
|
-1, 2, -1, -3, -1, 2, -1, 2, &
|
|
-2, 1, 1, -3, -2, 1, 1, 2, &
|
|
-1, -1, 2, -3, -1, -1, 2, 2, &
|
|
1, -2, 1, -3, 1, -2, 1, 2, &
|
|
2, -1, -1, -3, 2, -1, -1, 2 &
|
|
]),shape(system))
|
|
|
|
real, dimension(3,3), parameter :: &
|
|
I3 = real(reshape([&
|
|
1, 0, 0, &
|
|
0, 1, 0, &
|
|
0, 0, 1 &
|
|
],shape(I3))) !< 3x3 Identity
|
|
|
|
integer :: &
|
|
a, & !< index of active system
|
|
p, & !< index in potential system matrix
|
|
f, & !< index of my family
|
|
s, & !< index of my system in current family
|
|
f1, s1, e1, i, j, k !< indices for similar loops
|
|
|
|
|
|
!--------------------------------------------------------------------
|
|
!> Normal vector to twin plane and direction vector of the twin
|
|
!--------------------------------------------------------------------
|
|
|
|
a = 0
|
|
do f = 1, size(active,1) !< Active Twin Modes
|
|
do s = 1, active(f) !< Active twin systems
|
|
|
|
a = a + 1
|
|
p = sum(potential(1:f-1))+s
|
|
|
|
direction = [ system(1,p)*1.5, &
|
|
(system(1,p)+2.0*system(2,p))*sqrt(0.75), &
|
|
system(4,p)*cOverA ] ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(p/a)])
|
|
|
|
normal = [ system(5,p), &
|
|
(system(5,p)+2.0*system(6,p))/sqrt(3.0), &
|
|
system(8,p)/cOverA ] ! plane (hkil)->(h (h+2k)/sqrt(3) l/(p/a))
|
|
|
|
|
|
|
|
normal_vector(1:3,a) = normal /norm2(normal)
|
|
direction_vector(1:3,a) = direction / norm2(direction)
|
|
|
|
end do
|
|
end do
|
|
|
|
!--------------------------------------------------------------------
|
|
!> Magnitude of Characteristic shear for twinning modes
|
|
!--------------------------------------------------------------------
|
|
|
|
do f1 = 1,size(active,1) !< Active twin modes
|
|
s1 = sum(active(:f1-1)) + 1
|
|
e1 = sum(active(:f1))
|
|
select case(f1)
|
|
case (1)
|
|
characteristicShear(s1:e1) = (3.0-cOverA**2)/sqrt(3.0)/cOverA
|
|
case (2)
|
|
characteristicShear(s1:e1) = 1.0/cOverA
|
|
case (3)
|
|
characteristicShear(s1:e1) = (4.0*cOverA**2-9.0)/sqrt(48.0)/cOverA
|
|
case (4)
|
|
characteristicShear(s1:e1) = 2.0*(cOverA**2-2.0)/3.0/cOverA
|
|
end select
|
|
enddo
|
|
|
|
!> Write results for characteristic shear
|
|
|
|
write(6,*)'characteristic shear, for [1, 0, -1, 1],(-1, 0, 1, 2)'
|
|
write(6,*)characteristicShear(4)
|
|
write(6,*)'characteristic shear, for [-1, -1, 2, 6],(1, 1, -2, 1)'
|
|
write(6,*)characteristicShear(7)
|
|
write(6,*)'characteristic shear, for [1, 0, -1, -2],(1, 0, -1, 1)'
|
|
write(6,*)characteristicShear(13)
|
|
write(6,*)'characteristic shear, for [1, 1, -2, -3],(1, 1, -2, 2)'
|
|
write(6,*)characteristicShear(19)
|
|
|
|
!--------------------------------------------------------------------
|
|
!> SchmidMatrix = Outer product of direction and normal vectors.
|
|
!--------------------------------------------------------------------
|
|
|
|
do i = 1, sum(active)
|
|
forall(j=1:3, k=1:3) SchmidMatrix(j,k,i) = direction_vector(j,i) * normal_vector(k,i)
|
|
enddo
|
|
|
|
!--------------------------------------------------------------------
|
|
!> Correspondence Matrix = Reorientation * Shear
|
|
!--------------------------------------------------------------------
|
|
|
|
do i = 1, sum(active)
|
|
corresponcenceMatrix(1:3,1:3,i) = matmul(math_axisAngleToR(normal_vector(1:3,i),pi),&
|
|
I3+characteristicShear(i)*SchmidMatrix(1:3,1:3,i))
|
|
|
|
enddo
|
|
|
|
!> Write results for Correspondence Matrix
|
|
|
|
write(6,*)'correspondence matrix for [1, 0, -1, 1],(-1, 0, 1, 2)'
|
|
write(6,*)corresponcenceMatrix(1:3,1:3,4)
|
|
write(6,*)'oxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxox'
|
|
write(6,*)'correspondence matrix for [-1, -1, 2, 6],(1, 1, -2, 1)'
|
|
write(6,*)corresponcenceMatrix(1:3,1:3,7)
|
|
write(6,*)'oxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxox'
|
|
write(6,*)'correspondence matrix for [1, 0, -1, -2],(1, 0, -1, 1)'
|
|
write(6,*)corresponcenceMatrix(1:3,1:3,13)
|
|
write(6,*)'oxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxoxox'
|
|
write(6,*)'correspondence matrix for [1, 1, -2, -3],(1, 1, -2, 2)'
|
|
write(6,*)corresponcenceMatrix(1:3,1:3,19)
|
|
|
|
|
|
end program corresponcence_matrix
|
|
|