diff --git a/parallel_transport/notebook.ipynb b/parallel_transport/notebook.ipynb
new file mode 100644
index 0000000..88b3e24
--- /dev/null
+++ b/parallel_transport/notebook.ipynb
@@ -0,0 +1,649 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Parallel Transport Frames scratch work\n",
+ "\n",
+ "This is based on a couple things.\n",
+ "\n",
+ "It follows the paper, [\"Parallel Transport Approach to Curve Framing\" by Andrew J. Hanson & Hui Ma](https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.42.8103), which I first encountered via the [geom-types/src/ptf.org](https://github.com/thi-ng/geom/blob/master/geom-types/src/ptf.org) examples from Karsten Schmidt's work with [thi.ng/geom](https://github.com/thi-ng/geom). If I refer to any page number or section, I'm talking about this paper.\n",
+ "\n",
+ "It was also partly an experiment in learning a bit of [JupyterLab](https://github.com/jupyterlab/jupyterlab), which I only tried out after running into some [Julia](https://julialang.org/) problems on Jupyter. I ordinarily also work in [Maxima](http://maxima.sourceforge.net/), but decided to try out [SymPy](https://www.sympy.org/en/index.html). I was rather impressed with both: JupyterLab addresses many of my complaints with Jupyter, and I found SymPy quite capable and part of a more efficient workflow compared to manually transferring results from Maxima into Python.\n",
+ "\n",
+ "## TODO items\n",
+ "\n",
+ "- Don't require template to be a single closed path\n",
+ "- Figure out why I am leaving one lone vertex at the origin\n",
+ "- Figure out \"face_normals didn't match triangles, ignoring!\" (if it matters)\n",
+ "- Find a way to start out the trimesh view with things visible\n",
+ "- Figure out why this knot is not actually a [cinquefoil](https://en.wikipedia.org/wiki/Cinquefoil_knot)..."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import sympy\n",
+ "from sympy import symbols, cos, sin, diff, lambdify\n",
+ "from sympy.vector import CoordSys3D\n",
+ "import numpy\n",
+ "import quaternion\n",
+ "import stl.mesh\n",
+ "import trimesh\n",
+ "\n",
+ "import parallel_transport\n",
+ "import quat"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "(p,q) [Torus knot](https://en.wikipedia.org/wiki/Torus_knot) as a parametric formula (torus lies in XY plane):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "r,p,q,t = symbols(\"r,p,q,t\")\n",
+ "r = cos(q*t)+2\n",
+ "N = CoordSys3D('N')\n",
+ "knot_vec = (r * cos(p*t))*N.i + (r * sin(p*t))*N.j + (-sin(q*t))*N.k"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def vec2python(v):\n",
+ " return sympy.python(N.origin.locate_new('na', v).express_coordinates(N))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "q = Symbol('q')\n",
+ "t = Symbol('t')\n",
+ "p = Symbol('p')\n",
+ "e = ((cos(q*t) + 2)*cos(p*t), (cos(q*t) + 2)*sin(p*t), -sin(q*t))\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(vec2python(knot_vec))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "p = Symbol('p')\n",
+ "q = Symbol('q')\n",
+ "t = Symbol('t')\n",
+ "e = (-p*(cos(q*t) + 2)*sin(p*t) - q*sin(q*t)*cos(p*t), p*(cos(q*t) + 2)*cos(p*t) - q*sin(p*t)*sin(q*t), -q*cos(q*t))\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(vec2python(diff(knot_vec,t)))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[ 3. , 0. , -0. ],\n",
+ " [ 2.99960977, 0.04714521, -0.00628629],\n",
+ " [ 2.9984392 , 0.09427692, -0.01257233],\n",
+ " [ 2.99648866, 0.14138161, -0.01885787],\n",
+ " [ 2.99375876, 0.18844581, -0.02514266],\n",
+ " [ 2.99025034, 0.23545602, -0.03142647],\n",
+ " [ 2.9859645 , 0.28239879, -0.03770903],\n",
+ " [ 2.98090257, 0.32926067, -0.0439901 ],\n",
+ " [ 2.97506613, 0.37602826, -0.05026943],\n",
+ " [ 2.968457 , 0.42268818, -0.05654678]])"
+ ]
+ },
+ "execution_count": 6,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# can I do this without having to make it a Point?\n",
+ "torus52_expr = knot_vec.subs(p, 5).subs(q, 2)\n",
+ "torus52 = sympy.lambdify([t], N.origin.locate_new('na', torus52_expr).express_coordinates(N))\n",
+ "torus52_tan = sympy.lambdify([t], N.origin.locate_new('na', diff(torus52_expr, t)).express_coordinates(N))\n",
+ "#ts = numpy.linspace(0, 2*numpy.pi, 50000)\n",
+ "ts = numpy.linspace(0, 2*numpy.pi, 2000)\n",
+ "points = numpy.array(torus52(ts)).T\n",
+ "points[:10,:]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[-0. , 0.9912279 , -0.13216372],\n",
+ " [-0.01640804, 0.99109232, -0.13216196],\n",
+ " [-0.0328116 , 0.99068562, -0.1321567 ],\n",
+ " [-0.04920622, 0.99000792, -0.13214791],\n",
+ " [-0.06558742, 0.98905939, -0.13213562],\n",
+ " [-0.08195074, 0.98784029, -0.13211981],\n",
+ " [-0.09829171, 0.98635095, -0.13210049],\n",
+ " [-0.11460587, 0.98459179, -0.13207766],\n",
+ " [-0.13088879, 0.98256327, -0.1320513 ],\n",
+ " [-0.14713601, 0.98026595, -0.13202143]])"
+ ]
+ },
+ "execution_count": 7,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# analytical tangents:\n",
+ "tans = numpy.array(torus52_tan(ts)).T\n",
+ "tans = tans / numpy.linalg.norm(tans, axis=1)[:,numpy.newaxis]\n",
+ "tans[:10,:]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "SymPy is failing at the analytical integration, but also [numerical](https://docs.sympy.org/latest/modules/evalf.html#sums-and-integrals). I tried oscillatory quadrature, but it seems that despite the cos/sin everywhere, they're not in the right form."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# tor = parallel_transport.torsion_integral(torus52_expr, t, 0, 2*sympy.pi)\n",
+ "# Works (returns 1.408) but is horrendously slow\n",
+ "#tor_num = tor.evalf(4)\n",
+ "#tor_num"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "My own numerical approach seems to work:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "-1.4071908841197978"
+ ]
+ },
+ "execution_count": 9,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "parallel_transport.torsion_integral_approx(torus52_expr, t, 0, 2*numpy.pi, 0.001)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Do I want analytical tangents, or approximations? approximations may, oddly enough, be 'more correct' given that I am already approximating the curve."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "points = numpy.array(torus52(ts)).T\n",
+ "tan_approx = parallel_transport.approx_tangents(points)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "These are rather close regardless:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "9.38346943765395e-06"
+ ]
+ },
+ "execution_count": 11,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "numpy.mean(numpy.abs(tans - tan_approx))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "First, make some sort of template I can duplicate. Use a vertex ordering that can go directly into the STL library."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "template = numpy.array([\n",
+ " [0.0, 0.0, 0.0],\n",
+ " [0.5, 0.0, 0.0],\n",
+ " [0.5, 0.0, 0.5],\n",
+ "])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "and attempt the algorithm from \"Parallel Transport to Curve Framing\":"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Angular change: -1.4071908841197978, -80.62609862934731 deg\n",
+ "Approximate arc length: 64.25469846213409\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Precompute cross-products, angles, arc length, and angular difference:\n",
+ "t1 = tans\n",
+ "t2 = numpy.roll(tans, -1, axis=0)\n",
+ "cps = numpy.cross(t1, t2)\n",
+ "cps /= numpy.linalg.norm(cps, axis=1)[:,numpy.newaxis]\n",
+ "thetas = numpy.arccos(numpy.sum(t1*t2, axis=1))\n",
+ "ang = parallel_transport.torsion_integral_approx(torus52_expr, t, 0, 2*numpy.pi, 0.001)\n",
+ "print(\"Angular change: {}, {} deg\".format(ang, ang/numpy.pi*180))\n",
+ "total_length = parallel_transport.approx_arc_length(points)\n",
+ "print(\"Approximate arc length: {}\".format(total_length))\n",
+ "# One can add additional twists here if a multiple of a whole rotation:\n",
+ "# ang += 40 * 2*numpy.pi"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "done\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Build up the actual mesh, correcting spin as we go:\n",
+ "count = len(points)\n",
+ "data = numpy.zeros(count*2*3*len(template), dtype=stl.mesh.Mesh.dtype)\n",
+ "# Each step on the curve involves len(template) vertices which\n",
+ "# must connect to respective vertices of the next step. This\n",
+ "# requires two triangles for each vertex of the template, and\n",
+ "# each triangle requires 3 vertices.\n",
+ "vectors = data[\"vectors\"]\n",
+ "length = 0.0\n",
+ "# Use q for a running transformation of our template:\n",
+ "q = numpy.quaternion(1, 0, 0, 0)\n",
+ "for i in range(count):\n",
+ " t1 = tans[i,:]\n",
+ " b = cps[i]\n",
+ " theta = thetas[i]\n",
+ " # Add another transform (note we post-multiply):\n",
+ " q = quat.rotation_quaternion(b, theta) * q\n",
+ " length += numpy.linalg.norm(points[i,:] - points[(i+1) % count,:])\n",
+ " # Correction by arc length (see section 3.1):\n",
+ " theta2 = -ang * length / total_length\n",
+ " # Add a correction for twist. This doesn't update our\n",
+ " # running transform - it's only for final geometry:\n",
+ " q_corrected = quat.rotation_quaternion(t1, theta2) * q\n",
+ " # Connect every vertex in the frame (the template now rotated &\n",
+ " # shifted) to the 'last' frame location via 2 triangles\n",
+ " p_adj = quat.conjugate_by(template.copy(), q_corrected) + points[i,:]\n",
+ " if i == 0:\n",
+ " first_p_adj = p_adj\n",
+ " else:\n",
+ " for j,_ in enumerate(p_adj):\n",
+ " idx = len(template)*2*3*(i-1) + j*2\n",
+ " f1,f2 = parallel_transport.gen_faces(\n",
+ " p_adj[(j+1) % len(template),:],\n",
+ " p_adj[j,:],\n",
+ " last_p_adj[(j+1) % len(template),:],\n",
+ " last_p_adj[j,:],\n",
+ " )\n",
+ " vectors[idx] = f1\n",
+ " vectors[idx+1] = f2\n",
+ " last_p_adj = p_adj\n",
+ "\n",
+ "# Patch up first & last:\n",
+ "for j,_ in enumerate(last_p_adj):\n",
+ " idx = len(template)*2*3*(count-1) + j*2\n",
+ " f1,f2 = parallel_transport.gen_faces(\n",
+ " last_p_adj[j,:],\n",
+ " last_p_adj[(j+1) % len(template),:],\n",
+ " first_p_adj[j,:],\n",
+ " first_p_adj[(j+1) % len(template),:])\n",
+ " vectors[idx] = f1\n",
+ " vectors[idx+1] = f2\n",
+ "\n",
+ "mesh = stl.mesh.Mesh(data)\n",
+ "mesh_fname = \"test_2k.stl\"\n",
+ "mesh.save(mesh_fname)\n",
+ "print(\"done\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "face_normals didn't match triangles, ignoring!\n",
+ "/home/hodapp/.local/lib/python3.6/site-packages/IPython/core/display.py:694: UserWarning: Consider using IPython.display.IFrame instead\n",
+ " warnings.warn(\"Consider using IPython.display.IFrame instead\")\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 15,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "m = trimesh.load_mesh(mesh_fname)\n",
+ "m.show()"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.6.8"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/parallel_transport/parallel_transport.py b/parallel_transport/parallel_transport.py
new file mode 100644
index 0000000..45bd5b0
--- /dev/null
+++ b/parallel_transport/parallel_transport.py
@@ -0,0 +1,84 @@
+import numpy
+import sympy
+from sympy.vector import CoordSys3D
+
+def mtx_rotate_by_vector(b, theta):
+ """Returns 3x3 matrix for rotating around some 3D vector."""
+ # Source:
+ # Appendix A of "Parallel Transport to Curve Framing"
+ c = numpy.cos(theta)
+ s = numpy.sin(theta)
+ rot = numpy.array([
+ [c+(b[0]**2)*(1-c), b[0]*b[1]*(1-c)-s*b[2], b[2]*b[0]*(1-c)+s*b[1]],
+ [b[0]*b[1]*(1-c)+s*b[2], c+(b[1]**2)*(1-c), b[2]*b[1]*(1-c)-s*b[0]],
+ [b[0]*b[2]*(1-c)-s*b[1], b[1]*b[2]*(1-c)+s*b[0], c+(b[2]**2)*(1-c)],
+ ])
+ return rot
+
+def gen_faces(v1a, v1b, v2a, v2b):
+ """Returns faces (as arrays of vertices) connecting two pairs of
+ vertices."""
+ # Keep winding order consistent!
+ f1 = numpy.array([v1b, v1a, v2a])
+ f2 = numpy.array([v2b, v1b, v2a])
+ return f1, f2
+
+def approx_tangents(points):
+ """Returns an array of approximate tangent vectors. Assumes a
+ closed path, and approximates point I using neighbors I-1 and I+1 -
+ that is, treating the three points as a circle.
+
+ Input:
+ points -- Array of shape (N,3). points[0,:] is assumed to wrap around
+ to points[-1,:].
+
+ Output:
+ tangents -- Array of same shape as 'points'. Each row is normalized.
+ """
+ d = numpy.roll(points, -1, axis=0) - numpy.roll(points, +1, axis=0)
+ d = d/numpy.linalg.norm(d, axis=1)[:,numpy.newaxis]
+ return d
+
+def approx_arc_length(points):
+ p2 = numpy.roll(points, -1, axis=0)
+ return numpy.sum(numpy.linalg.norm(points - p2, axis=1))
+
+def torsion(v, arg):
+ """Returns an analytical SymPy expression for torsion of a 3D curve.
+
+ Inputs:
+ v -- SymPy expression returning a 3D vector
+ arg -- SymPy symbol for v's variable
+ """
+ # https://en.wikipedia.org/wiki/Torsion_of_a_curve#Alternative_description
+ dv1 = v.diff(arg)
+ dv2 = dv1.diff(arg)
+ dv3 = dv2.diff(arg)
+ v1_x_v2 = dv1.cross(dv2)
+ # This calls for the square of the norm in denominator - but that
+ # is just dot product with itself:
+ return v1_x_v2.dot(dv3) / (v1_x_v2.dot(v1_x_v2))
+
+def torsion_integral(curve_expr, var, a, b):
+ # The line integral from section 3.1 of "Parallel Transport to Curve
+ # Framing". This should work in theory, but with the functions I've
+ # actually tried, evalf() is ridiculously slow.
+ c = torsion(curve_expr, var)
+ return sympy.Integral(c * (sympy.diff(curve_expr, var).magnitude()), (var, a, b))
+
+def torsion_integral_approx(curve_expr, var, a, b, step):
+ # A numerical approximation of the line integral from section 3.1 of
+ # "Parallel Transport to Curve Framing"
+ N = CoordSys3D('N')
+ # Get a (callable) derivative function of the curve:
+ curve_diff = curve_expr.diff(var)
+ diff_fn = sympy.lambdify([var], N.origin.locate_new('na', curve_diff).express_coordinates(N))
+ # And a torsion function:
+ torsion_fn = sympy.lambdify([var], torsion(curve_expr, var))
+ # Generate values of 'var' to use:
+ vs = numpy.arange(a, b, step)
+ # Evaluate derivative function & torsion function over these:
+ d = numpy.linalg.norm(numpy.array(diff_fn(vs)).T, axis=1)
+ torsions = torsion_fn(vs)
+ # Turn this into basically a left Riemann sum (I'm lazy):
+ return -(d * torsions * step).sum()
diff --git a/parallel_transport/quat.py b/parallel_transport/quat.py
new file mode 100644
index 0000000..d8cbb90
--- /dev/null
+++ b/parallel_transport/quat.py
@@ -0,0 +1,24 @@
+import numpy
+import quaternion
+
+def conjugate_by(vec, quat):
+ """Turn 'vec' to a quaternion, conjugate it by 'quat', and return it."""
+ q2 = quat * vec2quat(vec) * quat.conjugate()
+ return quaternion.as_float_array(q2)[:,1:]
+
+def rotation_quaternion(axis, angle):
+ """Returns a quaternion for rotating by some axis and angle.
+
+ Inputs:
+ axis -- numpy array of shape (3,), with axis to rotate around
+ angle -- angle in radians by which to rotate
+ """
+ qc = numpy.cos(angle / 2)
+ qs = numpy.sin(angle / 2)
+ qv = qs * axis
+ return numpy.quaternion(qc, qv[0], qv[1], qv[2])
+
+def vec2quat(vs):
+ qs = numpy.zeros(vs.shape[0], dtype=numpy.quaternion)
+ quaternion.as_float_array(qs)[:,1:4] = vs
+ return qs
diff --git a/parallel_transport/spiral_parametric3.py b/parallel_transport/spiral_parametric3.py
new file mode 100644
index 0000000..2d88461
--- /dev/null
+++ b/parallel_transport/spiral_parametric3.py
@@ -0,0 +1,232 @@
+#!/usr/bin/env python3
+
+import sys
+import numpy
+import stl.mesh
+
+# TODO:
+# - This still has some strange errors around high curvature.
+# They are plainly visible in the cross-section.
+# - Check rotation direction
+# - Fix phase, which only works if 0 (due to how I work with y)
+# Things don't seem to line up right.
+# - Why is there still a gap when using Array modifier?
+# Check beginning and ending vertices maybe
+# - Organize this so that it generates both meshes when run
+
+# This is all rather tightly-coupled. Almost everything is specific
+# to the isosurface I was trying to generate. walk_curve may be able
+# to generalize to some other shapes.
+class ExplicitSurfaceThing(object):
+ def __init__(self, freq, phase, scale, inner, outer, rad, ext_phase):
+ self.freq = freq
+ self.phase = phase
+ self.scale = scale
+ self.inner = inner
+ self.outer = outer
+ self.rad = rad
+ self.ext_phase = ext_phase
+
+ def angle(self, z):
+ return self.freq*z + self.phase
+
+ def max_z(self):
+ # This value is the largest |z| for which 'radical' >= 0
+ # (thus, for x_cross to have a valid solution)
+ return (numpy.arcsin(self.rad / self.inner) - self.phase) / self.freq
+
+ def radical(self, z):
+ return self.rad*self.rad - self.inner*self.inner * (numpy.sin(self.angle(z)))**2
+
+ # Implicit curve function
+ def F(self, x, z):
+ return (self.outer*x - self.inner*numpy.cos(self.angle(z)))**2 + (self.inner*numpy.sin(self.angle(z)))**2 - self.rad**2
+
+ # Partial 1st derivatives of F:
+ def F_x(self, x, z):
+ return 2 * self.outer * self.outer * x - 2 * self.inner * self.outer * numpy.cos(self.angle(z))
+
+ def F_z(self, x, z):
+ return 2 * self.freq * self.inner * self.outer * numpy.sin(self.angle(z))
+
+ # Curvature:
+ def K(self, x, z):
+ a1 = self.outer**2
+ a2 = x**2
+ a3 = self.freq*z + self.phase
+ a4 = numpy.cos(a3)
+ a5 = self.inner**2
+ a6 = a4**2
+ a7 = self.freq**2
+ a8 = numpy.sin(a3)**2
+ a9 = self.outer**3
+ a10 = self.inner**3
+ return -((2*a7*a10*self.outer*x*a4 + 2*a7*a5*a1*a2)*a8 + (2*a7*self.inner*a9*x**3 + 2*a7*a10*self.outer*x)*a4 - 4*a7*a5*a1*a2) / ((a7*a5*a2*a8 + a5*a6 - 2*self.inner*self.outer*x*a4 + a1*a2) * numpy.sqrt(4*a7*a5*a1*a2*a8 + 4*a5*a1*a6 - 8*self.inner*a9*x*a4 + 4*a2*self.outer**4))
+
+ def walk_curve(self, x0, z0, eps, thresh = 1e-3, gd_thresh = 1e-7):
+ x, z = x0, z0
+ eps2 = eps*eps
+ verts = []
+ iters = 0
+ # Until we return to the same point at which we started...
+ while True:
+ iters += 1
+ verts.append([x, 0, z])
+ # ...walk around the curve by stepping perpendicular to the
+ # gradient by 'eps'. So, first find the gradient:
+ dx = self.F_x(x, z)
+ dz = self.F_z(x, z)
+ # Normalize it:
+ f = 1/numpy.sqrt(dx*dx + dz*dz)
+ nx, nz = dx*f, dz*f
+ # Find curvature at this point because it tells us a little
+ # about how far we can safely move:
+ K_val = abs(self.K(x, z))
+ eps_corr = 2 * numpy.sqrt(2*eps/K_val - eps*eps)
+ # Scale by 'eps' and use (-dz, dx) as perpendicular:
+ px, pz = -nz*eps_corr, nx*eps_corr
+ # Walk in that direction:
+ x += px
+ z += pz
+ # Moving in that direction is only good locally, and we may
+ # have deviated off the curve slightly. The implicit function
+ # tells us (sort of) how far away we are, and the gradient
+ # tells us how to minimize that:
+ #print("W: x={} z={} dx={} dz={} px={} pz={} K={} eps_corr={}".format(
+ # x, z, dx, dz, px, pz, K_val, eps_corr))
+ F_val = self.F(x, z)
+ count = 0
+ while abs(F_val) > gd_thresh:
+ count += 1
+ dx = self.F_x(x, z)
+ dz = self.F_z(x, z)
+ f = 1/numpy.sqrt(dx*dx + dz*dz)
+ nx, nz = dx*f, dz*f
+ # If F is negative, we want to increase it (thus, follow
+ # gradient). If F is positive, we want to decrease it
+ # (thus, opposite of gradient).
+ F_val = self.F(x, z)
+ x += -F_val*nx
+ z += -F_val*nz
+ # Yes, this is inefficient gradient-descent...
+ diff = numpy.sqrt((x-x0)**2 + (z-z0)**2)
+ #print("{} gradient-descent iters. diff = {}".format(count, diff))
+ if iters > 100 and diff < thresh:
+ #print("diff < eps, quitting")
+ #verts.append([x, 0, z])
+ break
+ data = numpy.array(verts)
+ return data
+
+ def x_cross(self, z, sign):
+ # Single cross-section point in XZ for y=0. Set sign for positive
+ # or negative solution.
+ n1 = numpy.sqrt(self.radical(z))
+ n2 = self.inner * numpy.cos(self.angle(z))
+ if sign > 0:
+ return (n2-n1) / self.outer
+ else:
+ return (n2+n1) / self.outer
+
+ def turn(self, points, dz):
+ # Note one full revolution is dz = 2*pi/freq
+ # How far to turn in radians (determined by dz):
+ rad = self.angle(dz)
+ c, s = numpy.cos(rad), numpy.sin(rad)
+ mtx = numpy.array([
+ [ c, s, 0],
+ [-s, c, 0],
+ [ 0, 0, 1],
+ ])
+ return points.dot(mtx) + [0, 0, dz]
+
+ def screw_360(self, z0_period_start, x_init, z_init, eps, dz, thresh, endcaps=False):
+ #z0 = -10 * 2*numpy.pi/freq / 2
+ z0 = z0_period_start * 2*numpy.pi/self.freq / 2
+ z1 = z0 + 2*numpy.pi/self.freq
+ #z1 = 5 * 2*numpy.pi/freq / 2
+ #z0 = 0
+ #z1 = 2*numpy.pi/freq
+ init_xsec = self.walk_curve(x_init, z_init, eps, thresh)
+ num_xsec_steps = init_xsec.shape[0]
+ zs = numpy.append(numpy.arange(z0, z1, dz), z1)
+ num_screw_steps = len(zs)
+ vecs = num_xsec_steps * num_screw_steps * 2
+ offset = 0
+ if endcaps:
+ offset = num_xsec_steps
+ vecs += 2*num_xsec_steps
+ print("Generating {} vertices...".format(vecs))
+ data = numpy.zeros(vecs, dtype=stl.mesh.Mesh.dtype)
+ v = data["vectors"]
+ # First endcap:
+ if endcaps:
+ center = init_xsec.mean(0)
+ for i in range(num_xsec_steps):
+ v[i][0,:] = init_xsec[(i + 1) % num_xsec_steps,:]
+ v[i][1,:] = init_xsec[i,:]
+ v[i][2,:] = center
+ # Body:
+ verts = init_xsec
+ for i,z in enumerate(zs):
+ verts_last = verts
+ verts = self.turn(init_xsec, z-z0)
+ if i > 0:
+ for j in range(num_xsec_steps):
+ # Vertex index:
+ vi = offset + (i-1)*num_xsec_steps*2 + j*2
+ v[vi][0,:] = verts[(j + 1) % num_xsec_steps,:]
+ v[vi][1,:] = verts[j,:]
+ v[vi][2,:] = verts_last[j,:]
+ #print("Write vertex {}".format(vi))
+ v[vi+1][0,:] = verts_last[(j + 1) % num_xsec_steps,:]
+ v[vi+1][1,:] = verts[(j + 1) % num_xsec_steps,:]
+ v[vi+1][2,:] = verts_last[j,:]
+ #print("Write vertex {} (2nd half)".format(vi+1))
+ # Second endcap:
+ if endcaps:
+ center = verts.mean(0)
+ for i in range(num_xsec_steps):
+ vi = num_xsec_steps * num_screw_steps * 2 + num_xsec_steps + i
+ v[vi][0,:] = center
+ v[vi][1,:] = verts[i,:]
+ v[vi][2,:] = verts[(i + 1) % num_xsec_steps,:]
+ v[:, :, 2] += z0 + self.ext_phase / self.freq
+ v[:, :, :] /= self.scale
+ mesh = stl.mesh.Mesh(data, remove_empty_areas=False)
+ print("Beginning z: {}".format(z0/self.scale))
+ print("Ending z: {}".format(z1/self.scale))
+ print("Period: {}".format((z1-z0)/self.scale))
+ return mesh
+
+surf1 = ExplicitSurfaceThing(
+ freq = 20,
+ phase = 0,
+ scale = 1/16, # from libfive
+ inner = 0.4 * 1/16,
+ outer = 2.0 * 1/16,
+ rad = 0.3 * 1/16,
+ ext_phase = 0)
+
+z_init = 0
+x_init = surf1.x_cross(z_init, 1)
+mesh1 = surf1.screw_360(-10, x_init, z_init, 0.000002, 0.001, 5e-4)
+fname = "spiral_inner0_one_period.stl"
+print("Writing {}...".format(fname))
+mesh1.save(fname)
+
+surf2 = ExplicitSurfaceThing(
+ freq = 10,
+ phase = 0,
+ scale = 1/16, # from libfive
+ inner = 0.9 * 1/16,
+ outer = 2.0 * 1/16,
+ rad = 0.3 * 1/16,
+ ext_phase = numpy.pi/2)
+
+z_init = 0
+x_init = surf2.x_cross(z_init, 1)
+mesh2 = surf2.screw_360(-5, x_init, z_init, 0.000002, 0.001, 5e-4)
+fname = "spiral_outer90_one_period.stl"
+print("Writing {}...".format(fname))
+mesh2.save(fname)