2018-12-08 12:32:55 +05:30
|
|
|
! ###################################################################
|
|
|
|
|
! 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
|
2018-12-08 12:32:55 +05:30
|
|
|
! All rights reserved.
|
|
|
|
|
!
|
|
|
|
|
! Redistribution and use in source and binary forms, with or without modification, are
|
|
|
|
|
! permitted provided that the following conditions are met:
|
|
|
|
|
!
|
|
|
|
|
! - Redistributions of source code must retain the above copyright notice, this list
|
|
|
|
|
! of conditions and the following disclaimer.
|
|
|
|
|
! - 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.
|
|
|
|
|
! - 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.
|
|
|
|
|
!
|
|
|
|
|
! 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
|
2019-09-20 21:06:16 +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_etal2015
|
|
|
|
|
! 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
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
|
|
2018-12-08 12:32:55 +05:30
|
|
|
module rotations
|
2019-05-17 02:26:48 +05:30
|
|
|
use prec
|
|
|
|
|
use IO
|
|
|
|
|
use math
|
|
|
|
|
use Lambert
|
2019-04-03 16:41:18 +05:30
|
|
|
use quaternions
|
|
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
private
|
2019-05-17 02:26:48 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
type, public :: rotation
|
|
|
|
|
type(quaternion), private :: q
|
|
|
|
|
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
|
|
|
!------------------------------------------
|
2019-09-20 11:20:30 +05:30
|
|
|
procedure, private :: rotRot__
|
|
|
|
|
generic, public :: operator(*) => rotRot__
|
2019-09-21 05:09:33 +05:30
|
|
|
generic, public :: rotate => rotVector,rotTensor2,rotTensor4
|
2019-09-20 11:20:30 +05:30
|
|
|
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
|
2019-09-21 06:06:37 +05:30
|
|
|
procedure, public :: rotTensor4sym
|
2019-09-20 11:20:30 +05:30
|
|
|
procedure, public :: misorientation
|
2020-01-14 16:03:18 +05:30
|
|
|
procedure, public :: standardize
|
2019-04-03 16:41:18 +05:30
|
|
|
end type rotation
|
2019-09-22 12:15:54 +05:30
|
|
|
|
2019-09-22 19:23:03 +05:30
|
|
|
public :: &
|
2020-01-29 13:55:39 +05:30
|
|
|
rotations_init, &
|
|
|
|
|
eu2om
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
contains
|
|
|
|
|
|
2019-09-22 19:23:03 +05:30
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
|
|
|
!> @brief doing self test
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
2019-09-23 00:40:39 +05:30
|
|
|
subroutine rotations_init
|
2019-09-22 19:23:03 +05:30
|
|
|
|
2019-09-23 18:07:36 +05:30
|
|
|
call quaternions_init
|
2020-01-14 16:03:18 +05:30
|
|
|
write(6,'(/,a)') ' <<<+- rotations init -+>>>'; flush(6)
|
2019-09-23 00:40:39 +05:30
|
|
|
call unitTest
|
2019-09-22 19:23:03 +05:30
|
|
|
|
|
|
|
|
end subroutine rotations_init
|
|
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-09-23 10:28:18 +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
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
class(rotation), intent(in) :: self
|
|
|
|
|
real(pReal), dimension(4) :: asQuaternion
|
|
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
asQuaternion = self%q%asArray()
|
2018-12-08 17:33:27 +05:30
|
|
|
|
|
|
|
|
end function asQuaternion
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-09-21 05:48:09 +05:30
|
|
|
pure function asEulers(self)
|
2019-02-01 14:31:54 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
class(rotation), intent(in) :: self
|
2019-09-21 05:48:09 +05:30
|
|
|
real(pReal), dimension(3) :: asEulers
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-21 05:48:09 +05:30
|
|
|
asEulers = qu2eu(self%q%asArray())
|
2018-12-08 12:32:55 +05:30
|
|
|
|
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
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
class(rotation), intent(in) :: self
|
2019-09-21 05:48:09 +05:30
|
|
|
real(pReal), dimension(4) :: asAxisAngle
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-21 05:48:09 +05:30
|
|
|
asAxisAngle = qu2ax(self%q%asArray())
|
2018-12-08 12:32:55 +05:30
|
|
|
|
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)
|
2019-02-01 14:31:54 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
class(rotation), intent(in) :: self
|
2019-09-21 05:48:09 +05:30
|
|
|
real(pReal), dimension(3,3) :: asMatrix
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-21 05:48:09 +05:30
|
|
|
asMatrix = qu2om(self%q%asArray())
|
2018-12-08 12:32:55 +05:30
|
|
|
|
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
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
class(rotation), intent(in) :: self
|
2019-09-21 05:48:09 +05:30
|
|
|
real(pReal), dimension(4) :: asRodrigues
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-21 05:48:09 +05:30
|
|
|
asRodrigues = qu2ro(self%q%asArray())
|
2019-02-01 14:31:54 +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
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
class(rotation), intent(in) :: self
|
|
|
|
|
real(pReal), dimension(3) :: asHomochoric
|
|
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
asHomochoric = qu2ho(self%q%asArray())
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
end function asHomochoric
|
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)
|
|
|
|
|
|
|
|
|
|
class(rotation), intent(out) :: self
|
|
|
|
|
real(pReal), dimension(4), intent(in) :: qu
|
|
|
|
|
|
2019-12-03 00:36:58 +05:30
|
|
|
if (dNeq(norm2(qu),1.0_pReal)) &
|
2019-12-02 20:52:27 +05:30
|
|
|
call IO_error(402,ext_msg='fromQuaternion')
|
|
|
|
|
|
|
|
|
|
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
|
|
|
|
|
|
|
|
class(rotation), intent(out) :: self
|
|
|
|
|
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
|
|
|
|
|
|
2019-09-20 19:27:39 +05:30
|
|
|
if (any(Eulers<0.0_pReal) .or. any(Eulers>2.0_pReal*PI) .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-20 11:20:30 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-09-21 05:48:09 +05:30
|
|
|
subroutine fromAxisAngle(self,ax,degrees,P)
|
2019-09-20 11:20:30 +05:30
|
|
|
|
|
|
|
|
class(rotation), 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
|
2019-09-20 11:20:30 +05:30
|
|
|
|
|
|
|
|
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
|
|
|
|
|
|
|
|
|
|
if (.not. present(P)) then
|
|
|
|
|
axis = ax(1:3)
|
2019-09-20 11:20:30 +05:30
|
|
|
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)')
|
2019-09-20 11:20:30 +05:30
|
|
|
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')
|
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-20 11:20:30 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-09-21 05:48:09 +05:30
|
|
|
subroutine fromMatrix(self,om)
|
2019-09-20 11:20:30 +05:30
|
|
|
|
|
|
|
|
class(rotation), 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
|
|
|
|
2019-09-20 11:20:30 +05:30
|
|
|
self%q = om2qu(om)
|
2018-12-08 20:14:00 +05:30
|
|
|
|
2019-09-21 05:48:09 +05:30
|
|
|
end subroutine fromMatrix
|
2019-06-07 18:00:16 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 20:14:00 +05:30
|
|
|
|
2019-09-20 11:20:30 +05:30
|
|
|
|
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
|
!> @brief: Rotate a rotation
|
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-09-21 05:09:33 +05:30
|
|
|
pure elemental function rotRot__(self,R) result(rRot)
|
2019-09-20 11:20:30 +05:30
|
|
|
|
|
|
|
|
type(rotation) :: rRot
|
2019-09-21 05:09:33 +05:30
|
|
|
class(rotation), intent(in) :: self,R
|
2019-09-20 11:20:30 +05:30
|
|
|
|
2019-09-21 05:09:33 +05:30
|
|
|
rRot = rotation(self%q*R%q)
|
2020-01-14 16:03:18 +05:30
|
|
|
call rRot%standardize()
|
2019-09-20 11:20:30 +05:30
|
|
|
|
|
|
|
|
end function rotRot__
|
|
|
|
|
|
|
|
|
|
|
2020-01-14 16:03:18 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
|
!> @brief quaternion representation with positive q
|
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
|
pure elemental subroutine standardize(self)
|
|
|
|
|
|
|
|
|
|
class(rotation), intent(inout) :: self
|
|
|
|
|
|
|
|
|
|
if (real(self%q) < 0.0_pReal) self%q = self%q%homomorphed()
|
|
|
|
|
|
|
|
|
|
end subroutine standardize
|
|
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-02-01 14:31:54 +05:30
|
|
|
!> @brief rotate a vector passively (default) or actively
|
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-09-21 05:09:33 +05:30
|
|
|
pure function rotVector(self,v,active) result(vRot)
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-20 11:20:30 +05:30
|
|
|
real(pReal), dimension(3) :: vRot
|
2019-04-03 16:41:18 +05:30
|
|
|
class(rotation), intent(in) :: self
|
|
|
|
|
real(pReal), intent(in), dimension(3) :: v
|
|
|
|
|
logical, intent(in), optional :: active
|
|
|
|
|
|
2019-09-20 21:15:23 +05:30
|
|
|
real(pReal), dimension(3) :: v_normed
|
|
|
|
|
type(quaternion) :: q
|
|
|
|
|
logical :: passive
|
2019-02-01 14:31:54 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
if (present(active)) then
|
|
|
|
|
passive = .not. active
|
|
|
|
|
else
|
|
|
|
|
passive = .true.
|
|
|
|
|
endif
|
|
|
|
|
|
2019-09-20 21:15:23 +05:30
|
|
|
if (dEq0(norm2(v))) then
|
|
|
|
|
vRot = v
|
2019-04-03 16:41:18 +05:30
|
|
|
else
|
2019-09-21 03:38:14 +05:30
|
|
|
v_normed = v/norm2(v)
|
2019-04-03 16:41:18 +05:30
|
|
|
if (passive) then
|
2019-09-20 21:15:23 +05:30
|
|
|
q = self%q * (quaternion([0.0_pReal, v_normed(1), v_normed(2), v_normed(3)]) * conjg(self%q) )
|
2019-04-03 16:41:18 +05:30
|
|
|
else
|
2019-09-20 21:15:23 +05:30
|
|
|
q = conjg(self%q) * (quaternion([0.0_pReal, v_normed(1), v_normed(2), v_normed(3)]) * self%q )
|
2019-04-03 16:41:18 +05:30
|
|
|
endif
|
2019-09-23 00:40:39 +05:30
|
|
|
vRot = q%aimag()*norm2(v)
|
2019-04-03 16:41:18 +05:30
|
|
|
endif
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
end function rotVector
|
|
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-09-21 05:09:33 +05:30
|
|
|
!> @brief rotate a rank-2 tensor passively (default) or actively
|
2019-02-01 14:31:54 +05:30
|
|
|
!> @details: rotation is based on rotation matrix
|
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-09-21 05:09:33 +05:30
|
|
|
pure function rotTensor2(self,T,active) result(tRot)
|
2019-02-01 14:31:54 +05:30
|
|
|
|
2019-09-21 05:09:33 +05:30
|
|
|
real(pReal), dimension(3,3) :: tRot
|
2019-04-03 16:41:18 +05:30
|
|
|
class(rotation), 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
|
|
|
|
|
|
|
|
|
|
logical :: passive
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
if (present(active)) then
|
|
|
|
|
passive = .not. active
|
|
|
|
|
else
|
|
|
|
|
passive = .true.
|
|
|
|
|
endif
|
|
|
|
|
|
|
|
|
|
if (passive) then
|
2019-09-21 05:48:09 +05:30
|
|
|
tRot = matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix()))
|
2019-04-03 16:41:18 +05:30
|
|
|
else
|
2019-09-21 05:48:09 +05:30
|
|
|
tRot = matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix())
|
2019-04-03 16:41:18 +05:30
|
|
|
endif
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-09-20 21:15:23 +05:30
|
|
|
end function rotTensor2
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-09-21 06:06:37 +05:30
|
|
|
|
2019-09-21 05:09:33 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
|
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
|
|
|
|
|
!> @brief rotate a rank-4 tensor passively (default) or actively
|
|
|
|
|
!> @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
|
|
|
|
|
class(rotation), intent(in) :: self
|
|
|
|
|
real(pReal), intent(in), dimension(3,3,3,3) :: T
|
|
|
|
|
logical, intent(in), optional :: active
|
|
|
|
|
|
|
|
|
|
real(pReal), dimension(3,3) :: R
|
|
|
|
|
integer :: i,j,k,l,m,n,o,p
|
|
|
|
|
|
|
|
|
|
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
|
|
|
|
|
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-09-21 06:06:37 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
|
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
|
|
|
|
|
!> @brief rotate a symmetric rank-4 tensor stored as (6,6) passively (default) or actively
|
|
|
|
|
!! ToDo: Need to check active/passive !!!
|
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
|
pure function rotTensor4sym(self,T,active) result(tRot)
|
|
|
|
|
|
|
|
|
|
real(pReal), dimension(6,6) :: tRot
|
|
|
|
|
class(rotation), intent(in) :: self
|
|
|
|
|
real(pReal), intent(in), dimension(6,6) :: T
|
|
|
|
|
logical, intent(in), optional :: active
|
|
|
|
|
|
|
|
|
|
if (present(active)) then
|
|
|
|
|
tRot = math_sym3333to66(rotTensor4(self,math_66toSym3333(T),active))
|
|
|
|
|
else
|
|
|
|
|
tRot = math_sym3333to66(rotTensor4(self,math_66toSym3333(T)))
|
|
|
|
|
endif
|
|
|
|
|
|
|
|
|
|
end function rotTensor4sym
|
|
|
|
|
|
|
|
|
|
|
2019-02-02 00:47:12 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
|
!> @brief misorientation
|
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-09-21 05:09:33 +05:30
|
|
|
pure elemental function misorientation(self,other)
|
2019-02-02 00:47:12 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
type(rotation) :: misorientation
|
|
|
|
|
class(rotation), intent(in) :: self, other
|
|
|
|
|
|
2020-01-14 16:03:18 +05:30
|
|
|
misorientation%q = other%q * conjg(self%q)
|
2019-02-02 00:47:12 +05:30
|
|
|
|
|
|
|
|
end function misorientation
|
|
|
|
|
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
|
|
|
|
!> @brief convert unit quaternion to rotation matrix
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function qu2om(qu) result(om)
|
|
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: qu
|
|
|
|
|
real(pReal), dimension(3,3) :: om
|
2019-04-17 01:47:56 +05:30
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal) :: qq
|
2019-04-17 01:47:56 +05:30
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
qq = qu(1)**2-sum(qu(2:4)**2)
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
|
2019-09-30 04:27:57 +05:30
|
|
|
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
|
2019-04-17 01:47:56 +05:30
|
|
|
|
2019-09-30 04:27:57 +05:30
|
|
|
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))
|
2019-04-17 01:47:56 +05:30
|
|
|
|
2019-09-30 04:27:57 +05:30
|
|
|
if (P < 0.0_pReal) om = transpose(om)
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
end function qu2om
|
2019-02-01 14:31:54 +05:30
|
|
|
|
|
|
|
|
|
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert unit quaternion to Euler angles
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function qu2eu(qu) result(eu)
|
|
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: qu
|
|
|
|
|
real(pReal), dimension(3) :: eu
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal) :: q12, q03, chi, chiInv
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
q03 = qu(1)**2+qu(4)**2
|
|
|
|
|
q12 = qu(2)**2+qu(3)**2
|
2019-04-17 01:47:56 +05:30
|
|
|
chi = sqrt(q03*q12)
|
|
|
|
|
|
|
|
|
|
degenerated: if (dEq0(chi)) then
|
2019-09-30 04:27:57 +05:30
|
|
|
eu = merge([atan2(-P*2.0_pReal*qu(1)*qu(4),qu(1)**2-qu(4)**2), 0.0_pReal, 0.0_pReal], &
|
|
|
|
|
[atan2( 2.0_pReal*qu(2)*qu(3),qu(2)**2-qu(3)**2), PI, 0.0_pReal], &
|
2019-04-17 01:47:56 +05:30
|
|
|
dEq0(q12))
|
|
|
|
|
else degenerated
|
2019-09-30 04:27:57 +05:30
|
|
|
chiInv = 1.0_pReal/chi
|
2019-09-20 21:06:16 +05:30
|
|
|
eu = [atan2((-P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)-qu(3)*qu(4))*chi ), &
|
2019-09-30 04:27:57 +05:30
|
|
|
atan2( 2.0_pReal*chi, q03-q12 ), &
|
2019-09-20 21:06:16 +05:30
|
|
|
atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )]
|
2019-04-17 01:47:56 +05:30
|
|
|
endif degenerated
|
|
|
|
|
where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function qu2eu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert unit quaternion to axis angle pair
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function qu2ax(qu) result(ax)
|
|
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: qu
|
|
|
|
|
real(pReal), dimension(4) :: ax
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal) :: omega, s
|
2019-04-17 01:47:56 +05:30
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
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]
|
2019-09-20 21:06:16 +05:30
|
|
|
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 ]
|
2019-04-17 01:47:56 +05:30
|
|
|
else
|
2019-09-20 21:06:16 +05:30
|
|
|
ax = [ qu(2), qu(3), qu(4), PI ]
|
2019-04-17 01:47:56 +05:30
|
|
|
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)
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: qu
|
|
|
|
|
real(pReal), dimension(4) :: ro
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal) :: s
|
|
|
|
|
real(pReal), parameter :: thr = 1.0e-8_pReal
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-23 02:47:17 +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
|
2019-09-20 21:06:16 +05:30
|
|
|
s = norm2(qu(2:4))
|
2019-04-17 01:47:56 +05:30
|
|
|
if (s < thr) then
|
2019-09-20 21:06:16 +05:30
|
|
|
ro = [0.0_pReal, 0.0_pReal, P, 0.0_pReal]
|
2019-04-17 01:47:56 +05:30
|
|
|
else
|
2019-09-20 21:06:16 +05:30
|
|
|
ro = [qu(2)/s,qu(3)/s,qu(4)/s, tan(acos(math_clip(qu(1),-1.0_pReal,1.0_pReal)))]
|
2019-04-17 01:47:56 +05:30
|
|
|
endif
|
|
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
end if
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
end function qu2ro
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert unit quaternion to homochoric
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function qu2ho(qu) result(ho)
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: qu
|
|
|
|
|
real(pReal), dimension(3) :: ho
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal) :: omega, f
|
2019-04-17 01:47:56 +05:30
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
omega = 2.0 * acos(math_clip(qu(1),-1.0_pReal,1.0_pReal))
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
if (dEq0(omega)) then
|
2019-09-30 04:27:57 +05:30
|
|
|
ho = [ 0.0_pReal, 0.0_pReal, 0.0_pReal ]
|
2019-04-03 16:41:18 +05:30
|
|
|
else
|
2019-09-20 21:06:16 +05:30
|
|
|
ho = qu(2:4)
|
2019-09-30 04:27:57 +05:30
|
|
|
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
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
end function qu2ho
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @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)
|
2019-04-17 01:47:56 +05:30
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: qu
|
|
|
|
|
real(pReal), dimension(3) :: cu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
cu = ho2cu(qu2ho(qu))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function qu2cu
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 16:21:00 +05:30
|
|
|
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
|
|
|
|
|
!> @brief convert rotation matrix to cubochoric
|
|
|
|
|
!> @details the original formulation (direct conversion) had (numerical?) issues
|
2019-04-17 01:47:56 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-09-20 11:20:30 +05:30
|
|
|
pure function om2qu(om) result(qu)
|
2019-04-17 01:47:56 +05:30
|
|
|
|
2019-04-17 16:21:00 +05:30
|
|
|
real(pReal), intent(in), dimension(3,3) :: om
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal), dimension(4) :: qu
|
2019-04-17 01:47:56 +05:30
|
|
|
|
2019-04-17 16:21:00 +05:30
|
|
|
qu = eu2qu(om2eu(om))
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
end function om2qu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-02-01 14:31:54 +05:30
|
|
|
!> @brief orientation matrix to Euler angles
|
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +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
|
|
|
|
|
|
2019-04-17 18:49:41 +05:30
|
|
|
if (abs(om(3,3)) < 1.0_pReal) then
|
2019-04-03 16:41:18 +05:30
|
|
|
zeta = 1.0_pReal/sqrt(1.0_pReal-om(3,3)**2.0_pReal)
|
|
|
|
|
eu = [atan2(om(3,1)*zeta,-om(3,2)*zeta), &
|
|
|
|
|
acos(om(3,3)), &
|
|
|
|
|
atan2(om(1,3)*zeta, om(2,3)*zeta)]
|
2019-04-17 18:49:41 +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
|
2019-04-17 18:49:41 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
end function om2eu
|
|
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert orientation matrix to axis angle pair
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
function om2ax(om) result(ax)
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-22 12:15:54 +05:30
|
|
|
real(pReal), intent(in), dimension(3,3) :: om
|
|
|
|
|
real(pReal), dimension(4) :: ax
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-22 12:15:54 +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
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-22 12:15:54 +05:30
|
|
|
external :: dgeev
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-22 12:15:54 +05:30
|
|
|
om_ = om
|
2019-04-03 16:41:18 +05:30
|
|
|
|
|
|
|
|
! first get the rotation angle
|
2019-09-22 12:15:54 +05:30
|
|
|
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))
|
|
|
|
|
|
|
|
|
|
if (dEq0(ax(4))) then
|
2019-09-30 04:27:57 +05:30
|
|
|
ax(1:3) = [ 0.0_pReal, 0.0_pReal, 1.0_pReal ]
|
2019-04-03 16:41:18 +05:30
|
|
|
else
|
2019-09-22 12:15:54 +05:30
|
|
|
call dgeev('N','V',3,om_,3,Wr,Wi,devNull,3,VR,3,work,size(work,1),ierr)
|
2020-01-30 03:23:19 +05:30
|
|
|
if (ierr /= 0) call IO_error(401,ext_msg='Error in om2ax: DGEEV return not zero')
|
2020-01-04 06:24:19 +05:30
|
|
|
#if defined(__GFORTRAN__) && __GNUC__<9 || defined(__INTEL_COMPILER) && INTEL_COMPILER<1800 || defined(__PGI)
|
2019-09-22 12:15:54 +05:30
|
|
|
i = maxloc(merge(1,0,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1)
|
|
|
|
|
#else
|
|
|
|
|
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)
|
|
|
|
|
#endif
|
2020-01-30 03:23:19 +05:30
|
|
|
if (i == 0) call IO_error(401,ext_msg='Error in om2ax Real: eigenvalue not found')
|
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
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
end function om2ax
|
|
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert rotation matrix to Rodrigues vector
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function om2ro(om) result(ro)
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(3,3) :: om
|
|
|
|
|
real(pReal), dimension(4) :: ro
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
ro = eu2ro(om2eu(om))
|
|
|
|
|
|
|
|
|
|
end function om2ro
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert rotation matrix to homochoric
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
function om2ho(om) result(ho)
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(3,3) :: om
|
|
|
|
|
real(pReal), dimension(3) :: ho
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
ho = ax2ho(om2ax(om))
|
|
|
|
|
|
|
|
|
|
end function om2ho
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert rotation matrix to cubochoric
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
function om2cu(om) result(cu)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(3,3) :: om
|
|
|
|
|
real(pReal), dimension(3) :: cu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
cu = ho2cu(om2ho(om))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function om2cu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief Euler angles to unit quaternion
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function eu2qu(eu) result(qu)
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal), intent(in), dimension(3) :: eu
|
|
|
|
|
real(pReal), dimension(4) :: qu
|
|
|
|
|
real(pReal), dimension(3) :: ee
|
|
|
|
|
real(pReal) :: cPhi, sPhi
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
ee = 0.5_pReal*eu
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
cPhi = cos(ee(2))
|
|
|
|
|
sPhi = sin(ee(2))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
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(qu(1) < 0.0_pReal) qu = qu * (-1.0_pReal)
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
end function eu2qu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief Euler angles to orientation matrix
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function eu2om(eu) result(om)
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(3) :: eu
|
|
|
|
|
real(pReal), dimension(3,3) :: om
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), dimension(3) :: c, s
|
|
|
|
|
|
|
|
|
|
c = cos(eu)
|
|
|
|
|
s = sin(eu)
|
|
|
|
|
|
|
|
|
|
om(1,1) = c(1)*c(3)-s(1)*s(3)*c(2)
|
|
|
|
|
om(1,2) = s(1)*c(3)+c(1)*s(3)*c(2)
|
|
|
|
|
om(1,3) = s(3)*s(2)
|
|
|
|
|
om(2,1) = -c(1)*s(3)-s(1)*c(3)*c(2)
|
|
|
|
|
om(2,2) = -s(1)*s(3)+c(1)*c(3)*c(2)
|
|
|
|
|
om(2,3) = c(3)*s(2)
|
|
|
|
|
om(3,1) = s(1)*s(2)
|
|
|
|
|
om(3,2) = -c(1)*s(2)
|
|
|
|
|
om(3,3) = c(2)
|
|
|
|
|
|
|
|
|
|
where(dEq0(om)) om = 0.0_pReal
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function eu2om
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert euler to axis angle
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function eu2ax(eu) result(ax)
|
|
|
|
|
|
|
|
|
|
real(pReal), intent(in), dimension(3) :: eu
|
|
|
|
|
real(pReal), dimension(4) :: ax
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal) :: t, delta, tau, alpha, sigma
|
|
|
|
|
|
2019-09-30 04:27:57 +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))
|
2019-04-17 01:47:56 +05:30
|
|
|
tau = sqrt(t**2+sin(sigma)**2)
|
|
|
|
|
|
2019-09-30 04:27:57 +05:30
|
|
|
alpha = merge(PI, 2.0_pReal*atan(tau/cos(sigma)), dEq(sigma,PI*0.5_pReal,tol=1.0e-15_pReal))
|
2019-04-17 01:47:56 +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
|
2019-04-17 01:47:56 +05:30
|
|
|
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
|
2019-09-30 04:27:57 +05:30
|
|
|
if (alpha < 0.0_pReal) ax = -ax ! ensure alpha is positive
|
2019-04-03 16:41:18 +05:30
|
|
|
end if
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
end function eu2ax
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief Euler angles to Rodrigues vector
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function eu2ro(eu) result(ro)
|
|
|
|
|
|
|
|
|
|
real(pReal), intent(in), dimension(3) :: eu
|
|
|
|
|
real(pReal), dimension(4) :: ro
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-04-17 01:47:56 +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
|
2019-09-30 04:27:57 +05:30
|
|
|
ro(4) = tan(ro(4)*0.5_pReal)
|
2019-04-03 16:41:18 +05:30
|
|
|
end if
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
end function eu2ro
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert Euler angles to homochoric
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function eu2ho(eu) result(ho)
|
2019-04-03 16:41:18 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(3) :: eu
|
|
|
|
|
real(pReal), dimension(3) :: ho
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
ho = ax2ho(eu2ax(eu))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function eu2ho
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert Euler angles to cubochoric
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
function eu2cu(eu) result(cu)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(3) :: eu
|
2019-04-03 16:41:18 +05:30
|
|
|
real(pReal), dimension(3) :: cu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
cu = ho2cu(eu2ho(eu))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function eu2cu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert axis angle pair to quaternion
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function ax2qu(ax) result(qu)
|
|
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: ax
|
|
|
|
|
real(pReal), dimension(4) :: qu
|
|
|
|
|
|
|
|
|
|
real(pReal) :: c, s
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
if (dEq0(ax(4))) then
|
2019-09-20 21:06:16 +05:30
|
|
|
qu = [ 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal ]
|
2019-04-17 01:47:56 +05:30
|
|
|
else
|
2019-09-30 04:27:57 +05:30
|
|
|
c = cos(ax(4)*0.5_pReal)
|
|
|
|
|
s = sin(ax(4)*0.5_pReal)
|
2019-09-20 21:06:16 +05:30
|
|
|
qu = [ c, ax(1)*s, ax(2)*s, ax(3)*s ]
|
2019-04-17 01:47:56 +05:30
|
|
|
end if
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function ax2qu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert axis angle pair to orientation matrix
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function ax2om(ax) result(om)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: ax
|
|
|
|
|
real(pReal), dimension(3,3) :: om
|
|
|
|
|
|
|
|
|
|
real(pReal) :: q, c, s, omc
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
c = cos(ax(4))
|
|
|
|
|
s = sin(ax(4))
|
2019-09-30 04:27:57 +05:30
|
|
|
omc = 1.0_pReal-c
|
2019-04-17 01:47:56 +05:30
|
|
|
|
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
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
q = omc*ax(1)*ax(2)
|
|
|
|
|
om(1,2) = q + s*ax(3)
|
|
|
|
|
om(2,1) = q - s*ax(3)
|
|
|
|
|
|
|
|
|
|
q = omc*ax(2)*ax(3)
|
|
|
|
|
om(2,3) = q + s*ax(1)
|
|
|
|
|
om(3,2) = q - s*ax(1)
|
|
|
|
|
|
|
|
|
|
q = omc*ax(3)*ax(1)
|
|
|
|
|
om(3,1) = q + s*ax(2)
|
|
|
|
|
om(1,3) = q - s*ax(2)
|
|
|
|
|
|
2019-09-30 04:27:57 +05:30
|
|
|
if (P > 0.0_pReal) om = transpose(om)
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
end function ax2om
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert axis angle pair to Euler angles
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function ax2eu(ax) result(eu)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: ax
|
|
|
|
|
real(pReal), dimension(3) :: eu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
eu = om2eu(ax2om(ax))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function ax2eu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert axis angle pair to Rodrigues vector
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function ax2ro(ax) result(ro)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: ax
|
|
|
|
|
real(pReal), dimension(4) :: ro
|
|
|
|
|
|
2019-09-30 04:27:57 +05:30
|
|
|
real(pReal), parameter :: thr = 1.0e-7_pReal
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
if (dEq0(ax(4))) then
|
|
|
|
|
ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal ]
|
|
|
|
|
else
|
|
|
|
|
ro(1:3) = ax(1:3)
|
|
|
|
|
! we need to deal with the 180 degree case
|
2019-09-30 04:27:57 +05:30
|
|
|
ro(4) = merge(IEEE_value(ro(4),IEEE_positive_inf),tan(ax(4)*0.5_pReal),abs(ax(4)-PI) < thr)
|
2019-04-17 01:47:56 +05:30
|
|
|
end if
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function ax2ro
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
|
!> @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
|
|
|
|
|
|
|
|
|
|
real(pReal) :: f
|
|
|
|
|
|
2019-09-30 04:27:57 +05:30
|
|
|
f = 0.75_pReal * ( ax(4) - sin(ax(4)) )
|
|
|
|
|
f = f**(1.0_pReal/3.0_pReal)
|
2019-04-17 01:47:56 +05:30
|
|
|
ho = ax(1:3) * f
|
|
|
|
|
|
|
|
|
|
end function ax2ho
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert axis angle pair to cubochoric
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
function ax2cu(ax) result(cu)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: ax
|
|
|
|
|
real(pReal), dimension(3) :: cu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
cu = ho2cu(ax2ho(ax))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function ax2cu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert Rodrigues vector to unit quaternion
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function ro2qu(ro) result(qu)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: ro
|
|
|
|
|
real(pReal), dimension(4) :: qu
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
qu = ax2qu(ro2ax(ro))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function ro2qu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-02-01 14:31:54 +05:30
|
|
|
!> @brief convert Rodrigues vector to rotation matrix
|
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
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
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
om = ax2om(ro2ax(ro))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
end function ro2om
|
|
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert Rodrigues vector to Euler angles
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function ro2eu(ro) result(eu)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: ro
|
|
|
|
|
real(pReal), dimension(3) :: eu
|
|
|
|
|
|
|
|
|
|
eu = om2eu(ro2om(ro))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function ro2eu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert Rodrigues vector to axis angle pair
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function ro2ax(ro) result(ax)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: ro
|
|
|
|
|
real(pReal), dimension(4) :: ax
|
|
|
|
|
|
|
|
|
|
real(pReal) :: ta, angle
|
|
|
|
|
|
|
|
|
|
ta = ro(4)
|
|
|
|
|
|
2019-09-23 03:16:27 +05:30
|
|
|
if (.not. IEEE_is_finite(ta)) then
|
2019-04-17 01:47:56 +05:30
|
|
|
ax = [ ro(1), ro(2), ro(3), PI ]
|
2019-09-23 03:16:27 +05:30
|
|
|
elseif (dEq0(ta)) then
|
2019-09-30 04:27:57 +05:30
|
|
|
ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ]
|
2019-04-17 01:47:56 +05:30
|
|
|
else
|
2019-09-30 04:27:57 +05:30
|
|
|
angle = 2.0_pReal*atan(ta)
|
|
|
|
|
ta = 1.0_pReal/norm2(ro(1:3))
|
2019-04-17 01:47:56 +05:30
|
|
|
ax = [ ro(1)/ta, ro(2)/ta, ro(3)/ta, angle ]
|
|
|
|
|
end if
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function ro2ax
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert Rodrigues vector to homochoric
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function ro2ho(ro) result(ho)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-09-22 12:15:54 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: ro
|
|
|
|
|
real(pReal), dimension(3) :: ho
|
2019-04-17 01:47:56 +05:30
|
|
|
|
2019-09-22 12:15:54 +05:30
|
|
|
real(pReal) :: f
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
if (dEq0(norm2(ro(1:3)))) then
|
2019-09-30 04:27:57 +05:30
|
|
|
ho = [ 0.0_pReal, 0.0_pReal, 0.0_pReal ]
|
2019-04-17 01:47:56 +05:30
|
|
|
else
|
2019-09-30 04:27:57 +05:30
|
|
|
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)
|
2019-04-17 01:47:56 +05:30
|
|
|
end if
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function ro2ho
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @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)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(4) :: ro
|
|
|
|
|
real(pReal), dimension(3) :: cu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
cu = ho2cu(ro2ho(ro))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function ro2cu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-02-01 14:31:54 +05:30
|
|
|
!> @brief convert homochoric to unit quaternion
|
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
pure function ho2qu(ho) result(qu)
|
|
|
|
|
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal), intent(in), dimension(3) :: ho
|
|
|
|
|
real(pReal), dimension(4) :: qu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
qu = ax2qu(ho2ax(ho))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
end function ho2qu
|
|
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert homochoric to rotation matrix
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function ho2om(ho) result(om)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(3) :: ho
|
|
|
|
|
real(pReal), dimension(3,3) :: om
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
om = ax2om(ho2ax(ho))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function ho2om
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert homochoric to Euler angles
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function ho2eu(ho) result(eu)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(3) :: ho
|
|
|
|
|
real(pReal), dimension(3) :: eu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
eu = ax2eu(ho2ax(ho))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function ho2eu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert homochoric to axis angle pair
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function ho2ax(ho) result(ax)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(3) :: ho
|
|
|
|
|
real(pReal), dimension(4) :: ax
|
|
|
|
|
|
|
|
|
|
integer :: i
|
|
|
|
|
real(pReal) :: hmag_squared, s, hm
|
|
|
|
|
real(pReal), parameter, dimension(16) :: &
|
|
|
|
|
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 ]
|
|
|
|
|
|
|
|
|
|
! normalize h and store the magnitude
|
|
|
|
|
hmag_squared = sum(ho**2.0_pReal)
|
|
|
|
|
if (dEq0(hmag_squared)) then
|
|
|
|
|
ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ]
|
|
|
|
|
else
|
|
|
|
|
hm = hmag_squared
|
|
|
|
|
|
|
|
|
|
! 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
|
|
|
|
|
end do
|
|
|
|
|
ax = [ho/sqrt(hmag_squared), 2.0_pReal*acos(s)]
|
|
|
|
|
end if
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function ho2ax
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert homochoric to Rodrigues vector
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-04-17 01:47:56 +05:30
|
|
|
pure function ho2ro(ho) result(ro)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(3) :: ho
|
|
|
|
|
real(pReal), dimension(4) :: ro
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
ro = ax2ro(ho2ax(ho))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function ho2ro
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert homochoric to cubochoric
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-09-21 05:09:33 +05:30
|
|
|
pure function ho2cu(ho) result(cu)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), intent(in), dimension(3) :: ho
|
|
|
|
|
real(pReal), dimension(3) :: cu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-05-17 02:26:48 +05:30
|
|
|
cu = Lambert_BallToCube(ho)
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
end function ho2cu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @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)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
real(pReal), intent(in), dimension(3) :: cu
|
2019-09-20 21:06:16 +05:30
|
|
|
real(pReal), dimension(4) :: qu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
qu = ho2qu(cu2ho(cu))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
end function cu2qu
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +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)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
real(pReal), intent(in), dimension(3) :: cu
|
|
|
|
|
real(pReal), dimension(3,3) :: om
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
om = ho2om(cu2ho(cu))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
end function cu2om
|
|
|
|
|
|
|
|
|
|
|
2019-04-17 01:47:56 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
|
!> @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)
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
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
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-02-01 14:31:54 +05:30
|
|
|
!> @brief convert cubochoric to axis angle pair
|
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
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
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
ax = ho2ax(cu2ho(cu))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
end function cu2ax
|
|
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +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)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
real(pReal), intent(in), dimension(3) :: cu
|
|
|
|
|
real(pReal), dimension(4) :: ro
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
ro = ho2ro(cu2ho(cu))
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
end function cu2ro
|
|
|
|
|
|
|
|
|
|
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2018-12-08 12:32:55 +05:30
|
|
|
!> @author Marc De Graef, Carnegie Mellon University
|
2019-04-17 01:47:56 +05:30
|
|
|
!> @brief convert cubochoric to homochoric
|
2019-02-01 14:31:54 +05:30
|
|
|
!---------------------------------------------------------------------------------------------------
|
2019-09-21 05:09:33 +05:30
|
|
|
pure function cu2ho(cu) result(ho)
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-04-03 16:41:18 +05:30
|
|
|
real(pReal), intent(in), dimension(3) :: cu
|
2019-04-17 01:47:56 +05:30
|
|
|
real(pReal), dimension(3) :: ho
|
2018-12-08 12:32:55 +05:30
|
|
|
|
2019-05-17 02:26:48 +05:30
|
|
|
ho = Lambert_CubeToBall(cu)
|
2019-04-17 01:47:56 +05:30
|
|
|
|
|
|
|
|
end function cu2ho
|
2018-12-08 12:32:55 +05:30
|
|
|
|
|
|
|
|
|
2019-09-22 19:23:03 +05:30
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
|
|
|
!> @brief check correctness of (some) rotations functions
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
|
|
|
subroutine unitTest
|
|
|
|
|
|
2019-09-23 10:28:18 +05:30
|
|
|
type(rotation) :: 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
|
|
|
|
|
character(len=pStringLen) :: msg
|
2019-09-22 19:23:03 +05:30
|
|
|
real :: A,B
|
|
|
|
|
integer :: i
|
|
|
|
|
|
2019-09-22 23:59:34 +05:30
|
|
|
do i=1,10
|
2019-09-22 19:23:03 +05:30
|
|
|
|
|
|
|
|
msg = ''
|
|
|
|
|
|
2019-09-23 03:16:27 +05:30
|
|
|
#if defined(__GFORTRAN__) && __GNUC__<9
|
2019-09-23 04:32:23 +05:30
|
|
|
if(i<7) cycle
|
2019-09-23 03:16:27 +05:30
|
|
|
#endif
|
|
|
|
|
|
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([2.0_pReal*PI,PI,2.0_pReal*PI])
|
|
|
|
|
elseif(i==4) then
|
|
|
|
|
qu = [0.0_pReal,0.0_pReal,1.0_pReal,0.0_pReal]
|
|
|
|
|
elseif(i==5) then
|
2019-09-23 03:16:27 +05:30
|
|
|
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
|
2019-09-23 10:28:18 +05:30
|
|
|
call random_number(x)
|
|
|
|
|
A = sqrt(x(3))
|
|
|
|
|
B = sqrt(1-0_pReal -x(3))
|
|
|
|
|
qu = [cos(2.0_pReal*PI*x(1))*A,&
|
|
|
|
|
sin(2.0_pReal*PI*x(2))*B,&
|
|
|
|
|
cos(2.0_pReal*PI*x(2))*B,&
|
|
|
|
|
sin(2.0_pReal*PI*x(1))*A]
|
2019-09-22 23:59:34 +05:30
|
|
|
if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
|
|
|
|
|
endif
|
2020-01-04 06:24:19 +05:30
|
|
|
#ifndef __PGI
|
2019-09-22 23:45:27 +05:30
|
|
|
if(dNeq0(norm2(om2qu(qu2om(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'om2qu/qu2om,'
|
|
|
|
|
if(dNeq0(norm2(eu2qu(qu2eu(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'eu2qu/qu2eu,'
|
|
|
|
|
if(dNeq0(norm2(ax2qu(qu2ax(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'ax2qu/qu2ax,'
|
|
|
|
|
if(dNeq0(norm2(ro2qu(qu2ro(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2qu/qu2ro,'
|
2019-09-23 10:28:18 +05:30
|
|
|
if(dNeq0(norm2(ho2qu(qu2ho(qu))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2qu/qu2ho,'
|
2019-09-22 23:45:27 +05:30
|
|
|
if(dNeq0(norm2(cu2qu(qu2cu(qu))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2qu/qu2cu,'
|
2020-01-04 06:24:19 +05:30
|
|
|
#endif
|
|
|
|
|
|
2019-09-22 23:45:27 +05:30
|
|
|
om = qu2om(qu)
|
2020-01-04 06:24:19 +05:30
|
|
|
#ifndef __PGI
|
2019-09-23 05:26:23 +05:30
|
|
|
if(dNeq0(norm2(om2qu(eu2om(om2eu(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'eu2om/om2eu,'
|
2019-09-22 23:45:27 +05:30
|
|
|
if(dNeq0(norm2(om2qu(ax2om(om2ax(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ax2om/om2ax,'
|
|
|
|
|
if(dNeq0(norm2(om2qu(ro2om(om2ro(om)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2om/om2ro,'
|
2019-09-23 10:28:18 +05:30
|
|
|
if(dNeq0(norm2(om2qu(ho2om(om2ho(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2om/om2ho,'
|
2019-09-22 23:45:27 +05:30
|
|
|
if(dNeq0(norm2(om2qu(cu2om(om2cu(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2om/om2cu,'
|
2020-01-04 06:24:19 +05:30
|
|
|
#endif
|
|
|
|
|
|
2019-09-22 23:45:27 +05:30
|
|
|
eu = qu2eu(qu)
|
2020-01-04 06:24:19 +05:30
|
|
|
#ifndef __PGI
|
2019-09-22 23:45:27 +05:30
|
|
|
if(dNeq0(norm2(eu2qu(ax2eu(eu2ax(eu)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ax2eu/eu2ax,'
|
|
|
|
|
if(dNeq0(norm2(eu2qu(ro2eu(eu2ro(eu)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2eu/eu2ro,'
|
2019-09-23 10:28:18 +05:30
|
|
|
if(dNeq0(norm2(eu2qu(ho2eu(eu2ho(eu)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2eu/eu2ho,'
|
2019-09-22 23:45:27 +05:30
|
|
|
if(dNeq0(norm2(eu2qu(cu2eu(eu2cu(eu)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2eu/eu2cu,'
|
2020-01-04 06:24:19 +05:30
|
|
|
#endif
|
|
|
|
|
|
2019-09-22 23:45:27 +05:30
|
|
|
ax = qu2ax(qu)
|
2020-01-04 06:24:19 +05:30
|
|
|
#ifndef __PGI
|
2019-09-22 23:45:27 +05:30
|
|
|
if(dNeq0(norm2(ax2qu(ro2ax(ax2ro(ax)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2ax/ax2ro,'
|
2019-09-23 10:28:18 +05:30
|
|
|
if(dNeq0(norm2(ax2qu(ho2ax(ax2ho(ax)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2ax/ax2ho,'
|
2019-09-22 23:45:27 +05:30
|
|
|
if(dNeq0(norm2(ax2qu(cu2ax(ax2cu(ax)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2ax/ax2cu,'
|
2020-01-04 06:24:19 +05:30
|
|
|
#endif
|
|
|
|
|
|
2019-09-22 23:45:27 +05:30
|
|
|
ro = qu2ro(qu)
|
2020-01-04 06:24:19 +05:30
|
|
|
#ifndef __PGI
|
2019-09-23 10:28:18 +05:30
|
|
|
if(dNeq0(norm2(ro2qu(ho2ro(ro2ho(ro)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2ro/ro2ho,'
|
2019-09-22 23:45:27 +05:30
|
|
|
if(dNeq0(norm2(ro2qu(cu2ro(ro2cu(ro)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2ro/ro2cu,'
|
2020-01-04 06:24:19 +05:30
|
|
|
#endif
|
|
|
|
|
|
2019-09-22 23:45:27 +05:30
|
|
|
ho = qu2ho(qu)
|
2020-01-04 06:24:19 +05:30
|
|
|
#ifndef __PGI
|
2019-09-22 23:45:27 +05:30
|
|
|
if(dNeq0(norm2(ho2qu(cu2ho(ho2cu(ho)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2ho/ho2cu,'
|
2020-01-04 06:24:19 +05:30
|
|
|
#endif
|
|
|
|
|
|
2019-09-23 10:28:18 +05:30
|
|
|
call R%fromMatrix(om)
|
|
|
|
|
|
|
|
|
|
call random_number(v3)
|
|
|
|
|
if(all(dNeq(R%rotVector(R%rotVector(v3),active=.true.),v3,1.0e-12_pReal))) &
|
|
|
|
|
msg = trim(msg)//'rotVector,'
|
|
|
|
|
|
|
|
|
|
call random_number(t33)
|
|
|
|
|
if(all(dNeq(R%rotTensor2(R%rotTensor2(t33),active=.true.),t33,1.0e-12_pReal))) &
|
|
|
|
|
msg = trim(msg)//'rotTensor2,'
|
|
|
|
|
|
|
|
|
|
call random_number(t3333)
|
|
|
|
|
if(all(dNeq(R%rotTensor4(R%rotTensor4(t3333),active=.true.),t3333,1.0e-12_pReal))) &
|
|
|
|
|
msg = trim(msg)//'rotTensor4,'
|
2019-09-22 23:45:27 +05:30
|
|
|
|
2020-01-27 01:23:13 +05:30
|
|
|
if(len_trim(msg) /= 0) call IO_error(0,ext_msg=msg)
|
2019-09-22 19:23:03 +05:30
|
|
|
|
|
|
|
|
enddo
|
2019-09-23 00:40:39 +05:30
|
|
|
|
2019-09-22 19:23:03 +05:30
|
|
|
end subroutine unitTest
|
|
|
|
|
|
|
|
|
|
|
2018-12-08 12:32:55 +05:30
|
|
|
end module rotations
|
