forked from susanhen/wave_tools
-
Notifications
You must be signed in to change notification settings - Fork 0
/
breaking_1d.py
71 lines (64 loc) · 2.39 KB
/
breaking_1d.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
from types import DynamicClassAttribute
import numpy as np
from numpy.lib import scimath as SM
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
from scipy.special import gamma
from scipy.interpolate import interp1d
from scipy.integrate import dblquad, quad, simps
from scipy import integrate
from matplotlib import cm
from help_tools import plotting_interface
from wave_tools import surface_core
from wave_tools import peak_tracking
import numpy.ma as ma
import math
bsurf = surface_core.spacetempSurface.surface_from_file('surfprofile')
bsurf.load_velocity('velprofile')
def breaking_tracking_old(surf, L, T):
pt = peak_tracking.get_PeakTracker(surf.x, surf.t, surf.eta, surf.vel)
pt.breaking_tracker()
msurf = np.zeros((np.size(surf.t), np.size(surf.x)))
xind = int(math.ceil(np.round(L/surf.dx)/2.)*2)
tind = int(np.round(T/surf.dt))
for i in range(0, pt.Nb):
tloc = pt.bindex[i,0]
xloc = pt.bindex[i,1]
dis = 0
speed = np.abs(pt.pc[i+1])
for j in range(0, tind):
if tloc + j >= np.size(surf.t):
break
for k in range(int(-xind/2), int(xind/2)):
if xloc - k < 0:
break
msurf[tloc+j, xloc-k] = 1
dis += surf.dt*(speed)
while dis >= surf.dx:
xloc -= 1
dis -= surf.dx
return msurf, pt
def breaking_tracking(surf, peakTracker, L, T):
mask = np.zeros(surf.eta.shape)
peak_dict = peakTracker.get_peak_dict()
ids_breaking_peaks = peakTracker.get_ids_breaking_peaks()
dt = peakTracker.dt
dx = peakTracker.dx
Nt = peakTracker.Nt
x_min = peakTracker.x[0]
N_Tb = int(T/dt) # number of points in breaking region in time
for id in ids_breaking_peaks:
peak = peak_dict[id]
t0_ind = peak.get_breaking_start_ind_t()
x0 = peak.get_breaking_start_x()
cb = peak.cb
#N_L = int(L/(np.abs(c)*dt)) # number of points in breaking region in space
N_L = int(L/dx)
xb = x0 + np.arange(0, N_L) * cb*dt
boundary_mask = np.where(xb>=0, 0, 1)
xb = ma.masked_array(xb, mask=boundary_mask).compressed()
xb_ind = ((xb-x_min)/dx).astype('int')
for i in range(0, N_L):
t_end_ind = np.min([t0_ind+i+N_Tb, Nt])
mask[t0_ind+i:t_end_ind, xb_ind[i]] = 1.0
return mask