Fig_08.ipynb

Kernel: Python 3 (system-wide)

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

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 "with_flow" | grep "sol3rd" | tee "data_directories.txt" > /dev/null

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

y_loc = "y=+00.00cm"
data_dirs = sorted([d for d in data_directories if y_loc in d])

in_to_mm = 1/25.4

linestyles = ['-','--','-.',':']
markstyles = ['o','^' ,'x' ,'s']
line_list = []

plt.figure(figsize=(4.5,1.5),layout="constrained")

for data_dir, linestyle, marker in zip(data_dirs,linestyles,markstyles):
    #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=',')
    c_t   = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',dtype=int,usecols=(0,))
    c_v   = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',usecols=(1,))
    
    # c_t is an integer number of dt=0.008 s
    c_t = c_t*0.008
    # space is an integer number of dx=0.3mm
    space = space*0.3
    
    name = os.path.basename(data_dir).split('_')[-1]
    amp_str = os.path.basename(data_dir).split('_')[1][len('a'):-len('cm')]
    amp = float(amp_str)
    
    mask = space < -200
    
    h = 50
    # wavemaker amplitude
    a = 10*amp/h
    g = 9810
    # theoretical speed at wavemaker amplitude
    F = 1 + 0.5*a - 3/20*(a**2) + (3/56)*(a**3)
    c0 = np.sqrt(g*h)
    theory_speed = c0*F
    
    # normalize
    space /= h
    space -= np.min(space)
    c_t   *= np.sqrt(g/h)
    #c_t   *= F
    c_t   -= np.min(c_t)
    
    result = stats.linregress(c_t[mask],space[mask])
    m, 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
    #print(c_conf)
    
    plt.plot(a,(m-F),'k',linestyle='',marker='.',label=fr'{a:.1f}')
    
    
plt.xlim(0,0.5)
plt.ylim(-0.1,0)

plt.xlabel(r'$a_0$')
plt.ylabel(r'$\Delta c$')

plt.yticks([0,-0.03,-0.06,-0.09])

# 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)

#plt.tight_layout()
plt.savefig("Fig_08.pdf",dpi=600)