213.sagews

r"""

q_lib - library of numerical quaternion functions

Authors:
    -- Doug Sweetser (2009): Initial version.
    
Examples:
    Parts
    sage: q1 = q(1, 2, 3, 4)
    sage: q1.scalar() == q(1, 0, 0, 0)
    True
    sage: q1.vector() == q(0, 2, 3, 4)
    True
    sage: q1.norm() == q(30, 0, 0, 0)
    True
     
    Add quaternions 3 times
    sage: q1 = q(1,1,1,1); q2 = q(1,2,3,4); q3 = q(3, 5, 7, 9)
    sage: test = add_n(q2, 3)
    sage: q3 == test
    True
    
    Test the various conjugates.
    sage: q1 = q(1, 2, 3, 4); q_time_rev = q(-1, 2, 3, 4); 
    sage: q_space_rev = q(1, -2, -3, -4); q_spacetime_rev = q(-1, -2, -3, -4)
    sage: q1.conj() == q_space_rev
    True
    sage: q1.conj().negate() == q_time_rev
    True
    sage: q1.negate() == q_spacetime_rev
    True
    
    Groups
    sage: q1 = q(2, 1, 2, 0); q2 = q(2/3, 1/3, 2/3, 0)
    sage: q1.u1() == q2
    True
    
    Test the group U(1)xSU(2).
    sage: q1 = q(0, 2, 1, 2); q2 = q(-sin(3), 2 * cos(3)/3, cos(3)/3, 2 * cos(3)/3)
    sage: q1.u1xsu2() == q2
    True
"""
#*****************************************************************************
#       Copyright (C) 2009 Doug Sweetser <[email protected]>
#
#  Distributed under the terms of the GNU General Public License (GPL)
#                  http://www.gnu.org/licenses/
#*****************************************************************************

class q(list):
    def __init__(self,t,x,y,z):
        self.t = t
        self.x = x
        self.y = y
        self.z = z
        self.qn = []
        list.__init__(self,[t,x,y,z])
    def __repr__(self):
        return "Quaternion tools."
    def show(self):
        print self.t, self.x, self.y, self.z
    def normalize(self):
        abs_value = sqrt(self.t^2 + self.x^2 + self.y^2 + self.z^2)
        return q(self.t/abs_value, self.x/abs_value, self.y/abs_value, self.z/abs_value)
        
    # Parts and sizes    
    def scalar(self):
        return q(self.t, 0, 0, 0)
    def vector(self):
        return q(0, self.x, self.y, self.z)
    def norm(self):
        return q(self.t * self.t + self.x * self.x + self.y * self.y + self.z * self.z, 0, 0, 0)

    # Conjugates
    def conj(self):
        return q(self.t, -self.x, -self.y, -self.z)
    def conj_1(self):
        return q(-self.t, self.x, -self.y, -self.z)   
    def conj_2(self):
        return q(-self.t, -self.x, self.y, -self.z)
    def conj_3(self):
        return q(-self.t, -self.x, -self.y, self.z)   
    def negate(self):
        return q(-self.t, -self.x, -self.y, -self.z) 

    def add(self, q1):
        return q(self.t + q1.t, self.x + q1.x, self.y + q1.y, self.z + q1.z)
    def add_n(self, q1, n):
        self.qn[:]=[]
        for i in range (0, n):
            self.t += q1.t
            self.x += q1.x
            self.y += q1.y
            self.z += q1.z
            self.qn.append(q(self.t, self.x, self.y, self.z))
        return q(self.t, self.x, self.y, self.z)
    def scalar_div(self, d):
        return(q(self.t/d, self.x/d, self.y/d, self.z/d));
    def scalar_x(self, d):
        return(q(self.t * d, self.x * d, self.y * d, self.z * d));

    def real(self):
        t = self.add(self.conj()).scalar_div(2)
        return q(t.t, t.x, t.y, t.z)
    def complex_x(self):
         tx =self.add(self).add(self.conj()).add(self.conj_1()).scalar_div(2)
         return q(tx.t, tx.x, tx.y, tx.z)
    def complex_y(self):
         ty = self.add(self).add(self.conj()).add(self.conj_2()).scalar_div(2)
         return q(ty.t, ty.x, ty.y, ty.z)
    def complex_z(self):
         tz = self.add(self).add(self.conj()).add(self.conj_3()).scalar_div(2)
         return q(tz.t, tz.x, tz.y, tz.z)
            
    def n_random(self, n):
        self.qn[:]=[]
        for i in range(n):
            self.qn.append(q(2 * random() - 1, 2 * random() - 1, 2 * random() - 1, 2 * random() - 1))
        return self
        
    def times(self, q2):
        return q(self.t * q2.t - self.x * q2.x - self.y * q2.y - self.z * q2.z, self.t * q2.x + self.x * q2.t + self.y * q2.z - self.z * q2.y, self.t * q2.y + self.y * q2.t + self.z * q2.x - self.x * q2.z, self.t * q2.z + self.z * q2.t + self.x * q2.y - self.y * q2.x)
        
# Exponentials and trig functions
#   exp(q) = (exp(t) cos(|R|), exp(t) sin(|R|) R/|R|)
    def exp(self):
        abs_vector = sqrt(self.x^2 + self.y^2 + self.z^2)
        if abs_vector == 0:
            return q(exp(self.t), 0, 0, 0)
        k = exp(self.t) * sin(abs_vector) / abs_vector
        return q(exp(self.t) * cos(abs_vector), k * self.x, k * self.y, k * self.z)

#   sin(q) = (sin(t) cosh(|R|), cos(t) sinh(|R|) R/|R|)
    def sin(self):
        abs_vector = sqrt(self.x^2 + self.y^2 + self.z^2)
        if abs_vector == 0:
            return q(sin(self.t), 0, 0, 0)
        k = cos(self.t) * sinh(abs_vector) / abs_vector
        return q(sin(self.t) * cosh(abs_vector), k * self.x, k * self.y, k * self.z)
        
#   cos(q) = (cos(t) cosh(|R|), sin(t) sinh(|R|) R/|R|)
    def cos(self):
        abs_vector = sqrt(self.x^2 + self.y^2 + self.z^2)
        if abs_vector == 0:
            return q(cos(self.t), 0, 0, 0)
        k = sin(self.t) * sinh(abs_vector) / abs_vector
        return q(cos(self.t) * cosh(abs_vector), k * self.x, k * self.y, k * self.z)
        
# Groups

    def u1(self):
        return self.normalize()
    def su2(self):
        return self.vector().exp()
    def u1xsu2(self):
        return self.su2().times(self.u1())

#auto

r"""Check all testing code"""

\text{Check all testing code}

sage: q1 = q(1, 2, 3, 4)
sage: q1.scalar() == q(1, 0, 0, 0)

{\rm True}

sage: q1.vector() == q(0, 2, 3, 4)

{\rm True}

sage: q1.norm() == q(30, 0, 0, 0)

{\rm True}

sage: q1.real()==q(1, 0, 0, 0)

{\rm True}

sage: q1.complex_x()==q(1, 2, 0, 0)

{\rm True}

sage: q1.complex_y()==q(1, 0, 3, 0)

{\rm True}

sage: q1.complex_z()==q(1, 0, 0, 4)

{\rm True}

sage: q1 = q(1, 2, 3, 4); q_time_rev = q(-1, 2, 3, 4); 
sage: q_space_rev = q(1, -2, -3, -4); q_spacetime_rev = q(-1, -2, -3, -4)
sage: q1.conj() == q_space_rev

{\rm True}

sage: q1.conj().negate() == q_time_rev

{\rm True}

sage: q1.negate() == q_spacetime_rev

{\rm True}

q1 = q(1,1,1,1); q2 = q(1,2,3,4); q3 = q(3, 5, 7, 9)
test = q1.add_n(q2, 2)
test == q3

{\rm True}

sage: q1 = q(2, 1, 2, 0); q2 = q(-79/81, -8/81, -16/81, 0)
sage: q1.group_u1(2) == q2

{\rm True}

q1 = q(2, 1, 2, 0).u1().show()

2/3 1/3 2/3 0