read("bonds.mpl"):
h1 := [0, 1, 0, 0]:
h2 := [0, 9*e, 1, -9*e]:
h3 := [0, -1/3-4*e, -2/3+4*e, 2/3+2*e]:
h4 := [0, 2/3+5*e, 1/3+4*e, 2/3-7*e]:
for i from 1 to 4 do
if Mul(h||i, Cj(h||i)) <> [1, 0, 0, 0] then
print("Error: Input quaternions are not normalized!"):
stop:
end if:
if h||i + Cj(h||i) <> [0, 0, 0, 0] then
print("Error: Input quaternions are not rotations quaternions!"):
stop:
end if:
od:
for i from 1 to 4 do
u||i := [t||i, 0, 0, 0] - h||i:
od:
eqs := RVec(Mulall(u1, u2, u3, u4)):
es := {op(eqs[2..-1]), eqs[1]*u-1}:
vars := [t1, t2, u, t3, t4]:
Gb1 := Groebner[Basis]([op(es)], tdeg(t1, t2, u, t3, t4)):
egb := select(T->not (u in indets(T)), Gb1):
egb := Groebner[Basis]({op(egb)}, plex(t4, t3, t2, t1)):
configcurve := solve(egb, [t1, t2, t3, t4]);
C := subs(t4 = -t, [rhs(configcurve[1][1]),
rhs(configcurve[1][2]),
rhs(configcurve[1][3]),
rhs(configcurve[1][4])]);
for i from 1 to 4 do
U||i := [C[i], 0, 0, 0] - h||i:
od:
factor(Mulall(U1, U2, U3, U4)[1]);
sol := [solve(%, t)];
bonds := [seq(subs(t = sol[i], C), i=1..4)];
u5 := u1: u6 := u2: u7 := u3: u8 := u4:
for j from 1 to nops(bonds) do
t1 := bonds[j][1]:
t2 := bonds[j][2]:
t3 := bonds[j][3]:
t4 := bonds[j][4]:
i := 1:
while (t||i^2 <> -1) do i := i + 1: od:
printf(" F_{%d,%d}(%a, %a, %a, %a) = %a\n\n", i, cc(i+2, 4), t1, t2, t3, t4, Mulall(seq(u||k, k=i..i+2))):
printf(" F_{%d,%d}(%a, %a, %a, %a) = %a\n\n", cc(i+2, 4), cc(i+4, 4), t1, t2, t3, t4, Mulall(seq(u||k, k=i+2..i+4))):
unassign('t1', 't2', 't3', 't4'):
od:
for j from 1 to nops(bonds) do
t1 := bonds[j][1]:
t2 := bonds[j][2]:
t3 := bonds[j][3]:
t4 := bonds[j][4]:
t0 := -t4:
printf("\n %a\n", map(simplify, subs(t = t0, C))):
for k from 1 to 4 do
for l from k+1 to 4 do
F||k||l := map(simplify, Mulall(seq([C[m], 0, 0, 0] - h||m, m=k+1..l))):
v||k||l := min(map(vanishingOrder, F||k||l, t, t0)):
d||k||l := 1/2*vanishingOrder(simplify(pNorm(F||k||l)), t, t0) - v||k||l:
printf(" d(%d,%d) = %a (v_{%d,%d}=%a)\n", k, l, d||k||l, k, l, v||k||l):
od:
od:
unassign('t0', 't1', 't2', 't3', 't4'):
od:
t0 := -I:
for i from 1 to 4 do
t||i := subs(t = t0, C[i]):
od:
F13 := map(factor, Mulall([C[2], 0, 0, 0] - h2, [C[3], 0, 0, 0] - h3));
map(simplify, subs(t = t0, F13));
factor(pNorm(F13));
v13 := min(map(vanishingOrder, F13, t, t0));
d13 := 1/2*vanishingOrder(simplify(pNorm(F13)), t, t0);
d13 := d13 - v13;