From 0cb87c50328356ddc367b9293707beb3c9a68fa2 Mon Sep 17 00:00:00 2001
From: "Edgar P. Burkhart" <git@edgarpierre.fr>
Date: Fri, 1 Apr 2022 13:07:37 +0200
Subject: [PATCH] Custom implementation of reflex3s (not working)

---
 swash/processing/reflex3s.py | 47 ++++++++++++++++++++++++++++++++++++
 1 file changed, 47 insertions(+)
 create mode 100644 swash/processing/reflex3s.py

diff --git a/swash/processing/reflex3s.py b/swash/processing/reflex3s.py
new file mode 100644
index 0000000..9b61030
--- /dev/null
+++ b/swash/processing/reflex3s.py
@@ -0,0 +1,47 @@
+import numpy as np
+from scipy import fft
+from scipy import signal as sgl
+from scipy import optimize as opti
+
+
+def reflex3s(eta, x, h, f):
+    #_welch = [sgl.welch(z, f) for z in eta]
+
+    #eta_psd = np.stack([_w[1] for _w in _welch])
+    #f_psd = _welch[0][0]
+    eta_psd = np.stack([fft.rfft(z) for z in eta])
+    f_psd = fft.rfftfreq(eta.shape[1])
+
+    eta_amp = np.abs(eta_psd) / eta_psd.shape[1]
+    eta_phase = np.angle(eta_psd)
+
+    g = 9.81
+    k = np.asarray([opti.root_scalar(
+        f=lambda k: k * np.tanh(k) - (2 * np.pi * f) ** 2 / g * h,
+        fprime=lambda k: np.tanh(k) + k * (1 - np.tanh(k)) ** 2,
+        x0=0.2,
+    ).root for f in f_psd])
+
+    dx = np.roll(x, 1) - x
+    dphi = np.roll(eta_phase, 1, axis=0) - eta_phase
+
+    ai = np.sqrt(
+        eta_amp**2
+        + np.roll(eta_amp, 1, axis=0)**2
+        - 2
+        * eta_amp
+        * np.roll(eta_amp, 1, axis=0)
+        * np.cos(dphi - k * dx[:, None])
+        / (2 * np.abs(np.sin(k * dx[:, None])))
+    )
+    ar = np.sqrt(
+        eta_amp**2
+        + np.roll(eta_amp, 1, axis=0)**2
+        - 2
+        * eta_amp
+        * np.roll(eta_amp, 1, axis=0)
+        * np.cos(dphi + k * dx[:, None])
+        / (2 * np.abs(np.sin(k * dx[:, None])))
+    )
+
+    return f_psd, ar / ai