From 7061abe2841e24af68ea9be741889d19e8cc6fc9 Mon Sep 17 00:00:00 2001
From: "Edgar P. Burkhart" <git@edgarpierre.fr>
Date: Thu, 3 Mar 2022 11:16:40 +0100
Subject: [PATCH] Post-process output images

---
 swash/config.ini         |  5 ++++-
 swash/processing/post.py | 33 +++++++++++++++++++++++----------
 2 files changed, 27 insertions(+), 11 deletions(-)

diff --git a/swash/config.ini b/swash/config.ini
index 4cfce0e..0310686 100644
--- a/swash/config.ini
+++ b/swash/config.ini
@@ -13,4 +13,7 @@ out=out_data
 [swash]
 input=sws/INPUT.sws
 swashrun=/opt/swash/swash/swashrun
-out=out_swash
+out=out
+
+[post]
+out=out_post
diff --git a/swash/processing/post.py b/swash/processing/post.py
index a1183eb..dafe5c1 100644
--- a/swash/processing/post.py
+++ b/swash/processing/post.py
@@ -25,32 +25,45 @@ vel = read_nohead_vect(root.joinpath("vel.dat"), n)
 n_t = botl.shape[0]
 
 # Cospectral calculations
-pos_x = n//10
+pos_x = n // 10
+f = 1 / 0.25
 
-eta = (dep - botl)[n_t//2:, pos_x]
-u = vel[n_t//2:, 0, pos_x]
+eta = (dep - botl)[n_t // 2 :, pos_x]
+u = vel[n_t // 2 :, 0, pos_x]
 
-phi_eta = np.abs(sgl.csd(eta, eta)[1])
-phi_u = np.abs(sgl.csd(u, u)[1])
-phi_eta_u = np.abs(sgl.csd(eta, u)[1])
+phi_eta = np.abs(sgl.csd(eta, eta, f))
+phi_u = np.abs(sgl.csd(u, u, f))
+phi_eta_u = np.abs(sgl.csd(eta, u, f))
 
-R = np.sqrt((phi_eta + phi_u - 2*phi_eta_u)/(phi_eta + phi_u + 2*phi_eta_u))
+R = np.sqrt(
+    (phi_eta[1] + phi_u[1] - 2 * phi_eta_u[1])
+    / (phi_eta[1] + phi_u[1] + 2 * phi_eta_u[1])
+)
 
 # Plotting
 fig, (ax_dep, ax_vel) = plt.subplots(2)
 
 ax_dep.plot((dep - botl)[:, pos_x], label="dep", color="#0066ff")
 ax_dep.set(xlabel="t (s)", ylabel="z (m)")
+ax_dep.autoscale(axis="x", tight=True)
+ax_dep.grid()
 
 ax_vel.plot(vel[:, 0, pos_x], label="vel")
 ax_vel.set(xlabel="t (s)", ylabel="U (m/s)")
+ax_vel.autoscale(axis="x", tight=True)
+ax_vel.grid()
 
 fig.tight_layout()
 
 fig_r, ax_r = plt.subplots()
 
-ax_r.plot(R)
-ax_r.set(ylim=(0,1))
+ax_r.plot(phi_eta[0], R)
+ax_r.autoscale(axis="x", tight=True)
+ax_r.set(ylim=(0, 1), xlabel="f (Hz)", ylabel="R")
 ax_r.grid()
 
-plt.show(block=True)
+out = pathlib.Path(config.get("post", "out"))
+out.mkdir(exist_ok=True)
+
+fig.savefig(out.joinpath(f"{pos_x}.png"))
+fig_r.savefig(out.joinpath(f"R{pos_x}.png"))