From b8da5508495755d42e0ed67ceafc242fda8cb71d Mon Sep 17 00:00:00 2001
From: "Edgar P. Burkhart" <git@edgarpierre.fr>
Date: Wed, 29 Jun 2022 09:31:07 +0200
Subject: [PATCH] Update zero_cross and add nansasena

---
 data/processing/nandasena.py  | 27 +++++++++++++++++++++++++++
 data/processing/zero_cross.py | 13 ++++++++++---
 2 files changed, 37 insertions(+), 3 deletions(-)
 create mode 100644 data/processing/nandasena.py

diff --git a/data/processing/nandasena.py b/data/processing/nandasena.py
new file mode 100644
index 0000000..298f20b
--- /dev/null
+++ b/data/processing/nandasena.py
@@ -0,0 +1,27 @@
+import numpy as np
+
+def u_nandasena(b, c, C_d, C_l, rho_s, rho_w=1e3, theta=0, g=9.81, **_):
+    return np.sqrt(
+        (2 * (rho_s/rho_w - 1) * g * c * np.cos(np.radians(theta)))
+        / C_l
+    )
+
+
+const = {
+    "a": 4,
+    "b": 2.5,
+    "c": 2,
+    "C_d": 1.95,
+    "C_l": 0.178,
+    "theta": 0,
+    "rho_s": 2.7e3,
+}
+
+u = u_nandasena(**const)
+u_ola = 14.5
+print(f"Nandasena: u={u:5.2f}m/s")
+print(f"Olaflow:   u={u_ola:5.2f}m/s")
+print(f"Gap:      du={abs(u_ola - u):5.2f}m/s")
+print(f"Gap:      du={abs(u_ola - u)/u:5.2%}")
+print(f"Gap:      du={abs(u_ola - u)/u_ola:5.2%}")
+
diff --git a/data/processing/zero_cross.py b/data/processing/zero_cross.py
index 8fd38a4..ce76b19 100644
--- a/data/processing/zero_cross.py
+++ b/data/processing/zero_cross.py
@@ -58,6 +58,7 @@ cr0 = np.where(np.diff(np.sign(z)))[0]
 wave = np.fromiter(
     (
         np.max(np.abs(z[cr0[i - 1] : cr0[i]])) + np.max(np.abs(z[cr0[i] : cr0[i + 1]]))
+        #(np.max(np.abs(raw_ts["z"][cr0[i - 1] : cr0[i]])) + np.max(np.abs(raw_ts["z"][cr0[i] : cr0[i + 1]]))) * 1e-2
         for i in range(1, len(cr0) - 1)
     ),
     dtype=np.single,
@@ -80,7 +81,7 @@ sj = s0 * 2 ** (j * dj)
 # sj = s0 * np.arange(1, 2 ** (J * dj))
 Mw = sj / dt
 Mlims = sj[[0, -1]]
-M = (np.abs(sgl.cwt(raw_ts["z"]*1e-2, sgl.morlet2, Mw))/np.var(raw_ts["z"]*1e-2))**2
+M = (np.abs(sgl.cwt(raw_ts["z"]*1e-2, sgl.morlet2, Mw))/np.std(raw_ts["z"]*1e-2))**2
 # M = np.abs(sgl.cwt(z, sgl.morlet, Mw))
 v = np.max(np.abs(M))
 
@@ -102,7 +103,8 @@ c = ax.contourf(t, sj, M, cmap="Greys", vmin=0, vmax=v)
 fig.colorbar(c)
 
 H13 = np.quantile(wave, 2 / 3)
-th = 10
+Hs = 4*np.std(raw_ts["z"]*1e-2)
+th = 2 * Hs
 log.info(f"Threshold: {th}m")
 bigw = np.where(wave > th)[0]
 ym = 1.1 * np.max(np.abs(z))
@@ -110,8 +112,11 @@ nw = wave.size / 2
 nlw = bigw.size
 log.info(f"Number of waves: {nw}")
 log.info(f"Number of waves >m: {nlw}")
-log.info(f"Proportion: {nlw/nw:e}")
+log.info(f"Frequency: {nlw/nw:e}")
+log.info(f"Frequency: 1/{nw/nlw:.0f}")
 log.info(f"H1/3: {H13}m")
+log.info(f"Hs: {Hs}")
+
 if bigw.size > 32:
     log.warning(f"Number of large waves: {bigw.size}")
     sys.exit()
@@ -122,6 +127,8 @@ for w in bigw:
     i1 = cr0[w + 2] + int(1200 / dt)
     # a = [t0[i0], t0[i1], *Mlims]
     # c = ax2.imshow(M[:, i0:i1], extent=a, aspect="auto", cmap="Spectral", vmin=-v, vmax=+v)
+    ws = np.ptp(raw_ts["z"][cr0[w]:cr0[w+2]]) * 1e-2
+    log.info(f"Wave [{w}] size: {ws:.2f}m")
     c = ax2.contourf(
         t[i0:i1],
         sj,