Fig_10.ipynb

Kernel: Python 3 (system-wide)

import os

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.transforms as mtransforms


import scipy.integrate as integrate
from scipy import signal

from scipy.signal import find_peaks
from matplotlib.lines import Line2D

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

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

def find_reflected_peaks(x,t,eta,a0,amin=0.5,amax=0.8,xmax=0,tmax=120):
    """
    Identify the reflected peaks using observed charactertics of the solitary wave data

    - reflected peak is always smaller than incident (amax = 0.8 a0)
    - reflected peak has a substantial value (amin = 0.5 a0)
    - reflected peaks do not occur in the channel (by definition)
    - reflected peaks occur shortly after the incident peak (t < 120)
    """
    T,X = np.meshgrid(t,x,indexing='ij')
    mask_loc = np.logical_and(X < xmax, T < tmax)
    mask_amp = np.logical_and(eta > amin*a0, eta < amax*a0)
    mask = np.logical_and(mask_loc,mask_amp)
    
    t_idx = np.zeros(x.shape,dtype=int)
    c_val = np.zeros(x.shape)
    for s_idx in range(len(space)):
        # find all reasonably large peaks
        peaks, _ = find_peaks(eta[:,s_idx],height=amin*a0)
        # find those within the mask
        valid_peaks = peaks[np.where(mask[peaks,s_idx])]
        # if any lie within the mask, store the values
        p = valid_peaks[0] if len(valid_peaks) > 0 else 0
        t_idx[s_idx] = p
        c_val[s_idx] = eta[p,s_idx] if p != 0 else np.nan

    t_val = t[t_idx]
    t_val[np.isnan(c_val)] = np.nan
    
    return t_val, c_val, t_idx

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

h0, g = 50, 9810
c0 = np.sqrt(g*h0)
t0 = h0/c0

fig, axes = plt.subplots(nrows=1,ncols=2,sharey='all',figsize=(5.0,1.85))

#amp_list = [ f'a{a/10:.1f}cm' for a in [5,10,15,20] ]
amp_list = [ 'a1.5cm' ]
linestyles = ['-','--','-.',':']

for amp,linestyle in zip(amp_list,linestyles):
    # choose only the same amplitude cases
    data_dirs = [d for d in data_directories if amp in d]
    # lexical sort will choose no_flow first and with_flow second

    line_kwargs = {
        'color':'k',
        #'linestyle':linestyle,
        'linewidth':1,
    }
    
    # choose the first axis
    ax = axes[0]
    data_dir = data_dirs[0]

    # load data
    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=',')
    t_val = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',usecols=(0,),dtype=int)
    c_val = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',usecols=(1,))
    t_idx = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',usecols=(2,),dtype=int)

    x   = (space*0.3)/h0
    t   = (time/125)/t0
    eta = elev/h0
    c   = c_val/h0
    a0  = float(amp[len('a'):-len('cm')])*10/h0

    ax.plot(x,c,**line_kwargs,label='incident')
    
    # find reflected peaks
    t_val, c_val, t_idx = find_reflected_peaks(x,t,eta,a0,amin=0.4)
    
    ax.plot(x,c_val,**line_kwargs,linestyle='--',label='reflected')

    # choose the second axis
    ax = axes[1]
    data_dir = data_dirs[1]

    # load data
    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=',')
    t_val = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',usecols=(0,),dtype=int)
    c_val = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',usecols=(1,))
    t_idx = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',usecols=(2,),dtype=int)

    x   = (space*0.3)/h0
    t   = (time/125)/t0
    eta = elev/h0
    c   = c_val/h0

    l1, = ax.plot(x,c,**line_kwargs,label='incident')
    
    # find reflected peaks
    t_val, c_val, t_idx = find_reflected_peaks(x,t,eta,a0,amin=0.3,amax=1.0)
    
    l2, = ax.plot(x,c_val,**line_kwargs,linestyle='--',label='reflected')


# combined paramaters
xlim = (-20,3)
ylim = (0,0.47)

#axes.flat[0].set_xlabel(r'$x/h_0$')
axes.flat[0].set_ylabel(r'$a$')

for ax in axes:
    ax.set_xlim(*xlim)
    ax.set_ylim(*ylim)
    ax.set_xlabel(r'$x$')
    
#for label, ax in zip(['a)', 'b)'], axes):
#    ax.set_title(label,loc='left',fontfamily='Times New Roman',fontsize='9.0')

plt.tight_layout()
plt.savefig('Fig_10.pdf')
plt.show()