quaternions.f90

Go to the documentation of this file.

# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/quaternions.f90"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/quaternions.f90"
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
module quaternions
  use prec
  use io
 
  implicit none
  private
 
  real(preal), parameter, public :: p = -1.0_preal                                                  
 
  type, public :: quaternion
    real(preal), private :: w = 0.0_preal
    real(preal), private :: x = 0.0_preal
    real(preal), private :: y = 0.0_preal
    real(preal), private :: z = 0.0_preal
 
 
  contains
    procedure, private :: add__
    procedure, private :: pos__
    generic,   public  :: operator(+) => add__,pos__
 
    procedure, private :: sub__
    procedure, private :: neg__
    generic,   public  :: operator(-) => sub__,neg__
 
    procedure, private :: mul_quat__
    procedure, private :: mul_scal__
    generic,   public  :: operator(*) => mul_quat__, mul_scal__
 
    procedure, private :: div_quat__
    procedure, private :: div_scal__
    generic,   public  :: operator(/) => div_quat__, div_scal__
 
    procedure, private :: eq__
    generic,   public  :: operator(==) => eq__
 
    procedure, private :: neq__
    generic,   public  :: operator(/=) => neq__
 
    procedure, private :: pow_quat__
    procedure, private :: pow_scal__
    generic,   public  :: operator(**) => pow_quat__, pow_scal__
 
    procedure, public  :: abs   => abs__
    procedure, public  :: conjg => conjg__
    procedure, public  :: real  => real__
    procedure, public  :: aimag => aimag__
 
    procedure, public  :: homomorphed
    procedure, public  :: asarray
    procedure, public  :: inverse
 
  end type
 
  interface assignment (=)
    module procedure assign_quat__
    module procedure assign_vec__
  end interface assignment (=)
  
  interface quaternion
    module procedure init__
  end interface quaternion
  
  interface abs
    procedure abs__
  end interface abs
  
  interface dot_product
    procedure dot_product__
  end interface dot_product
  
  interface conjg
    module procedure conjg__
  end interface conjg
  
  interface exp
    module procedure exp__
  end interface exp
  
  interface log
    module procedure log__
  end interface log
 
  interface real
    module procedure real__
  end interface real
 
  interface aimag
    module procedure aimag__
  end interface aimag
 
  public :: &
    quaternions_init, &
    assignment(=), &
    conjg, aimag, &
    log, exp, &
    real
 
contains
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine quaternions_init
 
  write(6,'(/,a)') ' <<<+-  quaternions init  -+>>>'; flush(6)
  call unittest
 
end subroutine quaternions_init
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) pure function init__(array)
 
  real(preal), intent(in), dimension(4) :: array
 
  init__%w = array(1)
  init__%x = array(2)
  init__%y = array(3)
  init__%z = array(4)
 
end function init__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
elemental pure subroutine assign_quat__(self,other)
 
  type(quaternion), intent(out) :: self
  type(quaternion), intent(in)  :: other
 
  self = [other%w,other%x,other%y,other%z]
  
end subroutine assign_quat__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
pure subroutine assign_vec__(self,other)
 
  type(quaternion), intent(out)                :: self
  real(preal),       intent(in), dimension(4)  :: other
 
  self%w = other(1)
  self%x = other(2)
  self%y = other(3)
  self%z = other(4)
 
end subroutine assign_vec__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function add__(self,other)
 
  class(quaternion), intent(in) :: self,other
 
  add__ = [ self%w,  self%x,  self%y ,self%z] &
        + [other%w, other%x, other%y,other%z]
  
end function add__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function pos__(self)
 
  class(quaternion), intent(in) :: self
 
  pos__ = self * (+1.0_preal)
  
end function pos__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function sub__(self,other)
 
  class(quaternion), intent(in) :: self,other
 
  sub__ = [ self%w,  self%x,  self%y ,self%z] &
        - [other%w, other%x, other%y,other%z]
  
end function sub__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function neg__(self)
 
  class(quaternion), intent(in) :: self
 
  neg__ = self * (-1.0_preal)
  
end function neg__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function mul_quat__(self,other)
 
  class(quaternion), intent(in) :: self, other
 
  mul_quat__%w = self%w*other%w - self%x*other%x -      self%y*other%y - self%z*other%z
  mul_quat__%x = self%w*other%x + self%x*other%w + p * (self%y*other%z - self%z*other%y)
  mul_quat__%y = self%w*other%y + self%y*other%w + p * (self%z*other%x - self%x*other%z)
  mul_quat__%z = self%w*other%z + self%z*other%w + p * (self%x*other%y - self%y*other%x)
 
end function mul_quat__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function mul_scal__(self,scal)
 
  class(quaternion), intent(in) :: self
  real(preal),       intent(in) :: scal
 
  mul_scal__ = [self%w,self%x,self%y,self%z]*scal
  
end function mul_scal__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function div_quat__(self,other)
 
  class(quaternion), intent(in) :: self, other
 
  div_quat__ = self * (conjg(other)/(abs(other)**2.0_preal))
 
end function div_quat__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function div_scal__(self,scal)
 
  class(quaternion), intent(in) :: self
  real(preal),       intent(in) :: scal
 
  div_scal__ = [self%w,self%x,self%y,self%z]/scal
 
end function div_scal__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
logical elemental pure function eq__(self,other)
 
  class(quaternion), intent(in) :: self,other
 
  eq__ = all(deq([ self%w, self%x, self%y, self%z], &
                 [other%w,other%x,other%y,other%z]))
 
end function eq__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
logical elemental pure function neq__(self,other)
 
  class(quaternion), intent(in) :: self,other
 
  neq__ = .not. self%eq__(other)
 
end function neq__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function pow_quat__(self,expon)
 
  class(quaternion), intent(in) :: self
  type(quaternion),  intent(in) :: expon
 
  pow_quat__ = exp(log(self)*expon)
 
end function pow_quat__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function pow_scal__(self,expon)
 
  class(quaternion), intent(in) :: self
  real(preal),       intent(in) :: expon
 
  pow_scal__ = exp(log(self)*expon)
 
end function pow_scal__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function exp__(a)
 
  class(quaternion), intent(in) :: a
  real(preal)                   :: absimag
 
  absimag = norm2(aimag(a))
 
  exp__ = merge(exp(a%w) * [               cos(absimag), &
                             a%x/absimag * sin(absimag), &
                             a%y/absimag * sin(absimag), &
                             a%z/absimag * sin(absimag)], &
                ieee_value(1.0_preal,ieee_signaling_nan), &
                dneq0(absimag))
 
end function exp__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function log__(a)
 
  class(quaternion), intent(in) :: a
  real(preal)                   :: absimag
 
  absimag = norm2(aimag(a))
 
  log__ = merge([log(abs(a)), &
                 a%x/absimag * acos(a%w/abs(a)), &
                 a%y/absimag * acos(a%w/abs(a)), &
                 a%z/absimag * acos(a%w/abs(a))], &
                ieee_value(1.0_preal,ieee_signaling_nan), &
                dneq0(absimag))
 
end function log__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
real(preal) elemental pure function abs__(self)
 
  class(quaternion), intent(in) :: self
 
  abs__ = norm2([self%w,self%x,self%y,self%z])
 
end function abs__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
real(preal) elemental pure function dot_product__(a,b)
 
  class(quaternion), intent(in) :: a,b
 
  dot_product__ = a%w*b%w + a%x*b%x + a%y*b%y + a%z*b%z
 
end function dot_product__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function conjg__(self)
 
  class(quaternion), intent(in) :: self
 
  conjg__ = [self%w,-self%x,-self%y,-self%z]
 
end function conjg__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function homomorphed(self)
 
  class(quaternion), intent(in) :: self
 
  homomorphed = - self
 
end function homomorphed
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
pure function asarray(self)
 
  real(preal), dimension(4)     :: asarray
  class(quaternion), intent(in) :: self
 
  asarray = [self%w,self%x,self%y,self%z]
 
end function asarray
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
pure function real__(self)
 
  real(preal)                   :: real__
  class(quaternion), intent(in) :: self
 
  real__ = self%w
 
end function real__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
pure function aimag__(self)
 
  real(preal), dimension(3)     :: aimag__
  class(quaternion), intent(in) :: self
 
  aimag__ = [self%x,self%y,self%z]
 
end function aimag__
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function inverse(self)
 
  class(quaternion), intent(in) :: self
 
  inverse = conjg(self)/abs(self)**2.0_preal
 
end function inverse
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine unittest
 
  real(pReal), dimension(4) :: qu
  type(quaternion)          :: q, q_2
 
  call random_number(qu)
  qu = (qu-0.5_preal) * 2.0_preal
  q  = quaternion(qu)
 
  q_2= qu
  if(any(dneq(q%asArray(),q_2%asArray())))             call io_error(0,ext_msg='assign_vec__')
 
  q_2 = q + q
  if(any(dneq(q_2%asArray(),2.0_preal*qu)))            call io_error(0,ext_msg='add__')
 
  q_2 = q - q
  if(any(dneq0(q_2%asArray())))                        call io_error(0,ext_msg='sub__')
 
  q_2 = q * 5.0_preal
  if(any(dneq(q_2%asArray(),5.0_preal*qu)))            call io_error(0,ext_msg='mul__')
 
  q_2 = q / 0.5_preal
  if(any(dneq(q_2%asArray(),2.0_preal*qu)))            call io_error(0,ext_msg='div__')
 
  q_2 = q * 0.3_preal
  if(dneq0(abs(q)) .and. q_2 == q)                     call io_error(0,ext_msg='eq__')
 
  q_2 = q
  if(q_2 /= q)                                         call io_error(0,ext_msg='neq__')
 
  if(dneq(abs(q),norm2(qu)))                           call io_error(0,ext_msg='abs__')
  if(dneq(abs(q)**2.0_preal, real(q*q%conjg()),1.0e-14_preal)) &
                                                       call io_error(0,ext_msg='abs__/*conjg')
 
  if(any(dneq(q%asArray(),qu)))                        call io_error(0,ext_msg='eq__')
  if(dneq(q%real(),       qu(1)))                      call io_error(0,ext_msg='real()')
  if(any(dneq(q%aimag(),  qu(2:4))))                   call io_error(0,ext_msg='aimag()')
 
  q_2 = q%homomorphed()
  if(q                 /= q_2*    (-1.0_preal))        call io_error(0,ext_msg='homomorphed')
  if(dneq(q_2%real(),     qu(1)*  (-1.0_preal)))       call io_error(0,ext_msg='homomorphed/real')
  if(any(dneq(q_2%aimag(),qu(2:4)*(-1.0_preal))))      call io_error(0,ext_msg='homomorphed/aimag')
 
  q_2 = conjg(q)
  if(dneq(abs(q),abs(q_2)))                            call io_error(0,ext_msg='conjg/abs')
  if(q /= conjg(q_2))                                  call io_error(0,ext_msg='conjg/involution')
  if(dneq(q_2%real(),     q%real()))                   call io_error(0,ext_msg='conjg/real')
  if(any(dneq(q_2%aimag(),q%aimag()*(-1.0_preal))))    call io_error(0,ext_msg='conjg/aimag')
 
  if(abs(q) > 0.0_preal) then
    q_2 = q * q%inverse()
    if(     dneq(real(q_2), 1.0_preal,1.0e-15_preal))  call io_error(0,ext_msg='inverse/real')
    if(any(dneq0(aimag(q_2),          1.0e-15_preal))) call io_error(0,ext_msg='inverse/aimag')
 
    q_2 = q/abs(q)
    q_2 = conjg(q_2) - inverse(q_2)
    if(any(dneq0(q_2%asArray(),1.0e-15_preal)))        call io_error(0,ext_msg='inverse/conjg')
  endif
  if(dneq(dot_product(qu,qu),dot_product(q,q)))        call io_error(0,ext_msg='dot_product')
 
 
 
 
 
 
 
 
end subroutine unittest
 
 
end module quaternions