From: Ian Jackson Date: Sun, 26 Nov 2017 10:41:29 +0000 (+0000) Subject: symbolic.py: turn into a module X-Git-Url: https://www.chiark.greenend.org.uk/ucgi/~ian/git?a=commitdiff_plain;h=ea18c8a3951b645f08c080cade9ab292ffb6c859;p=moebius3.git 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 --- diff --git a/Makefile b/Makefile index 11de207..5794c98 100644 --- 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 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 a42f3a6..b77d8cf 100644 --- 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 new mode 100644 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()