chiark / gitweb /
aa42d0766f58a9026420a4d5ebab1b63227b646c
[moebius3.git] / symbolic.py
1 #!/usr/bin/python
2
3 from sympy import *
4
5 import sys
6
7 import sys, codecs
8 if sys.stdout.encoding is None:
9   sys.stdout = codecs.open("/dev/stdout", "w", 'utf-8')
10
11 init_printing(use_unicode=True)
12
13 r, theta, s, la, mu = symbols('r theta s lambda mu')
14
15 # start      original formulation
16 # rightvars  replaces 
17
18 p_start = Matrix([
19   r * (1 - cos(theta)),
20   r * sin(theta),
21   mu * s,
22 ])
23
24 p_rightvars = p_start.subs( theta, s/r ).subs( r, 1/la )
25
26 def dbg(*args):
27   for vn in args:
28     print('\n' + vn + '\n')
29     pprint(eval(vn))
30
31 dbg('p_rightvars')
32
33 p_dirn_rightvars = diff(p_rightvars, s)
34
35 dbg('p_dirn_rightvars')
36
37 zeta = Wild('zeta')
38
39 p_nosing = (p_rightvars
40             .replace( 1-cos(zeta)  ,   2*sin(zeta/2)**2          )
41             .replace( sin(zeta)**2 ,   zeta*sinc(zeta)*sin(zeta) )
42             )
43 p_nosing[1] = (p_nosing[1]
44             .replace( sin(zeta) , zeta * sinc(zeta)                )
45                )
46
47 dbg('p_nosing')
48
49 t = symbols('t')
50
51 q_owncoords = p_nosing.replace(s,t).replace(la,-la)
52 q_dirn_owncoords = p_dirn_rightvars.replace(s,t).replace(la,-la)
53
54 q_dirn_owncoords_0 = q_dirn_owncoords.replace(t,0)
55
56 dbg('q_owncoords','q_dirn_owncoords','q_dirn_owncoords_0')
57
58 p2q_translate = p_nosing
59 #p2q_rotate_2d = Matrix([ p_dirn_rightvars[0:2],
60                          
61 p2q_rotate = eye(3)
62 p2q_rotate[0:2, 0] = Matrix([ -p_dirn_rightvars[1], p_dirn_rightvars[0] ])
63 p2q_rotate[0:2, 1] = p_dirn_rightvars[0:2]
64 #p2q_rotate.add_col([0,0])
65 #p2q_rotate.add_row([0,0,1])
66
67 dbg('p2q_rotate')
68
69 assert(Eq(p2q_rotate * Matrix([0,1,mu]), p_dirn_rightvars))
70
71 print('\n eye3 subs etc.\n')
72 dbg('''Eq(eye(3) * Matrix([1,0,mu]),
73      p_dirn_rightvars .cross(Matrix([0,0,1]) .subs(s,0)))''')
74
75 dbg('''Eq(p2q_rotate * Matrix([1,0,mu]),
76           p_dirn_rightvars .cross(Matrix([0,0,1])))''')
77
78 #eq = Eq(qmat * q_dirn_owncoords_0, p_dirn_rightvars)
79 #print
80 #pprint(eq)
81 #solve(eq, Q)