Symbolic dot and cross products with implicit rotation matrices

[Luke]

I'm trying to figure out how I could implement in Sympy a way to take
dot and cross products of unit vectors in different reference frames.
For example, suppose I have a reference frame 'A', along with 3
dextral orthonormal vectors: a1, a2, a3. I think I can figure out how
to code up something so that dot and cross products in the same
reference frame would give the following results:
In[1]: dot(a1,a2)
Out[1]: 0
In[2]: dot(a1,a1)
Out[2]: 1
In[3]: cross(a1,a2)
Out[3]: a3

I haven't done this yet, but I think it should be pretty
straightforward. What I would like to extend this to is to be able to
dot and cross multiply unit vectors from different reference frames.
For example, suppose I introduce a second reference frame 'B', along
with 3 dextral orthonormal vectors: b1, b2, b3, and I orient B by
aligning B with A then performing a simple rotation of B about the
a3=b3 axis by an angle q1. Then:
b1 = cos(q1)*a1 + sin(q1)*a2
b2 = -sin(q1)*a1 + cos(q1)*a2
b3 = a3

I would like to then be able to do:
In[4]: dot(a1,b1)
Out[4]: cos(q1)

etc...

I guess what I'm not sure about is how to structure all the objects so
that if a rotation matrix hasn't been defined between two reference
frames, an exception is raised and you can't take the dot or cross
product. It would be ideal to be able to handle long chains of simple
rotations so that every possible rotation matrix wouldn't need to be
defined explicitly, i.e., if I specify the orientation of B relative
to A (as above), and then the orientation of another reference frame C
relative to B, and then try to take a dot product between a1 and c1,
the two rotation matrices (C to B, B to A) get multplied and then used
to compute the dot product.

I'm familiar with classes and python fairly well, and I know the
kinematics well, but I'm by no means an experience object oriented
programmer, so I'm not sure about how the best way to structure things
would be.

Any suggestions on how I might start on something like this?

Thanks,
~Luke

[Alan Bromborsky]

Look at

http://wiki.sympy.org/wiki/Geometric_Algebra_Module

Then go to the section on rotations. For vectors the inner product is
the dot product. How do you want to parameterize each individual
rotation. Tell me that and I could show you how to do the coding.

[Alan Bromborsky]

Do you want to do this purely symbolically or are you going to put in
specific values? That is do you want sympy only to generate symbolic
formulas that you will later use in a program or do you want sympy to
also do the numerical evaluation?

[Luke]

I would like everything to be symbolic. If you have two reference
frames, say A and B, you can only take the dot product of a vector
expressed using the A coordinates with a vector expressed in the B
coordinates if you know the rotation matrix that relates the two
reference frames. For example:
Take A and B to be two reference frames that have a set of 3 dextral
orthonormal vectors ( i.e, dot(a1,a2)=dot(a2,a3)=dot(a1,a3)=0, cross
(a1,a2)=a3, cross(a2,a3)=a1, and cross(a3,a1)=a2). Define the
orientation of B relative to A by initially aligning a1 with b1, a2
with b2, a3 with b3, then performing a right handed rotation of B
relative to A by an angle q1 about an axis parallel to a3 (and b3).

Then:
b1 = cos(q1)*a1 + sin(q1)*a2
b2 = -sin(q1)*a1 + cos(q1)*a2
b3 = a3

Normally I would write a1,a2,a3,b1,b2,b3 in bold or with a carat above
them, the way you might see them in a physics class or a dynamics
class, but I think you get the idea.

Now, in order to take a dot product of two vectors that are expressed
in the same reference frame coordinates, we use the following rules:
dot(a1,a2) = 0
dot(a1,a3) = 0
dot(a2, a3) = 0
dot(a1,a1)=1
dot(a2,a2)=1
dot(a3,a3)=1

So any vector expression that only involves terms with a1,a2,a3, we
can easily a dot() function to implement the above rules for the dot
product. But what if we need to dot multiply a1 and b1? Then we need
to either 1) express a1 in terms of B coordinate basis vectors, or 2)
express b1 in terms of A coordinate basis vectors. Here is where we
need to rotation matrix defined above, and here is how I'm not sure of
how to best implement things. I'm looking for a good way to store the
rotation matrices with the reference frame objects.... If one defines
successive rotations, say, starting with A, then orienting B relative
to A, then C relative to B, then D relative to C, we want to have the
dot() function be able to resolve the dot product of say d1 and a1,
even though the rotation matrix between D and A wasn't explicitly
defined (instead is it is implicitly defined through a series of
matrix multiplications). One could imagine a long chain of rotations,
with various branches (as is common in the case of multi-body
systems). How to best store the and resolve the various rotations is
what I'm after.

The reason I would like to do this is to be ably to use the following
syntax to define the position of one point q relative to another p:
p_q_p = q2*a1 + q3*b2

and the be able to do:
dot(p_q_p,a1)
and get:
q2 - sin(q1)*q3

(recall, b2 = -sin(q1)*a1+cos(q1)*a2)

Hopefully this clarifies what I'm trying to do.

Thanks,
~Luke

[Alan Bromborsky]

Here is the code todo the frame rotation:

#!/usrlocal/bin/python
#rotate.py

import sympy.galgebra.GAsympy as GA

import sympy,numpy,sys

GA.set_main(sys.modules[__name__])

I = 0

def rotate_frame(frame_vec,rot_vec,rot_amp):
#rot_vec**2 = 1 and rot_amp in radians
rot_versor = I*rot_vec
theta = rot_amp/2
R = sympy.cos(theta)+sympy.sin(theta)*rot_versor
Rdag = R.rev()
new_frame_vec = []
for vec in frame_vec:
new_frame_vec.append(R*vec*Rdag)
return(new_frame_vec)

if __name__ == '__main__':
metric = '1 0 0,'+\
'0 1 0,'+\
'0 0 1'

GA.MV.setup('e_x e_y e_z',metric)
I = e_x*e_y*e_z #Calculate pseudo-scalar
GA.make_symbols('w_x w_y w_z theta') #create real scalar sympy symbols
old_frame = [e_x,e_y,e_z] #original orthogonal frame
rot_axis = w_x*e_x+w_y*e_y+w_z*e_z #rotation axis unit vector
new_frame = rotate_frame(old_frame,rot_axis,theta) #rotated
orthogonal frame
print 'rotation axis =',rot_axis
GA.MV.set_str_format(2) #print one vector component per line
print 'rotated frame:'
for (old_vec,new_vec) in zip(old_frame,new_frame):
print str(old_vec)[:-1]+"' =",new_vec

Here is the program output:

rotation axis =
{w_x}e_x+{w_y}e_y+{w_z}e_z

rotated frame:
e_x' = {cos(theta/2)**2 + w_x**2*sin(theta/2)**2 -
w_y**2*sin(theta/2)**2 - w_z**2*sin(theta/2)**2}e_x
+{-2*w_z*cos(theta/2)*sin(theta/2) + 2*w_x*w_y*sin(theta/2)**2}e_y
+{2*w_y*cos(theta/2)*sin(theta/2) + 2*w_x*w_z*sin(theta/2)**2}e_z

e_y' = {2*w_z*cos(theta/2)*sin(theta/2) + 2*w_x*w_y*sin(theta/2)**2}e_x
+{cos(theta/2)**2 + w_y**2*sin(theta/2)**2 - w_x**2*sin(theta/2)**2 -
w_z**2*sin(theta/2)**2}e_y
+{-2*w_x*cos(theta/2)*sin(theta/2) + 2*w_y*w_z*sin(theta/2)**2}e_z

e_z' = {-2*w_y*cos(theta/2)*sin(theta/2) + 2*w_x*w_z*sin(theta/2)**2}e_x
+{2*w_x*cos(theta/2)*sin(theta/2) + 2*w_y*w_z*sin(theta/2)**2}e_y
+{cos(theta/2)**2 + w_z**2*sin(theta/2)**2 - w_x**2*sin(theta/2)**2 -
w_y**2*sin(theta/2)**2}e_z


We assume that w_x**2+w_y**2+w_z**2 =1 and double angle substitution
could be used to modify the output. R in the function 'rotate_frame' is
a rotor which is a special case of a spinor. I is the pseudo-scalar for
3-dimensional euclidian space. w_x, w_y, and w_z could be expressed in
sperical coordinates.