223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

corrected math_misorientation

This commit is contained in:
Christoph Kords 2010-01-05 16:02:00 +00:00
parent ac080e0b52
commit aee0721ab2
1 changed files with 20 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
code/math.f90
View File

@ -302,37 +302,27 @@ subroutine math_misorientation(axis, angle, rot, ori1, ori2, symmetryType)
! apply symmetry operation to 1st orientation ! apply symmetry operation to 1st orientation
do s1 = 1,NsymOperations do s1 = 1,NsymOperations
ori1sym = math_mul33x33(mySymOperations(:,:,s1), ori1) ori1sym = math_mul33x33(mySymOperations(:,:,s1), ori1)
! calculate possible net rotation
this_rot = math_mul33x33(ori1sym, transpose(ori2))
! apply symmetry operation to 2nd orientation ! store smallest misorientation for an axis inside standard orientation triangle
do s2 = 1,NsymOperations ! calculate rotation axis
ori2sym = math_mul33x33(mySymOperations(:,:,s2), ori2) invNorm = ( (this_rot(1,2) - this_rot(2,1))**2.0_pReal &
+ (this_rot(2,3) - this_rot(3,2))**2.0_pReal &
! calculate 2 possible net rotations + (this_rot(3,1) - this_rot(1,3))**2.0_pReal ) ** (-0.5_pReal)
this_rot = math_mul33x33(ori1sym, transpose(ori2sym)) this_axis(1) = (this_rot(2,3) - this_rot(3,2)) * invNorm
this_axis(2) = (this_rot(3,1) - this_rot(1,3)) * invNorm
! store smallest misorientation for an axis inside standard orientation triangle this_axis(3) = (this_rot(1,2) - this_rot(2,1)) * invNorm
! calculate rotation axis
invNorm = ( (this_rot(1,2) - this_rot(2,1))**2.0_pReal & ! calculate rotation angle
+ (this_rot(2,3) - this_rot(3,2))**2.0_pReal & this_angle = abs(0.5_pReal * pi - asin(0.4999999_pReal * (this_rot(1,1) + this_rot(2,2) + this_rot(3,3) - 1.0_pReal)))
+ (this_rot(3,1) - this_rot(1,3))**2.0_pReal ) ** (-0.5_pReal)
this_axis(1) = (this_rot(2,3) - this_rot(3,2)) * invNorm if (abs(this_angle) < angle) then
this_axis(2) = (this_rot(3,1) - this_rot(1,3)) * invNorm angle = this_angle
this_axis(3) = (this_rot(1,2) - this_rot(2,1)) * invNorm rot = this_rot
axis = this_axis
if ( (this_axis(1) >= 0.0_pReal) .and. (this_axis(2) >= 0.0_pReal) .and. (this_axis(3) >= 0.0_pReal) & endif
.and. (this_axis(1) <= this_axis(2)) .and. (this_axis(2) <= this_axis(3)) ) then
! calculate rotation angle
this_angle = acos(0.5_pReal * (this_rot(1,1) + this_rot(2,2) + this_rot(3,3) - 1.0_pReal))
if (abs(this_angle) < angle) then
angle = this_angle
rot = this_rot
axis = this_axis
endif
endif
enddo
enddo enddo
! convert angle to degrees ! convert angle to degrees