Quaternion.h

//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2003-2009 by Willem van Straten
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/
 
// epsic/src/util/Quaternion.h
 
#ifndef __Quaternion_H
#define __Quaternion_H
 
#include "Traits.h"
#include "complex_promote.h"
#include "Vector.h"


enum QBasis { Hermitian, Unitary };

template<typename T, QBasis B = Unitary> 
class Quaternion {
  
public:
  T s0,s1,s2,s3;

  Quaternion (const T& a=0.0, const T& b=0.0, const T& c=0.0, const T& d=0.0)
    : s0(a), s1(b), s2(c), s3(d) { }

  template<typename U> Quaternion (T s, const Vector<3, U>& v)
    : s0(s), s1(v[0]), s2(v[1]), s3(v[2]) { }

  template<typename U> Quaternion (const Quaternion<U, B>& s)
    : s0(s.s0), s1(s.s1), s2(s.s2), s3(s.s3) { }

  template<typename U>
  const Quaternion& operator = (const Quaternion<U,B>& s)
    { s0=T(s.s0); s1=T(s.s1); s2=T(s.s2); s3=T(s.s3); return *this; }

  const Quaternion& operator += (const T& s)
    { s0+=s; return *this; }

  const Quaternion& operator -= (const T& s)
    { s0-=s; return *this; }

  const Quaternion& operator *= (const T& a)
    { s0*=a; s1*=a; s2*=a; s3*=a; return *this; }

  const Quaternion& operator /= (const T& a)
    { T d(1.0); d/=a; s0*=d; s1*=d; s2*=d; s3*=d; return *this; }

  template<class U>
  const Quaternion& operator += (const Quaternion<U,B>& s)
    { s0+=s.s0; s1+=s.s1; s2+=s.s2; s3+=s.s3; return *this; }

  template<class U>
  const Quaternion& operator -= (const Quaternion<U,B>& s)
    { s0-=s.s0; s1-=s.s1; s2-=s.s2; s3-=s.s3; return *this; }

  template<class U>
  const Quaternion& operator *= (const Quaternion<U,B>& s)
    { *this = *this * s; return *this; }

  bool operator == (const Quaternion& b) const
    { return  s0==b.s0 && s1==b.s1 && s2==b.s2 && s3==b.s3; }

  bool operator == (T scalar) const
    { return  s0==scalar && s1==0 && s2==0 && s3==0; }

  bool operator != (const Quaternion& b) const
    { return ! operator==(b); }

  const friend Quaternion operator - (Quaternion s)
    { s.s0=-s.s0; s.s1=-s.s1; s.s2=-s.s2; s.s3=-s.s3; return s; }

  T& operator [] (int n)
    { T* val = &s0; return val[n]; }
  
  T operator [] (int n) const
    { return *(&s0+n); }

  T get_scalar () const { return s0; }

  void set_scalar (T s) { s0 = s; }

  Vector<3,T> get_vector () const 
    { Vector<3,T> ret; ret[0]=s1; ret[1]=s2; ret[2]=s3; return ret; }

  template<typename U>
  void set_vector (const Vector<3,U>& v) { s1=v[0]; s2=v[1]; s3=v[2]; }

  static const Quaternion& identity();

  unsigned size () const { return 4; }
 
};

template<typename T, QBasis B, typename U>
const Quaternion<typename PromoteTraits<T,U>::promote_type,B>
operator + (const Quaternion<T,B>& a, const Quaternion<U,B>& b)
{
  Quaternion<typename PromoteTraits<T,U>::promote_type,B> ret(a);
  ret+=b;
  return ret;
}

template<typename T, QBasis B, typename U>
const Quaternion<typename PromoteTraits<T,U>::promote_type,B>
operator - (const Quaternion<T,B>& a, const Quaternion<U,B>& b)
{
  Quaternion<typename PromoteTraits<T,U>::promote_type,B> ret(a);
  ret-=b;
  return ret;
}


template<typename T, QBasis B, typename U>
const Quaternion<T,B>
operator * (const Quaternion<T,B>& a, const U& c)
{ 
  Quaternion<typename PromoteTraits<T,U>::promote_type,B> ret(a);
  ret*=c;
  return ret;
}


template<typename T, QBasis B, typename U>
const Quaternion<T,B>
operator * (const U& c, const Quaternion<T,B>& a)
{
  Quaternion<typename PromoteTraits<T,U>::promote_type,B> ret(a);
  ret*=c;
  return ret;
}

template<typename T, QBasis B, typename U>
const Quaternion<typename PromoteTraits<T,U>::promote_type,B>
operator / (const Quaternion<T,B>& a, const U& c)
{
  Quaternion<typename PromoteTraits<T,U>::promote_type,B> ret(a);
  ret/=c;
  return ret;
}

template<typename T, QBasis B>
const Quaternion<T,B>& Quaternion<T,B>::identity ()
{
  static Quaternion<T,B> I (1,0,0,0);
  return I;
}

template<typename T, QBasis B> struct DatumTraits< Quaternion<T,B> >
{
  ElementTraits<T> element_traits;
  static inline unsigned ndim () { return 4; }
  static inline T& element (Quaternion<T,B>& t, unsigned idim) 
  { return t[idim]; }
  static inline const T& element (const Quaternion<T,B>& t, unsigned idim)
  { return t[idim]; }
};

template<typename T, typename U>
const Quaternion<std::complex<typename PromoteTraits<T,U>::promote_type>,
                 Hermitian>
operator * (const Quaternion<std::complex<T>,Hermitian>& a,
            const Quaternion<std::complex<U>,Hermitian>& b)
{
  typedef std::complex<typename PromoteTraits<T,U>::promote_type> R;
  return Quaternion<R,Hermitian>
    ( a.s0*b.s0 + a.s1*b.s1 + a.s2*b.s2 + a.s3*b.s3 ,
      a.s0*b.s1 + a.s1*b.s0 + ci(a.s2*b.s3) - ci(a.s3*b.s2) ,
      a.s0*b.s2 - ci(a.s1*b.s3) + a.s2*b.s0 + ci(a.s3*b.s1) ,
      a.s0*b.s3 + ci(a.s1*b.s2) - ci(a.s2*b.s1) + a.s3*b.s0 );
}

template<typename T, typename U>
const Quaternion<typename PromoteTraits<T,U>::promote_type, Unitary>
operator * (const Quaternion<T,Unitary>& a, const Quaternion<U,Unitary>& b)
{
  return Quaternion<typename PromoteTraits<T,U>::promote_type, Unitary>
    (a.s0*b.s0 - a.s1*b.s1 - a.s2*b.s2 - a.s3*b.s3,
     a.s0*b.s1 + a.s1*b.s0 - a.s2*b.s3 + a.s3*b.s2,
     a.s0*b.s2 + a.s1*b.s3 + a.s2*b.s0 - a.s3*b.s1,
     a.s0*b.s3 - a.s1*b.s2 + a.s2*b.s1 + a.s3*b.s0);
}
 

template<typename T, QBasis B>
Quaternion<T,B> real (const Quaternion<std::complex<T>,B>& j)
{
  return Quaternion<T,B>
    (j.s0.real(), j.s1.real(), j.s2.real(), j.s3.real());
}

template<typename T, QBasis B>
Quaternion<T,B> imag (const Quaternion<std::complex<T>,B>& j)
{
  return Quaternion<T,B>
    (j.s0.imag(), j.s1.imag(), j.s2.imag(), j.s3.imag());
}

template<typename T>
Quaternion<T,Hermitian> conj (const Quaternion<T,Hermitian>& j)
{
  return Quaternion<T,Hermitian>
    (myconj(j.s0), myconj(j.s1), myconj(j.s2), -myconj(j.s3));
}

template<typename T>
Quaternion<T,Unitary> conj (const Quaternion<T,Unitary>& j)
{ 
  return Quaternion<T,Unitary>
    (myconj(j.s0), -myconj(j.s1), -myconj(j.s2), myconj(j.s3));
}
 

template<typename T>
Quaternion<T, Hermitian> herm (const Quaternion<T,Hermitian>& j)
{
  return Quaternion<T,Hermitian>
    (myconj(j.s0), myconj(j.s1), myconj(j.s2), myconj(j.s3));
}
 

template<typename T>
Quaternion<T, Unitary> herm (const Quaternion<T,Unitary>& j)
{
  return Quaternion<T,Unitary>
    (myconj(j.s0), -myconj(j.s1), -myconj(j.s2), -myconj(j.s3));
}
 

template<typename T, QBasis B>
Quaternion<T, B> inv (const Quaternion<T,B>& j) 
{
  T d (-1.0); d/=det(j);
  return Quaternion<T,B> (-d*j.s0, d*j.s1, d*j.s2, d*j.s3);
}
 

template<typename T>
T det (const Quaternion<T,Hermitian>& j)
{ return j.s0*j.s0 - j.s1*j.s1 - j.s2*j.s2 - j.s3*j.s3; }

template<typename T>
T det (const Quaternion<T,Unitary>& j)
{ return j.s0*j.s0 + j.s1*j.s1 + j.s2*j.s2 + j.s3*j.s3; }

template<typename T, QBasis B>
T trace (const Quaternion<T,B>& j)
{ return 2.0 * j.s0; }

template<typename T, QBasis B>
T norm (const Quaternion<std::complex<T>,B>& j)
{ return 2.0 * (norm(j.s0) + norm(j.s1) + norm(j.s2) + norm(j.s3)); }

template<typename T, QBasis B>
T norm (const Quaternion<T,B>& j)
{ return 2.0 * (j.s0*j.s0 + j.s1*j.s1 + j.s2*j.s2 + j.s3*j.s3); }
 
template<typename T, QBasis B>
T fabs (const Quaternion<T,B>& j)
{ return sqrt (norm(j)); }
 
// Return the positive definite square root of a Hermitian Quaternion
template<typename T, QBasis B>
const Quaternion<T,B> sqrt (const Quaternion<T,B>& h)
{
  T root_det = sqrt( det(h) );
  T scalar = sqrt( 0.5 * (h.s0 + root_det) );
 
  if (scalar == 0.0)
    return Quaternion<T,B> (0.0);
 
  return Quaternion<T,B> (scalar, h.get_vector()/(2*scalar));
}


 
template<typename T>
const Quaternion<T,Unitary> eigen (const Quaternion<T,Hermitian>& q)
{
  T p = norm( q.get_vector() );
 
  /*
    q.s0 == 0 is a special case used by calculate_Jacobi and is required
    for the Jacobi method to work on Hermitian matrices.  Unfortunately,
    at the moment I can't remember how I derived calculate_Jacobi for
    Hermitian matrices and I can't update it to work in the case of 
    q.s1 < 0 ... another time ...
 
    Willem van Straten - 30 November 2009
  */
  if (q.s1 < 0 && q.s0 != 0)
  {
    T m = 1.0 / sqrt( 2.0*p*(p-q.s1) );
    return Quaternion<T,Unitary> (m*q.s3, -m*q.s2, -m*(p-q.s1), 0.0);
  }
  else
  {
    T m = 1.0 / sqrt( 2.0*p*(p+q.s1) );
    return Quaternion<T,Unitary> (m*(p+q.s1), 0.0, -m*q.s3, m*q.s2);
  }
}

template<typename T>
std::ostream& operator<< (std::ostream& ostr, const Quaternion<T,Hermitian>& j)
{
  return ostr << "[h:" << j.s0 <<","<< j.s1 <<","<< j.s2 <<","<< j.s3 << "]";
}

template<typename T>
std::ostream& operator<< (std::ostream& ostr, const Quaternion<T,Unitary>& j)
{
  return ostr << "[u:" << j.s0 <<","<< j.s1 <<","<< j.s2 <<","<< j.s3 << "]";
}
 
#endif  /* not __Quaternion_H defined */