diff options
author | Karl Berry <karl@freefriends.org> | 2009-05-16 00:19:13 +0000 |
---|---|---|
committer | Karl Berry <karl@freefriends.org> | 2009-05-16 00:19:13 +0000 |
commit | bab45528d65eaafe68a705dbb2a57075c7b7cbd8 (patch) | |
tree | 10b4ae2b5195c8dede153ab89359ec00f55f325f /Build/source/utils/asymptote/knot.cc | |
parent | 8643d90372e9c31e0f461c15c596b60a545bd7d3 (diff) |
asymptote 1.72 sources (not integrated into build yet)
git-svn-id: svn://tug.org/texlive/trunk@13110 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/utils/asymptote/knot.cc')
-rw-r--r-- | Build/source/utils/asymptote/knot.cc | 840 |
1 files changed, 840 insertions, 0 deletions
diff --git a/Build/source/utils/asymptote/knot.cc b/Build/source/utils/asymptote/knot.cc new file mode 100644 index 00000000000..0180c4c6b2d --- /dev/null +++ b/Build/source/utils/asymptote/knot.cc @@ -0,0 +1,840 @@ +/***** + * knot.cc + * Andy Hammerlindl 2005/02/10 + * + * Describes a knot, a point and its neighbouring specifiers, used as an + * intermediate structure in solving paths. + *****/ + +#include "knot.h" +#include "util.h" + +#include "angle.h" +#include "settings.h" + +namespace camp { + +/***** Debugging *****/ +//bool tracing_solving=false; + +template <typename T> +ostream& info(ostream& o, string name, cvector<T>& v) +{ + if (settings::verbose > 3) { + o << name << ":\n\n"; + + for(Int i=0; i < (Int) v.size(); ++i) + o << v[i] << endl; + + o << endl; + } + return o; +} + +ostream& info(ostream& o, string name, knotlist& l) +{ + if (settings::verbose > 3) { + o << name << ":\n\n"; + + for(Int i=0; i < (Int) l.size(); ++i) + o << l[i] << endl; + + if (l.cyclic()) + o << "cyclic" << endl; + + o << endl; + } + return o; +} + +#define INFO(x) info(cerr,#x,x) + +/***** Auxillary computation functions *****/ + +// Computes the relative distance of a control point given the angles. +// The name is somewhat misleading as the velocity (with respect to the +// variable that parameterizes the path) is relative to the distance +// between the knots and even then, is actually three times this. +// The routine is based on the velocity function in Section 131 of the MetaPost +// code, but differs in that it automatically accounts for the tension and +// bounding with tension atleast. +double velocity(double theta, double phi, tension t) +{ + static const double VELOCITY_BOUND = 4.0; + static const double a = sqrt(2.0); + static const double b = 1.0/16.0; + static const double c = 1.5*(sqrt(5.0)-1.0); + static const double d = 1.5*(3.0-sqrt(5.0)); + + double st = sin(theta), ct = cos(theta), + sf = sin(phi), cf = cos(phi); + + double denom = t.val * (3.0 + c*ct + d*cf); + + double r = denom != 0.0 ? (2.0 + a*(st - b*sf)*(sf - b*st)*(ct-cf)) / denom + : VELOCITY_BOUND; + + //cerr << " velocity(" << theta << "," << phi <<")= " << r << endl; + if (r > VELOCITY_BOUND) + r = VELOCITY_BOUND; + + // Apply boundedness condition for tension atleast cases. + if (t.atleast) + { + double sine = sin(theta + phi); + if ((st >= 0.0 && sf >= 0.0 && sine > 0.0) || + (st <= 0.0 && sf <= 0.0 && sine < 0.0)) + { + double rmax = sf / sine; + if (r > rmax) + r = rmax; + } + } + + return r; +} + +double niceAngle(pair z) +{ + return z.gety() == 0 ? z.getx() >= 0 ? 0 : PI + : angle(z); +} + +// Ensures an angle is in the range between -PI and PI. +double reduceAngle(double angle) +{ + return angle > PI ? angle - 2.0*PI : + angle < -PI ? angle + 2.0*PI : + angle; +} + + +bool interesting(tension t) +{ + return t.val!=1.0 || t.atleast==true; +} + +bool interesting(spec *s) +{ + return !s->open(); +} + +ostream& operator<<(ostream& out, const knot& k) +{ + if (interesting(k.tin)) + out << k.tin << " "; + if (interesting(k.in)) + out << *k.in << " "; + out << k.z; + if (interesting(k.out)) + out << " " << *k.out; + if (interesting(k.tout)) + out << " " << k.tout; + return out; +} + + +eqn dirSpec::eqnOut(Int j, knotlist& l, cvector<double>&, cvector<double>&) +{ + // When choosing the control points, the path will come out the first knot + // going straight to the next knot rotated by the angle theta. + // Therefore, the angle theta we want is the difference between the + // specified heading given and the heading to the next knot. + double theta=reduceAngle(given-niceAngle(l[j+1].z-l[j].z)); + + // Give a simple linear equation to ensure this theta is picked. + return eqn(0.0,1.0,0.0,theta); +} + +eqn dirSpec::eqnIn(Int j, knotlist& l, cvector<double>&, cvector<double>&) +{ + double theta=reduceAngle(given-niceAngle(l[j].z-l[j-1].z)); + return eqn(0.0,1.0,0.0,theta); +} + +eqn curlSpec::eqnOut(Int j, knotlist& l, cvector<double>&, + cvector<double>& psi) +{ + double alpha=l[j].alpha(); + double beta=l[j+1].beta(); + + double chi=alpha*alpha*gamma/(beta*beta); + + double C=alpha*chi+3-beta; + double D=(3.0-alpha)*chi+beta; + + return eqn(0.0,C,D,-D*psi[j+1]); +} + +eqn curlSpec::eqnIn(Int j, knotlist& l, cvector<double>&, cvector<double>&) +{ + double alpha=l[j-1].alpha(); + double beta=l[j].beta(); + + double chi=beta*beta*gamma/(alpha*alpha); + + double A=(3-beta)*chi+alpha; + double B=beta*chi+3-alpha; + + return eqn(A,B,0.0,0.0); +} + +spec *controlSpec::outPartner(pair z) +{ + static curlSpec curl; + return cz==z ? (spec *)&curl : (spec *)new dirSpec(z-cz); +} + +spec *controlSpec::inPartner(pair z) +{ + static curlSpec curl; + return cz==z ? (spec *)&curl : (spec *)new dirSpec(cz-z); +} + +// Compute the displacement between points. The j-th result is the distance +// between knot j and knot j+1. +struct dzprop : public knotprop<pair> { + dzprop(knotlist& l) + : knotprop<pair>(l) {} + + pair solo(Int) { return pair(0,0); } + pair start(Int j) { return l[j+1].z - l[j].z; } + pair mid(Int j) { return l[j+1].z - l[j].z; } + pair end(Int) { return pair(0,0); } +}; + +// Compute the distance between points, using the already computed dz. This +// doesn't use the infomation in the knots, but the knotprop class is useful as +// it takes care of the iteration for us. +struct dprop : public knotprop<double> { + cvector<pair>& dz; + + dprop(knotlist &l, cvector<pair>& dz) + : knotprop<double>(l), dz(dz) {} + + double solo(Int j) { return length(dz[j]); } + double start(Int j) { return length(dz[j]); } + double mid(Int j) { return length(dz[j]); } + double end(Int j) { return length(dz[j]); } +}; + +// Compute the turning angles (psi) between points, using the already computed +// dz. +struct psiprop : public knotprop<double> { + cvector<pair>& dz; + + psiprop(knotlist &l, cvector<pair>& dz) + : knotprop<double>(l), dz(dz) {} + + double solo(Int) { return 0; } + + // We set the starting and ending psi to zero. + double start(Int) { return 0; } + double end(Int) { return 0; } + + double mid(Int j) { return niceAngle(dz[j]/dz[j-1]); } +}; + +struct eqnprop : public knotprop<eqn> { + cvector<double>& d; + cvector<double>& psi; + + eqnprop(knotlist &l, cvector<double>& d, cvector<double>& psi) + : knotprop<eqn>(l), d(d), psi(psi) {} + + eqn solo(Int) { + assert(False); + return eqn(0.0,1.0,0.0,0.0); + } + + eqn start(Int j) { + // Defer to the specifier, as it knows the specifics. + return dynamic_cast<endSpec *>(l[j].out)->eqnOut(j,l,d,psi); + } + + eqn end(Int j) { + return dynamic_cast<endSpec *>(l[j].in)->eqnIn(j,l,d,psi); + } + + eqn mid(Int j) { + double lastAlpha = l[j-1].alpha(); + double thisAlpha = l[j].alpha(); + double thisBeta = l[j].beta(); + double nextBeta = l[j+1].beta(); + + // Values based on the linear approximation of the curvature coming + // into the knot with respect to theta[j-1] and theta[j]. + double inFactor = 1.0/(thisBeta*thisBeta*d[j-1]); + double A = lastAlpha*inFactor; + double B = (3.0 - lastAlpha)*inFactor; + + // Values based on the linear approximation of the curvature going out of + // the knot with respect to theta[j] and theta[j+1]. + double outFactor = 1.0/(thisAlpha*thisAlpha*d[j]); + double C = (3.0 - nextBeta)*outFactor; + double D = nextBeta*outFactor; + + return eqn(A,B+C,D,-B*psi[j]-D*psi[j+1]); + } +}; + +// If the system of equations is homogeneous (ie. we are solving for x in +// Ax=0), there is no need to solve for theta; we can just use zeros for the +// thetas. In fact, our general solving method may not work in this case. +// A common example of this is +// +// a{curl 1}..{curl 1}b +// +// which arises when solving a one-length path a..b or in a larger path a +// section a--b. +bool homogeneous(cvector<eqn>& e) +{ + for(cvector<eqn>::iterator p=e.begin(); p!=e.end(); ++p) + if (p->aug != 0) + return false; + return true; +} + +// Checks whether the equation being solved will be solved as a straight +// path from the first point to the second. +bool straightSection(cvector<eqn>& e) +{ + return e.size()==2 && e.front().aug==0 && e.back().aug==0; +} + +struct weqn : public eqn { + double w; + weqn(double pre, double piv, double post, double aug, double w=0) + : eqn(pre,piv,post,aug), w(w) {} + + friend ostream& operator<< (ostream& out, const weqn we) + { + return out << (eqn &)we << " + " << we.w << " * theta[0]"; + } +}; + +weqn scale(weqn q) { + assert(q.pre == 0 && q.piv != 0); + return weqn(0,1,q.post/q.piv,q.aug/q.piv,q.w/q.piv); +} + +/* Recalculate the equations in the form: + * theta[j] + post * theta[j+1] = aug + w * theta[0] + * + * Used as the first step in solve cyclic equations. + */ +cvector<weqn> recalc(cvector<eqn>& e) +{ + Int n=(Int) e.size(); + cvector<weqn> we; + weqn lasteqn(0,1,0,0,1); + we.push_back(lasteqn); // As a placeholder. + for (Int j=1; j < n; j++) { + // Subtract a factor of the last equation so that the first entry is + // zero, then procede to scale it. + eqn& q=e[j]; + lasteqn=scale(weqn(0,q.piv-q.pre*lasteqn.post,q.post, + q.aug-q.pre*lasteqn.aug,-q.pre*lasteqn.w)); + we.push_back(lasteqn); + } + // To keep all of the infomation enocoded in the linear equations, we need + // to augment the computation to replace our trivial start weqn with a + // real one. To do this, we take one more step in the iteration and + // compute the weqn for j=n, since n=0 (mod n). + { + eqn& q=e[0]; + we.front()=scale(weqn(0,q.piv-q.pre*lasteqn.post,q.post, + q.aug-q.pre*lasteqn.aug,-q.pre*lasteqn.w)); + } + return we; +} + +double solveForTheta0(cvector<weqn>& we) +{ + // Solve for theta[0]=theta[n]. + // How we do this is essentially to write out the first equation as: + // + // theta[n] = aug[0] + w[0]*theta[0] - post[0]*theta[1] + // + // and then use the next equation to substitute in for theta[1]: + // + // theta[1] = aug[1] + w[1]*theta[0] - post[1]*theta[2] + // + // and so on until we have an equation just in terms of theta[0] and + // theta[n] (which are the same theta). + // + // The loop invariant maintained is that after j iterations, we have + // theta[n]= a + b*theta[0] + c*theta[j] + Int n=(Int) we.size(); + double a=0,b=0,c=1; + for (Int j=0;j<n;++j) { + weqn& q=we[j]; + a+=c*q.aug; + b+=c*q.w; + c=-c*q.post; + } + + // After the iteration we have + // + // theta[n] = a + b*theta[0] + c*theta[n] + // + // where theta[n]=theta[0], so + return a/(1.0-(b+c)); +} + +cvector<double> backsubCyclic(cvector<weqn>& we, double theta0) +{ + Int n=(Int) we.size(); + cvector<double> thetas; + double lastTheta=theta0; + for (Int j=1;j<=n;++j) + { + weqn& q=we[n-j]; + assert(q.pre == 0 && q.piv == 1); + double theta=-q.post*lastTheta+q.aug+q.w*theta0; + thetas.push_back(theta); + lastTheta=theta; + } + reverse(thetas.begin(),thetas.end()); + return thetas; +} + +// For the non-cyclic equations, do row operation to put the matrix into +// reduced echelon form, ie. calculates equivalent equations but with pre=0 and +// piv=1 for each eqn. +struct ref : public knotprop<eqn> { + cvector<eqn>& e; + eqn lasteqn; + + ref(knotlist& l, cvector<eqn>& e) + : knotprop<eqn>(l), e(e), lasteqn(0,1,0,0) {} + + // Scale the equation so that the pivot (diagonal) entry is one, and save + // the new equation as lasteqn. + eqn scale(eqn q) { + assert(q.pre == 0 && q.piv != 0); + return lasteqn = eqn(0,1,q.post/q.piv,q.aug/q.piv); + } + + eqn start(Int j) { + return scale(e[j]); + } + eqn mid(Int j) { + // Subtract a factor of the last equation so that the first entry is + // zero, then procede to scale it. + eqn& q=e[j]; + return scale(eqn(0,q.piv-q.pre*lasteqn.post,q.post, + q.aug-q.pre*lasteqn.aug)); + } + // The end case is the same as the middle case. +}; + +// Once the matrix is in reduced echelon form, we can solve for the values by +// back-substitution. This algorithm works from the bottom-up, so backCompute +// must be used to get the answer. +struct backsub : public knotprop<double> { + cvector<eqn>& e; + double lastTheta; + + backsub(knotlist& l, cvector<eqn>& e) + : knotprop<double>(l), e(e) {} + + double end(Int j) { + eqn& q=e[j]; + assert(q.pre == 0 && q.piv == 1 && q.post == 0); + double theta=q.aug; + lastTheta=theta; + return theta; + } + + double mid(Int j) { + eqn& q=e[j]; + assert(q.pre == 0 && q.piv == 1); + double theta=-q.post*lastTheta+q.aug; + lastTheta=theta; + return theta; + } + + // start is the same as mid. +}; + +// Once the equations have been determined, solve for the thetas. +cvector<double> solveThetas(knotlist& l, cvector<eqn>& e) +{ + if (homogeneous(e)) + // We are solving Ax=0, so a solution is zero for every theta. + return cvector<double>(e.size(),0); + else if (l.cyclic()) { + // The knotprop template is unusually unhelpful in this case, so I + // won't use it here. The algorithm breaks into three passes on the + // object. The old Asymptote code used a two-pass method, but I + // implemented this to stay closer to the MetaPost source code. + // This might be something to look at for optimization. + cvector<weqn> we=recalc(e); + INFO(we); + double theta0=solveForTheta0(we); + return backsubCyclic(we, theta0); + } + else { /* Non-cyclic case. */ + /* First do row operations to get it into reduced echelon form. */ + cvector<eqn> el=ref(l,e).compute(); + + /* Then, do back substitution. */ + return backsub(l,el).backCompute(); + } +} + +// Once thetas have been solved, determine the first control point of every +// join. +struct postcontrolprop : public knotprop<pair> { + cvector<pair>& dz; + cvector<double>& psi; + cvector<double>& theta; + + postcontrolprop(knotlist& l, cvector<pair>& dz, + cvector<double>& psi, cvector<double>& theta) + : knotprop<pair>(l), dz(dz), psi(psi), theta(theta) {} + + double phi(Int j) { + /* The third angle: psi + theta + phi = 0 */ + return -psi[j] - theta[j]; + } + + double vel(Int j) { + /* Use the standard velocity function. */ + return velocity(theta[j],phi(j+1),l[j].tout); + } + + // start is the same as mid. + + pair mid(Int j) { + // Put a control point at the relative distance determined by the velocity, + // and at an angle determined by theta. + return l[j].z + vel(j)*expi(theta[j])*dz[j]; + } + + // The end postcontrol is the same as the last knot. + pair end(Int j) { + return l[j].z; + } +}; + +// Determine the first control point of every join. +struct precontrolprop : public knotprop<pair> { + cvector<pair>& dz; + cvector<double>& psi; + cvector<double>& theta; + + precontrolprop(knotlist& l, cvector<pair>& dz, + cvector<double>& psi, cvector<double>& theta) + : knotprop<pair>(l), dz(dz), psi(psi), theta(theta) {} + + double phi(Int j) { + return -psi[j] - theta[j]; + } + + double vel(Int j) { + return velocity(phi(j),theta[j-1],l[j].tin); + } + + // The start precontrol is the same as the first knot. + pair start(Int j) { + return l[j].z; + } + pair mid(Int j) { + return l[j].z - vel(j)*expi(-phi(j))*dz[j-1]; + } + + // end is the same as mid. +}; + +// Puts solved controls into a protopath starting at the given index. +// By convention, the first knot is not coded, as it is assumed to be coded by +// the previous section (or it is the first breakpoint and encoded as a special +// case). +struct encodeControls : public knoteffect { + protopath& p; + Int k; + cvector<pair>& pre; + cvector<pair>& post; + + encodeControls(protopath& p, Int k, + cvector<pair>& pre, knotlist& l, cvector<pair>& post) + : knoteffect(l), p(p), k(k), pre(pre), post(post) {} + + void encodePre(Int j) { + p.pre(k+j)=pre[j]; + } + void encodePoint(Int j) { + p.point(k+j)=l[j].z; + } + void encodePost(Int j) { + p.post(k+j)=post[j]; + } + + void solo(Int) { +#if 0 + encodePoint(j); +#endif + } + void start(Int j) { +#if 0 + encodePoint(j); +#endif + encodePost(j); + } + void mid(Int j) { + encodePre(j); + encodePoint(j); + encodePost(j); + } + void end(Int j) { + encodePre(j); + encodePoint(j); + } +}; + +void encodeStraight(protopath& p, Int k, knotlist& l) +{ + pair a=l.front().z; + double at=l.front().tout.val; + pair b=l.back().z; + double bt=l.back().tin.val; + pair step=(b-a)/3.0; + + if (at==1.0 && bt==1.0) { + p.straight(k)=true; + p.post(k)=a+step; + p.pre(k+1)=b-step; + p.point(k+1)=b; + } + else { + p.post(k)=a+step/at; + p.pre(k+1)=b-step/bt; + p.point(k+1)=b; + } +} + +void solveSection(protopath& p, Int k, knotlist& l) +{ + if (l.length()>0) { + info(cerr, "solving section", l); + + // Calculate useful properties. + cvector<pair> dz = dzprop(l) .compute(); + cvector<double> d = dprop(l,dz).compute(); + cvector<double> psi = psiprop(l,dz).compute(); + + INFO(dz); INFO(d); INFO(psi); + + // Build and solve the linear equations for theta. + cvector<eqn> e = eqnprop(l,d,psi).compute(); + INFO(e); + + if (straightSection(e)) + // Handle straight section as special case. + encodeStraight(p,k,l); + else { + cvector<double> theta = solveThetas(l,e); + INFO(theta); + + // Calculate the control points. + cvector<pair> post = postcontrolprop(l,dz,psi,theta).compute(); + cvector<pair> pre = precontrolprop(l,dz,psi,theta).compute(); + + // Encode the results into the protopath. + encodeControls(p,k,pre,l,post).exec(); + } + } +} + +// Find the first breakpoint in the knotlist, ie. where we can start solving a +// non-cyclic section. If the knotlist is fully cyclic, then this returns +// NOBREAK. +// This must be called with a knot that has all of its implicit specifiers in +// place. +const Int NOBREAK=-1; +Int firstBreakpoint(knotlist& l) +{ + for (Int j=0;j<l.size();++j) + if (!l[j].out->open()) + return j; + return NOBREAK; +} + +// Once a breakpoint, a, is found, find where the next breakpoint after it is. +// This must be called with a knot that has all of its implicit specifiers in +// place, so that breakpoint can be identified by either an in or out specifier +// that is not open. +Int nextBreakpoint(knotlist& l, Int a) +{ + // This is guaranteed to terminate if a is the index of a breakpoint. If the + // path is non-cyclic it will stop at or before the last knot which must be a + // breakpoint. If the path is cyclic, it will stop at or before looping back + // around to a which is a breakpoint. + Int j=a+1; + while (l[j].in->open()) + ++j; + return j; +} + +// Write out the controls for section of the form +// a.. control b and c ..d +void writeControls(protopath& p, Int a, knotlist& l) +{ + // By convention, the first point will already be encoded. + p.straight(a)=dynamic_cast<controlSpec *>(l[a].out)->straight; + p.post(a)=dynamic_cast<controlSpec *>(l[a].out)->cz; + p.pre(a+1)=dynamic_cast<controlSpec *>(l[a+1].in)->cz; + p.point(a+1)=l[a+1].z; +} + +// Solves a path that has all of its specifiers laid out explicitly. +path solveSpecified(knotlist& l) +{ + protopath p(l.size(),l.cyclic()); + + Int first=firstBreakpoint(l); + if (first==NOBREAK) + /* We are solving a fully cyclic path, so do it in one swoop. */ + solveSection(p,0,l); + else { + // Encode the first point. + p.point(first)=l[first].z; + + // If the path is cyclic, we should stop where we started (modulo the + // length of the path); otherwise, just stop at the end. + Int last=l.cyclic() ? first+l.length() + : l.length(); + Int a=first; + while (a!=last) { + if (l[a].out->controlled()) { + assert(l[a+1].in->controlled()); + + // Controls are already picked, just write them out. + writeControls(p,a,l); + ++a; + } + else { + // Find the section a to b and solve it, putting the result (starting + // from index a into our protopath. + Int b=nextBreakpoint(l,a); + subknotlist section(l,a,b); + solveSection(p,a,section); + a=b; + } + } + + // For a non-cyclic path, the end control points need to be set. + p.controlEnds(); + } + + return p.fix(); +} + +/* If a knot is open on one side and restricted on the other, this replaces the + * open side with a restriction determined by the restriction on the other + * side. After this, any knot will either have two open specs or two + * restrictions. + */ +struct partnerUp : public knoteffect { + partnerUp(knotlist& l) + : knoteffect(l) {} + + void mid(Int j) { + knot& k=l[j]; + if (k.in->open() && !k.out->open()) + k.in=k.out->inPartner(k.z); + else if (!k.in->open() && k.out->open()) + k.out=k.in->outPartner(k.z); + } +}; + +/* Ensures a non-cyclic path has direction specifiers at the ends, adding curls + * if there are none. + */ +void curlEnds(knotlist& l) +{ + static curlSpec endSpec; + + if (!l.cyclic()) { + if (l.front().in->open()) + l.front().in=&endSpec; + if (l.back().out->open()) + l.back().out=&endSpec; + } +} + +/* If a point occurs twice in a row in a knotlist, write in controls + * between the two knots at that point (unless it already has controls). + */ +struct controlDuplicates : public knoteffect { + controlDuplicates(knotlist& l) + : knoteffect(l) {} + + void solo(Int) { /* One point ==> no duplicates */ } + // start is the same as mid. + void mid(Int j) { + knot &k1=l[j]; + knot &k2=l[j+1]; + if (!k1.out->controlled() && k1.z==k2.z) { + k1.out=k2.in=new controlSpec(k1.z,true); + } + } + void end(Int) { /* No next point to compare with. */ } +}; + +path solve(knotlist& l) +{ + if (l.empty()) + return path(); + else { + info(cerr, "input knotlist", l); + curlEnds(l); + controlDuplicates(l).exec(); + partnerUp(l).exec(); + info(cerr, "specified knotlist", l); + return solveSpecified(l); + } +} + +// Code for Testing +#if 0 +path solveSimple(cvector<pair>& z) +{ + // The two specifiers used: an open spec and a curl spec for the ends. + spec open; + +// curlSpec curl; +// curlSpec curly(2.0); +// dirSpec E(0); +// dirSpec N(PI/2.0); + + controlSpec here(pair(150,150)); + + // Encode the knots as open in the knotlist. + cvector<knot> nodes; + for (cvector<pair>::iterator p=z.begin(); p!=z.end(); ++p) { + knot k; + k.z=*p; + k.in=k.out=&open; + + nodes.push_back(k); + } + + // Substitute in a curl spec for the ends. + //nodes.front().out=nodes.back().in=&curl; + + // Test direction specifiers. + //nodes.front().tout=2; + //nodes.front().out=nodes.back().in=&curly; + + //nodes[0].out=nodes[0].in=&E; + nodes[1].out=nodes[2].in=&here; + + simpleknotlist l(nodes,false); + return solve(l); +} +#endif + +} // namespace camp |