initial comparison with katz & plotkin

data/debug2.py

@@ -0,0 +1,101 @@
+import numpy as np
+
+class DualQuaternion:
+    def __init__(self, real, dual):
+        self.real = real  # Quaternion
+        self.dual = dual  # Quaternion
+
+    @staticmethod
+    def from_translation(translation):
+        real = np.array([1, 0, 0, 0], dtype=translation.dtype)
+        translation = np.array(translation, dtype=translation.dtype)
+        dual = np.hstack(([0], 0.5 * translation))
+        return DualQuaternion(real, dual)
+
+    @staticmethod
+    def from_rotation(axis_origin, axis, angle):
+        axis = np.array(axis, dtype=angle.dtype)
+        axis /= np.linalg.norm(axis)
+        sin_half_angle = np.sin(angle / 2)
+        real = np.hstack(([np.cos(angle / 2)], sin_half_angle * axis)).astype(angle.dtype)
+
+        origin = np.array(axis_origin, dtype=angle.dtype)
+        dual_part = np.cross(-origin, sin_half_angle * axis)
+        dual = np.hstack(([0], dual_part)).astype(angle.dtype)
+
+        return DualQuaternion(real, dual)
+
+    def transform_point(self, point):
+        p = np.array([0, point[0], point[1], point[2]], dtype=self.real.dtype)
+        transformed_p = self._quaternion_mul(self._quaternion_mul(self.real, p), self._quaternion_conj(self.real)) + 2 * self._quaternion_mul(self.dual, self._quaternion_conj(self.real))
+        return transformed_p[1:4]
+
+    def __mul__(self, other):
+        real = self._quaternion_mul(self.real, other.real)
+        dual = self._quaternion_mul(self.real, other.dual) + self._quaternion_mul(self.dual, other.real)
+        return DualQuaternion(real, dual)
+
+    @staticmethod
+    def _quaternion_mul(q1, q2):
+        w1, x1, y1, z1 = q1
+        w2, x2, y2, z2 = q2
+        w = w1*w2 - x1*x2 - y1*y2 - z1*z2
+        x = w1*x2 + x1*w2 + y1*z2 - z1*y2
+        y = w1*y2 - x1*z2 + y1*w2 + z1*x2
+        z = w1*z2 + x1*y2 - y1*x2 + z1*w2
+        return np.array([w, x, y, z], dtype=q1.dtype)
+
+    @staticmethod
+    def _quaternion_conj(q):
+        return np.array([q[0], -q[1], -q[2], -q[3]], dtype=q.dtype)
+
+alpha = np.radians(5)  # Example angle in radians
+
+def freestream2(t):
+    return DualQuaternion.from_translation([-np.cos(alpha) * t, 0, -np.sin(alpha) * t])
+
+def pitching(t):
+    return DualQuaternion.from_rotation([0, 0, 0], [0, 1, 0], 0.25 * np.pi * np.sin(t))
+
+class Kinematics:
+    def __init__(self):
+        self.m_joints = []
+
+    def add_joint(self, joint):
+        self.m_joints.append(joint)
+
+    def displacement(self, t, n=0):
+        result = DualQuaternion(np.array([1, 0, 0, 0], dtype=t.dtype), np.array([0, 0, 0, 0], dtype=t.dtype))
+        end_joint = len(self.m_joints) if n == 0 else n
+        for i in range(end_joint):
+            result = result * self.m_joints[i](t)
+        return result
+
+    def relative_displacement(self, t0, t1, n=0):
+        displacement_t1 = self.displacement(t1, n)
+        displacement_t0 = self.displacement(t0, n)
+        displacement_t0_inv = DualQuaternion(displacement_t0.real, -displacement_t0.dual)
+        return displacement_t1 * displacement_t0_inv
+
+    def velocity(self, t, vertex, n=0):
+        EPS = np.sqrt(np.finfo(np.float32).eps)
+        complex_t = t + 1j * EPS
+        displacement_with_complex_step = self.relative_displacement(t, complex_t, n)
+        vertex_complex = np.array(vertex, dtype=np.complex64)
+        transformed_vertex = displacement_with_complex_step.transform_point(vertex_complex)
+        return transformed_vertex.imag / EPS
+
+    def velocity_magnitude(self, t, vertex):
+        return np.linalg.norm(self.velocity(t, vertex))
+
+# Example usage
+kinematics = Kinematics()
+kinematics.add_joint(freestream2)
+kinematics.add_joint(pitching)
+
+vertex = np.array([4.0, 0, 0], dtype=np.float32)  # Vertex without homogeneous coordinate
+t = 0.0
+velocity = kinematics.velocity(t, vertex)
+velocity_magnitude = kinematics.velocity_magnitude(t, vertex)
+print("Numerical vel: ", velocity)
+print("Velocity Magnitude:", velocity_magnitude)

data/katz/Makefile

@@ -0,0 +1,12 @@
+FC = gfortran
+executables := $(patsubst %.f,%.out,$(wildcard *.f))
+
+
+all:
+	make $(executables)
+
+%.out: %.f
+	$(FC) -O2 $*.f -o $*.out
+
+clean:
+	-rm *.out

data/katz/katz_cl.txt

@@ -0,0 +1,160 @@
+      0.00      -3.00
+      0.10      -0.56
+      0.20      -0.54
+      0.30      -0.53
+      0.40      -0.51
+      0.50      -0.49
+      0.60      -0.47
+      0.70      -0.45
+      0.80      -0.42
+      0.90      -0.40
+      1.00      -0.37
+      1.10      -0.34
+      1.20      -0.31
+      1.30      -0.28
+      1.40      -0.25
+      1.50      -0.22
+      1.60      -0.19
+      1.70      -0.16
+      1.80      -0.13
+      1.90      -0.11
+      2.00      -0.08
+      2.10      -0.06
+      2.20      -0.04
+      2.30      -0.02
+      2.40      -0.00
+      2.50       0.01
+      2.60       0.02
+      2.70       0.02
+      2.80       0.03
+      2.90       0.02
+      3.00       0.02
+      3.10       0.01
+      3.20       0.00
+      3.30      -0.01
+      3.40      -0.03
+      3.50      -0.05
+      3.60      -0.07
+      3.70      -0.09
+      3.80      -0.12
+      3.90      -0.15
+      4.00      -0.18
+      4.10      -0.21
+      4.20      -0.24
+      4.30      -0.27
+      4.40      -0.31
+      4.50      -0.34
+      4.60      -0.37
+      4.70      -0.41
+      4.80      -0.44
+      4.90      -0.47
+      5.00      -0.50
+      5.10      -0.53
+      5.20      -0.55
+      5.30      -0.58
+      5.40      -0.60
+      5.50      -0.61
+      5.60      -0.63
+      5.70      -0.64
+      5.80      -0.65
+      5.90      -0.65
+      6.00      -0.65
+      6.10      -0.65
+      6.20      -0.65
+      6.30      -0.64
+      6.40      -0.63
+      6.50      -0.61
+      6.60      -0.60
+      6.70      -0.58
+      6.80      -0.55
+      6.90      -0.53
+      7.00      -0.50
+      7.10      -0.47
+      7.20      -0.44
+      7.30      -0.41
+      7.40      -0.38
+      7.50      -0.35
+      7.60      -0.31
+      7.70      -0.28
+      7.80      -0.25
+      7.90      -0.21
+      8.00      -0.18
+      8.10      -0.15
+      8.20      -0.13
+      8.30      -0.10
+      8.40      -0.08
+      8.50      -0.06
+      8.60      -0.04
+      8.70      -0.02
+      8.80      -0.01
+      8.90       0.00
+      9.00       0.01
+      9.10       0.01
+      9.20       0.01
+      9.30       0.01
+      9.40      -0.00
+      9.50      -0.01
+      9.60      -0.02
+      9.70      -0.04
+      9.80      -0.06
+      9.90      -0.08
+     10.00      -0.11
+     10.10      -0.13
+     10.20      -0.16
+     10.30      -0.19
+     10.40      -0.22
+     10.50      -0.25
+     10.60      -0.29
+     10.70      -0.32
+     10.80      -0.35
+     10.90      -0.39
+     11.00      -0.42
+     11.10      -0.45
+     11.20      -0.48
+     11.30      -0.51
+     11.40      -0.54
+     11.50      -0.56
+     11.60      -0.58
+     11.70      -0.60
+     11.80      -0.62
+     11.90      -0.63
+     12.00      -0.64
+     12.10      -0.65
+     12.20      -0.66
+     12.30      -0.66
+     12.40      -0.66
+     12.50      -0.65
+     12.60      -0.64
+     12.70      -0.63
+     12.80      -0.61
+     12.90      -0.60
+     13.00      -0.58
+     13.10      -0.55
+     13.20      -0.53
+     13.30      -0.50
+     13.40      -0.47
+     13.50      -0.44
+     13.60      -0.41
+     13.70      -0.38
+     13.80      -0.34
+     13.90      -0.31
+     14.00      -0.28
+     14.10      -0.24
+     14.20      -0.21
+     14.30      -0.18
+     14.40      -0.15
+     14.50      -0.12
+     14.60      -0.10
+     14.70      -0.07
+     14.80      -0.05
+     14.90      -0.03
+     15.00      -0.02
+     15.10      -0.01
+     15.20       0.00
+     15.30       0.01
+     15.40       0.01
+     15.50       0.01
+     15.60       0.00
+     15.70      -0.00
+     15.80      -0.01
+     15.90      -0.03