not needed

This commit is contained in:
Martin Diehl 2022-06-11 23:04:21 +02:00
parent bbd4f01fb7
commit 9e0a0ee166
1 changed files with 10 additions and 666 deletions

View File

@ -61,7 +61,6 @@ module rotations
procedure, public :: asQuaternion procedure, public :: asQuaternion
procedure, public :: asEulers procedure, public :: asEulers
procedure, public :: asAxisAngle procedure, public :: asAxisAngle
procedure, public :: asRodrigues
procedure, public :: asMatrix procedure, public :: asMatrix
!------------------------------------------ !------------------------------------------
procedure, public :: fromQuaternion procedure, public :: fromQuaternion
@ -150,24 +149,7 @@ pure function asMatrix(self)
asMatrix = qu2om(self%q) asMatrix = qu2om(self%q)
end function asMatrix end function asMatrix
!---------------------------------------------------------------------------------------------------
pure function asRodrigues(self)
class(tRotation), intent(in) :: self
real(pReal), dimension(4) :: asRodrigues
asRodrigues = qu2ro(self%q)
end function asRodrigues
!---------------------------------------------------------------------------------------------------
pure function asHomochoric(self)
class(tRotation), intent(in) :: self
real(pReal), dimension(3) :: asHomochoric
asHomochoric = qu2ho(self%q)
end function asHomochoric
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
! Initialize rotation from different representations ! Initialize rotation from different representations
@ -508,71 +490,6 @@ pure function qu2ax(qu) result(ax)
end function qu2ax end function qu2ax
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert unit quaternion to Rodrigues vector
!---------------------------------------------------------------------------------------------------
pure function qu2ro(qu) result(ro)
real(pReal), intent(in), dimension(4) :: qu
real(pReal), dimension(4) :: ro
real(pReal) :: s
real(pReal), parameter :: thr = 1.0e-8_pReal
if (abs(qu(1)) < thr) then
ro = [qu(2), qu(3), qu(4), IEEE_value(1.0_pReal,IEEE_positive_inf)]
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
end if
end function qu2ro
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert unit quaternion to homochoric
!---------------------------------------------------------------------------------------------------
pure function qu2ho(qu) result(ho)
real(pReal), intent(in), dimension(4) :: qu
real(pReal), dimension(3) :: ho
real(pReal) :: omega, f
omega = 2.0 * acos(math_clip(qu(1),-1.0_pReal,1.0_pReal))
if (dEq0(omega,tol=1.e-5_pReal)) then
ho = [ 0.0_pReal, 0.0_pReal, 0.0_pReal ]
else
ho = qu(2:4)
f = 0.75_pReal * ( omega - sin(omega) )
ho = ho/norm2(ho)* f**(1.0_pReal/3.0_pReal)
end if
end function qu2ho
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert unit quaternion to cubochoric
!---------------------------------------------------------------------------------------------------
pure function qu2cu(qu) result(cu)
real(pReal), intent(in), dimension(4) :: qu
real(pReal), dimension(3) :: cu
cu = ho2cu(qu2ho(qu))
end function qu2cu
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief convert rotation matrix to unit quaternion !> @brief convert rotation matrix to unit quaternion
@ -669,48 +586,6 @@ function om2ax(om) result(ax)
end function om2ax end function om2ax
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert rotation matrix to Rodrigues vector
!---------------------------------------------------------------------------------------------------
pure function om2ro(om) result(ro)
real(pReal), intent(in), dimension(3,3) :: om
real(pReal), dimension(4) :: ro
ro = eu2ro(om2eu(om))
end function om2ro
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert rotation matrix to homochoric
!---------------------------------------------------------------------------------------------------
function om2ho(om) result(ho)
real(pReal), intent(in), dimension(3,3) :: om
real(pReal), dimension(3) :: ho
ho = ax2ho(om2ax(om))
end function om2ho
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert rotation matrix to cubochoric
!---------------------------------------------------------------------------------------------------
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
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University !> @author Marc De Graef, Carnegie Mellon University
!> @brief Euler angles to unit quaternion !> @brief Euler angles to unit quaternion
@ -794,55 +669,6 @@ pure function eu2ax(eu) result(ax)
end function eu2ax end function eu2ax
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief Euler angles to Rodrigues vector
!---------------------------------------------------------------------------------------------------
pure function eu2ro(eu) result(ro)
real(pReal), intent(in), dimension(3) :: eu
real(pReal), dimension(4) :: ro
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 ]
else
ro(4) = tan(ro(4)*0.5_pReal)
end if
end function eu2ro
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Euler angles to homochoric
!---------------------------------------------------------------------------------------------------
pure function eu2ho(eu) result(ho)
real(pReal), intent(in), dimension(3) :: eu
real(pReal), dimension(3) :: ho
ho = ax2ho(eu2ax(eu))
end function eu2ho
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Euler angles to cubochoric
!---------------------------------------------------------------------------------------------------
function eu2cu(eu) result(cu)
real(pReal), intent(in), dimension(3) :: eu
real(pReal), dimension(3) :: cu
cu = ho2cu(eu2ho(eu))
end function eu2cu
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University !> @author Marc De Graef, Carnegie Mellon University
!> @brief convert axis angle pair to quaternion !> @brief convert axis angle pair to quaternion
@ -916,465 +742,6 @@ pure function ax2eu(ax) result(eu)
end function ax2eu end function ax2eu
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert axis angle pair to Rodrigues vector
!---------------------------------------------------------------------------------------------------
pure function ax2ro(ax) result(ro)
real(pReal), intent(in), dimension(4) :: ax
real(pReal), dimension(4) :: ro
real(pReal), parameter :: thr = 1.0e-7_pReal
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
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
real(pReal) :: f
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
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert axis angle pair to cubochoric
!---------------------------------------------------------------------------------------------------
function ax2cu(ax) result(cu)
real(pReal), intent(in), dimension(4) :: ax
real(pReal), dimension(3) :: cu
cu = ho2cu(ax2ho(ax))
end function ax2cu
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Rodrigues vector to unit quaternion
!---------------------------------------------------------------------------------------------------
pure function ro2qu(ro) result(qu)
real(pReal), intent(in), dimension(4) :: ro
real(pReal), dimension(4) :: qu
qu = ax2qu(ro2ax(ro))
end function ro2qu
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Rodrigues vector to rotation matrix
!---------------------------------------------------------------------------------------------------
pure function ro2om(ro) result(om)
real(pReal), intent(in), dimension(4) :: ro
real(pReal), dimension(3,3) :: om
om = ax2om(ro2ax(ro))
end function ro2om
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Rodrigues vector to Euler angles
!---------------------------------------------------------------------------------------------------
pure function ro2eu(ro) result(eu)
real(pReal), intent(in), dimension(4) :: ro
real(pReal), dimension(3) :: eu
eu = om2eu(ro2om(ro))
end function ro2eu
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Rodrigues vector to axis angle pair
!---------------------------------------------------------------------------------------------------
pure function ro2ax(ro) result(ax)
real(pReal), intent(in), dimension(4) :: ro
real(pReal), dimension(4) :: ax
real(pReal) :: ta, angle
ta = ro(4)
if (.not. IEEE_is_finite(ta)) then
ax = [ ro(1), ro(2), ro(3), PI ]
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
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Rodrigues vector to homochoric
!---------------------------------------------------------------------------------------------------
pure function ro2ho(ro) result(ho)
real(pReal), intent(in), dimension(4) :: ro
real(pReal), dimension(3) :: ho
real(pReal) :: f
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
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Rodrigues vector to cubochoric
!---------------------------------------------------------------------------------------------------
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
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert homochoric to unit quaternion
!---------------------------------------------------------------------------------------------------
pure function ho2qu(ho) result(qu)
real(pReal), intent(in), dimension(3) :: ho
real(pReal), dimension(4) :: qu
qu = ax2qu(ho2ax(ho))
end function ho2qu
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert homochoric to rotation matrix
!---------------------------------------------------------------------------------------------------
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
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert homochoric to Euler angles
!---------------------------------------------------------------------------------------------------
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
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert homochoric to axis angle pair
!---------------------------------------------------------------------------------------------------
pure function ho2ax(ho) result(ax)
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)
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
enddo
ax = [ho/sqrt(hmag_squared), 2.0_pReal*acos(s)]
end if
end function ho2ax
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert homochoric to Rodrigues vector
!---------------------------------------------------------------------------------------------------
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
!--------------------------------------------------------------------------
!> @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)
q = (BETA/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2))
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
xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / PREF ]/SC
! reverse the coordinates back to order according to the original pyramid number
cu = xyz1(p(:,2))
endif center
end function ho2cu
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert cubochoric to unit quaternion
!---------------------------------------------------------------------------------------------------
pure function cu2qu(cu) result(qu)
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(4) :: qu
qu = ho2qu(cu2ho(cu))
end function cu2qu
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert cubochoric to rotation matrix
!---------------------------------------------------------------------------------------------------
pure function cu2om(cu) result(om)
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(3,3) :: om
om = ho2om(cu2ho(cu))
end function cu2om
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert cubochoric to Euler angles
!---------------------------------------------------------------------------------------------------
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
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert cubochoric to axis angle pair
!---------------------------------------------------------------------------------------------------
function cu2ax(cu) result(ax)
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(4) :: ax
ax = ho2ax(cu2ho(cu))
end function cu2ax
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert cubochoric to Rodrigues vector
!---------------------------------------------------------------------------------------------------
pure function cu2ro(cu) result(ro)
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(4) :: ro
ro = ho2ro(cu2ho(cu))
end function cu2ro
!--------------------------------------------------------------------------
!> @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
!--------------------------------------------------------------------------
pure function cu2ho(cu) result(ho)
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(3) :: ho, LamXYZ, XYZ
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
if (maxval(abs(cu)) > AP/2.0+eps) then
ho = IEEE_value(cu,IEEE_positive_inf)
return
end if
! transform to the sphere grid via the curved square, and intercept the zero point
center: if (all(dEq0(cu))) then
ho = 0.0_pReal
else center
! get pyramide and scale by grid parameter ratio
p = GetPyramidOrder(cu)
XYZ = cu(p(:,1)) * SC
! intercept all the points along the z-axis
special: if (all(dEq0(XYZ(1:2)))) then
LamXYZ = [ 0.0_pReal, 0.0_pReal, PREF * XYZ(3) ]
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)
q = PREK * XYZ(order(2))/ sqrt(R2-c)
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)
s = c * PI/(24.0*XYZ(3)**2)
c = c * sqrt(PI/24.0_pReal) / XYZ(3)
q = sqrt( 1.0 - s )
LamXYZ = [ T(order(2)) * q, T(order(1)) * q, PREF * XYZ(3) - c ]
end if special
! reverse the coordinates back to order according to the original pyramid number
ho = LamXYZ(p(:,2))
end if center
end function cu2ho
!--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Determine to which pyramid a point in a cubic grid belongs.
!--------------------------------------------------------------------------
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
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Multiply two quaternions. !> @brief Multiply two quaternions.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1425,17 +792,15 @@ subroutine selfTest()
do i = 1, 20 do i = 1, 20
if(i==1) then if(i==1) then
qu = om2qu(math_I3) qu = [1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal]
elseif(i==2) then elseif(i==2) then
qu = eu2qu([0.0_pReal,0.0_pReal,0.0_pReal]) qu = [1.0_pReal,-0.0_pReal,-0.0_pReal,-0.0_pReal]
elseif(i==3) then elseif(i==3) then
qu = eu2qu([TAU,PI,TAU]) qu = [0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal]
elseif(i==4) then elseif(i==4) then
qu = [0.0_pReal,0.0_pReal,1.0_pReal,0.0_pReal] qu = [0.0_pReal,0.0_pReal,1.0_pReal,0.0_pReal]
elseif(i==5) then elseif(i==5) then
qu = ro2qu([1.0_pReal,0.0_pReal,0.0_pReal,IEEE_value(1.0_pReal, IEEE_positive_inf)]) qu = [0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal]
elseif(i==6) then
qu = ax2qu([1.0_pReal,0.0_pReal,0.0_pReal,0.0_pReal])
else else
call random_number(x) call random_number(x)
A = sqrt(x(3)) A = sqrt(x(3))
@ -1448,37 +813,16 @@ subroutine selfTest()
endif endif
if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) error stop 'om2qu/qu2om' if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) error stop 'om2qu2om'
if(.not. quaternion_equal(eu2qu(qu2eu(qu)),qu)) error stop 'eu2qu/qu2eu' if(.not. quaternion_equal(eu2qu(qu2eu(qu)),qu)) error stop 'eu2qu2eu'
if(.not. quaternion_equal(ax2qu(qu2ax(qu)),qu)) error stop 'ax2qu/qu2ax' if(.not. quaternion_equal(ax2qu(qu2ax(qu)),qu)) error stop 'ax2qu2ax'
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) om = qu2om(qu)
if(.not. quaternion_equal(om2qu(eu2om(om2eu(om))),qu)) error stop 'eu2om/om2eu' if(.not. quaternion_equal(om2qu(eu2om(om2eu(om))),qu)) error stop 'eu2om2eu'
if(.not. quaternion_equal(om2qu(ax2om(om2ax(om))),qu)) error stop 'ax2om/om2ax' if(.not. quaternion_equal(om2qu(ax2om(om2ax(om))),qu)) error stop 'ax2om2ax'
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) eu = qu2eu(qu)
if(.not. quaternion_equal(eu2qu(ax2eu(eu2ax(eu))),qu)) error stop 'ax2eu/eu2ax' if(.not. quaternion_equal(eu2qu(ax2eu(eu2ax(eu))),qu)) error stop 'ax2eu2ax'
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) call R%fromMatrix(om)