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,