symbolic.py: turn into a module

this will let us use the sympy lambda facility in the actual scad generator

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

diff --git a/Makefile b/Makefile
index 11de2077598ea661e027f6e875be4e52ff6d636b..5794c98e4e6a9b145f87187d0738aceef37e4563 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -21,8 +21,8 @@ moebius-mesh.scad: meshscad $(PYLIBS)
 
 AUTO_TOPLEVELS := $(foreach m,$(USING_AUTOS),$(shell $(PLAY)/toplevel-find $m))
 
-symbolic.c:    symbolic.py
-               ./symbolic.py -q >$@.tmp && mv -f $@.tmp $@
+symbolic.c:    gensymbolic symbolic.py
+               ./gensymbolic -q >$@.tmp && mv -f $@.tmp $@
 
 findcurve.o:   findcurve.c symbolic.c
 

diff --git a/gensymbolic b/gensymbolic
new file mode 100755 (executable)
index 0000000..e0f08ee
--- /dev/null
+++ b/gensymbolic
@@ -0,0 +1,30 @@
+#!/usr/bin/python3
+
+import sys, codecs
+from optparse import OptionParser
+
+import sympy
+
+import symbolic
+from moedebug import dbg_enable, dbg
+
+opt_parser = OptionParser()
+opt_parser.add_option('-q',dest='quiet',action='store_true',
+                      default=False,help='suppress diagnostic output')
+opt_parser.add_option('-7',dest='ascii',action='store_true',
+                      default=False,help='use 7-bit output only')
+opt_parser.add_option('-w',dest='width',action='store',type='int',
+                      default=80,help='width for printing')
+(options, args) = opt_parser.parse_args()
+assert(not len(args))
+
+if (not options.ascii) and sys.stdout.encoding is None:
+  sys.stdout = codecs.open("/dev/stdout", "w", 'utf-8')
+
+dbg_enable = not options.quiet
+if dbg_enable:
+  sympy.init_printing(use_unicode=(not options.ascii),
+                      num_columns=options.width)
+
+symbolic.calculate()
+symbolic.gen_C()

diff --git a/moedebug.py b/moedebug.py
index a42f3a66b432251edf384ceabedca4467d5f7e3d..b77d8cff6f889af701703da7a298782c5ea2a768 100644 (file)
--- a/moedebug.py
+++ b/moedebug.py
@@ -3,5 +3,5 @@ dbg_enable = False
 
 def dbg(*args):
   if dbg_enable:
-    print(file=sys.stderr, 'D ', *args)
+    print('D ', *args, file=sys.stderr)
 

diff --git a/symbolic.py b/symbolic.py
old mode 100755 (executable)
new mode 100644 (file)
index 9d93f3b..8607034
--- a/symbolic.py
+++ b/symbolic.py
@@ -1,170 +1,156 @@
-#!/usr/bin/python3
 
 from sympy import *
 import itertools
 
-import sys, codecs
-from optparse import OptionParser
-
-opt_parser = OptionParser()
-opt_parser.add_option('-q',dest='quiet',action='store_true',
-                      default=False,help='suppress diagnostic output')
-opt_parser.add_option('-7',dest='ascii',action='store_true',
-                      default=False,help='use 7-bit output only')
-opt_parser.add_option('-w',dest='width',action='store',type='int',
-                      default=80,help='width for printing')
-(options, args) = opt_parser.parse_args()
-assert(not len(args))
-
-if (not options.ascii) and sys.stdout.encoding is None:
-  sys.stdout = codecs.open("/dev/stdout", "w", 'utf-8')
-
-init_printing(use_unicode=(not options.ascii), num_columns=options.width)
+from moedebug import dbg_enable
 
 r, theta, s, la, mu, kappa = symbols('r theta s lambda mu kappa')
 
 # start      original formulation
 # rightvars  replaces 
 
-p_start = Matrix([
-  r * (1 - cos(theta)),
-  r * sin(theta),
-  mu * s,
-])
-
-p_rightvars = p_start.subs( theta, s/r ).subs( r, 1/la )
-
 def dprint(*args):
-  if options.quiet: return
+  if not dbg_enable: return
   print(*args)
 
 def dbg(*args):
-  if options.quiet: return
+  if not dbg_enable: return
   for vn in args:
     print('\n    ' + vn + '\n')
     pprint(eval(vn))
     print('\n          =\n')
     pprint(cse(eval(vn)))
 
-dbg('p_rightvars')
+def calculate():
+  p_start = Matrix([
+    r * (1 - cos(theta)),
+    r * sin(theta),
+    mu * s,
+  ])
 
-p_dirn_rightvars = diff(p_rightvars, s)
+  p_rightvars = p_start.subs( theta, s/r ).subs( r, 1/la )
 
-dbg('p_dirn_rightvars')
+  dbg('p_rightvars')
 
-zeta = Wild('zeta')
+  p_dirn_rightvars = diff(p_rightvars, s)
 
-p_nosing = (p_rightvars
-            .replace( 1-cos(zeta)  ,   2*sin(zeta/2)**2          )
-            .replace( sin(zeta)**2 ,   zeta*sinc(zeta)*sin(zeta) )
-            )
-p_nosing[1] = (p_nosing[1]
-            .replace( sin(zeta) , zeta * sinc(zeta)                )
-               )
+  dbg('p_dirn_rightvars')
 
-dbg('p_nosing')
+  zeta = Wild('zeta')
 
-t = symbols('t')
+  p_nosing = (p_rightvars
+              .replace( 1-cos(zeta)  ,   2*sin(zeta/2)**2          )
+              .replace( sin(zeta)**2 ,   zeta*sinc(zeta)*sin(zeta) )
+              )
+  p_nosing[1] = (p_nosing[1]
+              .replace( sin(zeta) , zeta * sinc(zeta)                )
+                 )
 
-q_owncoords = p_nosing.replace(s,t).replace(la,-la)
-q_dirn_owncoords = p_dirn_rightvars.replace(s,t).replace(la,-la)
+  dbg('p_nosing')
 
-dbg('q_owncoords','q_dirn_owncoords')
-dbg('q_owncoords.replace(t,0)','q_dirn_owncoords.replace(t,0)')
+  t = symbols('t')
 
-p2q_translate = p_nosing
-#p2q_rotate_2d = Matrix([ p_dirn_rightvars[0:2],
-                         
-#p2q_rotate = eye(3)
-#p2q_rotate[0:2, 0] = Matrix([ p_dirn_rightvars[1], -p_dirn_rightvars[0] ])
-#p2q_rotate[0:2, 1] = p_dirn_rightvars[0:2]
+  q_owncoords = p_nosing.replace(s,t).replace(la,-la)
+  q_dirn_owncoords = p_dirn_rightvars.replace(s,t).replace(la,-la)
 
-p2q_rotate = Matrix([[  cos(theta), sin(theta), 0 ],
-                     [ -sin(theta), cos(theta), 0 ],
-                     [  0         , 0,          1 ]]).subs(theta,la*s)
-#p2q_rotate.add_col([0,0])
-#p2q_rotate.add_row([0,0,1])
+  dbg('q_owncoords','q_dirn_owncoords')
+  dbg('q_owncoords.replace(t,0)','q_dirn_owncoords.replace(t,0)')
 
-dbg('p2q_rotate')
+  p2q_translate = p_nosing
+  #p2q_rotate_2d = Matrix([ p_dirn_rightvars[0:2],
 
-q_dirn_maincoords = p2q_rotate * q_dirn_owncoords;
-q_maincoords = p2q_rotate * q_owncoords + p2q_translate
+  #p2q_rotate = eye(3)
+  #p2q_rotate[0:2, 0] = Matrix([ p_dirn_rightvars[1], -p_dirn_rightvars[0] ])
+  #p2q_rotate[0:2, 1] = p_dirn_rightvars[0:2]
 
-dbg('diff(p_dirn_rightvars,s)')
-dbg('diff(q_dirn_maincoords,t)')
-dbg('diff(q_dirn_maincoords,t).replace(t,0)')
+  p2q_rotate = Matrix([[  cos(theta), sin(theta), 0 ],
+                       [ -sin(theta), cos(theta), 0 ],
+                       [  0         , 0,          1 ]]).subs(theta,la*s)
+  #p2q_rotate.add_col([0,0])
+  #p2q_rotate.add_row([0,0,1])
 
-assert(Eq(p2q_rotate * Matrix([0,1,mu]), p_dirn_rightvars))
+  dbg('p2q_rotate')
 
-#for v in 's','t','la','mu':
-#  dbg('diff(q_maincoords,%s)' % v)
+  q_dirn_maincoords = p2q_rotate * q_dirn_owncoords;
+  q_maincoords = p2q_rotate * q_owncoords + p2q_translate
 
-#print('\n eye3 subs etc.\n')
-#dbg('''Eq(eye(3) * Matrix([1,0,mu]),
-#     p_dirn_rightvars .cross(Matrix([0,0,1]) .subs(s,0)))''')
+  dbg('diff(p_dirn_rightvars,s)')
+  dbg('diff(q_dirn_maincoords,t)')
+  dbg('diff(q_dirn_maincoords,t).replace(t,0)')
 
-#dbg('''Eq(p2q_rotate * Matrix([1,0,mu]),
-#          p_dirn_rightvars .cross(Matrix([0,0,1])))''')
+  assert(Eq(p2q_rotate * Matrix([0,1,mu]), p_dirn_rightvars))
 
-#eq = Eq(qmat * q_dirn_owncoords_0, p_dirn_rightvars)
-#print
-#pprint(eq)
-#solve(eq, Q)
+  #for v in 's','t','la','mu':
+  #  dbg('diff(q_maincoords,%s)' % v)
 
-dbg('q_maincoords.replace(t,0)','q_dirn_maincoords.replace(t,0)')
+  #print('\n eye3 subs etc.\n')
+  #dbg('''Eq(eye(3) * Matrix([1,0,mu]),
+  #     p_dirn_rightvars .cross(Matrix([0,0,1]) .subs(s,0)))''')
 
-dbg('q_maincoords','q_dirn_maincoords')
+  #dbg('''Eq(p2q_rotate * Matrix([1,0,mu]),
+  #          p_dirn_rightvars .cross(Matrix([0,0,1])))''')
 
-sinof_mu = sin(atan(mu))
-cosof_mu = cos(atan(mu))
+  #eq = Eq(qmat * q_dirn_owncoords_0, p_dirn_rightvars)
+  #print
+  #pprint(eq)
+  #solve(eq, Q)
 
-dbg('cosof_mu','sinof_mu')
+  dbg('q_maincoords.replace(t,0)','q_dirn_maincoords.replace(t,0)')
 
-o2p_rotate1 = Matrix([[ 1,  0,         0        ],
-                      [ 0,  cosof_mu, +sinof_mu ],
-                      [ 0, -sinof_mu,  cosof_mu ]])
+  dbg('q_maincoords','q_dirn_maincoords')
 
-check_dirn_p_s0 = o2p_rotate1 * p_dirn_rightvars.replace(s,0)
-check_dirn_p_s0.simplify()
-dbg('check_dirn_p_s0')
+  sinof_mu = sin(atan(mu))
+  cosof_mu = cos(atan(mu))
 
-o2p_rotate2 = Matrix([[  cos(kappa), 0, -sin(kappa) ],
-                      [  0,          1,  0          ],
-                      [ +sin(kappa), 0,  cos(kappa) ]])
+  dbg('cosof_mu','sinof_mu')
 
-p_dirn_orgcoords = o2p_rotate2 * o2p_rotate1 * p_dirn_rightvars
+  o2p_rotate1 = Matrix([[ 1,  0,         0        ],
+                        [ 0,  cosof_mu, +sinof_mu ],
+                        [ 0, -sinof_mu,  cosof_mu ]])
 
-check_dirn_p_s0 = p_dirn_orgcoords.replace(s,0)
-check_dirn_p_s0.simplify()
-dbg('check_dirn_p_s0')
+  check_dirn_p_s0 = o2p_rotate1 * p_dirn_rightvars.replace(s,0)
+  check_dirn_p_s0.simplify()
+  dbg('check_dirn_p_s0')
 
-check_accel_p_s0 = diff(p_dirn_orgcoords,s).replace(s,0)
-check_accel_p_s0.simplify()
-dbg('check_accel_p_s0')
+  o2p_rotate2 = Matrix([[  cos(kappa), 0, -sin(kappa) ],
+                        [  0,          1,  0          ],
+                        [ +sin(kappa), 0,  cos(kappa) ]])
 
-q_dirn_orgcoords = o2p_rotate2 * o2p_rotate1 * q_dirn_maincoords;
-q_orgcoords = o2p_rotate2 * o2p_rotate1 * q_maincoords;
-dbg('q_orgcoords','q_dirn_orgcoords')
+  p_dirn_orgcoords = o2p_rotate2 * o2p_rotate1 * p_dirn_rightvars
 
-sh, th = symbols('alpha beta')
+  check_dirn_p_s0 = p_dirn_orgcoords.replace(s,0)
+  check_dirn_p_s0.simplify()
+  dbg('check_dirn_p_s0')
 
-q_dirn_sqparm = q_dirn_orgcoords.replace(s, sh**2).replace(t, th**2)
-q_sqparm      = q_orgcoords     .replace(s, sh**2).replace(t, th**2)
+  check_accel_p_s0 = diff(p_dirn_orgcoords,s).replace(s,0)
+  check_accel_p_s0.simplify()
+  dbg('check_accel_p_s0')
 
-dprint('----------------------------------------')
-dbg('q_sqparm', 'q_dirn_sqparm')
-dprint('----------------------------------------')
-for v in 'sh','th','la','mu':
-  dbg('diff(q_sqparm,%s)' % v)
-  dbg('diff(q_dirn_sqparm,%s)' % v)
-dprint('----------------------------------------')
+  q_dirn_orgcoords = o2p_rotate2 * o2p_rotate1 * q_dirn_maincoords;
+  q_orgcoords = o2p_rotate2 * o2p_rotate1 * q_maincoords;
+  dbg('q_orgcoords','q_dirn_orgcoords')
 
-gamma = symbols('gamma')
+  global sh, th
+  sh, th = symbols('alpha beta')
 
-q_dirn_dirnscaled = q_dirn_sqparm * gamma
+  q_dirn_sqparm = q_dirn_orgcoords.replace(s, sh**2).replace(t, th**2)
+  q_sqparm      = q_orgcoords     .replace(s, sh**2).replace(t, th**2)
 
-result_dirnscaled = q_sqparm.col_join(q_dirn_dirnscaled)
-dbg('result_dirnscaled')
+  dprint('----------------------------------------')
+  dbg('q_sqparm', 'q_dirn_sqparm')
+  dprint('----------------------------------------')
+  for v in 'sh','th','la','mu':
+    dbg('diff(q_sqparm,%s)' % v)
+    dbg('diff(q_dirn_sqparm,%s)' % v)
+  dprint('----------------------------------------')
+
+  gamma = symbols('gamma')
+
+  q_dirn_dirnscaled = q_dirn_sqparm * gamma
+
+  global result_dirnscaled
+  result_dirnscaled = q_sqparm.col_join(q_dirn_dirnscaled)
+  dbg('result_dirnscaled')
 
 params = ('sh','th','la','mu','gamma','kappa')
 
@@ -245,5 +231,3 @@ def gen_C():
   gen_x_extract()
   gen_f_populate()
   gen_j_populate()
-
-gen_C()