Fig_22.ipynb

Kernel: Python 3 (system-wide)

import os
import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.transforms as mtransforms

from matplotlib import colors

from scipy import signal
from scipy.interpolate import InterpolatedUnivariateSpline

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

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

def t_idx_from_c_t(c_t,time):
    t_idx = np.zeros(c_t.shape,dtype=int)
    for idx, t_val in enumerate(c_t.flat):
        candidate = np.argwhere(time == t_val)
        assert(len(candidate) == 1)
        t_idx[idx] = candidate.flat[0]
        
    return t_idx

def find_idx_nearest(x,v):
    return np.argmin(np.abs(x-v))

def find_nearest_idx(x,val):
    return np.argmin(np.abs(x-val))

def solitary_amp(a0,h0,x,t,xoff=0,toff=0,g=9.81):
    '''Returns the amplitude of a 3rd order solitary wave on a meshgrid or as scalar'''
    assert isinstance(a0, (int, float)) and  isinstance(h0, (int, float)) \
        and isinstance(x, (int, float, np.ndarray)) and isinstance(t, (int, float, np.ndarray)) \
        and isinstance(toff, (int, float)) and isinstance(g, (int, float))
    
    from numpy import sqrt
    from numpy import cosh
    
    a = a0/h0
    c0 = sqrt(g*h0)
    F = 1 + (1/2)*a - (3/20)*(a**2) + (3/56)*(a**3)
    c = c0*F
    alpha = sqrt(3*a/4) * (1 - (5/8)*a + (71/128)*(a**2))
    xi = x/h0
    xioff = xoff/h0
    tau = c0*t/h0
    tauoff = c0*toff/h0
    X, T = np.meshgrid(xi,tau,indexing='ij')
    if X.size == 1 and T.size == 1:
        X = X.ravel()[0]
        T = T.ravel()[0]
    s = 1/(cosh(alpha*((X-xioff) - F*(T - tauoff))))
    return h0*(a*(s**2) - (3/4)*(a**2)*(s**2 - s**4) + (a**3)*( (5/8)*(s**2) - (151/80)*(s**4) + (101/80)*(s**6) ))


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

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

fig, axes = plt.subplots(nrows=2,ncols=1,sharex='all',
                         figsize=(4.0,2.4),layout='constrained')

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

for data_dir, ax in zip(data_dirs,axes):
    # load data
    c_t_kwargs = {'delimiter':',','dtype':int,'usecols':(0,)}
    c_v_kwargs = {'delimiter':',','dtype':float,'usecols':(1,)}
    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 ])

    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=',')
    
    x   = (space*0.3)/h0
    t   = (time/125)/t0
    eta = elev/h0

    T,X = np.meshgrid(t,x,indexing='ij')

    # :1 only compares the first two
    for cname1,cname2,rname1,rname2 in zip(crestnames[:1],crestnames[1:],troughnames[:-1],troughnames[1:]):
        c_v1 = np.loadtxt(cname1,**c_v_kwargs)/h0
        c_v2 = np.loadtxt(cname2,**c_v_kwargs)/h0
        
        r_v1 = np.loadtxt(rname1,**c_v_kwargs)/h0
        r_v1 = np.loadtxt(rname1,**c_v_kwargs)/h0

        c_t1_int = np.loadtxt(cname1,**c_t_kwargs)
        c_t2_int = np.loadtxt(cname2,**c_t_kwargs)
        
        r_t1_int = np.loadtxt(rname1,**c_t_kwargs)
        r_t2_int = np.loadtxt(rname2,**c_t_kwargs)

        c_t1 = c_t1_int/125
        c_t2 = c_t2_int/125

        c_t1 /= t0
        c_t2 /= t0

        crest_limit = 68 - x if "with_flow" in data_dir else 64 - x
        mask = c_t2 < crest_limit

        # find times when both crests are present
        # Return the sorted, unique values that are in both of the input arrays.
        match = np.intersect1d(c_t1_int[mask],c_t2_int[mask])

        if not match.size >= 1:
            print(f"no match between {os.path.basename(cname1)} and {os.path.basename(cname2)}")
            continue

        lengths = np.empty(match.shape)
        steepness = np.zeros(match.shape)
        t_idxs = np.empty(match.shape,dtype=int)
        loc1 = np.empty(match.shape)
        loc2 = np.empty(match.shape)
        loc3 = np.empty(match.shape)
        amp1 = np.empty(match.shape)
        amp2 = np.empty(match.shape)
        amp3 = np.empty(match.shape)
        for idx, t_int in enumerate(match):
            t_idxs[idx] = np.argwhere(time==t_int).flat[0]
            loc1[idx] = np.nanmedian(np.where(c_t1_int==t_int,x   ,np.nan))
            loc2[idx] = np.nanmedian(np.where(c_t2_int==t_int,x   ,np.nan))
            loc3[idx] = np.nanmedian(np.where(r_t1_int==t_int,x   ,np.nan))
            amp1[idx] = np.nanmedian(np.where(c_t2_int==t_int,c_v1,np.nan))
            amp2[idx] = np.nanmedian(np.where(c_t2_int==t_int,c_v2,np.nan))
            amp3[idx] = np.nanmedian(np.where(r_t1_int==t_int,r_v1,np.nan))
            
            lengths[idx] = loc1[idx] - loc2[idx]
            steepness[idx] = (np.nanmedian([amp1[idx],amp2[idx]]) - amp3[idx])/lengths[idx]
        
        
        #idx = find_nearest_idx(loc2,-18)
        #print(loc2[idx])
        #t_int = match[idx]
        
        t_int = 581
        t_idx = np.argwhere(t_int == time).flat[0]
        idx = np.argwhere(t_int == match).flat[0]
        
        print(t[t_idx])

        # where returns the same shape as c_t1_int so we want to reduce to a single value
        loc1 = loc1[idx]
        loc2 = loc2[idx]
        loc3 = loc3[idx]
        amp1 = np.amax(eta[t_idx,:])
        amp2 = amp2[idx]
        amp3 = amp3[idx]
        
        ax.plot(x,eta[t_idx,:],'k-')
        ampm = np.mean([amp1,amp2])
        linekwargs = {'color':'k','linestyle':'--','linewidth':0.9}
        ax.hlines(ampm,loc1,loc2,**linekwargs)
        ax.vlines(loc1,amp1,ampm,**linekwargs)
        ax.vlines(loc2,ampm,amp2,**linekwargs)
        #ax.vlines(loc3,amp3,ampm,**linekwargs)

        print(loc1,loc2,np.nanmean(lengths),np.nanmin(lengths),np.nanmax(lengths))
        print(amp1,np.amax(eta[t_idx,:]),c_v1[t_idx])

    sol3rd_profile = solitary_amp(amp1*h0,h0,x*h0,0,
                                  xoff=loc1*h0,toff=0,g=g)/h0
    ax.plot(x,sol3rd_profile,color='k',linestyle='-.',linewidth=0.95)

# If we want to show the fact that the offshore loc is aligned
#for ax in axes:
#    ax.axvline(-18,**linekwargs)

for ax in axes:
    ax.axhline(0,color='k',linewidth=0.5)
    ax.set_xlim(-20,0)
    ax.set_ylim(0,0.21)
    ax.set_ylabel(r'$\eta$')

axes[-1].set_xlabel(r'$x$')

#for label, ax in zip(['a)', 'b)'], axes):
#    trans = mtransforms.ScaledTranslation(-29/72, 7/72, fig.dpi_scale_trans)
#    ax.text(0.0, 1.0, label, transform=ax.transAxes + trans,
#            fontsize='9.0', verticalalignment='top', fontfamily='Times New Roman',
#            bbox=dict(facecolor='1.0', edgecolor='none', pad=2.0))

plt.savefig('Fig_22.pdf')

65.10519153493061
-9.402 -17.936999999999998 8.50009756097561 8.34 8.592
0.14640914490044735 0.14640914490044735 0.14718125118833145
65.10519153493061
-10.001999999999999 -19.185 9.214333333333334 9.113999999999997 9.326999999999998
0.19590700180798634 0.19590700180798634 0.19103806572821672