From 873cc73c8428690a3f705920a5ed21fb13d752cf Mon Sep 17 00:00:00 2001
From: "Edgar P. Burkhart" <git@edgarpierre.fr>
Date: Mon, 4 Apr 2022 10:26:58 +0200
Subject: [PATCH 1/3] PUV method test (with/out noise)

---
 swash/processing/r_test.py | 50 ++++++++++++++++++++++++++++++++++++++
 tasks.md                   |  2 +-
 2 files changed, 51 insertions(+), 1 deletion(-)
 create mode 100644 swash/processing/r_test.py

diff --git a/swash/processing/r_test.py b/swash/processing/r_test.py
new file mode 100644
index 0000000..c052c86
--- /dev/null
+++ b/swash/processing/r_test.py
@@ -0,0 +1,50 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from numpy import random
+import scipy.signal as sgl
+
+yi = random.normal(size=2**20)
+yr = 0.3 * np.roll(yi, -(2**10))
+
+figy, axy = plt.subplots()
+axy.plot(np.arange(2**10, 2**11), yi[2**10 : 2**11])
+axy.plot(np.arange(2**10, 2**11), yr[2**10 : 2**11])
+
+figf, axf = plt.subplots()
+axf.plot(*sgl.welch(yi))
+axf.plot(*sgl.welch(yr))
+
+eta = yi + yr
+u = -yi + yr
+
+
+def puv(eta, u):
+    f, phi_eta = sgl.welch(eta)
+    phi_u = sgl.welch(u)[1]
+    phi_eta_u = np.abs(sgl.csd(eta, u)[1].real)
+
+    return f, np.sqrt(
+        (phi_eta + phi_u - 2 * phi_eta_u) / (phi_eta + phi_u + 2 * phi_eta_u)
+    )
+
+
+figr, axr = plt.subplots()
+
+axr.plot(*puv(eta, u), label="Without noise")
+axr.plot(
+    *puv(
+        eta + 0.4 * random.normal(size=2**20), u + 0.4 * random.normal(size=2**20)
+    ),
+    label="With noise"
+)
+axr.grid()
+axr.autoscale(True, "x", tight=True)
+axr.set(ylim=(0, 1), ylabel="R", xlabel="f")
+axr.legend(loc="lower left")
+
+axpsd = axr.twinx()
+axpsd.plot(*sgl.welch(eta), label=r"PSD ($\eta$)", c="k", alpha=0.2, lw=1)
+axpsd.legend(loc="lower right")
+axpsd.set(ylabel="PSD")
+
+plt.show()
diff --git a/tasks.md b/tasks.md
index 9f2cee4..ad2c448 100644
--- a/tasks.md
+++ b/tasks.md
@@ -1,4 +1,4 @@
- Vérification méthodes calcul réflection avec données forcées en Python
+* Vérification méthodes calcul réflection avec données forcées en Python
 
 
 Vérification vague incidente possible avec digue et sans digue sur temps long (4h) avec spectre Jonswap

From 9e29a30b42c1e942cbd98b0ce2933c5bc54f2ccf Mon Sep 17 00:00:00 2001
From: "Edgar P. Burkhart" <git@edgarpierre.fr>
Date: Mon, 4 Apr 2022 10:59:17 +0200
Subject: [PATCH 2/3] PUV test with multiple reflection coefficients

---
 swash/processing/r_test.py | 31 +++++++++++++++----------------
 1 file changed, 15 insertions(+), 16 deletions(-)

diff --git a/swash/processing/r_test.py b/swash/processing/r_test.py
index c052c86..5a67f1c 100644
--- a/swash/processing/r_test.py
+++ b/swash/processing/r_test.py
@@ -4,7 +4,7 @@ from numpy import random
 import scipy.signal as sgl
 
 yi = random.normal(size=2**20)
-yr = 0.3 * np.roll(yi, -(2**10))
+yr = np.roll(yi, -(2**10))
 
 figy, axy = plt.subplots()
 axy.plot(np.arange(2**10, 2**11), yi[2**10 : 2**11])
@@ -14,8 +14,8 @@ figf, axf = plt.subplots()
 axf.plot(*sgl.welch(yi))
 axf.plot(*sgl.welch(yr))
 
-eta = yi + yr
-u = -yi + yr
+eta = lambda r: yi + r * yr
+u = lambda r: -yi + r * yr
 
 
 def puv(eta, u):
@@ -30,21 +30,20 @@ def puv(eta, u):
 
 figr, axr = plt.subplots()
 
-axr.plot(*puv(eta, u), label="Without noise")
-axr.plot(
-    *puv(
-        eta + 0.4 * random.normal(size=2**20), u + 0.4 * random.normal(size=2**20)
-    ),
-    label="With noise"
-)
+for r in np.arange(0, 1.1, 0.1):
+    axr.plot(*puv(eta(r), u(r)), c="k")
+    Rn = puv(
+        eta(r) + 0.4 * random.normal(size=2**20),
+        u(r) + 0.4 * random.normal(size=2**20),
+    )
+    axr.plot(
+        *Rn,
+        c="#ff6600",
+    )
+    axr.annotate(f"{r=:.1f}", (Rn[0][0], Rn[1][0]), bbox={"boxstyle": "square", "facecolor": "w"})
 axr.grid()
 axr.autoscale(True, "x", tight=True)
 axr.set(ylim=(0, 1), ylabel="R", xlabel="f")
-axr.legend(loc="lower left")
-
-axpsd = axr.twinx()
-axpsd.plot(*sgl.welch(eta), label=r"PSD ($\eta$)", c="k", alpha=0.2, lw=1)
-axpsd.legend(loc="lower right")
-axpsd.set(ylabel="PSD")
+axr.legend(("No noise", "40% noise"), loc="lower left")
 
 plt.show()

From 454655e430171b4fbfd41a01768e2091d1c5c56a Mon Sep 17 00:00:00 2001
From: "Edgar P. Burkhart" <git@edgarpierre.fr>
Date: Mon, 4 Apr 2022 11:55:23 +0200
Subject: [PATCH 3/3] Orbitals for huge wave

---
 data/config.ini             |  2 +-
 data/data/.gitignore        |  1 +
 data/processing/orbitals.py | 97 +++++++++++++++++++++++++++++++++++++
 tasks.md                    |  2 +-
 4 files changed, 100 insertions(+), 2 deletions(-)
 create mode 100644 data/processing/orbitals.py

diff --git a/data/config.ini b/data/config.ini
index 955650f..3daa05e 100644
--- a/data/config.ini
+++ b/data/config.ini
@@ -5,7 +5,7 @@ hires=bathyhires.dat
 hstru=Hstru.dat
 poro=Poro.dat
 psize=Psize.dat
-raw_ts=201702281700.raw,201702281730.raw
+raw_ts=cerema/raw/201702281700.raw,cerema/raw/201702281730.raw
 hires_step=0.5
 
 [out]
diff --git a/data/data/.gitignore b/data/data/.gitignore
index 47293d7..702fc79 100644
--- a/data/data/.gitignore
+++ b/data/data/.gitignore
@@ -1,2 +1,3 @@
 *.xyz
 *.raw
+cerema
\ No newline at end of file
diff --git a/data/processing/orbitals.py b/data/processing/orbitals.py
new file mode 100644
index 0000000..63eafb5
--- /dev/null
+++ b/data/processing/orbitals.py
@@ -0,0 +1,97 @@
+import argparse
+import configparser
+import logging
+import pathlib
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+parser = argparse.ArgumentParser(description="Plot orbitals")
+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("bathy")
+
+log.info("Starting time-series pre-processing")
+config = configparser.ConfigParser()
+config.read(args.config)
+
+inp_root = pathlib.Path(config.get("inp", "root"))
+out_root = pathlib.Path(config.get("out", "root"))
+
+out_ts = out_root.joinpath("ts.dat")
+
+raw_ts = []
+for tsi in config.get("inp", "raw_ts").split(","):
+    raw_ts.append(
+        np.loadtxt(
+            inp_root.joinpath(tsi),
+            dtype=[("state", int), ("z", float), ("y", float), ("x", float)],
+            delimiter=",",
+            max_rows=2304,
+        )
+    )
+n = len(raw_ts)
+raw_ts = np.concatenate(raw_ts)
+log.debug(f"{raw_ts=}")
+
+if (errs := np.count_nonzero(raw_ts["state"])) != 0:
+    log.warning(f"{errs} transmission errors!")
+    log.debug(f"{dict(zip(*np.unique(raw_ts['state'], return_counts=True)))}")
+
+t = np.linspace(0, 30 * 60 * n, 2304 * n + 1)[:-1]
+log.debug(f"{t=}")
+
+flt = (t > 1370) & (t < 1405)
+
+figt, axt = plt.subplots(3)
+axt[0].plot(t, raw_ts["x"])
+axt[1].plot(t, raw_ts["y"])
+axt[2].plot(t, raw_ts["z"])
+for ax in axt:
+    ax.axvline(t[flt].min(), c="k")
+    ax.axvline(t[flt].max(), c="k")
+    ax.grid()
+    ax.set(xlim=(t.min(), t.max()))
+
+ts_flt = raw_ts[flt]
+z0 = ts_flt["z"]
+figtz, axtz = plt.subplots(3)
+axtz[0].plot(t[flt], ts_flt["x"])
+axtz[1].plot(t[flt], ts_flt["y"])
+axtz[2].plot(t[flt], z0)
+for ax in axtz:
+    ax.grid()
+    ax.set(xlim=(t[flt].min(), t[flt].max()))
+
+fig3d = plt.figure()
+ax3d = fig3d.add_subplot(projection="3d")
+ax3d.plot(ts_flt["x"], ts_flt["y"], z0, c="#0066ff")
+ax3d.quiver3D(
+    ts_flt["x"][:-1],
+    ts_flt["y"][:-1],
+    z0[:-1],
+    np.diff(ts_flt["x"])[:],
+    np.diff(ts_flt["y"])[:],
+    np.diff(z0)[:],
+    color="#0066ff",
+)
+ax3d.set(xlabel="x (cm)", ylabel="y (cm)", zlabel="z (cm)")
+
+theta = np.angle(raw_ts["x"] + 1j * raw_ts["y"]).mean()
+fig2dv, ax2dv = plt.subplots()
+x0 = ts_flt["x"] * np.cos(theta) + ts_flt["y"] * np.sin(theta)
+ax2dv.plot(x0, z0, c="#0066ff", lw=1)
+ax2dv.quiver(
+    x0[:-1],
+    z0[:-1],
+    np.diff(x0)[:],
+    np.diff(z0)[:],
+    color="#0066ff",
+)
+ax2dv.grid()
+ax2dv.set(aspect="equal")
+
+plt.show()
diff --git a/tasks.md b/tasks.md
index ad2c448..900eaaa 100644
--- a/tasks.md
+++ b/tasks.md
@@ -4,7 +4,7 @@
 Vérification vague incidente possible avec digue et sans digue sur temps long (4h) avec spectre Jonswap
 
 
-Tracer trajectoires bouées (3d, 2d) autour de vague 15m
+* Tracer trajectoires bouées (3d, 2d) autour de vague 15m
 
 
 Format rapport: Journal of Geophysical Research