rupert-cube.sage.py

rupert-cube.sage.py

import re

def mat_of_quat(quat):
    [r, a, b, c] = quat
    H.<i,j,k> = QuaternionAlgebra(SR,-1,-1)
    N.<x,y,z> = QQ[]
    q = r + a * i + b * j + c * k
    qs = r - a * i - b * j - c * k
    m = q * (x * i + y * j + z * k) * qs
    return Matrix([[m[i+1].coefficient(v) for v in [x,y,z]] for i in range(3)])

cube_verts = [
 [1,1,-1],
 [-1,1,-1],
 [-1,1,1],
 [-1,-1,1],
 [1,-1,1],
 [1,-1,-1],
]
# print( [mat_of_quat([1,1,1,2]) * Matrix(v).transpose() for v in cube_verts])

# mkHullConstraint(u, v, w)
#
# Takes a triple of vertices u, v, w, and returns
# a pair of expressions e₁, e₂ such that
#    e₁ > 0 ∧ e₂ > 0
# is equivalent to say that, when
#        a rotation (1 + a î + b ĵ + c k̂) = M₁/n₁ is applied to u-->v and
#        a rotation (1 + a'î + b'ĵ + c'k̂) = M₂/n₂ is applied to w
#        we project from ℝ³ → ℝ² by dropping the z coordinate
#        we think of u-->v being on a counterclockwise traversal in an
#          x-right and y-up orientation of the plane
# then we have that w lies on the "inside" of the edge u-->v
#
# Those constraint are equivalent to saying:
#  n₂(M₁u × M₁v) + n₁(M₁v × M₂w) + n₁(M₂w × M₁u) > 0
def mkHullConstraint(u, v, w):
    N.<a1,b1,c1,a2,b2,c2> = QQ[]
    M1 = mat_of_quat([1,a1,b1,c1])
    M2 = mat_of_quat([1,a2,b2,c2])
    nn1 = 1 + a1^2 + b1^2 + c1^2
    nn2 = 1 + a2^2 + b2^2 + c2^2
    uu = vector(u)
    vv = vector(v)
    ww = vector(w)
    p1 = nn2 * (M1 * uu).cross_product(vector(M1 * vv))
    p2 = nn1 * (M1 * vv).cross_product(vector(M2 * ww))
    p3 = nn1 * (M2 * ww).cross_product(vector(M1 * uu))
    return [p1[0] + p2[0] + p3[0], p1[1] + p2[1] + p3[1]]

constraints = []

def mkRotConstraint(a,b,c):
    return [
         a * c - b,
         b * c + a,
         1 + c^2 - a^2 - b^2,
         ]

# expects a symbolic expression over the variables a1,b1,c1,a2,b2,c2
def addConstraint(c, x):
    s = str(x)
    s = re.sub("\\*", " ", s)
    s = s + " > 0"
    constraints.append(s)

for i in range(6):
    for j in range(6):
        k = mkHullConstraint(cube_verts[i], cube_verts[(i+1)%6], cube_verts[j])
        addConstraint(constraints, k[0].expand())
        addConstraint(constraints, k[1].expand())

N.<a1,b1,c1,a2,b2,c2> = QQ[]
for k in mkRotConstraint(a1,b1,c1):
    addConstraint(constraints, k)
for k in mkRotConstraint(a2,b2,c2):
    addConstraint(constraints, k)

body = " /\\\n".join(constraints)
script = f"""[ Rupert ]
(a1,b1,c1,a2,b2,c2)
0
(E a1)(E b1)(E c1) (E a2)(E b2)(E c2) [
{body}
].
finish
"""
print(f"cat <<\"EOF\" >input.txt\n{script}\nEOF\n")