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

Update zero_cross and add nansasena

Browse source
This commit is contained in:
Edgar P. Burkhart 2022-06-29 09:31:07 +02:00
parent f1bc55843a
commit b8da550849
Signed by: edpibu
GPG key ID: 9833D3C5A25BD227
2 changed files with 37 additions and 3 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

27
data/processing/nandasena.py Normal file
View file

@ -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%}")

13
data/processing/zero_cross.py
View file

@ -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,