8 if sys.stdout.encoding is None:
9 sys.stdout = codecs.open("/dev/stdout", "w", 'utf-8')
11 init_printing(use_unicode=True)
13 r, theta, s, la, mu = symbols('r theta s lambda mu')
15 # start original formulation
24 p_rightvars = p_start.subs( theta, s/r ).subs( r, 1/la )
28 print('\n ' + vn + '\n')
31 # pprint(cse(eval(vn)))
35 p_dirn_rightvars = diff(p_rightvars, s)
37 dbg('p_dirn_rightvars')
41 p_nosing = (p_rightvars
42 .replace( 1-cos(zeta) , 2*sin(zeta/2)**2 )
43 .replace( sin(zeta)**2 , zeta*sinc(zeta)*sin(zeta) )
45 p_nosing[1] = (p_nosing[1]
46 .replace( sin(zeta) , zeta * sinc(zeta) )
53 q_owncoords = p_nosing.replace(s,t).replace(la,-la)
54 q_dirn_owncoords = p_dirn_rightvars.replace(s,t).replace(la,-la)
56 dbg('q_owncoords','q_dirn_owncoords')
57 dbg('q_owncoords.replace(t,0)','q_dirn_owncoords.replace(t,0)')
59 p2q_translate = p_nosing
60 #p2q_rotate_2d = Matrix([ p_dirn_rightvars[0:2],
63 #p2q_rotate[0:2, 0] = Matrix([ p_dirn_rightvars[1], -p_dirn_rightvars[0] ])
64 #p2q_rotate[0:2, 1] = p_dirn_rightvars[0:2]
66 p2q_rotate = Matrix([[ cos(theta), sin(theta), 0 ],
67 [ -sin(theta), cos(theta), 0 ],
68 [ 0 , 0, 1 ]]).subs(theta,la*s)
69 #p2q_rotate.add_col([0,0])
70 #p2q_rotate.add_row([0,0,1])
74 q_dirn_maincoords = p2q_rotate * q_dirn_owncoords;
75 q_maincoords = p2q_rotate * q_owncoords + p2q_translate
77 dbg('diff(p_dirn_rightvars,s)')
78 dbg('diff(q_dirn_maincoords,t)')
79 dbg('diff(q_dirn_maincoords,t).replace(t,0)')
81 assert(Eq(p2q_rotate * Matrix([0,1,mu]), p_dirn_rightvars))
83 #for v in 's','t','la','mu':
84 # dbg('diff(q_maincoords,%s)' % v)
86 #print('\n eye3 subs etc.\n')
87 #dbg('''Eq(eye(3) * Matrix([1,0,mu]),
88 # p_dirn_rightvars .cross(Matrix([0,0,1]) .subs(s,0)))''')
90 #dbg('''Eq(p2q_rotate * Matrix([1,0,mu]),
91 # p_dirn_rightvars .cross(Matrix([0,0,1])))''')
93 #eq = Eq(qmat * q_dirn_owncoords_0, p_dirn_rightvars)
98 dbg('q_maincoords.replace(t,0)','q_dirn_maincoords.replace(t,0)')
100 dbg('q_maincoords','q_dirn_maincoords')