Custom implementation of reflex3s (fft, working, with plotting)

swash/processing/post_reflex.py

@@ -0,0 +1,68 @@
+import argparse
+import configparser
+import logging
+import pathlib
+
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy.signal as sgl
+
+from .reflex3s import reflex3s
+from .pa.reflex3S import reflex3S
+from .pa.PUV import PUV_dam
+
+parser = argparse.ArgumentParser(description="Post-process swash output")
+parser.add_argument("-v", "--verbose", action="count", default=0)
+parser.add_argument("-c", "--config", default="config.ini")
+args = parser.parse_args()
+
+logging.basicConfig(level=max((10, 20 - 10 * args.verbose)))
+log = logging.getLogger("post")
+
+log.info("Starting post-processing")
+config = configparser.ConfigParser()
+config.read(args.config)
+
+inp = pathlib.Path(config.get("post", "inp"))
+root = pathlib.Path(config.get("swash", "out"))
+
+log.info(f"Reading data from '{inp}'")
+x = np.load(inp.joinpath("xp.npy"))
+t = np.load(inp.joinpath("tsec.npy"))
+
+botl = np.load(inp.joinpath("botl.npy"))
+watl = np.load(inp.joinpath("watl.npy"))
+vel = np.load(inp.joinpath("vel.npy"))
+
+x0 = config.getint("post", "x0")
+arg_x0 = np.abs(x - x0).argmin()
+t0 = config.getfloat("post", "t0")
+arg_t0 = np.abs(t - t0).argmin()
+dt = np.diff(t).mean()
+f = 1 / dt
+
+idx = [arg_x0 - 5, arg_x0, arg_x0 + 7]
+wl = watl[arg_t0:, idx]
+
+res = reflex3S(*wl.T, *x[idx], botl[idx].mean(), f, 0.02, 0.2)
+f_, ai_, ar_ = reflex3s(wl.T, x[idx], botl[idx].mean(), f)
+ai, ar, _, _, _, _, fre = res
+
+out = pathlib.Path(config.get("post", "out")).joinpath(f"t{t0}x{x0}")
+log.info(f"Saving plots in '{out}'")
+out.mkdir(parents=True, exist_ok=True)
+
+plt.plot(fre, np.asarray(ar) / np.asarray(ai))
+plt.xlim((0, 0.3))
+plt.ylim((0, 1))
+plt.savefig(out.joinpath("reflex3s_pa.pdf"))
+plt.clf()
+plt.plot(f_, ar_ / ai_)
+plt.xlim((0, 0.3))
+plt.ylim((0, 1))
+plt.savefig(out.joinpath("reflex3s.pdf"))
+
+# pressk = np.load(inp.joinpath("pressk.npy"))
+# press = pressk[:, 1, :]
+# res = PUV_dam(vel[arg_t0:, 0, 500], press[arg_t0:, 500], botl[500], f, botl[500]/2, botl[500]/2)
+# print(res)

swash/processing/reflex3s.py

@@ -5,43 +5,41 @@ 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])
+    f_psd = fft.rfftfreq(eta.shape[1], 1 / f)
 
-    eta_amp = np.abs(eta_psd) / eta_psd.shape[1]
+    eta_amp = np.abs(eta_psd) / (f_psd.size - 1)
     eta_phase = np.angle(eta_psd)
 
     g = 9.81
-    k = np.asarray([opti.root_scalar(
+    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,
+                fprime=lambda k: np.tanh(k) + k * (1 - np.tanh(k) ** 2),
-        x0=0.2,
+                x0=0.5,
-    ).root for f in f_psd])
+            ).root
-
+            / h
-    dx = np.roll(x, 1) - x
+            for f in f_psd
-    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
+    s1 = 1 + np.exp(2j * k * (x[1] - x[0])) + np.exp(2j * k * (x[2] - x[0]))
+    s2 = 1 + np.exp(-2j * k * (x[1] - x[0])) + np.exp(-2j * k * (x[2] - x[0]))
+    s12 = 3
+    s3 = (
+        eta_amp[0] * np.exp(-1j * (eta_phase[0]))
+        + eta_amp[1] * np.exp(-1j * (eta_phase[1] - k * (x[1] - x[0])))
+        + eta_amp[2] * np.exp(-1j * (eta_phase[2] - k * (x[2] - x[0])))
+    )
+    s4 = (
+        eta_amp[0] * np.exp(-1j * (eta_phase[0]))
+        + eta_amp[1] * np.exp(-1j * (eta_phase[1] + k * (x[1] - x[0])))
+        + eta_amp[2] * np.exp(-1j * (eta_phase[2] + k * (x[2] - x[0])))
+    )
+    s5 = s1 * s2 - s12**2
+
+    ai = np.abs((s2 * s3 - s12 * s4) / s5)
+    ar = np.abs((s1 * s4 - s12 * s3) / s5)
+
+    return f_psd, ai, ar