1 files changed, 533 insertions, 0 deletions
diff --git a/Build/source/utils/asymptote/knot.h b/Build/source/utils/asymptote/knot.h
new file mode 100644
index 00000000000..e087f271e61
--- /dev/null
+++ b/Build/source/utils/asymptote/knot.h
@@ -0,0 +1,533 @@
+/*****
+ * knot.h
+ * Andy Hammerlindl 200/02/10
+ *
+ * Describes a knot, a point and its neighbouring specifiers, used as an
+ * intermediate structure in solving paths.
+ *****/
+
+#ifndef KNOT_H
+#define KNOT_H
+
+#include <iostream>
+#include <vector>
+#include <algorithm>
+
+#include "mod.h"
+#include "pair.h"
+#include "path.h"
+
+namespace camp {
+
+using mem::vector;
+
+// The choice of branch cuts of atan2 disguishes between y=+0.0 and y=-0.0 in
+// the case where x<0.  This can lead to strange looking paths being
+// calculated from guides of the form a..b..cycle.  To avoid these degenerate
+// cases, the niceAngle routine moves the branch cut so that the sign of a
+// zero won't matter.
+double niceAngle(pair z);
+
+// A cyclic vector: ie. a vector where the index is taken mod the size of the
+// vector.
+template <typename T>
+class cvector : public vector<T> {
+public:
+  cvector() {}
+  cvector(size_t n) : vector<T>(n) {}
+  cvector(size_t n, const T& t) : vector<T>(n,t) {}
+  cvector(const vector<T>& v) : vector<T>(v) {}
+
+  T& operator[](Int j) {
+    return vector<T>::operator[](imod(j,(Int) this->size()));
+  }
+  const T& operator[](Int j) const {
+    return vector<T>::operator[](imod(j,(Int) this->size()));
+  }
+};
+
+// Forward declaration.
+class knotlist;
+
+/* A linear equation (one of a set of equations to solve for direction through
+   knots in a path).  The i-th equation is:
+
+   pre*theta[i-1] + piv*theta[i] + post*theta[i+1] = aug
+
+   where indices are taken mod n.
+*/
+struct eqn {
+  double pre,piv,post,aug;
+  eqn(double pre, double piv, double post, double aug)
+    : pre(pre), piv(piv), post(post), aug(aug) {}
+
+  friend ostream& operator<< (ostream& out, const eqn& e)
+  {
+    return out << e.pre << " * pre + "
+               << e.piv << " * piv + "
+               << e.post << " * post = "
+               << e.aug;
+  }
+};
+
+
+// A direction specifier, telling how the path behaves coming in or out of a
+// point.  The base class represents the "open" specifier.
+class spec : public gc {
+public:
+  virtual ~spec() {}
+
+  // If the knot is open, it gives no restriction on the behavior of the
+  // path.
+  virtual bool open() { return true; }
+  virtual bool controlled() { return false; }
+  virtual pair control() {return pair(0.0,0.0);}
+  virtual double curl() { return -1.0; }
+  virtual pair dir() { return pair(0.0,0.0); }
+
+  // When a knot has a restriction on one side but is open on the other, the
+  // restriction implies a restriction on the other side.  This is the partner
+  // restriction defined here, where the pair argument is for the location of
+  // the knot.
+  virtual spec *outPartner(pair) { return this; }
+  virtual spec *inPartner(pair) { return this; }
+
+  virtual void print(ostream&) const {}
+};
+
+inline ostream& operator<< (ostream& out, spec& s)
+{
+  s.print(out);
+  return out;
+}
+
+
+// Specifier used at an endpoint.
+class endSpec : public spec {
+public:
+  bool open() { return false; }
+
+  // Returns an equation used to solve for the thetas along the knot.  These are
+  // called by eqnprop in the non-cyclic case for the first and last equations.
+  virtual eqn eqnOut(Int j, knotlist& l,
+                     cvector<double>& d, cvector<double>& psi) = 0;
+  virtual eqn eqnIn (Int j, knotlist& l,
+                     cvector<double>& d, cvector<double>& psi) = 0;
+};
+
+// A specifier with a given direction (in radians).
+class dirSpec : public endSpec {
+  double given;
+public:
+  // Direction should be given in the range [-PI,PI]
+  dirSpec(double given)
+    : given(given) {}
+  dirSpec(pair z)
+    : given(niceAngle(z)) {}
+
+  pair dir() { return expi(given); }
+
+  eqn eqnOut(Int j, knotlist& l, cvector<double>& d, cvector<double>& psi);
+  eqn eqnIn (Int j, knotlist& l, cvector<double>& d, cvector<double>& psi);
+
+  void print(ostream& out) const {
+    out << "{dir(" << degrees(given) << ")}";
+  }
+};
+
+// A curl specifier.  The curvature at the end knot should be gamma times the
+// curvature at the neighbouring knot.
+class curlSpec : public endSpec {
+  double gamma;
+public:
+  // Gamma should be non-negative.
+  curlSpec(double gamma=1.0)
+    : gamma(gamma) {
+    if(gamma < 0)
+      reportError("curl cannot be less than 0");
+  }
+
+  double curl() { return gamma; }
+
+  eqn eqnOut(Int j, knotlist& l, cvector<double>& d, cvector<double>& psi);
+  eqn eqnIn (Int j, knotlist& l, cvector<double>& d, cvector<double>& psi);
+
+  void print(ostream& out) const {
+    out << "{curl " << gamma << "}";
+  }
+};
+
+
+// A specifier with a control point.  All information for this portion of the
+// curve has been determined.
+class controlSpec : public spec {
+public:
+  pair cz;
+  bool straight;
+
+  controlSpec(pair cz, bool straight=false)
+    : cz(cz), straight(straight) {}
+
+  bool open() { return false; }
+  bool controlled() { return true; }
+  pair control() { return cz; }
+
+  // The partner spec will be a dirSpec in the same direction the specifier
+  // takes the path, unless the velocity is zero, then it uses a curl
+  // specifier.
+  spec *outPartner(pair);
+  spec *inPartner(pair);
+
+  void print(ostream& out) const {
+    // NOTE: This format cannot be read back in.
+    out << "{control " << cz << "}";
+  }
+};
+
+
+// The tension information for one side of a knot.
+struct tension {
+  double val;
+  bool atleast;
+
+  tension(double val=1.0, bool atleast=false)
+    : val(val), atleast(atleast) {
+    if(val < 0.75)
+      reportError("tension cannot be less than 3/4");
+  }
+};
+
+inline
+ostream& operator<<(ostream& out, tension t)
+{
+  return out << "tension" << (t.atleast ? " atleast " : " ") << t.val;
+}
+
+// A knot, a point with specifiers to double the path coming in and going out
+// of the knot.
+struct knot {
+  pair z;
+  spec *in;
+  spec *out;
+  tension tin, tout;
+
+  knot() {}
+  knot(pair z, spec *in, spec *out,
+       tension tin=tension(), tension tout=tension())
+    : z(z), in(in), out(out), tin(tin), tout(tout) {}
+
+  double alpha() { return 1.0/tout.val; }
+  double beta() { return 1.0/tin.val; }
+};
+
+ostream& operator<<(ostream& out, const knot& k);
+
+// Abstract base class for a section of a guide.
+class knotlist {
+public:
+  virtual ~knotlist() {}
+
+  virtual Int length() = 0;
+  virtual bool cyclic() = 0;
+
+  // Returns the number of knots.
+  Int size() {
+    return cyclic() ? length() : length() + 1;
+  }
+  bool empty() {
+    return size()==0;
+  }
+
+  virtual knot& cell(Int) = 0;
+  virtual knot& operator[] (Int i) {
+#if 0
+    assert(cyclic() || (0 <= i && i <= length()));  // Bounds check.
+#endif
+    return cell(i);
+  }
+  knot& front() {
+    return (*this)[0];
+  }
+  knot& back() {
+    return (*this)[length()];
+  }
+};
+
+
+// Defines a knotlist as a piece of another knotlist.
+class subknotlist : public knotlist {
+  knotlist& l;
+  Int a,b;
+public:
+  subknotlist(knotlist& l, Int a, Int b)
+    : l(l), a(a), b(b) {}
+
+  Int length() { return b-a; }
+  bool cyclic() { return false; }
+  knot& cell(Int i) { return l[a+i]; }
+};
+
+struct simpleknotlist : public knotlist {
+  cvector<knot> nodes;
+  bool cycles;
+
+  simpleknotlist(cvector<knot> nodes, bool cycles=false)
+    : nodes(nodes), cycles(cycles) {}
+
+  Int length() { return cycles ? (Int) nodes.size() : (Int) nodes.size() - 1; }
+  bool cyclic() { return cycles; }
+  knot& cell(Int j) { return nodes[j]; }
+};
+
+// A protopath is a path being made.
+struct protopath {
+  bool cycles;
+  Int n;
+  mem::vector<solvedKnot> nodes;
+
+  protopath(Int n, bool cycles)
+    : cycles(cycles), n(n), nodes(n) {}
+
+  solvedKnot& operator[](Int j) {
+    return nodes[imod(j,n)];
+  }
+
+  bool& straight(Int j) {
+    return (*this)[j].straight;
+  }
+  pair& pre(Int j) {
+    return (*this)[j].pre;
+  }
+  pair& point(Int j) {
+    return (*this)[j].point;
+  }
+  pair& post(Int j) {
+    return (*this)[j].post;
+  }
+
+  void controlEnds() {
+    if (!cycles) {
+      solvedKnot& start=(*this)[0];
+      solvedKnot& end=(*this)[n-1];
+
+      start.pre=start.point;
+      end.post=end.point;
+    }
+  }
+  // Once all the controls are set, return the final (constant) path.
+  path fix() {
+    return path(nodes,n,cycles);
+  }
+};
+
+
+// Represents properties that can be computed along a knotlist.
+// Examples include distances (d), turning angles (psi), and the linear
+// equations used to solve for the thetas.
+template <typename T>
+class knotprop {
+protected:
+  knotlist& l;
+
+  // Calculate the property for the usual case in the iteration (and for a
+  // cyclic knot, the only case), at the index given.
+  virtual T mid(Int) = 0;
+
+  // The special cases, these default to the usual case: mid.
+  virtual T solo(Int j)     // Calculates the property for a list of length 0.
+  {
+    return mid(j);
+  }
+  virtual T start(Int j)    // Calculates it at the start of the list.
+  {
+    return mid(j);
+  }
+  virtual T end(Int j)      // Calculate it at the end.
+  {
+    return mid(j);
+  }
+
+  virtual cvector<T> linearCompute()
+  {
+    Int n=l.length();
+    cvector<T> v;
+    if (n==0)
+      v.push_back(solo(0));
+    else {
+      v.push_back(start(0));
+      for (Int j=1; j<n; ++j)
+        v.push_back(mid(j));
+      v.push_back(end(n));
+    }
+    return v;
+  }
+
+  virtual cvector<T> cyclicCompute()
+  {
+    Int n=l.length();
+    cvector<T> v;
+    for (Int j=0; j<n; ++j)
+      v.push_back(mid(j));
+    return v;
+  }
+
+  virtual cvector<T> linearBackCompute()
+  {
+    Int n=l.length();
+    cvector<T> v;
+    if (n==0)
+      v.push_back(solo(0));
+    else {
+      v.push_back(end(n));
+      for (Int j=1; j<n; ++j)
+        v.push_back(mid(n-j));
+      v.push_back(start(0));
+    }
+    return v;
+  }
+
+  virtual cvector<T> cyclicBackCompute()
+  {
+    Int n=l.length();
+    cvector<T> v;
+    for (Int j=1; j<=n; ++j)
+      v.push_back(mid(n-j));
+    return v;
+  }
+
+public:
+  virtual ~knotprop() {}
+
+  virtual cvector<T> compute() {
+    return l.cyclic() ? cyclicCompute() : linearCompute();
+  }
+
+  // Compute the values in the opposite order.  This is needed for instance if
+  // the i-th calculation needed a result computed in the i+1-th, such as in the
+  // back substitution for solving thetas.
+  virtual cvector<T> backCompute() {
+    cvector<T> v=l.cyclic() ? cyclicBackCompute() : linearBackCompute();
+
+    // Even though they are computed in the backwards order, return them in the
+    // standard order.
+    reverse(v.begin(),v.end());
+    return v;
+  }
+
+  knotprop(knotlist& l)
+    : l(l) {}
+};
+
+// A knot transforms, it takes in one knotlist and transforms it knot for knot
+// into a new one.
+class knottrans : public knotprop<knot> {
+protected:
+  virtual knot mid(Int j) {
+    /* By default, just copy the knot. */
+    return l[j];
+  }
+
+public:
+  virtual ~knottrans() {}
+
+  knottrans(knotlist& l)
+    : knotprop<knot>(l) {}
+
+  virtual simpleknotlist trans() {
+    return simpleknotlist(compute(),l.cyclic());
+  }
+};
+
+// Like a knotprop, but it doesn't compute a vector of values for the knot.  It
+// iterates over the knotlist calling method for side-effect.  For instance,
+// this is used to plug control points into protopaths.
+class knoteffect {
+protected:
+  knotlist& l;
+
+  virtual void mid(Int) = 0;
+
+  // The special cases, these default to the usual case: mid.
+  virtual void solo(Int j) {
+    mid(j);
+  }
+  virtual void start(Int j) {
+    mid(j);
+  }
+  virtual void end(Int j) {
+    mid(j);
+  }
+
+  virtual void linearExec()
+  {
+    Int n=l.length();
+    if (n==0)
+      solo(0);
+    else {
+      start(0);
+      for (Int j=1; j<n; ++j)
+        mid(j);
+      end(n);
+    }
+  }
+
+  virtual void cyclicExec()
+  {
+    Int n=l.length();
+    for (Int j=0; j<n; ++j)
+      mid(j);
+  }
+
+  virtual void linearBackExec()
+  {
+    Int n=l.length();
+    if (n==0)
+      solo(0);
+    else {
+      end(n);
+      for (Int j=1; j<n; ++j)
+        mid(n-j);
+      start(0);
+    }
+  }
+
+  virtual void cyclicBackExec()
+  {
+    Int n=l.length();
+    for (Int j=1; j<=n; ++j)
+      mid(n-j);
+  }
+
+public:
+  virtual ~knoteffect() {}
+
+  virtual void exec() {
+    if (l.cyclic())
+      cyclicExec();
+    else
+      linearExec();
+  }
+
+  virtual void backCompute() {
+    if (l.cyclic())
+      cyclicBackExec();
+    else
+      linearBackExec();
+  }
+
+  knoteffect(knotlist& l)
+    : l(l) {}
+};
+
+path solve(knotlist& l);
+
+path solveSimple(cvector<pair>& z);
+
+double velocity(double theta, double phi, tension t);
+
+} // namespace camp
+
+GC_DECLARE_PTRFREE(camp::eqn);
+GC_DECLARE_PTRFREE(camp::tension);
+
+#endif // KNOT_H