add misc files

misc/eigenvalues.py

@@ -0,0 +1,40 @@
+#! /usr/bin/env python3
+import numpy as np
+
+A = np.array([[0, 0, 0], [0, 0, -1], [0, -1, 0]])
+B = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]])
+
+
+def similarity_transform(matrix):
+    L, S = np.linalg.eig(matrix)
+    L = np.diag(L)
+    S = S.transpose()
+    assert np.allclose(np.matmul(S.transpose(), np.matmul(L, S)), matrix)
+    return L, S
+
+
+def plusminus(matrix):
+    L, S = similarity_transform(matrix)
+
+    def signed(op):
+        return 0.5 * np.matmul(S.transpose(), np.matmul(op(L, np.abs(L)), S))
+
+    plus = signed(np.add)
+    minus = signed(np.subtract)
+
+    assert np.allclose(matrix, plus + minus)
+    return plus, minus
+
+
+Aplus, Aminus = plusminus(A)
+Bplus, Bminus = plusminus(B)
+
+print("A+")
+print(Aplus)
+print("A-")
+print(Aminus)
+print()
+print("B+")
+print(Bplus)
+print("B-")
+print(Bminus)

misc/eigenvalues_rotated.py

@@ -0,0 +1,55 @@
+#! /usr/bin/env python3
+import sympy as sp
+
+
+A = sp.Matrix([[0, 0, 0], [0, 0, -1], [0, -1, 0]])
+B = sp.Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 0]])
+
+kx, ky = sp.symbols("kx ky")
+
+A = kx * A + ky * B
+
+print("flux")
+sp.pprint(A)
+
+eigenvalues = A.eigenvals()
+
+# sp.pprint(eigenvalues)
+
+eigenvectors = A.eigenvects()
+
+# sp.pprint(eigenvectors)
+
+S, L = A.diagonalize()
+Labs = (L*L)**0.5
+
+plus = S*(L + Labs)*S.inv()/2
+minus = S*(L - Labs)*S.inv()/2
+
+Ahat = plus + minus
+
+
+plus.simplify()
+minus.simplify()
+
+print("Positive flux")
+r = sp.symbols("r")
+sp.pprint(plus.subs((kx*kx + ky*ky)**0.5, r))
+print("Negative flux")
+sp.pprint(minus.subs((kx*kx + ky*ky)**0.5, r))
+print("r = ")
+sp.pprint((kx*kx + ky*ky)**0.5)
+
+
+print("Evaluations:")
+print("plus (1, 0)")
+sp.pprint(plus.subs(kx, 1).subs(ky, 0))
+print("plus (0, 1)")
+sp.pprint(plus.subs(kx, 0).subs(ky, 1))
+print("plus (1/sqrt(2), 1/sqrt(2)")
+sp.pprint(plus.subs(kx, 1/2**0.5).subs(ky, 1/2**0.5))
+
+print("positive flux - negative flux")
+S = plus - minus
+S.simplify
+sp.pprint(S.subs((kx*kx + ky*ky)**0.5, r))

misc/euler.py

@@ -0,0 +1,33 @@
+import numpy as np
+
+rstar = 0.5  # Fixed radius
+eps = 1.0
+M = 0.1  # Indirectly sets p_inf
+gamma = 1.4  # Solid
+
+p_inf = 1.0 / (gamma * M * M)
+print(f"p_inf: {p_inf}")
+
+dx = 10000.0
+dy = 100.0
+f = (1 - (dx*dx + dy*dy))/(rstar*rstar)
+print(f"f: {f}")
+
+# print(eps*dy/(2*np.pi*np.sqrt(p_inf)*rstar * rstar) * np.exp(f / 2))
+
+u = 1.0 - eps*dy/(2*np.pi*np.sqrt(p_inf)*rstar * rstar) * np.exp(f / 2)
+v = 0.0 + eps*dx/(2*np.pi*np.sqrt(p_inf)*rstar * rstar) * np.exp(f / 2)
+
+print(f"sub p: {eps*eps*(gamma - 1)*M*M / (8*np.pi*np.pi*p_inf*rstar*rstar)*np.exp(f)}")
+
+
+rho = np.power(1.0 - eps*eps*(gamma - 1)*M*M / (
+    8*np.pi*np.pi*p_inf*rstar*rstar)*np.exp(f), 1.0/(gamma - 1))
+p = (rho**gamma)*p_inf
+print(f"p: {p}")
+e = p / (gamma - 1) + rho*(u**2 + v**2) / 2
+
+print(f"rho: {v}")
+print(f"u: {u}")
+print(f"v: {v}")
+print(f"e: {e}")