Fig_23.ipynb

Kernel: Python 3 (system-wide)

import os
import numpy as np
import matplotlib.pyplot as plt

from matplotlib import colors

from scipy import signal, stats
from scipy.interpolate import InterpolatedUnivariateSpline
from scipy.stats import t as t_dist

plt.style.use('../jfm.mplstyle')

!find ../Experimental_Dataset -type f -name "elevation.csv.gz" | xargs dirname | sort | grep "undular" | grep "y=+00.00cm" | tee "data_directories.txt" > /dev/null

with open("data_directories.txt","r") as data_directories_file:
    data_directories = data_directories_file.read().splitlines()

speeds = {}
speeds_conf = {}

for data_dir in data_directories:
    crestnames = sorted([ f"{data_dir}/{name}" for name in os.listdir(data_dir) if "crests" in name ])
    troughnames = sorted([ f"{data_dir}/{name}" for name in os.listdir(data_dir) if "troughs" in name ])

    c_t_kwargs = {'delimiter':',','dtype':int,'usecols':(0,)}
    c_v_kwargs = {'delimiter':',','dtype':float,'usecols':(1,)}

    time  = np.loadtxt(f"{data_dir}/time.csv",delimiter=',',dtype=int)
    space = np.loadtxt(f"{data_dir}/space.csv",delimiter=',',dtype=int)
    elev  = np.loadtxt(f"{data_dir}/elevation.csv",delimiter=',')
    
    h0, a0 = 50, 5
    g = 9810
    t0 = np.sqrt(h0/g)

    x   = space*0.3/h0
    t   = (time/125)/t0
    eta = elev/h0
    
    T,X = np.meshgrid(t,x,indexing='ij')
    
    c_t_kwargs = {'delimiter':',','dtype':int,'usecols':(0,)}
    c_v_kwargs = {'delimiter':',','dtype':float,'usecols':(1,)}

    linestyles = ['-','--','-.',':']

    F = 1.075
    ub = 0.0977

    speed = np.empty((4,))
    speed_conf = np.empty((4,))
    
    for idx,name in enumerate(crestnames):
        c_t = (np.loadtxt(name,**c_t_kwargs)/125)/t0

        crest_limit = 68 - x if "with_flow" in data_dir else 63 - x
        mask = c_t < crest_limit

        w = signal.get_window('boxcar',30)
        w /= np.sum(w)
        c_t_smooth = np.convolve(c_t,w,'same')

        # speed
        result = stats.linregress(c_t[mask],x[mask])
        c, b = result.slope, result.intercept
        tinv = lambda alpha, df: abs(t_dist.ppf(alpha/2, df))
        # Let M=1 such that N-M-1 = N-2
        c_conf = tinv(0.05,len(c_t[mask])-2)*result.stderr

        speed[idx] = c
        speed_conf[idx] = c_conf
        
    condition = "with-flow" if "with_flow" in data_dir else "no-flow"
    speeds[condition] = speed.copy()
    speeds_conf[condition] = speed_conf.copy()
    
fig, ax = plt.subplots(figsize=(4.5,1.5),layout='constrained')

rng = np.arange(4)

ax.errorbar(rng,speeds["with-flow"]-speeds["no-flow"],yerr=speeds_conf["no-flow"]+speeds_conf["with-flow"],
           capsize=2,barsabove=True,linestyle='',color='k')

# flow velocity at water surface at x=-20
us = -0.0834

plt.axhline(us*7/8,color='k',linestyle='--',linewidth=0.7)
plt.axhline(us*(4*np.sqrt(2)/9),color='k',linestyle='-',linewidth=0.7)

ax.set_xticks(rng, ("1st","2nd","4rd","5th"))

ax.set_ylim(-0.1,0)

ax.set_xlabel('crest')
ax.set_ylabel(r'$\Delta c$')

plt.savefig('Fig_23.pdf',dpi=600)
plt.show()