DAMASK_EICMD/src/rotations.f90

1522 lines
54 KiB
Fortran
Raw Normal View History

! ###################################################################
! Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University
2020-01-29 13:55:39 +05:30
! Modified 2017-2020, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
! All rights reserved.
!
2020-04-22 18:56:29 +05:30
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
2020-04-22 18:56:29 +05:30
! - Redistributions of source code must retain the above copyright notice, this list
! of conditions and the following disclaimer.
2020-04-22 18:56:29 +05:30
! - Redistributions in binary form must reproduce the above copyright notice, this
! list of conditions and the following disclaimer in the documentation and/or
! other materials provided with the distribution.
2020-04-22 18:56:29 +05:30
! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
! of its contributors may be used to endorse or promote products derived from
! this software without specific prior written permission.
!
2020-04-22 18:56:29 +05:30
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ###################################################################
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief rotation storage and conversion
2020-04-22 18:56:29 +05:30
!> @details: rotation is internally stored as quaternion. It can be inialized from different
!> representations and also returns itself in different representations.
2019-02-01 14:47:20 +05:30
!
! All methods and naming conventions based on Rowenhorst et al. 2015
2019-02-01 14:47:20 +05:30
! Convention 1: coordinate frames are right-handed
! Convention 2: a rotation angle ω is taken to be positive for a counterclockwise rotation
! when viewing from the end point of the rotation axis towards the origin
! Convention 3: rotations will be interpreted in the passive sense
! Convention 4: Euler angle triplets are implemented using the Bunge convention,
! with the angular ranges as [0, 2π],[0, π],[0, 2π]
! Convention 5: the rotation angle ω is limited to the interval [0, π]
2020-01-14 16:03:18 +05:30
! Convention 6: the real part of a quaternion is positive, Re(q) > 0
! Convention 7: P = -1
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
module rotations
use IO
use math
2020-04-22 18:56:29 +05:30
2019-04-03 16:41:18 +05:30
implicit none
private
2021-01-02 22:27:11 +05:30
real(pReal), parameter :: P = -1.0_pReal !< parameter for orientation conversion.
2022-01-29 20:00:59 +05:30
type, public :: tRotation
2021-01-02 22:27:11 +05:30
real(pReal), dimension(4) :: q
2019-04-03 16:41:18 +05:30
contains
procedure, public :: asQuaternion
2019-09-21 05:48:09 +05:30
procedure, public :: asEulers
procedure, public :: asAxisAngle
procedure, public :: asRodrigues
procedure, public :: asMatrix
2019-04-03 16:41:18 +05:30
!------------------------------------------
2019-12-02 20:52:27 +05:30
procedure, public :: fromQuaternion
2019-09-21 05:48:09 +05:30
procedure, public :: fromEulers
procedure, public :: fromAxisAngle
procedure, public :: fromMatrix
2019-04-03 16:41:18 +05:30
!------------------------------------------
procedure, private :: rotRot__
generic, public :: operator(*) => rotRot__
2019-09-21 05:09:33 +05:30
generic, public :: rotate => rotVector,rotTensor2,rotTensor4
procedure, public :: rotVector
2019-09-20 21:15:23 +05:30
procedure, public :: rotTensor2
2019-09-21 05:09:33 +05:30
procedure, public :: rotTensor4
2021-11-20 19:45:59 +05:30
procedure, public :: rotStiffness
procedure, public :: misorientation
2020-01-14 16:03:18 +05:30
procedure, public :: standardize
2022-01-29 20:00:59 +05:30
end type tRotation
2020-04-22 18:56:29 +05:30
real(pReal), parameter :: &
PREF = sqrt(6.0_pReal/PI), &
A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), &
AP = PI**(2.0_pReal/3.0_pReal), &
SC = A/AP, &
BETA = A/2.0_pReal, &
R1 = (3.0_pReal*PI/4.0_pReal)**(1.0_pReal/3.0_pReal), &
R2 = sqrt(2.0_pReal), &
PI12 = PI/12.0_pReal, &
PREK = R1 * 2.0_pReal**(1.0_pReal/4.0_pReal)/BETA
2019-09-22 19:23:03 +05:30
public :: &
2020-01-29 13:55:39 +05:30
rotations_init, &
eu2om
contains
2019-09-22 19:23:03 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-14 00:58:53 +05:30
!> @brief Do self test.
2019-09-22 19:23:03 +05:30
!--------------------------------------------------------------------------------------------------
subroutine rotations_init
2019-09-22 19:23:03 +05:30
print'(/,1x,a)', '<<<+- rotations init -+>>>'; flush(IO_STDOUT)
2020-05-21 14:15:52 +05:30
print'(/,1x,a)', 'D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015'
print'( 1x,a)', 'https://doi.org/10.1088/0965-0393/23/8/083501'
2020-05-21 14:15:52 +05:30
2020-05-16 20:35:03 +05:30
call selfTest
2019-09-22 19:23:03 +05:30
end subroutine rotations_init
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
! Return rotation in different representations
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
2019-09-20 19:23:49 +05:30
pure function asQuaternion(self)
2019-02-01 14:31:54 +05:30
2022-01-29 20:00:59 +05:30
class(tRotation), intent(in) :: self
real(pReal), dimension(4) :: asQuaternion
2020-04-22 18:56:29 +05:30
2021-01-02 22:27:11 +05:30
asQuaternion = self%q
end function asQuaternion
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
2019-09-21 05:48:09 +05:30
pure function asEulers(self)
2020-04-22 18:56:29 +05:30
2022-01-29 20:00:59 +05:30
class(tRotation), intent(in) :: self
real(pReal), dimension(3) :: asEulers
2020-04-22 18:56:29 +05:30
2021-01-02 22:27:11 +05:30
asEulers = qu2eu(self%q)
2019-09-21 05:48:09 +05:30
end function asEulers
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
2019-09-21 05:48:09 +05:30
pure function asAxisAngle(self)
2019-02-01 14:31:54 +05:30
2022-01-29 20:00:59 +05:30
class(tRotation), intent(in) :: self
real(pReal), dimension(4) :: asAxisAngle
2020-04-22 18:56:29 +05:30
2021-01-02 22:27:11 +05:30
asAxisAngle = qu2ax(self%q)
2019-09-21 05:48:09 +05:30
end function asAxisAngle
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
2019-09-21 05:48:09 +05:30
pure function asMatrix(self)
2020-04-22 18:56:29 +05:30
2022-01-29 20:00:59 +05:30
class(tRotation), intent(in) :: self
real(pReal), dimension(3,3) :: asMatrix
2020-04-22 18:56:29 +05:30
2021-01-02 22:27:11 +05:30
asMatrix = qu2om(self%q)
2019-09-21 05:48:09 +05:30
end function asMatrix
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
2019-09-21 05:48:09 +05:30
pure function asRodrigues(self)
2019-02-01 14:31:54 +05:30
2022-01-29 20:00:59 +05:30
class(tRotation), intent(in) :: self
real(pReal), dimension(4) :: asRodrigues
2020-04-22 18:56:29 +05:30
2021-01-02 22:27:11 +05:30
asRodrigues = qu2ro(self%q)
2020-04-22 18:56:29 +05:30
2019-09-21 05:48:09 +05:30
end function asRodrigues
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
2019-09-20 19:23:49 +05:30
pure function asHomochoric(self)
2019-02-01 14:31:54 +05:30
2022-01-29 20:00:59 +05:30
class(tRotation), intent(in) :: self
real(pReal), dimension(3) :: asHomochoric
2020-04-22 18:56:29 +05:30
2021-01-02 22:27:11 +05:30
asHomochoric = qu2ho(self%q)
end function asHomochoric
2020-04-22 18:56:29 +05:30
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
! Initialize rotation from different representations
!---------------------------------------------------------------------------------------------------
2019-12-02 20:52:27 +05:30
subroutine fromQuaternion(self,qu)
2022-01-29 20:00:59 +05:30
class(tRotation), intent(out) :: self
2019-12-02 20:52:27 +05:30
real(pReal), dimension(4), intent(in) :: qu
if (dNeq(norm2(qu),1.0_pReal,1.0e-8_pReal)) call IO_error(402,ext_msg='fromQuaternion')
2019-12-02 20:52:27 +05:30
self%q = qu
end subroutine fromQuaternion
!---------------------------------------------------------------------------------------------------
2019-09-21 05:48:09 +05:30
subroutine fromEulers(self,eu,degrees)
2019-06-07 18:00:16 +05:30
2022-01-29 20:00:59 +05:30
class(tRotation), intent(out) :: self
2019-06-07 18:00:16 +05:30
real(pReal), dimension(3), intent(in) :: eu
2019-09-20 07:44:37 +05:30
logical, intent(in), optional :: degrees
2019-09-20 19:23:49 +05:30
real(pReal), dimension(3) :: Eulers
2019-09-20 07:44:37 +05:30
if (.not. present(degrees)) then
2019-09-20 19:23:49 +05:30
Eulers = eu
2019-09-20 07:44:37 +05:30
else
2019-09-20 19:23:49 +05:30
Eulers = merge(eu*INRAD,eu,degrees)
2019-09-20 07:44:37 +05:30
endif
if (any(Eulers<0.0_pReal) .or. any(Eulers>TAU) .or. Eulers(2) > PI) &
2019-09-21 05:48:09 +05:30
call IO_error(402,ext_msg='fromEulers')
2019-09-20 19:23:49 +05:30
self%q = eu2qu(Eulers)
2019-09-21 05:48:09 +05:30
end subroutine fromEulers
!---------------------------------------------------------------------------------------------------
2019-09-21 05:48:09 +05:30
subroutine fromAxisAngle(self,ax,degrees,P)
2022-01-29 20:00:59 +05:30
class(tRotation), intent(out) :: self
real(pReal), dimension(4), intent(in) :: ax
logical, intent(in), optional :: degrees
2019-09-20 19:23:49 +05:30
integer, intent(in), optional :: P
2019-09-20 19:27:39 +05:30
real(pReal) :: angle
real(pReal),dimension(3) :: axis
if (.not. present(degrees)) then
2019-09-20 19:23:49 +05:30
angle = ax(4)
else
angle = merge(ax(4)*INRAD,ax(4),degrees)
endif
2020-04-22 18:56:29 +05:30
2019-09-20 19:23:49 +05:30
if (.not. present(P)) then
axis = ax(1:3)
else
2019-09-20 19:27:39 +05:30
axis = ax(1:3) * merge(-1.0_pReal,1.0_pReal,P == 1)
2019-09-21 05:48:09 +05:30
if(abs(P) /= 1) call IO_error(402,ext_msg='fromAxisAngle (P)')
endif
2019-09-20 19:27:39 +05:30
if(dNeq(norm2(axis),1.0_pReal) .or. angle < 0.0_pReal .or. angle > PI) &
2019-09-21 05:48:09 +05:30
call IO_error(402,ext_msg='fromAxisAngle')
2020-04-22 18:56:29 +05:30
2019-09-20 19:23:49 +05:30
self%q = ax2qu([axis,angle])
2019-09-21 05:48:09 +05:30
end subroutine fromAxisAngle
!---------------------------------------------------------------------------------------------------
2019-09-21 05:48:09 +05:30
subroutine fromMatrix(self,om)
2022-01-29 20:00:59 +05:30
class(tRotation), intent(out) :: self
real(pReal), dimension(3,3), intent(in) :: om
2019-09-20 19:23:49 +05:30
2019-09-20 19:27:39 +05:30
if (dNeq(math_det33(om),1.0_pReal,tol=1.0e-5_pReal)) &
2019-09-21 05:48:09 +05:30
call IO_error(402,ext_msg='fromMatrix')
2019-09-20 19:23:49 +05:30
self%q = om2qu(om)
2019-09-21 05:48:09 +05:30
end subroutine fromMatrix
2019-06-07 18:00:16 +05:30
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
!> @brief: Rotate a rotation
!---------------------------------------------------------------------------------------------------
2019-09-21 05:09:33 +05:30
pure elemental function rotRot__(self,R) result(rRot)
2020-04-22 18:56:29 +05:30
2022-01-29 20:00:59 +05:30
type(tRotation) :: rRot
class(tRotation), intent(in) :: self,R
2020-04-22 18:56:29 +05:30
2022-01-29 20:00:59 +05:30
rRot = tRotation(multiply_quaternion(self%q,R%q))
2020-01-14 16:03:18 +05:30
call rRot%standardize()
end function rotRot__
2020-01-14 16:03:18 +05:30
!---------------------------------------------------------------------------------------------------
2020-04-22 18:56:29 +05:30
!> @brief quaternion representation with positive q
2020-01-14 16:03:18 +05:30
!---------------------------------------------------------------------------------------------------
pure elemental subroutine standardize(self)
2022-01-29 20:00:59 +05:30
class(tRotation), intent(inout) :: self
2020-04-22 18:56:29 +05:30
if (sign(1.0_pReal,self%q(1)) < 0.0_pReal) self%q = - self%q
2020-01-14 16:03:18 +05:30
end subroutine standardize
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
2021-01-02 22:27:11 +05:30
!> @brief Rotate a vector passively (default) or actively.
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
2019-09-21 05:09:33 +05:30
pure function rotVector(self,v,active) result(vRot)
2020-04-22 18:56:29 +05:30
real(pReal), dimension(3) :: vRot
2022-01-29 20:00:59 +05:30
class(tRotation), intent(in) :: self
2019-04-03 16:41:18 +05:30
real(pReal), intent(in), dimension(3) :: v
logical, intent(in), optional :: active
2020-04-22 18:56:29 +05:30
2021-01-02 22:27:11 +05:30
real(pReal), dimension(4) :: v_normed, q
logical :: passive
2020-04-22 18:56:29 +05:30
2019-04-03 16:41:18 +05:30
if (present(active)) then
passive = .not. active
else
passive = .true.
endif
2020-04-22 18:56:29 +05:30
2019-09-20 21:15:23 +05:30
if (dEq0(norm2(v))) then
vRot = v
2019-04-03 16:41:18 +05:30
else
2021-01-02 22:27:11 +05:30
v_normed = [0.0_pReal,v]/norm2(v)
2019-04-03 16:41:18 +05:30
if (passive) then
2021-01-02 22:27:11 +05:30
q = multiply_quaternion(self%q, multiply_quaternion(v_normed, conjugate_quaternion(self%q)))
2019-04-03 16:41:18 +05:30
else
2021-01-02 22:27:11 +05:30
q = multiply_quaternion(conjugate_quaternion(self%q), multiply_quaternion(v_normed, self%q))
2019-04-03 16:41:18 +05:30
endif
2021-01-02 22:27:11 +05:30
vRot = q(2:4)*norm2(v)
2019-04-03 16:41:18 +05:30
endif
end function rotVector
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
2021-01-02 22:27:11 +05:30
!> @brief Rotate a rank-2 tensor passively (default) or actively.
!> @details: Rotation is based on rotation matrix
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
2019-09-21 05:09:33 +05:30
pure function rotTensor2(self,T,active) result(tRot)
2020-04-22 18:56:29 +05:30
2019-09-21 05:09:33 +05:30
real(pReal), dimension(3,3) :: tRot
2022-01-29 20:00:59 +05:30
class(tRotation), intent(in) :: self
2019-09-21 05:09:33 +05:30
real(pReal), intent(in), dimension(3,3) :: T
2019-04-03 16:41:18 +05:30
logical, intent(in), optional :: active
2020-04-22 18:56:29 +05:30
2019-04-03 16:41:18 +05:30
logical :: passive
2019-04-03 16:41:18 +05:30
if (present(active)) then
passive = .not. active
else
passive = .true.
endif
2020-04-22 18:56:29 +05:30
tRot = merge(matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())), &
matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix()), &
passive)
2019-09-20 21:15:23 +05:30
end function rotTensor2
2019-09-21 05:09:33 +05:30
!---------------------------------------------------------------------------------------------------
2021-11-20 19:45:59 +05:30
!> @brief Rotate a rank-4 tensor passively (default) or actively.
2019-09-21 05:09:33 +05:30
!> @details: rotation is based on rotation matrix
!! ToDo: Need to check active/passive !!!
!---------------------------------------------------------------------------------------------------
pure function rotTensor4(self,T,active) result(tRot)
real(pReal), dimension(3,3,3,3) :: tRot
2022-01-29 20:00:59 +05:30
class(tRotation), intent(in) :: self
2019-09-21 05:09:33 +05:30
real(pReal), intent(in), dimension(3,3,3,3) :: T
logical, intent(in), optional :: active
2020-04-22 18:56:29 +05:30
2019-09-21 05:09:33 +05:30
real(pReal), dimension(3,3) :: R
integer :: i,j,k,l,m,n,o,p
2021-11-20 19:45:59 +05:30
2019-09-21 05:09:33 +05:30
if (present(active)) then
2019-09-21 05:48:09 +05:30
R = merge(transpose(self%asMatrix()),self%asMatrix(),active)
2019-09-21 05:09:33 +05:30
else
2019-09-21 05:48:09 +05:30
R = self%asMatrix()
2019-09-21 05:09:33 +05:30
endif
tRot = 0.0_pReal
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
tRot(i,j,k,l) = tRot(i,j,k,l) &
+ R(i,m) * R(j,n) * R(k,o) * R(l,p) * T(m,n,o,p)
enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo
end function rotTensor4
!---------------------------------------------------------------------------------------------------
!> @brief Rotate a rank-4 stiffness tensor in Voigt 6x6 notation passively (default) or actively.
2021-11-20 19:45:59 +05:30
!> @details: https://scicomp.stackexchange.com/questions/35600
!! ToDo: Need to check active/passive !!!
!---------------------------------------------------------------------------------------------------
pure function rotStiffness(self,C,active) result(cRot)
real(pReal), dimension(6,6) :: cRot
2022-01-29 20:00:59 +05:30
class(tRotation), intent(in) :: self
2021-11-20 19:45:59 +05:30
real(pReal), intent(in), dimension(6,6) :: C
logical, intent(in), optional :: active
real(pReal), dimension(3,3) :: R
real(pReal), dimension(6,6) :: M
if (present(active)) then
R = merge(transpose(self%asMatrix()),self%asMatrix(),active)
else
R = self%asMatrix()
endif
M = reshape([R(1,1)**2, R(2,1)**2, R(3,1)**2, &
2021-11-20 19:45:59 +05:30
R(2,1)*R(3,1), R(1,1)*R(3,1), R(1,1)*R(2,1), &
R(1,2)**2, R(2,2)**2, R(3,2)**2, &
2021-11-20 19:45:59 +05:30
R(2,2)*R(3,2), R(1,2)*R(3,2), R(1,2)*R(2,2), &
R(1,3)**2, R(2,3)**2, R(3,3)**2, &
2021-11-20 19:45:59 +05:30
R(2,3)*R(3,3), R(1,3)*R(3,3), R(1,3)*R(2,3), &
2.0_pReal*R(1,2)*R(1,3), 2.0_pReal*R(2,2)*R(2,3), 2.0_pReal*R(3,2)*R(3,3), &
R(2,2)*R(3,3)+R(2,3)*R(3,2), R(1,2)*R(3,3)+R(1,3)*R(3,2), R(1,2)*R(2,3)+R(1,3)*R(2,2), &
2.0_pReal*R(1,3)*R(1,1), 2.0_pReal*R(2,3)*R(2,1), 2.0_pReal*R(3,3)*R(3,1), &
R(2,3)*R(3,1)+R(2,1)*R(3,3), R(1,3)*R(3,1)+R(1,1)*R(3,3), R(1,3)*R(2,1)+R(1,1)*R(2,3), &
2.0_pReal*R(1,1)*R(1,2), 2.0_pReal*R(2,1)*R(2,2), 2.0_pReal*R(3,1)*R(3,2), &
R(2,1)*R(3,2)+R(2,2)*R(3,1), R(1,1)*R(3,2)+R(1,2)*R(3,1), R(1,1)*R(2,2)+R(1,2)*R(2,1)],[6,6])
cRot = matmul(M,matmul(C,transpose(M)))
end function rotStiffness
!---------------------------------------------------------------------------------------------------
!> @brief Misorientation.
!---------------------------------------------------------------------------------------------------
2019-09-21 05:09:33 +05:30
pure elemental function misorientation(self,other)
2020-04-22 18:56:29 +05:30
2022-01-29 20:00:59 +05:30
type(tRotation) :: misorientation
class(tRotation), intent(in) :: self, other
2020-04-22 18:56:29 +05:30
2021-11-18 20:55:51 +05:30
2021-01-06 18:07:45 +05:30
misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q))
end function misorientation
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
2021-11-20 19:45:59 +05:30
!> @brief Convert unit quaternion to rotation matrix.
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function qu2om(qu) result(om)
real(pReal), intent(in), dimension(4) :: qu
real(pReal), dimension(3,3) :: om
2020-04-22 18:56:29 +05:30
real(pReal) :: qq
2021-11-20 19:45:59 +05:30
qq = qu(1)**2-sum(qu(2:4)**2)
om(1,1) = qq+2.0_pReal*qu(2)**2
om(2,2) = qq+2.0_pReal*qu(3)**2
om(3,3) = qq+2.0_pReal*qu(4)**2
om(1,2) = 2.0_pReal*(qu(2)*qu(3)-qu(1)*qu(4))
om(2,3) = 2.0_pReal*(qu(3)*qu(4)-qu(1)*qu(2))
om(3,1) = 2.0_pReal*(qu(4)*qu(2)-qu(1)*qu(3))
om(2,1) = 2.0_pReal*(qu(3)*qu(2)+qu(1)*qu(4))
om(3,2) = 2.0_pReal*(qu(4)*qu(3)+qu(1)*qu(2))
om(1,3) = 2.0_pReal*(qu(2)*qu(4)+qu(1)*qu(3))
if (sign(1.0_pReal,P) < 0.0_pReal) om = transpose(om)
end function qu2om
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
2021-11-20 19:45:59 +05:30
!> @brief Convert unit quaternion to Euler angles.
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function qu2eu(qu) result(eu)
real(pReal), intent(in), dimension(4) :: qu
real(pReal), dimension(3) :: eu
2020-04-22 18:56:29 +05:30
2020-04-11 15:06:37 +05:30
real(pReal) :: q12, q03, chi
2020-04-22 18:56:29 +05:30
2021-11-20 19:45:59 +05:30
q03 = qu(1)**2+qu(4)**2
q12 = qu(2)**2+qu(3)**2
chi = sqrt(q03*q12)
2020-04-22 18:56:29 +05:30
2020-04-11 15:06:37 +05:30
degenerated: if (dEq0(q12)) then
eu = [atan2(-P*2.0_pReal*qu(1)*qu(4),qu(1)**2-qu(4)**2), 0.0_pReal, 0.0_pReal]
elseif (dEq0(q03)) then
eu = [atan2( 2.0_pReal*qu(2)*qu(3),qu(2)**2-qu(3)**2), PI, 0.0_pReal]
else degenerated
eu = [atan2((-P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)-qu(3)*qu(4))*chi ), &
atan2( 2.0_pReal*chi, q03-q12 ), &
atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )]
endif degenerated
where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+TAU,[TAU,PI,TAU])
end function qu2eu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert unit quaternion to axis angle pair
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function qu2ax(qu) result(ax)
real(pReal), intent(in), dimension(4) :: qu
real(pReal), dimension(4) :: ax
2020-04-22 18:56:29 +05:30
real(pReal) :: omega, s
if (dEq0(sum(qu(2:4)**2))) then
2019-04-17 16:21:00 +05:30
ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] ! axis = [001]
elseif (dNeq0(qu(1))) then
s = sign(1.0_pReal,qu(1))/norm2(qu(2:4))
omega = 2.0_pReal * acos(math_clip(qu(1),-1.0_pReal,1.0_pReal))
ax = [ qu(2)*s, qu(3)*s, qu(4)*s, omega ]
else
ax = [ qu(2), qu(3), qu(4), PI ]
end if
end function qu2ax
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert unit quaternion to Rodrigues vector
!---------------------------------------------------------------------------------------------------
pure function qu2ro(qu) result(ro)
2020-04-22 18:56:29 +05:30
real(pReal), intent(in), dimension(4) :: qu
real(pReal), dimension(4) :: ro
2020-04-22 18:56:29 +05:30
real(pReal) :: s
real(pReal), parameter :: thr = 1.0e-8_pReal
2020-04-22 18:56:29 +05:30
if (abs(qu(1)) < thr) then
ro = [qu(2), qu(3), qu(4), IEEE_value(1.0_pReal,IEEE_positive_inf)]
2019-04-03 16:41:18 +05:30
else
s = norm2(qu(2:4))
if (s < thr) then
ro = [0.0_pReal, 0.0_pReal, P, 0.0_pReal]
else
ro = [qu(2)/s,qu(3)/s,qu(4)/s, tan(acos(math_clip(qu(1),-1.0_pReal,1.0_pReal)))]
endif
2020-04-22 18:56:29 +05:30
2019-04-03 16:41:18 +05:30
end if
end function qu2ro
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert unit quaternion to homochoric
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function qu2ho(qu) result(ho)
2019-04-03 16:41:18 +05:30
real(pReal), intent(in), dimension(4) :: qu
real(pReal), dimension(3) :: ho
2020-04-22 18:56:29 +05:30
real(pReal) :: omega, f
omega = 2.0 * acos(math_clip(qu(1),-1.0_pReal,1.0_pReal))
2020-04-22 18:56:29 +05:30
2020-04-26 18:28:07 +05:30
if (dEq0(omega,tol=1.e-5_pReal)) then
ho = [ 0.0_pReal, 0.0_pReal, 0.0_pReal ]
2019-04-03 16:41:18 +05:30
else
ho = qu(2:4)
f = 0.75_pReal * ( omega - sin(omega) )
ho = ho/norm2(ho)* f**(1.0_pReal/3.0_pReal)
2019-04-03 16:41:18 +05:30
end if
end function qu2ho
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert unit quaternion to cubochoric
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
2019-09-21 05:09:33 +05:30
pure function qu2cu(qu) result(cu)
2020-04-22 18:56:29 +05:30
real(pReal), intent(in), dimension(4) :: qu
real(pReal), dimension(3) :: cu
cu = ho2cu(qu2ho(qu))
end function qu2cu
!---------------------------------------------------------------------------------------------------
2019-04-17 16:21:00 +05:30
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
2022-05-20 10:00:07 +05:30
!> @brief convert rotation matrix to unit quaternion
2019-04-17 16:21:00 +05:30
!> @details the original formulation (direct conversion) had (numerical?) issues
!---------------------------------------------------------------------------------------------------
pure function om2qu(om) result(qu)
2019-04-17 16:21:00 +05:30
real(pReal), intent(in), dimension(3,3) :: om
real(pReal), dimension(4) :: qu
real(pReal) :: trace,s
trace = math_trace33(om)
if(trace > 0.0_pReal) then
s = 0.5_pReal / sqrt(trace+1.0_pReal)
qu = [0.25_pReal/s, (om(3,2)-om(2,3))*s,(om(1,3)-om(3,1))*s,(om(2,1)-om(1,2))*s]
else
if( om(1,1) > om(2,2) .and. om(1,1) > om(3,3) ) then
s = 2.0_pReal * sqrt( 1.0_pReal + om(1,1) - om(2,2) - om(3,3))
qu = [ (om(3,2) - om(2,3)) /s,0.25_pReal * s,(om(1,2) + om(2,1)) / s,(om(1,3) + om(3,1)) / s]
elseif (om(2,2) > om(3,3)) then
s = 2.0_pReal * sqrt( 1.0_pReal + om(2,2) - om(1,1) - om(3,3))
qu = [ (om(1,3) - om(3,1)) /s,(om(1,2) + om(2,1)) / s,0.25_pReal * s,(om(2,3) + om(3,2)) / s]
else
s = 2.0_pReal * sqrt( 1.0_pReal + om(3,3) - om(1,1) - om(2,2) )
qu = [ (om(2,1) - om(1,2)) /s,(om(1,3) + om(3,1)) / s,(om(2,3) + om(3,2)) / s,0.25_pReal * s]
endif
endif
if(sign(1.0_pReal,qu(1))<0.0_pReal) qu =-1.0_pReal * qu
2022-05-20 10:00:07 +05:30
qu(2:4) = merge(qu(2:4),qu(2:4)*P,dEq0(qu(2:4)))
end function om2qu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
2022-05-20 10:00:07 +05:30
!> @brief convert orientation matrix to Euler angles
!> @details Two step check for special cases to avoid invalid operations (not needed for python)
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function om2eu(om) result(eu)
2019-04-03 16:41:18 +05:30
real(pReal), intent(in), dimension(3,3) :: om
real(pReal), dimension(3) :: eu
real(pReal) :: zeta
2020-04-22 18:56:29 +05:30
if (dNeq(abs(om(3,3)),1.0_pReal,1.e-8_pReal)) then
zeta = 1.0_pReal/sqrt(math_clip(1.0_pReal-om(3,3)**2,1e-64_pReal,1.0_pReal))
2019-04-03 16:41:18 +05:30
eu = [atan2(om(3,1)*zeta,-om(3,2)*zeta), &
acos(math_clip(om(3,3),-1.0_pReal,1.0_pReal)), &
2019-04-03 16:41:18 +05:30
atan2(om(1,3)*zeta, om(2,3)*zeta)]
2020-04-22 18:56:29 +05:30
else
2019-09-23 05:26:23 +05:30
eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ]
2019-04-03 16:41:18 +05:30
end if
where(abs(eu) < 1.e-8_pReal) eu = 0.0_pReal
where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+TAU,[TAU,PI,TAU])
2020-04-22 18:56:29 +05:30
end function om2eu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert orientation matrix to axis angle pair
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
function om2ax(om) result(ax)
2019-04-03 16:41:18 +05:30
real(pReal), intent(in), dimension(3,3) :: om
real(pReal), dimension(4) :: ax
2020-04-22 18:56:29 +05:30
real(pReal) :: t
real(pReal), dimension(3) :: Wr, Wi
real(pReal), dimension((64+2)*3) :: work
real(pReal), dimension(3,3) :: VR, devNull, om_
2019-09-22 19:23:03 +05:30
integer :: ierr, i
2020-04-22 18:56:29 +05:30
om_ = om
2020-04-22 18:56:29 +05:30
2019-04-03 16:41:18 +05:30
! first get the rotation angle
t = 0.5_pReal * (math_trace33(om) - 1.0_pReal)
2019-04-03 16:41:18 +05:30
ax(4) = acos(math_clip(t,-1.0_pReal,1.0_pReal))
2020-04-22 18:56:29 +05:30
2019-04-03 16:41:18 +05:30
if (dEq0(ax(4))) then
ax(1:3) = [ 0.0_pReal, 0.0_pReal, 1.0_pReal ]
2019-04-03 16:41:18 +05:30
else
call dgeev('N','V',3,om_,3,Wr,Wi,devNull,3,VR,3,work,size(work,1),ierr)
if (ierr /= 0) error stop 'LAPACK error'
i = findloc(cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal),.true.,dim=1) !find eigenvalue (1,0)
if (i == 0) error stop 'om2ax conversion failed'
2019-04-03 16:41:18 +05:30
ax(1:3) = VR(1:3,i)
where ( dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) &
ax(1:3) = sign(ax(1:3),-P *[om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])
endif
end function om2ax
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert rotation matrix to Rodrigues vector
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function om2ro(om) result(ro)
2019-04-03 16:41:18 +05:30
real(pReal), intent(in), dimension(3,3) :: om
real(pReal), dimension(4) :: ro
ro = eu2ro(om2eu(om))
end function om2ro
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert rotation matrix to homochoric
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
function om2ho(om) result(ho)
2019-04-03 16:41:18 +05:30
real(pReal), intent(in), dimension(3,3) :: om
real(pReal), dimension(3) :: ho
ho = ax2ho(om2ax(om))
end function om2ho
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert rotation matrix to cubochoric
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
function om2cu(om) result(cu)
real(pReal), intent(in), dimension(3,3) :: om
real(pReal), dimension(3) :: cu
cu = ho2cu(om2ho(om))
end function om2cu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief Euler angles to unit quaternion
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function eu2qu(eu) result(qu)
2019-04-03 16:41:18 +05:30
real(pReal), intent(in), dimension(3) :: eu
real(pReal), dimension(4) :: qu
real(pReal), dimension(3) :: ee
real(pReal) :: cPhi, sPhi
ee = 0.5_pReal*eu
2020-04-22 18:56:29 +05:30
cPhi = cos(ee(2))
sPhi = sin(ee(2))
qu = [ cPhi*cos(ee(1)+ee(3)), &
-P*sPhi*cos(ee(1)-ee(3)), &
-P*sPhi*sin(ee(1)-ee(3)), &
-P*cPhi*sin(ee(1)+ee(3))]
if(sign(1.0_pReal,qu(1)) < 0.0_pReal) qu = qu * (-1.0_pReal)
end function eu2qu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief Euler angles to orientation matrix
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function eu2om(eu) result(om)
2020-04-22 18:56:29 +05:30
real(pReal), intent(in), dimension(3) :: eu
real(pReal), dimension(3,3) :: om
2020-04-22 18:56:29 +05:30
real(pReal), dimension(3) :: c, s
c = cos(eu)
s = sin(eu)
2020-04-22 18:56:29 +05:30
om(1,1) = c(1)*c(3)-s(1)*s(3)*c(2)
om(2,1) = -c(1)*s(3)-s(1)*c(3)*c(2)
om(3,1) = s(1)*s(2)
2020-10-12 08:59:48 +05:30
om(1,2) = s(1)*c(3)+c(1)*s(3)*c(2)
om(2,2) = -s(1)*s(3)+c(1)*c(3)*c(2)
om(3,2) = -c(1)*s(2)
2020-10-12 08:59:48 +05:30
om(1,3) = s(3)*s(2)
om(2,3) = c(3)*s(2)
om(3,3) = c(2)
2020-04-22 18:56:29 +05:30
where(abs(om)<1.0e-12_pReal) om = 0.0_pReal
end function eu2om
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert euler to axis angle
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function eu2ax(eu) result(ax)
2020-04-22 18:56:29 +05:30
real(pReal), intent(in), dimension(3) :: eu
real(pReal), dimension(4) :: ax
2020-04-22 18:56:29 +05:30
real(pReal) :: t, delta, tau, alpha, sigma
2020-04-22 18:56:29 +05:30
t = tan(eu(2)*0.5_pReal)
sigma = 0.5_pReal*(eu(1)+eu(3))
delta = 0.5_pReal*(eu(1)-eu(3))
tau = sqrt(t**2+sin(sigma)**2)
2020-04-22 18:56:29 +05:30
alpha = merge(PI, 2.0_pReal*atan(tau/cos(sigma)), dEq(sigma,PI*0.5_pReal,tol=1.0e-15_pReal))
2020-04-22 18:56:29 +05:30
if (dEq0(alpha)) then ! return a default identity axis-angle pair
ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ]
2019-04-03 16:41:18 +05:30
else
ax(1:3) = -P/tau * [ t*cos(delta), t*sin(delta), sin(sigma) ] ! passive axis-angle pair so a minus sign in front
ax(4) = alpha
if (sign(1.0_pReal,alpha) < 0.0_pReal) ax = -ax ! ensure alpha is positive
2019-04-03 16:41:18 +05:30
end if
2020-04-22 18:56:29 +05:30
end function eu2ax
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief Euler angles to Rodrigues vector
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function eu2ro(eu) result(ro)
real(pReal), intent(in), dimension(3) :: eu
real(pReal), dimension(4) :: ro
2020-04-22 18:56:29 +05:30
ro = eu2ax(eu)
if (ro(4) >= PI) then
ro(4) = IEEE_value(ro(4),IEEE_positive_inf)
elseif(dEq0(ro(4))) then
ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal ]
2019-04-03 16:41:18 +05:30
else
ro(4) = tan(ro(4)*0.5_pReal)
2019-04-03 16:41:18 +05:30
end if
2020-04-22 18:56:29 +05:30
end function eu2ro
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Euler angles to homochoric
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function eu2ho(eu) result(ho)
2019-04-03 16:41:18 +05:30
real(pReal), intent(in), dimension(3) :: eu
real(pReal), dimension(3) :: ho
ho = ax2ho(eu2ax(eu))
end function eu2ho
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Euler angles to cubochoric
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
function eu2cu(eu) result(cu)
real(pReal), intent(in), dimension(3) :: eu
2019-04-03 16:41:18 +05:30
real(pReal), dimension(3) :: cu
cu = ho2cu(eu2ho(eu))
end function eu2cu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert axis angle pair to quaternion
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function ax2qu(ax) result(qu)
2020-04-22 18:56:29 +05:30
real(pReal), intent(in), dimension(4) :: ax
real(pReal), dimension(4) :: qu
real(pReal) :: c, s
if (dEq0(ax(4))) then
qu = [ 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal ]
else
c = cos(ax(4)*0.5_pReal)
s = sin(ax(4)*0.5_pReal)
qu = [ c, ax(1)*s, ax(2)*s, ax(3)*s ]
end if
end function ax2qu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert axis angle pair to orientation matrix
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function ax2om(ax) result(om)
real(pReal), intent(in), dimension(4) :: ax
real(pReal), dimension(3,3) :: om
2020-04-22 18:56:29 +05:30
real(pReal) :: q, c, s, omc
c = cos(ax(4))
s = sin(ax(4))
omc = 1.0_pReal-c
2019-09-20 19:23:49 +05:30
om(1,1) = ax(1)**2*omc + c
om(2,2) = ax(2)**2*omc + c
om(3,3) = ax(3)**2*omc + c
q = omc*ax(1)*ax(2)
om(1,2) = q + s*ax(3)
om(2,1) = q - s*ax(3)
2020-04-22 18:56:29 +05:30
q = omc*ax(2)*ax(3)
om(2,3) = q + s*ax(1)
om(3,2) = q - s*ax(1)
2020-04-22 18:56:29 +05:30
q = omc*ax(3)*ax(1)
om(3,1) = q + s*ax(2)
om(1,3) = q - s*ax(2)
if (P > 0.0_pReal) om = transpose(om)
end function ax2om
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert axis angle pair to Euler angles
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function ax2eu(ax) result(eu)
real(pReal), intent(in), dimension(4) :: ax
real(pReal), dimension(3) :: eu
eu = om2eu(ax2om(ax))
end function ax2eu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert axis angle pair to Rodrigues vector
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function ax2ro(ax) result(ro)
real(pReal), intent(in), dimension(4) :: ax
real(pReal), dimension(4) :: ro
2020-04-22 18:56:29 +05:30
real(pReal), parameter :: thr = 1.0e-7_pReal
2020-04-22 18:56:29 +05:30
if (dEq0(ax(4))) then
ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal ]
2020-04-22 18:56:29 +05:30
else
ro(1:3) = ax(1:3)
! we need to deal with the 180 degree case
ro(4) = merge(IEEE_value(ro(4),IEEE_positive_inf),tan(ax(4)*0.5_pReal),abs(ax(4)-PI) < thr)
end if
end function ax2ro
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert axis angle pair to homochoric
!---------------------------------------------------------------------------------------------------
pure function ax2ho(ax) result(ho)
real(pReal), intent(in), dimension(4) :: ax
real(pReal), dimension(3) :: ho
2020-04-22 18:56:29 +05:30
real(pReal) :: f
2020-04-22 18:56:29 +05:30
f = 0.75_pReal * ( ax(4) - sin(ax(4)) )
f = f**(1.0_pReal/3.0_pReal)
ho = ax(1:3) * f
end function ax2ho
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert axis angle pair to cubochoric
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
function ax2cu(ax) result(cu)
real(pReal), intent(in), dimension(4) :: ax
real(pReal), dimension(3) :: cu
2020-04-29 18:22:18 +05:30
cu = ho2cu(ax2ho(ax))
end function ax2cu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Rodrigues vector to unit quaternion
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function ro2qu(ro) result(qu)
real(pReal), intent(in), dimension(4) :: ro
real(pReal), dimension(4) :: qu
2020-04-22 18:56:29 +05:30
qu = ax2qu(ro2ax(ro))
end function ro2qu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
2019-02-01 14:31:54 +05:30
!> @brief convert Rodrigues vector to rotation matrix
!---------------------------------------------------------------------------------------------------
pure function ro2om(ro) result(om)
2019-04-03 16:41:18 +05:30
real(pReal), intent(in), dimension(4) :: ro
real(pReal), dimension(3,3) :: om
2019-04-03 16:41:18 +05:30
om = ax2om(ro2ax(ro))
end function ro2om
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Rodrigues vector to Euler angles
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function ro2eu(ro) result(eu)
real(pReal), intent(in), dimension(4) :: ro
real(pReal), dimension(3) :: eu
2020-04-22 18:56:29 +05:30
eu = om2eu(ro2om(ro))
end function ro2eu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Rodrigues vector to axis angle pair
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function ro2ax(ro) result(ax)
real(pReal), intent(in), dimension(4) :: ro
real(pReal), dimension(4) :: ax
2020-04-22 18:56:29 +05:30
real(pReal) :: ta, angle
2020-04-22 18:56:29 +05:30
ta = ro(4)
2020-04-22 18:56:29 +05:30
if (.not. IEEE_is_finite(ta)) then
ax = [ ro(1), ro(2), ro(3), PI ]
2020-04-22 18:56:29 +05:30
elseif (dEq0(ta)) then
ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ]
else
angle = 2.0_pReal*atan(ta)
ta = 1.0_pReal/norm2(ro(1:3))
ax = [ ro(1)/ta, ro(2)/ta, ro(3)/ta, angle ]
end if
end function ro2ax
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Rodrigues vector to homochoric
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function ro2ho(ro) result(ho)
real(pReal), intent(in), dimension(4) :: ro
real(pReal), dimension(3) :: ho
2020-04-22 18:56:29 +05:30
real(pReal) :: f
2020-04-22 18:56:29 +05:30
if (dEq0(norm2(ro(1:3)))) then
ho = [ 0.0_pReal, 0.0_pReal, 0.0_pReal ]
else
f = merge(2.0_pReal*atan(ro(4)) - sin(2.0_pReal*atan(ro(4))),PI, IEEE_is_finite(ro(4)))
ho = ro(1:3) * (0.75_pReal*f)**(1.0_pReal/3.0_pReal)
end if
end function ro2ho
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Rodrigues vector to cubochoric
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
2019-09-21 05:09:33 +05:30
pure function ro2cu(ro) result(cu)
real(pReal), intent(in), dimension(4) :: ro
real(pReal), dimension(3) :: cu
cu = ho2cu(ro2ho(ro))
end function ro2cu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
2019-02-01 14:31:54 +05:30
!> @brief convert homochoric to unit quaternion
!---------------------------------------------------------------------------------------------------
pure function ho2qu(ho) result(qu)
real(pReal), intent(in), dimension(3) :: ho
real(pReal), dimension(4) :: qu
2019-04-03 16:41:18 +05:30
qu = ax2qu(ho2ax(ho))
end function ho2qu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert homochoric to rotation matrix
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function ho2om(ho) result(om)
real(pReal), intent(in), dimension(3) :: ho
real(pReal), dimension(3,3) :: om
om = ax2om(ho2ax(ho))
end function ho2om
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert homochoric to Euler angles
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function ho2eu(ho) result(eu)
real(pReal), intent(in), dimension(3) :: ho
real(pReal), dimension(3) :: eu
eu = ax2eu(ho2ax(ho))
end function ho2eu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert homochoric to axis angle pair
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function ho2ax(ho) result(ax)
real(pReal), intent(in), dimension(3) :: ho
real(pReal), dimension(4) :: ax
2020-04-22 18:56:29 +05:30
integer :: i
real(pReal) :: hmag_squared, s, hm
real(pReal), parameter, dimension(16) :: &
2020-04-22 18:56:29 +05:30
tfit = [ 1.0000000000018852_pReal, -0.5000000002194847_pReal, &
-0.024999992127593126_pReal, -0.003928701544781374_pReal, &
-0.0008152701535450438_pReal, -0.0002009500426119712_pReal, &
-0.00002397986776071756_pReal, -0.00008202868926605841_pReal, &
+0.00012448715042090092_pReal, -0.0001749114214822577_pReal, &
+0.0001703481934140054_pReal, -0.00012062065004116828_pReal, &
+0.000059719705868660826_pReal, -0.00001980756723965647_pReal, &
+0.000003953714684212874_pReal, -0.00000036555001439719544_pReal ]
2020-04-22 18:56:29 +05:30
! normalize h and store the magnitude
2021-11-26 01:22:22 +05:30
hmag_squared = sum(ho**2)
if (dEq0(hmag_squared)) then
ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ]
else
hm = hmag_squared
2020-04-22 18:56:29 +05:30
! convert the magnitude to the rotation angle
s = tfit(1) + tfit(2) * hmag_squared
do i=3,16
hm = hm*hmag_squared
s = s + tfit(i) * hm
2021-07-16 14:00:45 +05:30
enddo
ax = [ho/sqrt(hmag_squared), 2.0_pReal*acos(s)]
end if
end function ho2ax
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert homochoric to Rodrigues vector
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
pure function ho2ro(ho) result(ro)
real(pReal), intent(in), dimension(3) :: ho
real(pReal), dimension(4) :: ro
ro = ax2ro(ho2ax(ho))
end function ho2ro
2020-04-29 18:23:27 +05:30
!--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief convert homochoric to cubochoric
!--------------------------------------------------------------------------
pure function ho2cu(ho) result(cu)
real(pReal), intent(in), dimension(3) :: ho
real(pReal), dimension(3) :: cu, xyz1, xyz3
real(pReal), dimension(2) :: Tinv, xyz2
real(pReal) :: rs, qxy, q2, sq2, q, tt
integer, dimension(3,2) :: p
rs = norm2(ho)
if (rs > R1+1.e-6_pReal) then
cu = IEEE_value(cu,IEEE_positive_inf)
return
endif
center: if (all(dEq0(ho))) then
cu = 0.0_pReal
else center
p = GetPyramidOrder(ho)
xyz3 = ho(p(:,1))
! inverse M_3
xyz2 = xyz3(1:2) * sqrt( 2.0*rs/(rs+abs(xyz3(3))) )
! inverse M_2
qxy = sum(xyz2**2)
special: if (dEq0(qxy)) then
Tinv = 0.0_pReal
else special
q2 = qxy + maxval(abs(xyz2))**2
sq2 = sqrt(q2)
2022-01-29 19:19:45 +05:30
q = (BETA/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2))
2020-04-29 18:23:27 +05:30
tt = (minval(abs(xyz2))**2+maxval(abs(xyz2))*sq2)/R2/qxy
Tinv = q * sign(1.0_pReal,xyz2) * merge([ 1.0_pReal, acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12], &
[ acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12, 1.0_pReal], &
abs(xyz2(2)) <= abs(xyz2(1)))
endif special
! inverse M_1
2022-01-29 19:19:45 +05:30
xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / PREF ]/SC
2020-04-29 18:23:27 +05:30
! reverse the coordinates back to order according to the original pyramid number
cu = xyz1(p(:,2))
endif center
end function ho2cu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert cubochoric to unit quaternion
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
2019-09-21 05:09:33 +05:30
pure function cu2qu(cu) result(qu)
2019-04-03 16:41:18 +05:30
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(4) :: qu
qu = ho2qu(cu2ho(cu))
end function cu2qu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
2019-02-01 14:31:54 +05:30
!> @brief convert cubochoric to rotation matrix
!---------------------------------------------------------------------------------------------------
2019-09-21 05:09:33 +05:30
pure function cu2om(cu) result(om)
2019-04-03 16:41:18 +05:30
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(3,3) :: om
2019-04-03 16:41:18 +05:30
om = ho2om(cu2ho(cu))
end function cu2om
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert cubochoric to Euler angles
!---------------------------------------------------------------------------------------------------
2019-09-21 05:09:33 +05:30
pure function cu2eu(cu) result(eu)
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(3) :: eu
eu = ho2eu(cu2ho(cu))
end function cu2eu
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
2019-02-01 14:31:54 +05:30
!> @brief convert cubochoric to axis angle pair
!---------------------------------------------------------------------------------------------------
function cu2ax(cu) result(ax)
2019-04-03 16:41:18 +05:30
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(4) :: ax
2019-04-03 16:41:18 +05:30
ax = ho2ax(cu2ho(cu))
end function cu2ax
2019-02-01 14:31:54 +05:30
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
2019-02-01 14:31:54 +05:30
!> @brief convert cubochoric to Rodrigues vector
!---------------------------------------------------------------------------------------------------
2019-09-21 05:09:33 +05:30
pure function cu2ro(cu) result(ro)
2019-04-03 16:41:18 +05:30
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(4) :: ro
2019-04-03 16:41:18 +05:30
ro = ho2ro(cu2ho(cu))
end function cu2ro
2020-04-22 18:56:29 +05:30
!--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief map from 3D cubic grid to 3D ball
!--------------------------------------------------------------------------
2020-04-29 18:22:18 +05:30
pure function cu2ho(cu) result(ho)
2020-04-22 18:56:29 +05:30
2020-04-29 18:12:21 +05:30
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(3) :: ho, LamXYZ, XYZ
2020-04-22 18:56:29 +05:30
real(pReal), dimension(2) :: T
real(pReal) :: c, s, q
real(pReal), parameter :: eps = 1.0e-8_pReal
integer, dimension(3,2) :: p
integer, dimension(2) :: order
2020-04-29 18:12:21 +05:30
if (maxval(abs(cu)) > AP/2.0+eps) then
ho = IEEE_value(cu,IEEE_positive_inf)
2020-04-22 18:56:29 +05:30
return
end if
! transform to the sphere grid via the curved square, and intercept the zero point
2020-04-29 18:12:21 +05:30
center: if (all(dEq0(cu))) then
ho = 0.0_pReal
2020-04-22 18:56:29 +05:30
else center
! get pyramide and scale by grid parameter ratio
2020-04-29 18:12:21 +05:30
p = GetPyramidOrder(cu)
2022-01-29 19:19:45 +05:30
XYZ = cu(p(:,1)) * SC
2020-04-22 18:56:29 +05:30
! intercept all the points along the z-axis
special: if (all(dEq0(XYZ(1:2)))) then
2022-01-29 19:19:45 +05:30
LamXYZ = [ 0.0_pReal, 0.0_pReal, PREF * XYZ(3) ]
2020-04-22 18:56:29 +05:30
else special
order = merge( [2,1], [1,2], abs(XYZ(2)) <= abs(XYZ(1))) ! order of absolute values of XYZ
q = PI12 * XYZ(order(1))/XYZ(order(2)) ! smaller by larger
c = cos(q)
s = sin(q)
2022-01-29 19:19:45 +05:30
q = PREK * XYZ(order(2))/ sqrt(R2-c)
2020-04-22 18:56:29 +05:30
T = [ (R2*c - 1.0), R2 * s] * q
! transform to sphere grid (inverse Lambert)
! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero]
c = sum(T**2)
2022-05-20 10:00:07 +05:30
s = c * PI/(24.0*XYZ(3)**2)
c = c * sqrt(PI/24.0_pReal) / XYZ(3)
2020-04-22 18:56:29 +05:30
q = sqrt( 1.0 - s )
2022-01-29 19:19:45 +05:30
LamXYZ = [ T(order(2)) * q, T(order(1)) * q, PREF * XYZ(3) - c ]
end if special
2020-04-22 18:56:29 +05:30
! reverse the coordinates back to order according to the original pyramid number
2020-04-29 18:12:21 +05:30
ho = LamXYZ(p(:,2))
2020-04-22 18:56:29 +05:30
end if center
2020-04-22 18:56:29 +05:30
2020-04-29 18:22:18 +05:30
end function cu2ho
2020-04-22 18:56:29 +05:30
!--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
2021-01-02 22:27:11 +05:30
!> @brief Determine to which pyramid a point in a cubic grid belongs.
2020-04-22 18:56:29 +05:30
!--------------------------------------------------------------------------
pure function GetPyramidOrder(xyz)
real(pReal),intent(in),dimension(3) :: xyz
integer, dimension(3,2) :: GetPyramidOrder
if (((abs(xyz(1)) <= xyz(3)).and.(abs(xyz(2)) <= xyz(3))) .or. &
((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then
GetPyramidOrder = reshape([[1,2,3],[1,2,3]],[3,2])
else if (((abs(xyz(3)) <= xyz(1)).and.(abs(xyz(2)) <= xyz(1))) .or. &
((abs(xyz(3)) <= -xyz(1)).and.(abs(xyz(2)) <= -xyz(1)))) then
GetPyramidOrder = reshape([[2,3,1],[3,1,2]],[3,2])
else if (((abs(xyz(1)) <= xyz(2)).and.(abs(xyz(3)) <= xyz(2))) .or. &
((abs(xyz(1)) <= -xyz(2)).and.(abs(xyz(3)) <= -xyz(2)))) then
GetPyramidOrder = reshape([[3,1,2],[2,3,1]],[3,2])
else
GetPyramidOrder = -1 ! should be impossible, but might simplify debugging
end if
end function GetPyramidOrder
2019-09-22 19:23:03 +05:30
!--------------------------------------------------------------------------------------------------
2021-01-02 22:27:11 +05:30
!> @brief Multiply two quaternions.
!--------------------------------------------------------------------------------------------------
pure function multiply_quaternion(qu1,qu2)
real(pReal), dimension(4), intent(in) :: qu1, qu2
real(pReal), dimension(4) :: multiply_quaternion
multiply_quaternion(1) = qu1(1)*qu2(1) - qu1(2)*qu2(2) - qu1(3)*qu2(3) - qu1(4)*qu2(4)
multiply_quaternion(2) = qu1(1)*qu2(2) + qu1(2)*qu2(1) + P * (qu1(3)*qu2(4) - qu1(4)*qu2(3))
multiply_quaternion(3) = qu1(1)*qu2(3) + qu1(3)*qu2(1) + P * (qu1(4)*qu2(2) - qu1(2)*qu2(4))
multiply_quaternion(4) = qu1(1)*qu2(4) + qu1(4)*qu2(1) + P * (qu1(2)*qu2(3) - qu1(3)*qu2(2))
end function multiply_quaternion
!--------------------------------------------------------------------------------------------------
!> @brief Calculate conjugate complex of a quaternion.
!--------------------------------------------------------------------------------------------------
pure function conjugate_quaternion(qu)
real(pReal), dimension(4), intent(in) :: qu
real(pReal), dimension(4) :: conjugate_quaternion
conjugate_quaternion = [qu(1), -qu(2), -qu(3), -qu(4)]
end function conjugate_quaternion
!--------------------------------------------------------------------------------------------------
!> @brief Check correctness of some rotations functions.
2019-09-22 19:23:03 +05:30
!--------------------------------------------------------------------------------------------------
subroutine selfTest()
2020-04-22 18:56:29 +05:30
2022-01-29 20:00:59 +05:30
type(tRotation) :: R
real(pReal), dimension(4) :: qu, ax, ro
real(pReal), dimension(3) :: x, eu, ho, v3
real(pReal), dimension(3,3) :: om, t33
real(pReal), dimension(3,3,3,3) :: t3333
2021-11-20 19:45:59 +05:30
real(pReal), dimension(6,6) :: C
2019-09-22 19:23:03 +05:30
real :: A,B
integer :: i
2021-01-02 22:27:11 +05:30
do i = 1, 20
2020-04-22 18:56:29 +05:30
2019-09-22 23:59:34 +05:30
if(i==1) then
qu = om2qu(math_I3)
elseif(i==2) then
qu = eu2qu([0.0_pReal,0.0_pReal,0.0_pReal])
elseif(i==3) then
qu = eu2qu([TAU,PI,TAU])
2019-09-22 23:59:34 +05:30
elseif(i==4) then
qu = [0.0_pReal,0.0_pReal,1.0_pReal,0.0_pReal]
elseif(i==5) then
qu = ro2qu([1.0_pReal,0.0_pReal,0.0_pReal,IEEE_value(1.0_pReal, IEEE_positive_inf)])
2019-09-22 23:59:34 +05:30
elseif(i==6) then
qu = ax2qu([1.0_pReal,0.0_pReal,0.0_pReal,0.0_pReal])
else
call random_number(x)
A = sqrt(x(3))
B = sqrt(1-0_pReal -x(3))
qu = [cos(TAU*x(1))*A,&
sin(TAU*x(2))*B,&
cos(TAU*x(2))*B,&
sin(TAU*x(1))*A]
2019-09-22 23:59:34 +05:30
if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
endif
2021-11-20 19:45:59 +05:30
if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) error stop 'om2qu/qu2om'
if(.not. quaternion_equal(eu2qu(qu2eu(qu)),qu)) error stop 'eu2qu/qu2eu'
if(.not. quaternion_equal(ax2qu(qu2ax(qu)),qu)) error stop 'ax2qu/qu2ax'
if(.not. quaternion_equal(ro2qu(qu2ro(qu)),qu)) error stop 'ro2qu/qu2ro'
if(.not. quaternion_equal(ho2qu(qu2ho(qu)),qu)) error stop 'ho2qu/qu2ho'
if(.not. quaternion_equal(cu2qu(qu2cu(qu)),qu)) error stop 'cu2qu/qu2cu'
om = qu2om(qu)
if(.not. quaternion_equal(om2qu(eu2om(om2eu(om))),qu)) error stop 'eu2om/om2eu'
if(.not. quaternion_equal(om2qu(ax2om(om2ax(om))),qu)) error stop 'ax2om/om2ax'
if(.not. quaternion_equal(om2qu(ro2om(om2ro(om))),qu)) error stop 'ro2om/om2ro'
if(.not. quaternion_equal(om2qu(ho2om(om2ho(om))),qu)) error stop 'ho2om/om2ho'
if(.not. quaternion_equal(om2qu(cu2om(om2cu(om))),qu)) error stop 'cu2om/om2cu'
eu = qu2eu(qu)
if(.not. quaternion_equal(eu2qu(ax2eu(eu2ax(eu))),qu)) error stop 'ax2eu/eu2ax'
if(.not. quaternion_equal(eu2qu(ro2eu(eu2ro(eu))),qu)) error stop 'ro2eu/eu2ro'
if(.not. quaternion_equal(eu2qu(ho2eu(eu2ho(eu))),qu)) error stop 'ho2eu/eu2ho'
if(.not. quaternion_equal(eu2qu(cu2eu(eu2cu(eu))),qu)) error stop 'cu2eu/eu2cu'
ax = qu2ax(qu)
if(.not. quaternion_equal(ax2qu(ro2ax(ax2ro(ax))),qu)) error stop 'ro2ax/ax2ro'
if(.not. quaternion_equal(ax2qu(ho2ax(ax2ho(ax))),qu)) error stop 'ho2ax/ax2ho'
if(.not. quaternion_equal(ax2qu(cu2ax(ax2cu(ax))),qu)) error stop 'cu2ax/ax2cu'
ro = qu2ro(qu)
if(.not. quaternion_equal(ro2qu(ho2ro(ro2ho(ro))),qu)) error stop 'ho2ro/ro2ho'
if(.not. quaternion_equal(ro2qu(cu2ro(ro2cu(ro))),qu)) error stop 'cu2ro/ro2cu'
ho = qu2ho(qu)
if(.not. quaternion_equal(ho2qu(cu2ho(ho2cu(ho))),qu)) error stop 'cu2ho/ho2cu'
call R%fromMatrix(om)
2020-04-22 18:56:29 +05:30
call random_number(v3)
if (any(dNeq(R%rotVector(R%rotVector(v3),active=.true.),v3,1.0e-12_pReal))) &
error stop 'rotVector'
2020-04-22 18:56:29 +05:30
call random_number(t33)
if (any(dNeq(R%rotTensor2(R%rotTensor2(t33),active=.true.),t33,1.0e-12_pReal))) &
error stop 'rotTensor2'
2020-04-22 18:56:29 +05:30
call random_number(t3333)
if (any(dNeq(R%rotTensor4(R%rotTensor4(t3333),active=.true.),t3333,1.0e-12_pReal))) &
error stop 'rotTensor4'
2019-09-22 19:23:03 +05:30
2021-11-20 19:45:59 +05:30
call random_number(C)
C = C+transpose(C)
if (any(dNeq(R%rotStiffness(C), &
math_3333toVoigt66_stiffness(R%rotate(math_Voigt66to3333_stiffness(C))),1.0e-12_pReal))) &
2021-11-20 19:45:59 +05:30
error stop 'rotStiffness'
call R%fromQuaternion(qu * (1.0_pReal + merge(+5.e-9_pReal,-5.e-9_pReal, mod(i,2) == 0))) ! allow reasonable tolerance for ASCII/YAML
2021-11-20 19:45:59 +05:30
end do
contains
pure recursive function quaternion_equal(qu1,qu2) result(ok)
real(pReal), intent(in), dimension(4) :: qu1,qu2
logical :: ok
ok = all(dEq(qu1,qu2,1.0e-7_pReal))
if(dEq0(qu1(1),1.0e-12_pReal)) &
ok = ok .or. all(dEq(-1.0_pReal*qu1,qu2,1.0e-7_pReal))
end function quaternion_equal
2020-05-16 20:35:03 +05:30
end subroutine selfTest
2019-09-22 19:23:03 +05:30
end module rotations