helixish: bugfixes

Signed-off-by: Ian Jackson <ijackson@chiark.greenend.org.uk>

diff --git a/helixish.py b/helixish.py
index 850b225985874f5835489b0f6cb35ab6a8cb502e..d93d75f516f59d786acdd017c0149de675746b02 100644 (file)
--- a/helixish.py
+++ b/helixish.py
@@ -5,15 +5,17 @@ import numpy as np
 from numpy import cos, sin
 
 import sys
+import subprocess
+
 from moedebug import dbg
 from moenp import *
 
-from math import atan2
+from math import atan2, atan, sqrt
 
 import symbolic
 
-def augment(v): return np.append(v, 1)
-def augment0(v): return np.append(v, 0)
+def augment(v, augwith=1): return np.append(v, 1)
+def augment0(v): return augment(v, 0)
 def unaugment(v): return v[0:3]
 
 def matmultiply(mat,vect):
@@ -22,8 +24,14 @@ def matmultiply(mat,vect):
   # but that doesn't work in Python 2
   return np.array((vect * np.matrix(mat).T))[0,:]
 
-def augmatmultiply(mat,unaugvect):
-  return unaugment(matmultiply(mat, augment(unaugvect)))
+def matmatmultiply(mat1,mat2):
+  # both are "array"s
+  # we would prefer to write   mat1 @ mat2
+  # but that doesn't work in Python 2
+  return np.array((np.matrix(mat1) * np.matrix(mat2)))
+
+def augmatmultiply(mat,unaugvect, augwith=1):
+  return unaugment(matmultiply(mat, augment(unaugvect, augwith)))
 
 findcurve_subproc = None
 
@@ -63,8 +71,8 @@ class HelixishCurve():
     dbg(dPQplane_basis)
     dPQplane_into = np.linalg.inv(dPQplane_basis)
 
-    dp_plane = augmatmultiply(dPQplane_into, dp)
-    dq_plane = augmatmultiply(dPQplane_into, dq)
+    dp_plane = augmatmultiply(dPQplane_into, dp, augwith=0)
+    dq_plane = augmatmultiply(dPQplane_into, dq, augwith=0)
     q_plane  = augmatmultiply(dPQplane_into, q)
     dist_pq_plane = np.linalg.norm(q_plane)
 
@@ -113,33 +121,35 @@ class HelixishCurve():
       start_t = dist_pq_plane * .35
       start_mu = 0.05
 
-    bodge = max( q_plane[2] * mu,
+    bodge = max( q_plane[2] * start_mu,
                  (start_s + start_t) * 0.1 )
     start_s += 0.5 * bodge
     start_t += 0.5 * bodge
     start_kappa = 0
     start_gamma = 1
 
-    tilt = atan(mu)
+    tilt = atan(start_mu)
     tilt_basis = np.array([
-      1,     0,           0,         0,
-      0,   cos(tilt), -sin(tilt),    0,
-      0,   sin(tilt),  cos(tilt),    0,
-      0,     0,           0,         1,
+      [ 1,     0,           0,         0 ],
+      [ 0,   cos(tilt), -sin(tilt),    0 ],
+      [ 0,   sin(tilt),  cos(tilt),    0 ],
+      [ 0,     0,           0,         1 ],
     ])
-    findcurve_basis = augmatmultiply(dPQplane_basis, tilt_basis)
+    findcurve_basis = matmatmultiply(dPQplane_basis, tilt_basis)
     findcurve_into = np.linalg.inv(findcurve_basis)
 
-    q_findcurve = unaugment(findcurve_into, augment(q))
-    dq_findcurve = unaugment(findcurve_into, augment0(dq))
+    q_findcurve = augmatmultiply(findcurve_into, q)
+    dq_findcurve = augmatmultiply(findcurve_into, dq, augwith=0)
 
-    findcurve_target = np.concatenate(q_findcurve, dq_findcurve)
+    findcurve_target = np.hstack((q_findcurve, dq_findcurve))
     findcurve_start = (sqrt(start_s), sqrt(start_t), start_la,
                        start_mu, start_gamma, start_kappa)
     
     findcurve_epsilon = dist_pq_plane * 0.01
 
+    global findcurve_subproc
     if findcurve_subproc is None:
+      dbg('STARTING FINDCURVE')
       findcurve_subproc = subprocess.Popen(
         ['./findcurve'],
         bufsize=1,
@@ -155,7 +165,7 @@ class HelixishCurve():
                                  findcurve_start,
                                  [findcurve_epsilon]))
     dbg('RUNNING FINDCURVE', *findcurve_input)
-    print(findcurve_subproc.stdin, *findcurve_input)
+    print(*findcurve_input, file=findcurve_subproc.stdin)
     findcurve_subproc.stdin.flush()
 
     while True: