Add euler symbolic algebra helpers

euler/misc/eigenvalues_euler.py

@@ -0,0 +1,47 @@
+#! /usr/bin/env python3
+import sympy as sp
+
+sp.init_printing(use_unicode=True, use_latex=True)
+
+rho = sp.symbols("rho", real=True, positive=True, nonzero=True)
+ru, rv, e = sp.symbols("rhou rhov e", real=True)
+gamma = sp.symbols("gamma", real=True, positive=True, nonzero=True)
+p = (gamma - 1) * (e - (ru ** 2 + rv ** 2) / (2 * rho))
+c = sp.sqrt(gamma * p / rho)
+
+u = ru / rho
+v = rv / rho
+
+E = sp.Matrix([rho * u, rho * u * u + p, rho * u * v, u * (e + p)])
+F = sp.Matrix([rho * v, rho * u * v, rho * v * v + p, v * (e + p)])
+
+A = E.jacobian([rho, ru, rv, e])
+B = F.jacobian([rho, ru, rv, e])
+
+print()
+for key in A.eigenvals():
+    sp.pprint(key)
+print()
+for key in B.eigenvals():
+    sp.pprint(key)
+
+# Eigenvalues are u, u +- c and v, v +- c
+# sp.pprint((key - c).simplify())  # Remainder after subtracting c
+
+# kx, ky = sp.symbols("kx ky")
+# Arot = kx*A + ky*B
+
+S, L = A.diagonalize()
+
+sp.pprint(S)
+sp.pprint(L)
+
+Labs = (L * L) ** 0.5
+
+plus = S * (L + Labs) * S.inv() / 2
+sp.pprint(plus)
+minus = S * (L - Labs) * S.inv() / 2
+sp.pprint(minus)
+
+S = plus - minus
+sp.pprint(S)

euler/misc/euler_vortex_evaluate.py

@@ -1,3 +1,4 @@
+#! /usr/bin/env python3
 import numpy as np
 
 rstar = 0.5  # Fixed radius