From b85ee9d64a536937912544c7bbd5b98b635b7e8d Mon Sep 17 00:00:00 2001
From: Christian C <cc@localhost>
Date: Mon, 11 Nov 2024 12:29:32 -0800
Subject: Initial commit

---
 code/sunlab/common/mathlib/base.py | 57 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 57 insertions(+)
 create mode 100644 code/sunlab/common/mathlib/base.py

(limited to 'code/sunlab/common/mathlib/base.py')

diff --git a/code/sunlab/common/mathlib/base.py b/code/sunlab/common/mathlib/base.py
new file mode 100644
index 0000000..38ab14c
--- /dev/null
+++ b/code/sunlab/common/mathlib/base.py
@@ -0,0 +1,57 @@
+import numpy as np
+
+
+def angle(a, b):
+    """# Get Angle Between Row Vectors"""
+    from numpy import arctan2, pi
+
+    theta_a = arctan2(a[:, 1], a[:, 0])
+    theta_b = arctan2(b[:, 1], b[:, 0])
+    d_theta = theta_b - theta_a
+    assert (-pi <= d_theta) and (d_theta <= pi), "Theta difference outside of [-π,π]"
+    return d_theta
+
+
+def normalize(column):
+    """# Normalize Column Vector"""
+    from numpy.linalg import norm
+
+    return column / norm(column, axis=0)
+
+
+def winding(xy_grid, trajectory_start, trajectory_end):
+    """# Get Winding Number on Grid"""
+    from numpy import zeros, cross, clip, arcsin
+
+    trajectories = trajectory_end - trajectory_start
+    winding = zeros((xy_grid.shape[0]))
+    for idx, trajectory in enumerate(trajectories):
+        r = xy_grid - trajectory_start[idx]
+        cross = cross(normalize(trajectory), normalize(r))
+        cross = clip(cross, -1, 1)
+        theta = arcsin(cross)
+        winding += theta
+    return winding
+
+
+def vorticity(xy_grid, trajectory_start, trajectory_end):
+    """# Get Vorticity Number on Grid"""
+    from numpy import zeros, cross
+
+    trajectories = trajectory_end - trajectory_start
+    vorticity = zeros((xy_grid.shape[0]))
+    for idx, trajectory in enumerate(trajectories):
+        r = xy_grid - trajectory_start[idx]
+        vorticity += cross(normalize(trajectory), normalize(r))
+    return vorticity
+
+
+def data_range(data):
+    """# Get the range of values for each row"""
+    from numpy import min, max
+
+    return min(data, axis=0), max(data, axis=0)
+
+
+np.normalize = normalize
+np.range = data_range
-- 
cgit v1.2.1