2 from __future__ import print_function
4 from sympy.vector.vector import *
7 from sympy.utilities.lambdify import lambdify, implemented_function
10 from moedebug import *
13 if not dbg_enabled(): return
17 if not dbg_enabled(): return
19 print('\n ' + vn + '\n')
24 N = CoordSysCartesian('N')
28 def vector_symbols(vnames):
30 for vname in vnames.split(' '):
32 for cname in 'i j k'.split(' '):
33 v += getattr(N, cname) * symbols(vname + '_' + cname)
37 A, B, C, D = vector_symbols('A B C D')
47 cost_ABCD = (QR & QR) / (BC & BC)
54 def ourccode(*a, **kw):
55 return ccode(*a, user_functions={'sinc':'sinc'}, **kw)
61 for l in s.split('\n'):
64 def cse_prep_cprint(v, tmp_prefix):
65 # => v, but also having cprint'd the common subexpression assignments
66 sym_iter = map((lambda i: symbols('%s%d' % (tmp_prefix,i))),
68 (defs, vs) = cse(v, symbols=sym_iter)
69 for defname, defval in defs:
70 cprint('double '+ourccode(defval, assign_to=defname))
73 def cassign(v, assign_to, tmp_prefix):
74 v = cse_prep_cprint(v, tmp_prefix)
75 cprint(ourccode(v, assign_to=assign_to))
77 def gen_diff(current, smalls):
80 j = zeros(len(params),0)
84 d = diff(current, paramv)
88 j = cse_prep_cprint(j, 'jtmp')
89 for ix in range(0, j.cols):
90 cprint(ourccode(j.col(ix), 'J_COL'))
91 cprint('J_END_COL(%d)' % ix)
95 cprint('if (!IS_SMALL(' + ourccode(small) + ')) {')
96 gen_diff(current, smalls)
97 cprint('} else { /* %s small */' % small)
98 gen_diff(current.replace(
100 1 - small*small/factorial(3) - small**4/factorial(5),
104 cprint('} /* %s small */' % small)
107 cprintraw('// AUTOGENERATED - DO NOT EDIT\n')
112 def gen_f_populate():
113 cprint('#define F_POPULATE')
114 cassign(cost_ABCD,'F','ftmp')
117 def gen_j_populate():
118 cprint('#define J_POPULATE')
119 gen_diff(result_dirnscaled, (sh*sh*la, th*th*la))
129 # https://github.com/sympy/sympy/issues/13642
130 # "lambdify sinc gives wrong answer!"
132 sinc_fixed = Function('sinc_fixed')
133 implemented_function(sinc_fixed, lambda x: np.sinc(x/np.pi))
134 out = out.subs(sinc,sinc_fixed)
135 p = list(map(eval,params))
136 return lambdify(p, out)