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