Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Avatar for JFM-2024-1227.
 Download
3482 views
unlisted
ubuntu2204
Tweet
Share
Share
Kernel: Python 3 (system-wide)
import os
import numpy as np
import matplotlib.pyplot as plt

from scipy import signal
from scipy.interpolate import InterpolatedUnivariateSpline

plt.style.use('thesis')

!find . -maxdepth 3 -type d -name "prepared_data" | sort \
    | xargs -I {} find  {} -mindepth 1 -type d | grep sol3rd | grep -v ".ipynb" | tee "data_directories.txt"
./x=+00.00cm/prepared_data/sol3rd_a1.5cm_x=+00.00cm ./x=+00.00cm/prepared_data/sol3rd_a0.5cm_x=+00.00cm ./x=+00.00cm/prepared_data/sol3rd_a2.0cm_x=+00.00cm ./x=+00.00cm/prepared_data/sol3rd_a1.0cm_x=+00.00cm ./x=-02.54cm/prepared_data/sol3rd_a2.0cm_x=-02.54cm ./x=-02.54cm/prepared_data/sol3rd_a0.5cm_x=-02.54cm ./x=-02.54cm/prepared_data/sol3rd_a1.5cm_x=-02.54cm ./x=-02.54cm/prepared_data/sol3rd_a1.0cm_x=-02.54cm ./x=+02.54cm/prepared_data/sol3rd_a2.0cm_x=+02.54cm ./x=+02.54cm/prepared_data/sol3rd_a1.5cm_x=+02.54cm ./x=+02.54cm/prepared_data/sol3rd_a1.0cm_x=+02.54cm ./x=+02.54cm/prepared_data/sol3rd_a0.5cm_x=+02.54cm ./x=-07.62cm/prepared_data/sol3rd_a2.0cm_x=-07.62cm ./x=-07.62cm/prepared_data/sol3rd_a0.5cm_x=-07.62cm ./x=-07.62cm/prepared_data/sol3rd_a1.0cm_x=-07.62cm ./x=-07.62cm/prepared_data/sol3rd_a1.5cm_x=-07.62cm ./x=+07.62cm/prepared_data/sol3rd_a1.5cm_x=+07.62cm ./x=+07.62cm/prepared_data/sol3rd_a2.0cm_x=+07.62cm ./x=+07.62cm/prepared_data/sol3rd_a1.0cm_x=+07.62cm ./x=+07.62cm/prepared_data/sol3rd_a0.5cm_x=+07.62cm ./x=-100.0cm/prepared_data/sol3rd_a2.0cm_x=-100.0cm ./x=-100.0cm/prepared_data/sol3rd_a1.5cm_x=-100.0cm ./x=-100.0cm/prepared_data/sol3rd_a1.0cm_x=-100.0cm ./x=-100.0cm/prepared_data/sol3rd_a0.5cm_x=-100.0cm ./x=-15.24cm/prepared_data/sol3rd_a1.5cm_x=-15.24cm ./x=-15.24cm/prepared_data/sol3rd_a2.0cm_x=-15.24cm ./x=-15.24cm/prepared_data/sol3rd_a0.5cm_x=-15.24cm ./x=-15.24cm/prepared_data/sol3rd_a1.0cm_x=-15.24cm ./x=+15.24cm/prepared_data/sol3rd_a0.5cm_x=+15.24cm ./x=+15.24cm/prepared_data/sol3rd_a1.0cm_x=+15.24cm ./x=+15.24cm/prepared_data/sol3rd_a2.0cm_x=+15.24cm ./x=+15.24cm/prepared_data/sol3rd_a1.5cm_x=+15.24cm ./x=-50.00cm/prepared_data/sol3rd_a2.0cm_x=-50.00cm ./x=-50.00cm/prepared_data/sol3rd_a0.5cm_x=-50.00cm ./x=-50.00cm/prepared_data/sol3rd_a1.5cm_x=-50.00cm ./x=-50.00cm/prepared_data/sol3rd_a1.0cm_x=-50.00cm ./x=+50.00cm/prepared_data/sol3rd_a2.0cm_x=+50.00cm ./x=+50.00cm/prepared_data/sol3rd_a0.5cm_x=+50.00cm ./x=+50.00cm/prepared_data/sol3rd_a1.5cm_x=+50.00cm ./x=+50.00cm/prepared_data/sol3rd_a1.0cm_x=+50.00cm ./y=+00.00cm/prepared_data/sol3rd_a2.0cm_y=+00.00cm ./y=+00.00cm/prepared_data/sol3rd_a0.5cm_y=+00.00cm ./y=+00.00cm/prepared_data/sol3rd_a1.0cm_y=+00.00cm ./y=+00.00cm/prepared_data/sol3rd_a1.5cm_y=+00.00cm 
with open("data_directories.txt","r") as data_directories_file:
    data_directories = data_directories_file.read().splitlines()

x_loc = "-100.0cm"
#x_loc = "-50.00cm"
data_dirs = sorted([d for d in data_directories if x_loc in d])

in_to_mm = 1/25.4

fig, axes = plt.subplots(nrows=3,ncols=1,sharex='all',
                         figsize=(5.9, 5.0))

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

axes.flat[0].axhline(0,color='k',linewidth=0.5)

for data_dir, linestyle in zip(data_dirs,linestyles):
    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,))
    
    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)
    
    
    lag = 0.008*c_t
    lag -= min(lag)
    
    print(space[np.argmin(lag)])
    
    c_idx = c_t - time[0]
    
    # interpolate surface for a "better" phase estimate
    rng = range(len(space))
    u = elev[c_idx-1,rng]
    c = elev[c_idx  ,rng]
    l = elev[c_idx+1,rng]
    n = u-l
    d = 2*u + 2*l - 4*c
    s = n/d
    s[np.where(np.abs(s)>0.5)] = 0.0
    s /= 125

    # normalization
    amp_str = os.path.basename(data_dir).split('_')[1]
    a0 = int(amp_str[len('a'):-len('cm')][::2])
    h0 = 50
    b0 = 229
    g  = 9810
    
    F     = 1 + (1/2)*(a0/h0) - (3/20)*((a0/h0)**2) + (3/56)*((a0/h0)**3)
    alpha = np.sqrt(3*a0/(4*h0))*(1 - (5/8)*(a0/h0) + (71/128)*((a0/h0)**2))
    c0    = np.sqrt(g*h0)
    
    #t0 = h0/(F*c0)
    #t0 = h0/(alpha*F*c0)
    #t0 = h0/np.sqrt(9810*(h0+a0))
    #t0 = h0/np.sqrt(9810*h0)
    #pos = 0.3 * space / b0
    pos = 0.3 * space / h0
    
    # low pass filtering
    w = signal.get_window('bartlett',30)
    w /= np.sum(w)
    c_v = signal.convolve(c_v,w,'same')
    lag = signal.convolve(lag+s,w,'same')
    for _ in range(3):
        lag = signal.convolve(lag,w,'same')
    
    kwargs = {'color':'k','linestyle':linestyle}
    
    # plotting
    ax = axes.flat[0]
    
    x_idx = np.argmin(np.abs(pos-(12)))
    print(pos[x_idx])
    print(np.amax(c_v/a0-1),c_v[x_idx]/a0-1)
    line, = ax.plot(pos,c_v/a0-1,**kwargs)
    
    line_list.append(line)
    
    ax = axes.flat[1]
    
    t0 = h0/c0
    
    ax.plot(pos,lag/t0,**kwargs)
    
    ax = axes.flat[2]
    
    t0 = h0/(alpha*F*c0)
    
    ax.plot(pos,lag/t0,**kwargs)


axes.flat[0].set_yticks([-0.2,0,0.2])
axes.flat[1].set_yticks([0,0.5,1.0,1.5])
axes.flat[2].set_yticks([0,0.2,0.4,0.6])
    
#plt.xlim(-0.55,2.7)
#plt.xlabel(r'$y/b_0$')

plt.xlim(-2.5,12)
plt.xlabel(r'$y/h_0$')
    
axes.flat[0].set_ylim(-0.2,0.3)
axes.flat[1].set_ylim(0,1.8)
axes.flat[2].set_ylim(0,0.7)

axes.flat[0].set_ylabel(r'$a/a_0-1$')
axes.flat[1].set_ylabel(r'$\Delta t / t_0$')
axes.flat[2].set_ylabel(r'$\alpha F_s \times \Delta t / t_0$')

axes.flat[1].legend(line_list, [ fr"0.{n}" for n in range(1,5) ],
                    loc='center left', bbox_to_anchor=(1,0.5),title=r'$a_0/h_0$')

plt.tight_layout()
plt.savefig('crest_amp_arrival_time_x=-20.png')
plt.show()
2012 12.0 0.19543728882891687 -0.1314860458966205 1994 12.0 0.23987049063114685 -0.12199618269831425 1920 12.0 0.2325078300211192 -0.13292997386700833 1925 12.0 0.23754417223798474 -0.13862347317123136
Image in a Jupyter notebook
#x_loc = "-100.0cm"
x_loc = "-50.00cm"
data_dirs = sorted([d for d in data_directories if x_loc in d])

in_to_mm = 1/25.4

fig, axes = plt.subplots(nrows=2,ncols=1,sharex='all',
                         figsize=(5.9, 3.0))

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

axes.flat[0].axhline(0,color='k',linewidth=0.5)

for data_dir, linestyle in zip(data_dirs,linestyles):
    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,))
    
    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)
    
    
    lag = 0.008*c_t
    lag -= min(lag)
    
    c_idx = c_t - time[0]
    
    # interpolate surface for a "better" phase estimate
    rng = range(len(space))
    u = elev[c_idx-1,rng]
    c = elev[c_idx  ,rng]
    l = elev[c_idx+1,rng]
    n = u-l
    d = 2*u + 2*l - 4*c
    s = n/d
    s[np.where(np.abs(s)>0.5)] = 0.0
    s /= 125

    # normalization
    amp_str = os.path.basename(data_dir).split('_')[1]
    a0 = int(amp_str[len('a'):-len('cm')][::2])
    h0 = 50
    b0 = 229
    g  = 9810
    
    F     = 1 + (1/2)*(a0/h0) - (3/20)*((a0/h0)**2) + (3/56)*((a0/h0)**3)
    alpha = np.sqrt(3*a0/(4*h0))*(1 - (5/8)*(a0/h0) + (71/128)*((a0/h0)**2))
    c0    = np.sqrt(g*h0)
    t0 = h0/c0
    #t0 = h0/(F*c0)
    #t0 = h0/(alpha*F*c0)
    #t0 = h0/np.sqrt(9810*(h0+a0))
    #t0 = h0/np.sqrt(9810*h0)
    #pos = 0.3 * space / b0
    pos = 0.3 * space / h0
    
    # low pass filtering
    w = signal.get_window('bartlett',30)
    w /= np.sum(w)
    c_v = signal.convolve(c_v,w,'same')
    lag = signal.convolve(lag+s,w,'same')
    for _ in range(3):
        lag = signal.convolve(lag,w,'same')
    
    kwargs = {'color':'k','linestyle':linestyle}
    
    # plotting
    ax = axes.flat[0]
    
    x_idx = np.argmin(np.abs(pos-(12)))
    print(pos[x_idx])
    print(np.amax(c_v/a0-1),c_v[x_idx]/a0-1)
    
    line, = ax.plot(pos,c_v/a0-1,**kwargs)
    
    line_list.append(line)
    
    ax = axes.flat[1]
    
    ax.plot(pos,lag/t0,**kwargs)


axes.flat[0].set_yticks([-0.2,0,0.2])
#axes.flat[1].set_yticks([0,0.2,0.4,0.6])
    
#plt.xlim(-0.55,2.7)
#plt.xlabel(r'$y/b_0$')

plt.xlim(-2.5,12)
plt.xlabel(r'$y/h_0$')
    
axes.flat[0].set_ylim(-0.2,0.3)
#axes.flat[1].set_ylim(0,0.7)

axes.flat[0].set_ylabel(r'$a/a_0-1$')
axes.flat[1].set_ylabel(r'$\alpha F_s \times \Delta t / t_0$')

axes.flat[0].legend(line_list, [ fr"0.{n}" for n in range(1,5) ],
                    loc='center left', bbox_to_anchor=(1,0.5),title=r'$a_0/h_0$')

plt.tight_layout()
plt.savefig('crest_amp_arrival_time_linear.png')
plt.show()
12.0 0.22112945563099085 -0.1582698348515621 12.0 0.26004500465027003 -0.15565166961859234 12.0 0.2537532609634183 -0.1597648740256965 12.0 0.259366411367874 -0.1637003037837571
Image in a Jupyter notebook