symbolic.py: wip differentiation

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

diff --git a/symbolic.py b/symbolic.py
index 6cc5c2b00b28b68dc648cf57a8979f94c4ee7c11..d0bcdc24f41393d66b1168c7062c8e2384728070 100755 (executable)
--- a/symbolic.py
+++ b/symbolic.py
@@ -8,7 +8,7 @@ import sys, codecs
 if sys.stdout.encoding is None:
   sys.stdout = codecs.open("/dev/stdout", "w", 'utf-8')
 
-init_printing(use_unicode=True)
+init_printing(use_unicode=True, num_columns=280)
 
 r, theta, s, la, mu, kappa = symbols('r theta s lambda mu kappa')
 
@@ -27,8 +27,8 @@ def dbg(*args):
   for vn in args:
     print('\n    ' + vn + '\n')
     pprint(eval(vn))
-#    print('\n          =\n')
-#    pprint(cse(eval(vn)))
+    print('\n          =\n')
+    pprint(cse(eval(vn)))
 
 dbg('p_rightvars')
 
@@ -125,3 +125,44 @@ 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')
+
+q_dirn_orgcoords = o2p_rotate2 * o2p_rotate1 * q_dirn_maincoords;
+q_orgcoords = o2p_rotate2 * o2p_rotate1 * q_maincoords;
+dbg('q_orgcoords','q_dirn_orgcoords')
+
+sh, th = symbols('alpha beta')
+
+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)
+
+print '----------------------------------------'
+dbg('q_sqparm', 'q_dirn_sqparm')
+print '----------------------------------------'
+for v in 'sh','th','la','mu':
+  dbg('diff(q_sqparm,%s)' % v)
+  dbg('diff(q_dirn_sqparm,%s)' % v)
+print '----------------------------------------'
+
+gamma = symbols('gamma')
+
+q_dirn_dirnscaled = q_dirn_sqparm * gamma
+
+result_dirnscaled = q_sqparm.col_join(q_dirn_dirnscaled)
+dbg('result_dirnscaled')
+
+for sm_sh in ((), (sh*sh*la,)):
+  for sm_th in ((), (th*th*la,)):
+    smalls = sm_sh + sm_th
+    result_fordiff = result_dirnscaled
+    for small in smalls:
+      result_fordiff = result_fordiff.replace(
+        sinc(small),
+        1 - small*small/factorial(3) - small**4/factorial(5)
+      )
+    print smalls
+    dbg('result_fordiff');
+
+#bad = q_orgcoords[0]
+#badd = diff(bad, la)
+
+#dbg('bad','badd')