m2cce/internship
m2cce/internship
1
Fork 0
Code Issues Activity

Processing reflection coefficient with PUV method

Browse source
This commit is contained in:
Edgar P. Burkhart 2022-03-03 10:54:09 +01:00
parent e2b495c5a9
commit 588318e4ce
Signed by: edpibu
GPG key ID: 9833D3C5A25BD227
1 changed files with 25 additions and 5 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

30
swash/processing/post.py
View file

@ -4,6 +4,7 @@ import pathlib
import numpy as np import numpy as np
import pandas as pd import pandas as pd
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import scipy.signal as sgl
from .read_swash import * from .read_swash import *
@ -17,20 +18,39 @@ root = pathlib.Path(config.get("swash", "out"))
bathy = pd.read_hdf(cache.joinpath("bathy.h5"), "bathy") bathy = pd.read_hdf(cache.joinpath("bathy.h5"), "bathy")
n = bathy.index.size n = bathy.index.size
# xp = read_nohead_scalar(root.joinpath("xp.dat"), n)
botl = read_nohead_scalar(root.joinpath("botl.dat"), n) botl = read_nohead_scalar(root.joinpath("botl.dat"), n)
dep = np.maximum(0, read_nohead_scalar(root.joinpath("dep.dat"), n)) dep = np.maximum(0, read_nohead_scalar(root.joinpath("dep.dat"), n))
vel = read_nohead_vect(root.joinpath("vel.dat"), n) vel = read_nohead_vect(root.joinpath("vel.dat"), n)
# watl = read_nohead_scalar(root.joinpath("watl.dat"), n)
n_t = botl.shape[0]
# Cospectral calculations
pos_x = n//10
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])
R = np.sqrt((phi_eta + phi_u - 2*phi_eta_u)/(phi_eta + phi_u + 2*phi_eta_u))
# Plotting
fig, (ax_dep, ax_vel) = plt.subplots(2) fig, (ax_dep, ax_vel) = plt.subplots(2)
ax_dep.plot((dep - botl)[:, n//2], label="dep", color="#0066ff") ax_dep.plot((dep - botl)[:, pos_x], label="dep", color="#0066ff")
ax_dep.set(xlabel="t (s)", ylabel="z (m)") ax_dep.set(xlabel="t (s)", ylabel="z (m)")
ax_vel.plot(vel[:, 0, n//2], label="vel") ax_vel.plot(vel[:, 0, pos_x], label="vel")
ax_vel.set(xlabel="t (s)", ylabel="U (m/s)") ax_vel.set(xlabel="t (s)", ylabel="U (m/s)")
fig.tight_layout() fig.tight_layout()
#fig.legend()
fig_r, ax_r = plt.subplots()
ax_r.plot(R)
ax_r.set(ylim=(0,1))
ax_r.grid()
plt.show(block=True) plt.show(block=True)