#!/usr/bin/env python3
# quaternion.py
# Martin Miller
# Created: 2015/10/16
"""A quaternion class. Quaternions are represented in the form:
    a+bi+cj+dk"""
import numpy

class Quaternion(object):
    def __init__(self,a,b=None,c=None,d=None):
        """either the tuple a,b,c,d is provided where a is the scalar part and
        b,c,d is the vector part, or the quaternion is given as a 2x2 complex
        matrix 'a'"""
        if b is None:
            qy=a
            a=qy[0,0].real
            b=qy[0,0].imag
            c=qy[0,1].real
            d=qy[0,1].imag
        self.q=numpy.matrix([a,b,c,d]).T

    def column(self):
        """Returns in column representation"""
        return self.q.reshape((4,1))

    def complexMatrix(self):
        """Returns in 2x2 complex matrix representation. Suitable for unit
        quaternions only."""
        a=self.q[0,0]
        b=self.q[1,0]
        c=self.q[2,0]
        d=self.q[3,0]
        compMat=numpy.matrix([[complex(a,b),complex(c,d)],
                              [complex(-c,d),complex(a,-b)]])
        return compMat

    def scalar(self):
        return self.q[0,0]

    def vector(self):
        return self.q[1:4,0]

    def conj(self):
        """Returns the conjugate of the quaterion. For the conjugation that
        transforms a point see Quaternion.conjugate()"""
        s=self.scalar()
        v=-self.vector()
        return self.__class__(s,*v.A1)

    def conjugate(self,x):
        """Conjugates the vector or pure quaternion x with the quaternion.
        This is the operation x'=qxq^-1. For the conjugate of a quaternion see
        Quaternion.conj()"""
        if type(x)!=type(self):
            x=self.__class__(0,*x.A1)
        y=self*x*self.inverse()
        return y

    def reflect(self,x):
        """Reflects a vector or pure quaternion with self, a pure quaternion
        representing the plane of reflection.
        Based on:
            http://www.euclideanspace.com/maths/geometry/affine/reflection/quaternion/index.htm"""
        if type(x)!=type(self):
            x=self.__class__(0,*x.A1)
        y=self*x*self
        return y


    def __str__(self):
        return self.q.__str__()

    def __mul__(self,other):
        """Performs scalar or quaternion multiplication depending on the type of
        other"""
        if type(other) is not type(self): # scalar
            scaled=self.q*other
            return self.__class__(*scaled.A1)
        a1=self.q[0,0]
        b1=self.q[1,0]
        c1=self.q[2,0]
        d1=self.q[3,0]

        a2=other.q[0,0]
        b2=other.q[1,0]
        c2=other.q[2,0]
        d2=other.q[3,0]

# From https://en.wikipedia.org/wiki/Quaternion#Hamilton_product
        a3=a1*a2-b1*b2-c1*c2-d1*d2
        b3=a1*b2+b1*a2+c1*d2-d1*c2
        c3=a1*c2-b1*d2+c1*a2+d1*b2
        d3=a1*d2+b1*c2-c1*b2+d1*a2
        return self.__class__(a3,b3,c3,d3)

    def __add__(self,other):
        """Performs quaternion addition"""
        q1=self.column()
        q2=other.column()
        qy=q1+q2
        return self.__class__(*qy.A1)

    def inverse(self):
        normSquared=(self*self.conj()).q[0,0]
        return self.conj()*(1/normSquared)


def main():
    '''
    scalar=0.800103145191266
    vector=[-0.191341716182545,0.331413574035592,0.461939766255643]
    '''
    vector=[0.000000000000000,0.000000000000000,0.382683432365090]
    scalar=0.923879532511287
    q=Quaternion(scalar,*vector)
    print(q)
    print(q.conj())
    print(q.inverse())
    print(q.column())
    print(q.complexMatrix())
    print(q*q.conj())
    print((q+q.conj())*0.5)
    x=numpy.matrix("1;0;5")
    print(q.conjugate(x).vector())
    ground=Quaternion(0,0,0,1)
    print(ground.reflect(x).vector())

if __name__=='__main__':
    main()