From ecea6c665a171162baeddd60caad00bc5d72391b Mon Sep 17 00:00:00 2001
From: "Edgar P. Burkhart" <git@edgarpierre.fr>
Date: Thu, 7 Apr 2022 14:47:42 +0200
Subject: [PATCH] Zero-crossings method to find large waves

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

diff --git a/swash/processing/zero_cross.py b/swash/processing/zero_cross.py
new file mode 100644
index 0000000..dcfb46a
--- /dev/null
+++ b/swash/processing/zero_cross.py
@@ -0,0 +1,60 @@
+import argparse
+import configparser
+import logging
+import pathlib
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+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("x.npy"))
+t = np.load(inp.joinpath("t.npy"))
+
+watl = np.load(inp.joinpath("watl.npy"))
+
+# Cospectral calculations
+x0 = config.getint("post", "x0")
+arg_x0 = np.abs(x - x0).argmin()
+
+w0 = watl[:, arg_x0]
+cr0 = np.where(np.diff(np.sign(w0)))[0]
+
+wave = np.fromiter(
+    (
+        np.abs(
+            np.max(np.abs(w0[cr0[i - 1] : cr0[i]]))
+            + np.max(np.abs(w0[cr0[i] : cr0[i + 1]]))
+        )
+        for i in range(1, len(cr0) - 1)
+    ),
+    dtype=np.single,
+)
+
+i0 = np.argmax(wave)
+
+out = pathlib.Path(config.get("post", "out")).joinpath(f"x{x0}")
+log.info(f"Saving plots in '{out}'")
+out.mkdir(parents=True, exist_ok=True)
+
+fig, ax = plt.subplots()
+ax.plot(t[cr0[1:-1]] * 1e-3, wave)
+fig.savefig(out.joinpath("wsize.pdf"))
+
+fig2, ax2 = plt.subplots()
+ax2.plot(t[cr0[i0 - 5] : cr0[i0 + 7]], w0[cr0[i0 - 5] : cr0[i0 + 7]])
+fig2.savefig(out.joinpath("maxw.pdf"))