def rot(u, a):
    a *= pi / 180
    return cos(a) * matrix.identity(len(u)) + sin(a) * u.cross_product_matrix() + (1 - cos(a)) * u.tensor_product(u)



a = vector([1, 0, 0])
b = vector([0, 1, 0])
c = vector([0, 0, 1])

d = vector([1, 1, 1])

Rc = rot(c, 45)
R = Rc

a = (Rc * a).normalized()
b = (Rc * b).normalized()

Ra = rot(a, 45)
R = Ra * R

b = (Ra * b).normalized()
c = (Ra * c).normalized()

Rb = rot(b, 45)
R = Rb * R

c = (Rb * c).normalized()
a = (Rb * a).normalized()


d = R * d


plt = plot(a, color="red")+plot(b, color="blue")+plot(c, color="yellow") + plot(d, color="black")


########## two rots #######
x = vector([1, 0, 0])
y = vector([0, 1, 0])
z = vector([0, 0, 1])

d = vector([1, 1, 1])

d = rot(x, 180 / pi * asin(1/sqrt(3))) * rot(y, -45) * d

a = rot(x, 180 / pi * asin(1/sqrt(3))) * rot(y, -45) * x
b = rot(x, 180 / pi * asin(1/sqrt(3))) * rot(y, -45) * y
c = rot(x, 180 / pi * asin(1/sqrt(3))) * rot(y, -45) * z



plt2 = plot(a, color="red")+plot(b, color="blue")+plot(c, color="yellow") + plot(d, color="black")


print pi.n()/2
print 180 / pi.n() * acos(d.dot_product(a)).n()
print 180 / pi.n() * acos(d.dot_product(b)).n()
print 180 / pi.n() * acos(d.dot_product(c)).n()
print a.N()
print d.N()
show(plt)
show(plt2)


show(cube(color=['red', 'blue', 'green'], frame_thickness=2,
         frame_color='brown', opacity=0.8), frame=True)