# ======================================================================== #
#
# 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 a spherical four-bar linkage
# ===================================
h1 := [0, 1, 0, 0]:
h2 := [0, 0, 1, 0]:
h3 := [0, 0, 0, 1]:
h4 := [0, 3/5, 4/5, 0]:

# 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]);

# Compute bonds
# =============
vars:=[t1, t2, t3, t4]:
eqss := [op(egb), eqs[1]]:
bonds := solve(eqss, vars):
rr := simplify(convert(bonds, radical)):
for r in rr do
  print(r):
od:
stop:

# Verify vanishing of certain coupler maps:
cc := proc(k, n)
  if (k > n) then: return cc(k-n, n): end if:
  if (k < 1) then: return cc(k+n, n): end if:
  return k:
end proc:
pNorm := proc(k)
  return Mul(k, Cj(k))[1]:
end proc:

for j from 1 to nops(rr) do
  assign(rr[j]):
  i := 1:
  while (t||i^2 <> -1) do i := i + 1: od:
  u5 := u1: u6 := u2: u7 := u3: u8 := u4:
  print(t1, t2, t3, t4):
  print(seq(k, k=i..i+2), Mulall(seq(u||k, k=i..i+2)), seq(cc(k, 4), k=i+2..i+4), Mulall(seq(u||k, k=i+2..i+4)));
  unassign('t1', 't2', 't3', 't4'):
od: