# ======================================================================== #
#
# Copyright 2012, Gabor Hegedus, Josef Schicho, and Hans-Peter Schroecker
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# ======================================================================== #

read("bonds.mpl"):

# define planar four-bar linkage
# ==============================
h1 := [0, e,   0, 1]:
h2 := [0, 0,   e, 1]:
h3 := [0, 0,   0, 1]:
h4 := [0, e, 2*e, 1]:

# sanity checks of input data
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:

# Ideal of configuration curve
# ============================

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 := allvalues(solve(egb, [t1, t2, t3, t4])):
simplify(configcurve[1]);
simplify(configcurve[2]);

# parametrized equation of (parts of) configuration curve
C1 := map(simplify, subs(t4 = t, [rhs(configcurve[1][1][1]),
                                  rhs(configcurve[1][1][2]),
                                  rhs(configcurve[1][1][3]),
                                  rhs(configcurve[1][1][4])]));
C2 := map(simplify, subs(t4 = t, [rhs(configcurve[2][1][1]),
                                  rhs(configcurve[2][1][2]),
                                  rhs(configcurve[2][1][3]),
                                  rhs(configcurve[2][1][4])]));
vars:=[t1, t2, t3, t4]:
eqss := [op(egb), eqs[1]]:
bonds := solve(eqss, vars);

# The next step looses conjugate complex bonds. This is no problem as
# the bond structure of conjugate complex bonds is equal.
bonds := simplify(convert(bonds, radical)):

for b in bonds do
  print(b):
od:

# Verify vanishing of certain coupler polynomials:
u5 := u1: u6 := u2: u7 := u3: u8 := u4:

for j from 1 to nops(bonds) do
  assign(bonds[j]):
  printf("	beta = (%a, %a, %a, %a)\n\n", t1, t2, t3, t4):
  for i from 1 to 4 do
    for k from i+1 to i+3 do
      F||i||k := Mulall(seq(u||l, l=i+1..k)):
      if F||i||k = [0, 0, 0, 0] then
        printf("	F_{%d,%d}(beta) = %a\n\n", i, cc(k, 4), F||i||k):
      end if:
    od:
  od:
  unassign('t1', 't2', 't3', 't4'):
od:

C1 := subs(t4 = t, [rhs(configcurve[1][1][1]),
                    rhs(configcurve[1][1][2]),
                    rhs(configcurve[1][1][3]),
                    rhs(configcurve[1][1][4])]):
C2 := subs(t4 = t, [rhs(configcurve[2][1][1]),
                    rhs(configcurve[2][1][2]),
                    rhs(configcurve[2][1][3]),
                    rhs(configcurve[2][1][4])]):
C1 := map(simplify, C1);
C2 := map(simplify, C2);

# Compute all local distances
# ===========================
for j from 1 to nops(bonds) do
  assign(bonds[j]):
  t0 := t4:
  C := subs(t = t0, C1):
  if (map(simplify, C - [t1, t2, t3, t4]) <> [0, 0, 0, 0]) then
    C := C2:
  else
    C := C1:
  end if:
  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 := 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(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:

# Details of computation of d14 and bond (-I, I, -I, I)
# =====================================================
t0 := I:
tmp := simplify(subs(t = t0, C1)):
if (tmp[1]^2 + 1 = 0) and (tmp[2]^2 + 1 = 0) and (tmp[2]^2 + 1 = 0) and (tmp[3]^2 + 1 = 0) then
  C := C1:
else
  C := C2:
end if:
C0 := simplify(subs(t = t0, C));
for i from 1 to 4 do
  t||i := subs(t = t0, C0[i]):
od:

F14 := map(factor, Mulall([C[2], 0, 0, 0] - h2, [C[3], 0, 0, 0] - h3, [C[4], 0, 0, 0] - h4)):
v14 := min(map(vanishingOrder, F14, t, t0));
d14 := 1/2*vanishingOrder(simplify(pNorm(F14)), t, t0);
d14 := d14 - v14;