From cb1a94c1ecc1b263e8c379d394a79b247bf1b1d3 Mon Sep 17 00:00:00 2001
From: shfe <shfe@dtu.dk>
Date: Wed, 27 Nov 2019 07:42:15 +0100
Subject: [PATCH] some functions for wave modeling and soil modeling
---
wetb/soil/KlinkvortModel.py | 121 ++++++++
wetb/soil/LeblancModel.py | 189 +++++++++++
wetb/soil/__init__.py | 0
wetb/soil/hyperbolic_model.py | 48 +++
wetb/soil/py_model.py | 167 ++++++++++
wetb/wave/__init__.py | 0
wetb/wave/dispersion.py | 157 ++++++++++
wetb/wave/numba_misc.py | 569 ++++++++++++++++++++++++++++++++++
wetb/wave/spectrum.py | 531 +++++++++++++++++++++++++++++++
wetb/wave/zerocrossing.py | 281 +++++++++++++++++
10 files changed, 2063 insertions(+)
create mode 100644 wetb/soil/KlinkvortModel.py
create mode 100644 wetb/soil/LeblancModel.py
create mode 100644 wetb/soil/__init__.py
create mode 100644 wetb/soil/hyperbolic_model.py
create mode 100644 wetb/soil/py_model.py
create mode 100644 wetb/wave/__init__.py
create mode 100644 wetb/wave/dispersion.py
create mode 100644 wetb/wave/numba_misc.py
create mode 100644 wetb/wave/spectrum.py
create mode 100644 wetb/wave/zerocrossing.py
diff --git a/wetb/soil/KlinkvortModel.py b/wetb/soil/KlinkvortModel.py
new file mode 100644
index 00000000..bf9b7f58
--- /dev/null
+++ b/wetb/soil/KlinkvortModel.py
@@ -0,0 +1,121 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Jun 18 13:34:12 2018
+
+@author: shfe
+
+Note:
+This script is used for Klinkvort model used for cyclic accumulation calculation
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+plt.close('all')
+
+class KlinkvortModel(object):
+
+
+ def __init__(self, xi_b, xi_c, N):
+ """
+ Initialize the input parameters
+ """
+
+ self.xi_b = xi_b
+ self.xi_c = xi_c
+ self.N = N
+
+ def alpha(self):
+ """
+ Calculate the parameter alpha
+ """
+ # Calculate the Tb, Tc
+ if not 0 <= self.xi_b <= 1:
+ raise ValueError('The maximum load ratio should between 0 to 1')
+
+ if not -1 <= self.xi_c <= 1:
+ raise ValueError('The cyclic load ratio should between -1 to 1')
+
+ if self.xi_b <= 0.022:
+ Tb = 0
+ else:
+ Tb = 0.61*self.xi_b - 0.013
+
+ Tc = (self.xi_c + 0.63)*(self.xi_c - 1)*(self.xi_c - 1.64)
+ alpha = Tb * Tc
+
+ return Tb, Tc, alpha
+
+ def evoluation(self, y1):
+ """
+ Evoluation model for cyclic accumulated rotation
+ """
+
+ _, _, alpha = self.alpha()
+ yn = y1 * self.N**alpha
+
+ return yn
+
+ def equivalentnum(self, y1, yn, alpha):
+ """
+ Equivalent number for cylces
+ """
+ num = (yn/y1)**(1/alpha)
+
+ return num
+
+def superposition(xi_b, xi_c, N_cases, y1_cases, s0 = 0):
+
+
+ # Check the array equal number of items
+ if not (np.size(xi_b) == np.size(xi_c) and np.size(xi_b) == np.size(N_cases)):
+ raise ValueError('Size for input should be identical')
+ print(s0)
+ yn = np.zeros(np.size(xi_b))
+ for i, (b, c, N, y1) in enumerate(zip(xi_b, xi_c, N_cases, y1_cases)):
+
+ if i==0:
+
+ em = KlinkvortModel(b, c, N)
+ _, _, alpha = em.alpha()
+ if s0 == 0:
+ yn[i] = em.evoluation(y1)
+ else:
+ N_eq = (s0/y1)**(1/alpha)
+ if N_eq < 1e8:
+
+# N_eq = (s0/y1)**(1/alpha)
+ yn[i] = y1*(N+N_eq)**alpha
+ else:
+ yn[i] = s0
+
+
+
+ else:
+ em = KlinkvortModel(b, c, N)
+ _, _, alpha = em.alpha()
+ N_eq = (yn[i-1]/y1)**(1/alpha)
+ if N_eq < 1e8:
+# N_eq = (yn[i-1]/y1)**(1/alpha)
+ yn[i] = y1*(N_eq + N)**alpha
+ else:
+ yn[i] = yn[i-1]
+
+ return yn
+
+def alpha_super(xi_b, xi_c, N_cases):
+ # Check the alpha value
+ if not (np.size(xi_b) == np.size(xi_c) and np.size(xi_b) == np.size(N_cases)):
+ raise ValueError('Size for input should be identical')
+
+ alpha_list = np.zeros(np.size(xi_b))
+ for i, (b, c, N) in enumerate(zip(xi_b, xi_c, N_cases)):
+
+ em = KlinkvortModel(b, c, N)
+ _, _, alpha = em.alpha()
+ alpha_list[i] = alpha
+
+ return alpha_list
+
+
+
+
\ No newline at end of file
diff --git a/wetb/soil/LeblancModel.py b/wetb/soil/LeblancModel.py
new file mode 100644
index 00000000..a6dbb51c
--- /dev/null
+++ b/wetb/soil/LeblancModel.py
@@ -0,0 +1,189 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Mar 25 16:15:25 2018
+
+@author: shfe@dtu.dk
+
+The model proposed by Leblanc
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+plt.close('all')
+
+class LeblancModel(object):
+
+
+ def __init__(self, xi_b, xi_c, N):
+ """
+ Initialize the input parameters
+ """
+
+ self.xi_b = xi_b
+ self.xi_c = xi_c
+ self.N = N
+ self.alpha = 0.31
+
+ def coeff(self, rd = 0.38):
+ """
+ Calculate the parameter alpha
+ """
+ # Calculate the Tb, Tc
+ if not 0 <= self.xi_b <= 1:
+ raise ValueError('The maximum load ratio should between 0 to 1')
+
+ if not -1 <= self.xi_c <= 1:
+ raise ValueError('The cyclic load ratio should between -1 to 1')
+
+ if rd == 0.38:
+ if self.xi_b <= 0.051:
+ Tb = 0
+ else:
+ Tb = 0.4238*self.xi_b - 0.0217
+
+ elif rd == 0.04:
+ if self.xi_b <= 0.1461:
+ Tb = 0
+ else:
+ Tb = 0.3087*self.xi_b - 0.0451
+
+ else:
+ raise ValueError('The relative density value is not validated yet')
+
+
+ if -1 <= self.xi_c < -0.65:
+ Tc = 13.71*self.xi_c + 13.71
+
+ elif -0.65 <= self.xi_c < 0:
+ Tc = -5.54*self.xi_c + 1.2
+
+ else:
+ Tc = -1.2*self.xi_c + 1.2
+
+ return Tb, Tc
+
+ def evoluation(self, y1):
+ """
+ Evoluation model for cyclic accumulated rotation
+ """
+
+ Tb, Tc = self.coeff()
+ yn = y1 * Tb * Tc * self.N**0.31
+
+ return yn
+
+ def equivalentnum(self, y1, yn):
+ """
+ Equivalent number for cylces
+ """
+ Tb, Tc = self.coeff()
+ num = (yn/(y1*Tb*Tc))**(1/self.alpha)
+
+ return num
+
+def superposition(xi_b, xi_c, N_cases, y1_cases, s0 = 0, y0 = 0):
+
+
+ # Check the array equal number of items
+ if not (np.size(xi_b) == np.size(xi_c) and np.size(xi_b) == np.size(N_cases)):
+ raise ValueError('Size for input should be identical')
+# print(s0)
+ dyn = np.zeros(np.size(xi_b))
+ yn = np.zeros(np.size(xi_b))
+ y_static = np.zeros(np.size(xi_b))
+ for i, (b, c, N, y1) in enumerate(zip(xi_b, xi_c, N_cases, y1_cases)):
+
+ if i==0:
+
+ em = LeblancModel(b, c, N)
+ Tb, Tc = em.coeff()
+ if s0 == 0:
+ dyn[i] = em.evoluation(y1)
+ else:
+ N_eq = (s0/(y1*Tb*Tc))**(1/em.alpha)
+ if N_eq < 1e8:
+ dyn[i] = (y1*Tb*Tc)*(N+N_eq)**em.alpha
+ else:
+ dyn[i] = s0
+
+ y_static[i] = max(y1, y0)
+
+ else:
+ em = LeblancModel(b, c, N)
+ Tb, Tc = em.coeff()
+ N_eq = (dyn[i-1]/(y1*Tb*Tc))**(1/em.alpha)
+ if N_eq < 1e8:
+ dyn[i] = (y1*Tb*Tc)*(N_eq + N)**em.alpha
+ else:
+ dyn[i] = dyn[i-1]
+
+ y_static[i] = max(y_static[i-1], y1)
+
+ yn[i] = dyn[i] + y_static[i]
+
+ return dyn, y_static, yn
+
+def coeff_super(xi_b, xi_c, N_cases):
+ # Check the alpha value
+ if not (np.size(xi_b) == np.size(xi_c) and np.size(xi_b) == np.size(N_cases)):
+ raise ValueError('Size for input should be identical')
+
+ Tb_list = np.zeros(np.size(xi_b))
+ Tc_list = np.zeros(np.size(xi_c))
+ for i, (b, c, N) in enumerate(zip(xi_b, xi_c, N_cases)):
+
+ em = LeblancModel(b, c, N)
+ Tb, Tc = em.coeff()
+ Tb_list[i] = Tb
+ Tc_list[i] = Tc
+
+ return Tb_list, Tc_list
+
+
+def moment_ult(L, D, gamma):
+ """
+ Ultimate static moment
+
+ Parameters
+ ----------------
+ L -- Pile length
+ D -- Pile diameter
+ gamma -- Submerged unit weight
+
+ Output
+ ----------------
+ m_ult -- Ultimate static moment
+ """
+
+
+ m_ult = L**3*D*gamma*0.6
+
+ return m_ult
+
+def rot_moment(m, L, D, gamma, pa = 101, rd = 0.04):
+ """
+ Ultimate static moment
+
+ Parameters
+ ----------------
+ L -- Pile length
+ D -- Pile diameter
+ gamma -- Submerged unit weight
+ m -- Moment
+
+ Output
+ ----------------
+ theta - Rotational angle in radians
+ """
+
+ if rd == 0.04:
+ m_ult = L**3*D*gamma*0.6/1e3 # unit: MN
+ theta_norm = 0.038 * (m/m_ult)**2.33
+ theta = theta_norm*np.sqrt((L*gamma)/pa)
+
+ else:
+ m_ult = L**3*D*gamma*1.2/1e3 # unit: MN
+ theta_norm = 0.042 * (m/m_ult)**1.92
+ theta = theta_norm*np.sqrt((L*gamma)/pa)
+
+ return theta
+
diff --git a/wetb/soil/__init__.py b/wetb/soil/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/wetb/soil/hyperbolic_model.py b/wetb/soil/hyperbolic_model.py
new file mode 100644
index 00000000..f08e8203
--- /dev/null
+++ b/wetb/soil/hyperbolic_model.py
@@ -0,0 +1,48 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Apr 22 17:06:06 2018
+
+@author: shfe
+
+Hyperbolic model used for soil load displacement behaviour
+"""
+
+import numpy as np
+
+def hyperbolic_load(y, pu, a, b):
+ """
+ Static hyperbolic model developed by Lee
+
+ Parameters
+ ----------------
+ y -- displacement
+ pu -- ultimate load
+ a, b -- coefficients
+
+ Output
+ ----------------
+ p -- soil resistence pressure
+ """
+
+ p = y/(a + b*y) * pu
+
+ return p
+
+def hyperbolic_disp(p, pu, a, b):
+ """
+ Static hyperbolic model developed by Lee
+
+ Parameters
+ ----------------
+ p -- load
+ pu -- ultimate load
+ a, b -- coefficients
+
+ Output
+ ----------------
+ y -- displacement
+ """
+
+ y = p * a/(pu - p*b)
+
+ return y
\ No newline at end of file
diff --git a/wetb/soil/py_model.py b/wetb/soil/py_model.py
new file mode 100644
index 00000000..9d51b53c
--- /dev/null
+++ b/wetb/soil/py_model.py
@@ -0,0 +1,167 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Mar 25 16:15:25 2018
+
+@author: shfe@dtu.dk
+
+Static p-y soil model and cyclic degradation model
+"""
+
+import numpy as np
+
+def py_static(y, pu, yc):
+ """
+ Static p-y model developed by Matlock
+
+ Parameters
+ ----------------
+ y -- displacement
+ pu -- ultimate stress
+ yc -- ultimate displacement
+
+ Output
+ ----------------
+ p -- soil resistence pressure
+ """
+
+ if y/yc <= 8:
+ p = 0.5*pu*(y/yc)**(1/3)
+
+ else:
+ p = pu
+
+ return p
+
+def strain_static(p, pu, yc):
+ """
+ Static strain using p-y model developed by Matlock
+
+ Parameters
+ ----------------
+ p -- stress
+ pu -- ultimate stress
+ yc -- ultimate displacement
+
+ Output
+ ----------------
+ y -- soil strain
+ """
+ if p <= pu:
+ y = (p/(0.5*pu))**3*yc
+
+ else:
+ raise ValueError("The applied stress is LARGER than ultimate strength")
+
+ return y
+
+def load_unload_py(p, pu, yc):
+ """
+ load/unload module using p-y model developed by Matlock
+
+ Parameters
+ ----------------
+ p -- stress
+ pu -- ultimate stress
+ yc -- ultimate displacement
+
+ Output
+ ----------------
+ y -- soil strain
+ y_cum -- cumulative soil displacement
+ """
+
+ y = strain_static(p, pu, yc)
+
+ if p <= 0.5*pu: # elastic
+ y_cum = 0
+
+ else:
+ y_cum = y - p/(0.5*pu/yc)
+
+ return y, y_cum
+
+
+
+
+def degra_coeff(N, y1, b):
+ """
+ Degradation factor for p-y model
+
+ Parameters
+ ----------------
+ N -- number of cycle
+ y1 -- displacement predicted by static p-y curve
+ b -- pile diameter
+
+ Output
+ ----------------
+ lambda -- Degradation factor
+ """
+ # check
+
+ lambda_N = np.abs(y1/(0.2*b)*np.log(N))
+
+ try:
+ lambda_N <= 1
+ except:
+ print("soil must be degraded!!!")
+
+ return lambda_N
+
+def py_cyclic(y, pu, yc, N, y1, b):
+ """
+ Degradation p-y model with cyclic load
+
+ Parameters
+ ----------------
+ y -- displacement
+ pu -- ultimate stress
+ yc -- ultimate displacement
+ N -- number of cycle
+ y1 -- displacement predicted by static p-y curve
+ b -- pile diameter
+
+ Output
+ ----------------
+ p -- degraded soil resistence pressure
+ """
+
+ lambda_N = degra_coeff(N, y1, b)
+ p = py_static(y, pu, yc) * (1 - lambda_N)
+
+ return p
+
+def py_parcel(load_parcel, cycle_parcel, pu, yc, b):
+ """
+ Degradation p-y model after a parcel of cyclic load
+
+ Parameters
+ ----------------
+ load_parcel -- load parcel
+ cycle_parcel -- number of cycle of load parcel
+ yc -- ultimate displacement
+ b -- pile diameter
+
+ Output
+ ----------------
+ p -- degraded soil resistence pressure
+ """
+
+ pu_N = np.zeros(np.size(load_parcel)+1)
+ pu_N[0] = pu
+ lambda_N = np.zeros(np.size(load_parcel))
+ y_N = np.zeros(np.size(load_parcel))
+
+ for i in np.arange(np.size(load_parcel)):
+ load_i = load_parcel[i]
+ cycle_i = cycle_parcel[i]
+ if np.isnan(load_i):
+ y_i = 0
+ else:
+ y_i = strain_static(load_i, pu_N[i], yc)
+ y_N[i] = y_i
+ lambda_i = degra_coeff(cycle_i, y_i, b)
+ lambda_N[i] = lambda_i
+ pu_N[i+1] = pu_N[i] * (1-lambda_i)
+
+ return y_N, lambda_N, pu_N
diff --git a/wetb/wave/__init__.py b/wetb/wave/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/wetb/wave/dispersion.py b/wetb/wave/dispersion.py
new file mode 100644
index 00000000..467f7063
--- /dev/null
+++ b/wetb/wave/dispersion.py
@@ -0,0 +1,157 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Nov 6 21:06:05 2017
+
+@author: shfe
+
+Note: The module contains an assortment of tools for working with linear wave theory calculations.
+"""
+# py2/3 compatible:
+from __future__ import absolute_import, division, print_function, unicode_literals
+import numpy as np
+import matplotlib.pyplot as plt
+plt.close('all')
+
+g = 9.806 # gravitational acceleration in m^2/s
+
+def wave_number(g, omega, h):
+ """
+ Computes the wave number for given frequency and water depth
+ (linear dispersion relationship)
+ :param omega: -- wave frequency
+ :param g: -- gravitational acceleration (defaults to 9.806 m/s^2)
+ :param h: -- the water depth
+ :returns k: the wave number
+ """
+
+ p = omega**2 * h / g
+ q = dispersion(p)
+ k = q * omega**2 / g
+
+ return k
+
+
+def frequency(g, k, h):
+ """
+ Computes the frequency for a given wave number and water depth
+ (linear dispersion relationship)
+ :param k: the wave number
+ :param g: -- gravitational acceleration
+ :param h: -- the water depth
+ :returns omega: -- wave frequency
+ """
+
+ return np.sqrt(g * k * np.tanh(k * h))
+
+
+def dispersion(p, tol=1e-14, max_iter=100):
+ """
+ The linear dispersion relation in non-dimensional form:
+ finds q, given p
+ q = gk/omega^2 non-d wave number
+ p = omega^2 h / g non-d water depth
+ Starts with the Fenton and McKee approximation, then iterates with Newton's
+ method until accurate to within tol.
+ :param p: non-dimensional water depth
+ :param tol=1e-14: acceptable tolerance
+ :param max_iter: maximum number of iterations to accept
+ (it SHOULD converge in well less than 100)
+ """
+ if p <= 0.0:
+ raise ValueError("Non dimensional water depth d must be >= 0.0")
+ # First guess (from Fenton and McKee):
+ q = np.tanh(p ** 0.75) ** (-2.0 / 3.0)
+
+ iter = 0
+ f = q * np.tanh(q * p) - 1
+ while abs(f) > tol:
+ qp = q * p
+ fp = qp / (np.cosh(qp) ** 2) + np.tanh(qp)
+ q = q - f / fp
+ f = q * np.tanh(q * p) - 1
+ iter += 1
+ if iter > max_iter:
+ raise RuntimeError("Maximum number of iterations reached in dispersion()")
+ return q
+
+
+def max_u(a, omega, g, h, z=None):
+ """
+ Compute the maximum Eulerian horizontal velocity at a given depth
+ :param a: -- wave amplitude (1/2 the height)
+ :param omega: -- wave frequency
+ :param g: -- gravitational acceleration (defaults to 9.806 m/s^2)
+ :param h: -- the water depth
+ :param z=None: -- the depth at which to compute the maximum velocity.
+ if None, then h is used (the bottom)
+ """
+
+ z = h if z is None else z
+ k = wave_number(g, omega, h)
+ u = a * omega * (np.cosh(k * (h + z)) / np.sinh(k * h))
+
+ return u
+
+
+def amp_scale_at_depth(g, omega, h, z):
+ """
+ Compute the scale factor of the orbital amplitude at the given depth
+ :param g: -- gravitational acceleration
+ :param omega: -- the wave frequency
+ :param h: -- the water depth
+ :param z: -- the depth at which to compute the scale factor
+ """
+
+ k = wave_number(g, omega, h)
+
+ return np.cosh(k * (h + z)) / np.cosh(k * (h))
+
+
+def celerity(k, h, g=g):
+ """
+ Compute the celerity (wave speed, phase speed) for a given wave number and depth
+ :param k: -- the wave number
+ :param h: -- the water depth
+ :param g=g: -- gravitational acceleration (defaults to 9.806 m/s^2)
+ """
+
+ C = np.sqrt(g / k * np.tanh(k * h))
+
+ return C
+
+
+def group_speed(k, h, g=g):
+ """
+ Compute the group speed for a given wave number and depth
+ :param k: -- the wave number
+ :param h: -- the water depth
+ :param g=g: -- gravitational acceleration (defaults to 9.806 m/s^2)
+ """
+
+ n = 1.0 / 2 * (1 + (2 * k * h / np.sinh(2 * k * h)))
+ Cg = n * celerity(k, h, g)
+
+ return Cg
+
+
+def shoaling_coeff(omega, g, h0, h2):
+ """
+ Compute the shoaling coeff for two depths: ho and h2.
+ The shoaling coeff is the ratio of wave height (H2) at a particular
+ point of interest to the original or deep water wave height (H0).
+ Pass in h0 = None for deep water
+ :param omega: -- the wave frequency
+ :param g: -- gravitational acceleration
+ :param h0: -- the initial water depth
+ :param h2: -- the depth at which to compute the shoaling coeff
+ """
+
+ k2 = wave_number(g, omega, h2)
+ Cg2 = group_speed(k2, h2, g)
+ if h0 is not None:
+ k0 = wave_number(g, omega, h0)
+ Cg0 = group_speed(k0, h0, g)
+ Ks = np.sqrt(Cg0 / Cg2)
+ return Ks
+ else: # Deep water
+ return np.sqrt((g / (2 * omega)) / Cg2)
diff --git a/wetb/wave/numba_misc.py b/wetb/wave/numba_misc.py
new file mode 100644
index 00000000..0e6872dc
--- /dev/null
+++ b/wetb/wave/numba_misc.py
@@ -0,0 +1,569 @@
+"""
+Created on 6. okt. 2016
+
+@author: pab
+"""
+from __future__ import absolute_import, division
+from numba import jit, float64, int64, int32, int8, void
+import numpy as np
+
+
+@jit(int64(int64[:], int8[:]))
+def _findcross(ind, y):
+ ix, dcross, start, v = 0, 0, 0, 0
+ n = len(y)
+ if y[0] < v:
+ dcross = -1 # first is a up-crossing
+ elif y[0] > v:
+ dcross = 1 # first is a down-crossing
+ elif y[0] == v:
+ # Find out what type of crossing we have next time..
+ for i in range(1, n):
+ start = i
+ if y[i] < v:
+ ind[ix] = i - 1 # first crossing is a down crossing
+ ix += 1
+ dcross = -1 # The next crossing is a up-crossing
+ break
+ elif y[i] > v:
+ ind[ix] = i - 1 # first crossing is a up-crossing
+ ix += 1
+ dcross = 1 # The next crossing is a down-crossing
+ break
+
+ for i in range(start, n - 1):
+ if ((dcross == -1 and y[i] <= v and v < y[i + 1]) or
+ (dcross == 1 and v <= y[i] and y[i + 1] < v)):
+
+ ind[ix] = i
+ ix += 1
+ dcross = -dcross
+ return ix
+
+
+def findcross(xn):
+ """Return indices to zero up and downcrossings of a vector
+ """
+ ind = np.empty(len(xn), dtype=np.int64)
+ m = _findcross(ind, xn)
+ return ind[:m]
+
+
+def _make_findrfc(cmp1, cmp2):
+
+ @jit(int64(int64[:], float64[:], float64), nopython=True)
+ def local_findrfc(t, y, h):
+ # cmp1, cmp2 = (a_le_b, a_lt_b) if method==0 else (a_lt_b, a_le_b)
+
+ n = len(y)
+ j, t0, z0 = 0, 0, 0
+ y0 = y[t0]
+ # The rainflow filter
+ for ti in range(1, n):
+ fpi = y0 + h
+ fmi = y0 - h
+ yi = y[ti]
+
+ if z0 == 0:
+ if cmp1(yi, fmi):
+ z1 = -1
+ elif cmp1(fpi, yi):
+ z1 = +1
+ else:
+ z1 = 0
+ t1, y1 = (t0, y0) if z1 == 0 else (ti, yi)
+ else:
+ if (((z0 == +1) and cmp1(yi, fmi)) or
+ ((z0 == -1) and cmp2(yi, fpi))):
+ z1 = -1
+ elif (((z0 == +1) and cmp2(fmi, yi)) or
+ ((z0 == -1) and cmp1(fpi, yi))):
+ z1 = +1
+ else:
+ raise ValueError
+ # warnings.warn('Something wrong, i={}'.format(tim1))
+
+ # Update y1
+ if z1 != z0:
+ t1, y1 = ti, yi
+ elif z1 == -1:
+ t1, y1 = (t0, y0) if y0 < yi else (ti, yi)
+ elif z1 == +1:
+ t1, y1 = (t0, y0) if y0 > yi else (ti, yi)
+
+ # Update y if y0 is a turning point
+ if abs(z0 - z1) == 2:
+ j += 1
+ t[j] = t0
+
+ # Update t0, y0, z0
+ t0, y0, z0 = t1, y1, z1
+ # end
+
+ # Update y if last y0 is greater than (or equal) threshold
+ if cmp2(h, abs(y0 - y[t[j]])):
+ j += 1
+ t[j] = t0
+ return j + 1
+ return local_findrfc
+
+
+@jit(int32(float64, float64), nopython=True)
+def a_le_b(a, b):
+ return a <= b
+
+
+@jit(int32(float64, float64), nopython=True)
+def a_lt_b(a, b):
+ return a < b
+
+
+_findrfc_le = _make_findrfc(a_le_b, a_lt_b)
+_findrfc_lt = _make_findrfc(a_lt_b, a_le_b)
+
+
+@jit(int64(int64[:], float64[:], float64), nopython=True)
+def _findrfc(ind, y, h):
+ n = len(y)
+ t_start = 0
+ nc = n // 2 - 1
+ ix = 0
+ for i in range(nc):
+ Tmi = t_start + 2 * i
+ Tpl = t_start + 2 * i + 2
+ xminus = y[2 * i]
+ xplus = y[2 * i + 2]
+
+ if(i != 0):
+ j = i - 1
+ while ((j >= 0) and (y[2 * j + 1] <= y[2 * i + 1])):
+ if (y[2 * j] < xminus):
+ xminus = y[2 * j]
+ Tmi = t_start + 2 * j
+ j -= 1
+ if (xminus >= xplus):
+ if (y[2 * i + 1] - xminus >= h):
+ ind[ix] = Tmi
+ ix += 1
+ ind[ix] = (t_start + 2 * i + 1)
+ ix += 1
+ # goto L180 continue
+ else:
+ j = i + 1
+ while (j < nc):
+ if (y[2 * j + 1] >= y[2 * i + 1]):
+ break # goto L170
+ if((y[2 * j + 2] <= xplus)):
+ xplus = y[2 * j + 2]
+ Tpl = (t_start + 2 * j + 2)
+ j += 1
+ else:
+ if ((y[2 * i + 1] - xminus) >= h):
+ ind[ix] = Tmi
+ ix += 1
+ ind[ix] = (t_start + 2 * i + 1)
+ ix += 1
+ # iy = i
+ continue
+
+ # goto L180
+ # L170:
+ if (xplus <= xminus):
+ if ((y[2 * i + 1] - xminus) >= h):
+ ind[ix] = Tmi
+ ix += 1
+ ind[ix] = (t_start + 2 * i + 1)
+ ix += 1
+ elif ((y[2 * i + 1] - xplus) >= h):
+ ind[ix] = (t_start + 2 * i + 1)
+ ix += 1
+ ind[ix] = Tpl
+ ix += 1
+
+ # L180:
+ # iy=i
+ # /* for i */
+ return ix
+
+
+def findrfc(y, h, method=0):
+ n = len(y)
+ t = np.zeros(n, dtype=np.int64)
+ findrfc_ = [_findrfc_le, _findrfc_lt, _findrfc][method]
+ m = findrfc_(t, y, h)
+ return t[:m]
+
+
+@jit(void(float64[:], float64[:], float64[:], float64[:],
+ float64[:], float64[:], float64, float64,
+ int32, int32, int32, int32), nopython=True)
+def _finite_water_disufq(rvec, ivec, rA, iA, w, kw, h, g, nmin, nmax, m, n):
+ # kfact is set to 2 in order to exploit the symmetry.
+ # If you set kfact to 1, you must uncomment all statements
+ # including the expressions: rvec[iz2], rvec[iv2], ivec[iz2] and ivec[iv2].
+
+ kfact = 2.0
+ for ix in range(nmin - 1, nmax):
+ # for (ix = nmin-1;ix<nmax;ix++) {
+ kw1 = kw[ix]
+ w1 = w[ix]
+ tmp1 = np.tanh(kw1 * h)
+ # Cg, wave group velocity
+ Cg = 0.5 * g * (tmp1 + kw1 * h * (1.0 - tmp1 * tmp1)) / w1 # OK
+ tmp1 = 0.5 * g * (kw1 / w1) * (kw1 / w1)
+ tmp2 = 0.5 * w1 * w1 / g
+ tmp3 = g * kw1 / (w1 * Cg)
+ tmp4 = kw1 / np.sinh(2.0 * kw1 * h) if kw1 * h < 300.0 else 0.0
+
+ # Difference frequency effects finite water depth
+ Edij = (tmp1 - tmp2 + tmp3) / (1.0 - g * h / (Cg * Cg)) - tmp4 # OK
+
+ # Sum frequency effects finite water depth
+ Epij = (3.0 * (tmp1 - tmp2) /
+ (1.0 - tmp1 / kw1 * np.tanh(2.0 * kw1 * h)) +
+ 3.0 * tmp2 - tmp1) # OK
+ # printf("Edij = %f Epij = %f \n", Edij,Epij);
+
+ ixi = ix * m
+ iz1 = 2 * ixi
+ # iz2 = n*m-ixi;
+ for i in range(m):
+ rrA = rA[ixi] * rA[ixi]
+ iiA = iA[ixi] * iA[ixi]
+ riA = rA[ixi] * iA[ixi]
+
+ # Sum frequency effects along the diagonal
+ rvec[iz1] += kfact * (rrA - iiA) * Epij
+ ivec[iz1] += kfact * 2.0 * riA * Epij
+ # rvec[iz2] += kfact*(rrA-iiA)*Epij;
+ # ivec[iz2] -= kfact*2.0*riA*Epij;
+ # iz2++;
+
+ # Difference frequency effects along the diagonal
+ # are only contributing to the mean
+ rvec[i] += 2.0 * (rrA + iiA) * Edij
+ ixi += 1
+ iz1 += 1
+ # }
+ for jy in range(ix + 1, nmax):
+ # w1 = w[ix];
+ # kw1 = kw[ix];
+ w2 = w[jy]
+ kw2 = kw[jy]
+ tmp1 = g * (kw1 / w1) * (kw2 / w2)
+ tmp2 = 0.5 / g * (w1 * w1 + w2 * w2 + w1 * w2)
+ tmp3 = 0.5 * g * \
+ (w1 * kw2 * kw2 + w2 * kw1 * kw1) / (w1 * w2 * (w1 + w2))
+ tmp4 = (1 - g * (kw1 + kw2) / (w1 + w2) / (w1 + w2) *
+ np.tanh((kw1 + kw2) * h))
+ Epij = (tmp1 - tmp2 + tmp3) / tmp4 + tmp2 - 0.5 * tmp1 # OK */
+
+ tmp2 = 0.5 / g * (w1 * w1 + w2 * w2 - w1 * w2) # OK*/
+ tmp3 = -0.5 * g * \
+ (w1 * kw2 * kw2 - w2 * kw1 * kw1) / (w1 * w2 * (w1 - w2))
+ tmp4 = (1.0 - g * (kw1 - kw2) / (w1 - w2) /
+ (w1 - w2) * np.tanh((kw1 - kw2) * h))
+ Edij = (tmp1 - tmp2 + tmp3) / tmp4 + tmp2 - 0.5 * tmp1 # OK */
+ # printf("Edij = %f Epij = %f \n", Edij,Epij);
+
+ ixi = ix * m
+ jyi = jy * m
+ iz1 = ixi + jyi
+ iv1 = jyi - ixi
+ # iz2 = (n*m-iz1)
+ # iv2 = n*m-iv1
+ for i in range(m):
+ # for (i=0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
+ rrA = rA[ixi] * rA[jyi] # rrA = rA[i][ix]*rA[i][jy];
+ iiA = iA[ixi] * iA[jyi] # iiA = iA[i][ix]*iA[i][jy];
+ riA = rA[ixi] * iA[jyi] # riA = rA[i][ix]*iA[i][jy];
+ irA = iA[ixi] * rA[jyi] # irA = iA[i][ix]*rA[i][jy];
+
+ # Sum frequency effects */
+ tmp1 = kfact * 2.0 * (rrA - iiA) * Epij
+ tmp2 = kfact * 2.0 * (riA + irA) * Epij
+ rvec[iz1] += tmp1 # rvec[i][jy+ix] += tmp1
+ ivec[iz1] += tmp2 # ivec[i][jy+ix] += tmp2
+ # rvec[iz2] += tmp1 # rvec[i][n*m-(jy+ix)] += tmp1
+ # ivec[iz2] -= tmp2 # ivec[i][n*m-(jy+ix)] -= tmp2
+ # iz2++
+
+ # Difference frequency effects */
+ tmp1 = kfact * 2.0 * (rrA + iiA) * Edij
+ tmp2 = kfact * 2.0 * (riA - irA) * Edij
+ rvec[iv1] += tmp1 # rvec[i][jy-ix] += tmp1
+ ivec[iv1] += tmp2 # ivec[i][jy-ix] -= tmp2
+
+ # rvec[iv2] += tmp1
+ # ivec[iv2] -= tmp2
+ # iv2 += 1
+ ixi += 1
+ jyi += 1
+ iz1 += 1
+ iv1 += 1
+
+
+@jit(void(float64[:], float64[:], float64[:], float64[:],
+ float64[:], float64[:], float64, float64,
+ int32, int32, int32, int32), nopython=True)
+def _deep_water_disufq(rvec, ivec, rA, iA, w, kw, h, g, nmin, nmax, m, n):
+ # kfact is set to 2 in order to exploit the symmetry.
+ # If you set kfact to 1, you must uncomment all statements
+ # including the expressions: rvec[iz2], rvec[iv2], ivec[iz2] and ivec[iv2].
+
+ kfact = 2.0
+ for ix in range(nmin - 1, nmax):
+ ixi = ix * m
+ iz1 = 2 * ixi
+ iz2 = n * m - ixi
+ kw1 = kw[ix]
+ Epij = kw1
+ for _i in range(m):
+ rrA = rA[ixi] * rA[ixi]
+ iiA = iA[ixi] * iA[ixi]
+ riA = rA[ixi] * iA[ixi]
+
+ # Sum frequency effects along the diagonal
+ tmp1 = kfact * (rrA - iiA) * Epij
+ tmp2 = kfact * 2.0 * riA * Epij
+ rvec[iz1] += tmp1
+ ivec[iz1] += tmp2
+ ixi += 1
+ iz1 += 1
+ # rvec[iz2] += tmp1
+ # ivec[iz2] -= tmp2
+ iz2 += 1
+
+ # Difference frequency effects are zero along the diagonal
+ # and are thus not contributing to the mean.
+ for jy in range(ix + 1, nmax):
+ kw2 = kw[jy]
+ Epij = 0.5 * (kw2 + kw1)
+ Edij = -0.5 * (kw2 - kw1)
+ # printf("Edij = %f Epij = %f \n", Edij,Epij);
+
+ ixi = ix * m
+ jyi = jy * m
+ iz1 = ixi + jyi
+ iv1 = jyi - ixi
+ iz2 = (n * m - iz1)
+ iv2 = (n * m - iv1)
+ for _i in range(m):
+ rrA = rA[ixi] * rA[jyi] # rrA = rA[i][ix]*rA[i][jy]
+ iiA = iA[ixi] * iA[jyi] # iiA = iA[i][ix]*iA[i][jy]
+ riA = rA[ixi] * iA[jyi] # riA = rA[i][ix]*iA[i][jy]
+ irA = iA[ixi] * rA[jyi] # irA = iA[i][ix]*rA[i][jy]
+
+ # Sum frequency effects
+ tmp1 = kfact * 2.0 * (rrA - iiA) * Epij
+ tmp2 = kfact * 2.0 * (riA + irA) * Epij
+ rvec[iz1] += tmp1 # rvec[i][ix+jy] += tmp1
+ ivec[iz1] += tmp2 # ivec[i][ix+jy] += tmp2
+ # rvec[iz2] += tmp1 # rvec[i][n*m-(ix+jy)] += tmp1
+ # ivec[iz2] -= tmp2 # ivec[i][n*m-(ix+jy)] -= tmp2
+ iz2 += 1
+
+ # Difference frequency effects */
+ tmp1 = kfact * 2.0 * (rrA + iiA) * Edij
+ tmp2 = kfact * 2.0 * (riA - irA) * Edij
+
+ rvec[iv1] += tmp1 # rvec[i][jy-ix] += tmp1
+ ivec[iv1] += tmp2 # ivec[i][jy-ix] += tmp2
+
+ # rvec[iv2] += tmp1 # rvec[i][n*m-(jy-ix)] += tmp1
+ # ivec[iv2] -= tmp2 # ivec[i][n*m-(jy-ix)] -= tmp2
+ iv2 += 1
+ ixi += 1
+ jyi += 1
+ iz1 += 1
+ iv1 += 1
+
+
+def disufq(rA, iA, w, kw, h, g, nmin, nmax, m, n):
+ """
+ DISUFQ Is an internal function to spec2nlsdat
+
+ Parameters
+ ----------
+ rA, iA = real and imaginary parts of the amplitudes (size m X n).
+ w = vector with angular frequencies (w>=0)
+ kw = vector with wavenumbers (kw>=0)
+ h = water depth (h >=0)
+ g = constant acceleration of gravity
+ nmin = minimum index where rA(:,nmin) and iA(:,nmin) is
+ greater than zero.
+ nmax = maximum index where rA(:,nmax) and iA(:,nmax) is
+ greater than zero.
+ m = size(rA,1),size(iA,1)
+ n = size(rA,2),size(iA,2), or size(rvec,2),size(ivec,2)
+
+ returns
+ -------
+ rvec, ivec = real and imaginary parts of the resultant (size m X n).
+
+ DISUFQ returns the summation of difference frequency and sum
+ frequency effects in the vector vec = rvec +sqrt(-1)*ivec.
+ The 2'nd order contribution to the Stokes wave is then calculated by
+ a simple 1D Fourier transform, real(FFT(vec)).
+
+ """
+ rvec = np.zeros(n * m)
+ ivec = np.zeros(n * m)
+
+ if h > 10000: # { /* deep water /Inifinite water depth */
+ _deep_water_disufq(rvec, ivec, rA, iA, w, kw, h, g, nmin, nmax, m, n)
+ else:
+ _finite_water_disufq(rvec, ivec, rA, iA, w, kw, h, g, nmin, nmax, m, n)
+ return rvec, ivec
+
+
+@jit(int32[:](float64[:], float64[:], float64[:, :]))
+def _findrfc3_astm(array_ext, a, array_out):
+ """
+ Rain flow without time analysis
+
+ Return [ampl ampl_mean nr_of_cycle]
+
+ By Adam Nieslony
+ Visit the MATLAB Central File Exchange for latest version
+ http://www.mathworks.com/matlabcentral/fileexchange/3026
+ """
+ n = len(array_ext)
+ po = 0
+ # The original rainflow counting by Nieslony, unchanged
+ j = -1
+ c_nr1 = 1
+ for i in range(n):
+ j += 1
+ a[j] = array_ext[i]
+ while j >= 2 and abs(a[j - 1] - a[j - 2]) <= abs(a[j] - a[j - 1]):
+ ampl = abs((a[j - 1] - a[j - 2]) / 2)
+ mean = (a[j - 1] + a[j - 2]) / 2
+ if j == 2:
+ a[0] = a[1]
+ a[1] = a[2]
+ j = 1
+ if (ampl > 0):
+ array_out[po, :] = (ampl, mean, 0.5)
+ po += 1
+ else:
+ a[j - 2] = a[j]
+ j = j - 2
+ if (ampl > 0):
+ array_out[po, :] = (ampl, mean, 1.0)
+ po += 1
+ c_nr1 += 1
+
+ c_nr2 = 1
+ for i in range(j):
+ ampl = abs(a[i] - a[i + 1]) / 2
+ mean = (a[i] + a[i + 1]) / 2
+ if (ampl > 0):
+ array_out[po, :] = (ampl, mean, 0.5)
+ po += 1
+ c_nr2 += 1
+ return c_nr1, c_nr2
+
+
+@jit(int32[:](float64[:], float64[:], float64[:], float64[:], float64[:, :]))
+def _findrfc5_astm(array_ext, array_t, a, t, array_out):
+ """
+ Rain flow with time analysis
+
+ returns
+ [ampl ampl_mean nr_of_cycle cycle_begin_time cycle_period_time]
+
+ By Adam Nieslony
+ Visit the MATLAB Central File Exchange for latest version
+ http://www.mathworks.com/matlabcentral/fileexchange/3026
+ """
+ n = len(array_ext)
+ po = 0
+ # The original rainflow counting by Nieslony, unchanged
+ j = -1
+ c_nr1 = 1
+ for i in range(n):
+ j += 1
+ a[j] = array_ext[i]
+ t[j] = array_t[i]
+ while (j >= 2) and (abs(a[j - 1] - a[j - 2]) <= abs(a[j] - a[j - 1])):
+ ampl = abs((a[j - 1] - a[j - 2]) / 2)
+ mean = (a[j - 1] + a[j - 2]) / 2
+ period = (t[j - 1] - t[j - 2]) * 2
+ atime = t[j - 2]
+ if j == 2:
+ a[0] = a[1]
+ a[1] = a[2]
+ t[0] = t[1]
+ t[1] = t[2]
+ j = 1
+ if (ampl > 0):
+ array_out[po, :] = (ampl, mean, 0.5, atime, period)
+ po += 1
+ else:
+ a[j - 2] = a[j]
+ t[j - 2] = t[j]
+ j = j - 2
+ if (ampl > 0):
+ array_out[po, :] = (ampl, mean, 1.0, atime, period)
+ po += 1
+ c_nr1 += 1
+
+ c_nr2 = 1
+ for i in range(j):
+ # for (i=0; i<j; i++) {
+ ampl = abs(a[i] - a[i + 1]) / 2
+ mean = (a[i] + a[i + 1]) / 2
+ period = (t[i + 1] - t[i]) * 2
+ atime = t[i]
+ if (ampl > 0):
+ array_out[po, :] = (ampl, mean, 0.5, atime, period)
+ po += 1
+ c_nr2 += 1
+ return c_nr1, c_nr2
+
+
+def findrfc_astm(tp, t=None):
+ """
+ Return rainflow counted cycles
+
+ Nieslony's Matlab implementation of the ASTM standard practice for rainflow
+ counting ported to a Python C module.
+
+ Parameters
+ ----------
+ tp : array-like
+ vector of turningpoints (NB! Only values, not sampled times)
+ t : array-like
+ vector of sampled times
+
+ Returns
+ -------
+ sig_rfc : array-like
+ array of shape (n,3) or (n, 5) with:
+ sig_rfc[:,0] Cycles amplitude
+ sig_rfc[:,1] Cycles mean value
+ sig_rfc[:,2] Cycle type, half (=0.5) or full (=1.0)
+ sig_rfc[:,3] cycle_begin_time (only if t is given)
+ sig_rfc[:,4] cycle_period_time (only if t is given)
+ """
+
+ y1 = np.atleast_1d(tp).ravel()
+ n = len(y1)
+ a = np.zeros(n)
+ if t is None:
+ sig_rfc = np.zeros((n, 3))
+ cnr = _findrfc3_astm(y1, a, sig_rfc)
+ else:
+ t1 = np.atleast_1d(t).ravel()
+ sig_rfc = np.zeros((n, 5))
+ t2 = np.zeros(n)
+ cnr = _findrfc5_astm(y1, t1, a, t2, sig_rfc)
+ # the sig_rfc was constructed too big in rainflow.rf3, so
+ # reduce the sig_rfc array as done originally by a matlab mex c function
+ # n = len(sig_rfc)
+ return sig_rfc[:n - cnr[0]]
+
+
+if __name__ == '__main__':
+ pass
diff --git a/wetb/wave/spectrum.py b/wetb/wave/spectrum.py
new file mode 100644
index 00000000..128d2f9f
--- /dev/null
+++ b/wetb/wave/spectrum.py
@@ -0,0 +1,531 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Mar 2 13:33:21 2018
+
+@author: shfe
+"""
+
+from __future__ import absolute_import, division
+
+import warnings
+
+from scipy.stats import moment
+import scipy.integrate as integrate
+import scipy.special as sp
+import numpy as np
+from numpy import (atleast_1d, minimum, maximum, ones_like,
+ exp, log, sqrt, where, pi, isfinite, cosh, zeros_like, flatnonzero)
+
+import matplotlib.pyplot as plt
+
+
+def sech(x):
+ return 1.0 / cosh(x)
+
+class ModelSpectrum(object):
+ type = 'ModelSpectrum'
+
+ def __init__(self, Hm0=7.0, Tp=11.0, **kwds):
+ self.Hm0 = Hm0
+ self.Tp = Tp
+
+ def chk_seastate(self):
+ """ Check if seastate is valid
+ """
+
+ if self.Hm0 < 0:
+ raise ValueError('Hm0 can not be negative!')
+
+ if self.Tp <= 0:
+ raise ValueError('Tp must be positve!')
+
+ if self.Hm0 == 0.0:
+ warnings.warn('Hm0 is zero!')
+
+ self._chk_extra_param()
+
+ def _chk_extra_param(self):
+ pass
+
+def jonswap_peakfact(Hm0, Tp):
+ """ Jonswap peakedness factor, gamma, given Hm0 and Tp
+
+ Parameters
+ ----------
+ Hm0 : significant wave height [m].
+ Tp : peak period [s]
+
+ Returns
+ -------
+ gamma : Peakedness parameter of the JONSWAP spectrum
+
+ Details
+ -------
+ A standard value for GAMMA is 3.3. However, a more correct approach is
+ to relate GAMMA to Hm0 and Tp:
+ D = 0.036-0.0056*Tp/sqrt(Hm0)
+ gamma = exp(3.484*(1-0.1975*D*Tp**4/(Hm0**2)))
+ This parameterization is based on qualitative considerations of deep water
+ wave data from the North Sea, see Torsethaugen et. al. (1984)
+ Here GAMMA is limited to 1..7.
+
+ NOTE: The size of GAMMA is the common shape of Hm0 and Tp.
+
+ Examples
+ --------
+ >>> import wafo.spectrum.models as wsm
+ >>> import pylab as plb
+ >>> Tp,Hs = plb.meshgrid(range(4,8),range(2,6))
+ >>> gam = wsm.jonswap_peakfact(Hs,Tp)
+
+ >>> Hm0 = plb.linspace(1,20)
+ >>> Tp = Hm0
+ >>> [T,H] = plb.meshgrid(Tp,Hm0)
+ >>> gam = wsm.jonswap_peakfact(H,T)
+ >>> v = plb.arange(0,8)
+
+ >>> Hm0 = plb.arange(1,11)
+ >>> Tp = plb.linspace(2,16)
+ >>> T,H = plb.meshgrid(Tp,Hm0)
+ >>> gam = wsm.jonswap_peakfact(H,T)
+
+ h = plb.contourf(Tp,Hm0,gam,v);h=plb.colorbar()
+ h = plb.plot(Tp,gam.T)
+ h = plb.xlabel('Tp [s]')
+ h = plb.ylabel('Peakedness parameter')
+
+ plb.close('all')
+
+ See also
+ --------
+ jonswap
+ """
+ Hm0, Tp = atleast_1d(Hm0, Tp)
+
+ x = Tp / sqrt(Hm0)
+
+ gam = ones_like(x)
+
+ k1 = flatnonzero(x <= 5.14285714285714)
+ if k1.size > 0: # limiting gamma to [1 7]
+ xk = x.take(k1)
+ D = 0.036 - 0.0056 * xk # approx 5.061*Hm0**2/Tp**4*(1-0.287*log(gam))
+ # gamma
+ gam.put(k1, minimum(exp(3.484 * (1.0 - 0.1975 * D * xk ** 4.0)), 7.0))
+
+ return gam
+
+
+def jonswap_seastate(u10, fetch=150000., method='lewis', g=9.81,
+ output='dict'):
+ """
+ Return Jonswap seastate from windspeed and fetch
+
+ Parameters
+ ----------
+ U10 : real scalar
+ windspeed at 10 m above mean water surface [m/s]
+ fetch : real scalar
+ fetch [m]
+ method : 'hasselman73' seastate according to Hasselman et. al. 1973
+ 'hasselman76' seastate according to Hasselman et. al. 1976
+ 'lewis' seastate according to Lewis and Allos 1990
+ g : real scalar
+ accelaration of gravity [m/s**2]
+ output : 'dict' or 'list'
+
+ Returns
+ -------
+ seastate: dict where
+ Hm0 : significant wave height [m]
+ Tp : peak period [s]
+ gamma : jonswap peak enhancement factor.
+ sigmaA,
+ sigmaB : jonswap spectral width parameters.
+ Ag : jonswap alpha, normalization factor.
+
+ Example
+ --------
+ >>> import wafo.spectrum.models as wsm
+ >>> fetch = 10000; u10 = 10
+ >>> ss = wsm.jonswap_seastate(u10, fetch, output='dict')
+ >>> for key in sorted(ss.keys()): key, ss[key]
+ ('Ag', 0.016257903375341734)
+ ('Hm0', 0.51083679198275533)
+ ('Tp', 2.7727680999585265)
+ ('gamma', 2.4824142635861119)
+ ('sigmaA', 0.07531733139517202)
+ ('sigmaB', 0.09191208451225134)
+ >>> S = wsm.Jonswap(**ss)
+ >>> S.Hm0
+ 0.51083679198275533
+
+ # Alternatively
+ >>> ss1 = wsm.jonswap_seastate(u10, fetch, output='list')
+ >>> S1 = wsm.Jonswap(*ss1)
+ >>> S1.Hm0
+ 0.51083679198275533
+
+ See also
+ --------
+ Jonswap
+
+
+ References
+ ----------
+ Lewis, A. W. and Allos, R.N. (1990)
+ JONSWAP's parameters: sorting out the inconscistencies.
+ Ocean Engng, Vol 17, No 4, pp 409-415
+
+ Hasselmann et al. (1973)
+ Measurements of Wind-Wave Growth and Swell Decay during the Joint
+ North Sea Project (JONSWAP).
+ Ergansungsheft, Reihe A(8), Nr. 12, Deutschen Hydrografischen Zeitschrift.
+
+ Hasselmann et al. (1976)
+ A parametric wave prediction model.
+ J. phys. oceanogr. Vol 6, pp 200-228
+
+ """
+
+ # The following formulas are from Lewis and Allos 1990:
+ zeta = g * fetch / (u10 ** 2) # dimensionless fetch, Table 1
+ # zeta = min(zeta, 2.414655013429281e+004)
+ if method.startswith('h'):
+ if method[-1] == '3': # Hasselman et.al (1973)
+ A = 0.076 * zeta ** (-0.22)
+ # dimensionless peakfrequency, Table 1
+ ny = 3.5 * zeta ** (-0.33)
+ # dimensionless surface variance, Table 1
+ epsilon1 = 9.91e-8 * zeta ** 1.1
+ else: # Hasselman et.al (1976)
+ A = 0.0662 * zeta ** (-0.2)
+ ny = 2.84 * zeta ** (-0.3) # dimensionless peakfrequency, Table 1
+ # dimensionless surface variance, Eq.4
+ epsilon1 = 1.6e-7 * zeta
+
+ sa = 0.07
+ sb = 0.09
+ gam = 3.3
+ else:
+ A = 0.074 * zeta ** (-0.22) # Eq. 10
+ ny = 3.57 * zeta ** (-0.33) # dimensionless peakfrequency, Eq. 11
+ # dimensionless surface variance, Eq.12
+ epsilon1 = 3.512e-4 * A * ny ** (-4.) * zeta ** (-0.1)
+ sa = 0.05468 * ny ** (-0.32) # Eq. 13
+ sb = 0.078314 * ny ** (-0.16) # Eq. 14
+ gam = maximum(17.54 * zeta ** (-0.28384), 1) # Eq. 15
+
+ Tp = u10 / (ny * g) # Table 1
+ Hm0 = 4 * sqrt(epsilon1) * u10 ** 2. / g # Table 1
+ if output[0] == 'l':
+ return Hm0, Tp, gam, sa, sb, A
+ else:
+ return dict(Hm0=Hm0, Tp=Tp, gamma=gam, sigmaA=sa, sigmaB=sb, Ag=A)
+
+
+def _gengamspec(wn, N=5, M=4):
+ """ Return Generalized gamma spectrum in dimensionless form
+
+ Parameters
+ ----------
+ wn : arraylike
+ normalized frequencies, w/wp.
+ N : scalar
+ defining the decay of the high frequency part.
+ M : scalar
+ defining the spectral width around the peak.
+
+ Returns
+ -------
+ S : arraylike
+ spectral values, same size as wn.
+
+ The generalized gamma spectrum in non-
+ dimensional form is defined as:
+
+ S = G0.*wn.**(-N).*exp(-B*wn.**(-M)) for wn > 0
+ = 0 otherwise
+ where
+ B = N/M
+ C = (N-1)/M
+ G0 = B**C*M/gamma(C), Normalizing factor related to Bretschneider form
+
+ Note that N = 5, M = 4 corresponds to a normalized
+ Bretschneider spectrum.
+
+ Examples
+ --------
+ >>> import wafo.spectrum.models as wsm
+ >>> import numpy as np
+ >>> wn = np.linspace(0,4,5)
+ >>> wsm._gengamspec(wn, N=6, M=2)
+ array([ 0. , 1.16765216, 0.17309961, 0.02305179, 0.00474686])
+
+ See also
+ --------
+ Bretschneider
+ Jonswap,
+ Torsethaugen
+
+
+ References
+ ----------
+ Torsethaugen, K. (2004)
+ "Simplified Double Peak Spectral Model for Ocean Waves"
+ In Proc. 14th ISOPE
+ """
+ w = atleast_1d(wn)
+ S = zeros_like(w)
+
+ k = flatnonzero(w > 0.0)
+ if k.size > 0:
+ B = N / M
+ C = (N - 1.0) / M
+
+ # A = Normalizing factor related to Bretschneider form
+ # A = B**C*M/gamma(C)
+ # S[k] = A*wn[k]**(-N)*exp(-B*wn[k]**(-M))
+ logwn = log(w.take(k))
+ logA = (C * log(B) + log(M) - sp.gammaln(C))
+ S.put(k, exp(logA - N * logwn - B * exp(-M * logwn)))
+ return S
+
+class Jonswap(ModelSpectrum):
+
+ """
+ Jonswap spectral density model
+
+ Member variables
+ ----------------
+ Hm0 : significant wave height (default 7 (m))
+ Tp : peak period (default 11 (sec))
+ gamma : peakedness factor determines the concentraton
+ of the spectrum on the peak frequency.
+ Usually in the range 1 <= gamma <= 7.
+ default depending on Hm0, Tp, see jonswap_peakedness)
+ sigmaA : spectral width parameter for w<wp (default 0.07)
+ sigmaB : spectral width parameter for w<wp (default 0.09)
+ Ag : normalization factor used when gamma>1:
+ N : scalar defining decay of high frequency part. (default 5)
+ M : scalar defining spectral width around the peak. (default 4)
+ method : String defining method used to estimate Ag when gamma>1
+ 'integration': Ag = 1/gaussq(Gf*ggamspec(wn,N,M),0,wnc) (default)
+ 'parametric' : Ag = (1+f1(N,M)*log(gamma)**f2(N,M))/gamma
+ 'custom' : Ag = Ag
+ wnc : wc/wp normalized cut off frequency used when calculating Ag
+ by integration (default 6)
+ Parameters
+ ----------
+ w : array-like
+ angular frequencies [rad/s]
+
+ Description
+ -----------
+ The JONSWAP spectrum is defined as
+
+ S(w) = A * Gf * G0 * wn**(-N)*exp(-N/(M*wn**M))
+ where
+ G0 = Normalizing factor related to Bretschneider form
+ A = Ag * (Hm0/4)**2 / wp (Normalization factor)
+ Gf = j**exp(-.5*((wn-1)/s)**2) (Peak enhancement factor)
+ wn = w/wp
+ wp = angular peak frequency
+ s = sigmaA for wn <= 1
+ sigmaB for 1 < wn
+ j = gamma, (j=1, => Bretschneider spectrum)
+
+ The JONSWAP spectrum is assumed to be especially suitable for the North
+ Sea, and does not represent a fully developed sea. It is a reasonable model
+ for wind generated sea when the seastate is in the so called JONSWAP range,
+ i.e., 3.6*sqrt(Hm0) < Tp < 5*sqrt(Hm0)
+
+ The relation between the peak period and mean zero-upcrossing period
+ may be approximated by
+ Tz = Tp/(1.30301-0.01698*gamma+0.12102/gamma)
+
+ Examples
+ ---------
+ >>> import pylab as plb
+ >>> import wafo.spectrum.models as wsm
+ >>> S = wsm.Jonswap(Hm0=7, Tp=11,gamma=1)
+ >>> S2 = wsm.Bretschneider(Hm0=7, Tp=11)
+ >>> w = plb.linspace(0,5)
+ >>> all(np.abs(S(w)-S2(w))<1.e-7)
+ True
+
+ h = plb.plot(w,S(w))
+ plb.close('all')
+
+ See also
+ --------
+ Bretschneider
+ Tmaspec
+ Torsethaugen
+
+ References
+ -----------
+ Torsethaugen et al. (1984)
+ Characteristica for extreme Sea States on the Norwegian continental shelf.
+ Report No. STF60 A84123. Norwegian Hydrodyn. Lab., Trondheim
+
+ Hasselmann et al. (1973)
+ Measurements of Wind-Wave Growth and Swell Decay during the Joint
+ North Sea Project (JONSWAP).
+ Ergansungsheft, Reihe A(8), Nr. 12, Deutschen Hydrografischen Zeitschrift.
+ """
+
+ type = 'Jonswap'
+
+ def __init__(self, Hm0=7.0, Tp=11.0, gamma=None, sigmaA=0.07, sigmaB=0.09,
+ Ag=None, N=5, M=4, method='integration', wnc=6.0,
+ chk_seastate=True):
+ super(Jonswap, self).__init__(Hm0, Tp)
+ self.N = N
+ self.M = M
+ self.sigmaA = sigmaA
+ self.sigmaB = sigmaB
+ self.gamma = gamma
+ self.Ag = Ag
+ self.method = method
+ self.wnc = wnc
+
+ if self.gamma is None or not isfinite(self.gamma) or self.gamma < 1:
+ self.gamma = jonswap_peakfact(Hm0, Tp)
+
+ self._pre_calculate_ag()
+
+ if chk_seastate:
+ self.chk_seastate()
+
+ def _chk_extra_param(self):
+ Tp = self.Tp
+ Hm0 = self.Hm0
+ gam = self.gamma
+ outsideJonswapRange = Tp > 5 * sqrt(Hm0) or Tp < 3.6 * sqrt(Hm0)
+ if outsideJonswapRange:
+ txt0 = """
+ Hm0=%g,Tp=%g is outside the JONSWAP range.
+ The validity of the spectral density is questionable.
+ """ % (Hm0, Tp)
+ warnings.warn(txt0)
+
+ if gam < 1 or 7 < gam:
+ txt = """
+ The peakedness factor, gamma, is possibly too large.
+ The validity of the spectral density is questionable.
+ """
+ warnings.warn(txt)
+
+ def _localspec(self, wn):
+ Gf = self.peak_e_factor(wn)
+ return Gf * _gengamspec(wn, self.N, self.M)
+
+ def _check_parametric_ag(self, N, M, gammai):
+ parameters_ok = 3 <= N <= 50 or 2 <= M <= 9.5 and 1 <= gammai <= 20
+ if not parameters_ok:
+ raise ValueError('Not knowing the normalization because N, ' +
+ 'M or peakedness parameter is out of bounds!')
+ if self.sigmaA != 0.07 or self.sigmaB != 0.09:
+ warnings.warn('Use integration to calculate Ag when ' + 'sigmaA!=0.07 or sigmaB!=0.09')
+
+ def _parametric_ag(self):
+ """
+ Original normalization
+
+ NOTE: that Hm0**2/16 generally is not equal to intS(w)dw
+ with this definition of Ag if sa or sb are changed from the
+ default values
+ """
+ self.method = 'parametric'
+
+ N = self.N
+ M = self.M
+ gammai = self.gamma
+ f1NM = 4.1 * (N - 2 * M ** 0.28 + 5.3) ** (-1.45 * M ** 0.1 + 0.96)
+ f2NM = ((2.2 * M ** (-3.3) + 0.57) * N ** (-0.58 * M ** 0.37 + 0.53) -
+ 1.04 * M ** (-1.9) + 0.94)
+ self.Ag = (1 + f1NM * log(gammai) ** f2NM) / gammai
+ # if N == 5 && M == 4,
+ # options.Ag = (1+1.0*log(gammai).**1.16)/gammai
+ # options.Ag = (1-0.287*log(gammai))
+ # options.normalizeMethod = 'Three'
+ # elseif N == 4 && M == 4,
+ # options.Ag = (1+1.1*log(gammai).**1.19)/gammai
+
+ self._check_parametric_ag(N, M, gammai)
+
+ def _custom_ag(self):
+ self.method = 'custom'
+ if self.Ag <= 0:
+ raise ValueError('Ag must be larger than 0!')
+
+ def _integrate_ag(self):
+ # normalizing by integration
+ self.method = 'integration'
+ if self.wnc < 1.0:
+ raise ValueError('Normalized cutoff frequency, wnc, ' +
+ 'must be larger than one!')
+ area1, unused_err1 = integrate.quad(self._localspec, 0, 1)
+ area2, unused_err2 = integrate.quad(self._localspec, 1, self.wnc)
+ area = area1 + area2
+ self.Ag = 1.0 / area
+
+ def _pre_calculate_ag(self):
+ """ PRECALCULATEAG Precalculate normalization.
+ """
+ if self.gamma == 1:
+ self.Ag = 1.0
+ self.method = 'parametric'
+ elif self.Ag is not None:
+ self._custom_ag()
+ else:
+ norm_ag = dict(i=self._integrate_ag,
+ p=self._parametric_ag,
+ c=self._custom_ag)[self.method[0]]
+ norm_ag()
+
+ def peak_e_factor(self, wn):
+ """ PEAKENHANCEMENTFACTOR
+ """
+ w = maximum(atleast_1d(wn), 0.0)
+ sab = where(w > 1, self.sigmaB, self.sigmaA)
+
+ wnm12 = 0.5 * ((w - 1.0) / sab) ** 2.0
+ Gf = self.gamma ** (exp(-wnm12))
+ return Gf
+
+ def __call__(self, wi):
+ """ JONSWAP spectral density
+ """
+ w = atleast_1d(wi)
+ if (self.Hm0 > 0.0):
+
+ N = self.N
+ M = self.M
+ wp = 2 * pi / self.Tp
+ wn = w / wp
+ Ag = self.Ag
+ Hm0 = self.Hm0
+ Gf = self.peak_e_factor(wn)
+ S = ((Hm0 / 4.0) ** 2 / wp * Ag) * Gf * _gengamspec(wn, N, M)
+ else:
+ S = zeros_like(w)
+ return S
+
+if __name__ == '__main__':
+ plt.close('all')
+ Hm0 = 0.8
+ Tp = 4.0
+ S = Jonswap(Hm0=Hm0, Tp=Tp, gamma = 3.3)
+ w = np.arange(0, 50, 0.005)
+ f = w / (2 * pi)
+ Sw = S(w)
+ Sf = Sw * 2 * pi
+ plt.figure(figsize = (6, 6))
+ plt.plot(f, Sf, lw = 2) # change to Hz
+ plt.xlabel(r'$f$ $[Hz]$', fontsize = 20)
+ plt.ylabel(r'$S$', fontsize = 20)
+ plt.xlim([0, 1])
+ plt.grid()
+ plt.tight_layout()
\ No newline at end of file
diff --git a/wetb/wave/zerocrossing.py b/wetb/wave/zerocrossing.py
new file mode 100644
index 00000000..9eafd9ae
--- /dev/null
+++ b/wetb/wave/zerocrossing.py
@@ -0,0 +1,281 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 19 16:19:17 2018
+
+@author: shfe
+
+This script is used to plot the probability of exceedence of wave elevation for Linear and Nonlinear Wave
+"""
+import numpy as np
+from numpy import atleast_1d, sign, mod, zeros
+from wetb.wave import numba_misc
+from wetb.signal.spectrum import psd
+import warnings
+
+#try:
+# from wafo import c_library as clib # @UnresolvedImport
+#except ImportError:
+# warnings.warn('c_library not found. Check its compilation.')
+# clib = None
+
+def xor(a, b):
+ """
+ Return True only when inputs differ.
+ """
+ return a ^ b
+
+def _findcross(x, method='numba'):
+ '''Return indices to zero up and downcrossings of a vector
+ '''
+
+ return numba_misc.findcross(x)
+
+def findcross(x, v=0.0, kind=None, method='clib'):
+ '''
+ Return indices to level v up and/or downcrossings of a vector
+
+ Parameters
+ ----------
+ x : array_like
+ vector with sampled values.
+ v : scalar, real
+ level v.
+ kind : string
+ defines type of wave or crossing returned. Possible options are
+ 'dw' : downcrossing wave
+ 'uw' : upcrossing wave
+ 'cw' : crest wave
+ 'tw' : trough wave
+ 'd' : downcrossings only
+ 'u' : upcrossings only
+ None : All crossings will be returned
+
+ Returns
+ -------
+ ind : array-like
+ indices to the crossings in the original sequence x.
+
+ Example
+ -------
+ >>> from matplotlib import pyplot as plt
+ >>> import wafo.misc as wm
+ >>> ones = np.ones
+ >>> np.allclose(findcross([0, 1, -1, 1], 0), [0, 1, 2])
+ True
+ >>> v = 0.75
+ >>> t = np.linspace(0,7*np.pi,250)
+ >>> x = np.sin(t)
+ >>> ind = wm.findcross(x,v) # all crossings
+ >>> np.allclose(ind, [ 9, 25, 80, 97, 151, 168, 223, 239])
+ True
+ >>> ind2 = wm.findcross(x,v,'u')
+ >>> np.allclose(ind2, [ 9, 80, 151, 223])
+ True
+ >>> ind3 = wm.findcross(x,v,'d')
+ >>> np.allclose(ind3, [ 25, 97, 168, 239])
+ True
+ >>> ind4 = wm.findcross(x,v,'d', method='2')
+ >>> np.allclose(ind4, [ 25, 97, 168, 239])
+ True
+
+ t0 = plt.plot(t,x,'.',t[ind],x[ind],'r.', t, ones(t.shape)*v)
+ t0 = plt.plot(t[ind2],x[ind2],'o')
+ plt.close('all')
+
+ See also
+ --------
+ crossdef
+ wavedef
+ '''
+ xn = np.int8(sign(atleast_1d(x).ravel() - v)) # @UndefinedVariable
+ ind = _findcross(xn, method)
+ if ind.size == 0:
+ warnings.warn('No level v = %0.5g crossings found in x' % v)
+ return ind
+
+ if kind not in ('du', 'all', None):
+ if kind == 'd': # downcrossings only
+ t_0 = int(xn[ind[0] + 1] > 0)
+ ind = ind[t_0::2]
+ elif kind == 'u': # upcrossings only
+ t_0 = int(xn[ind[0] + 1] < 0)
+ ind = ind[t_0::2]
+ elif kind in ('dw', 'uw', 'tw', 'cw'):
+ # make sure that the first is a level v down-crossing
+ # if kind=='dw' or kind=='tw'
+ # or that the first is a level v up-crossing
+ # if kind=='uw' or kind=='cw'
+
+ first_is_down_crossing = int(xn[ind[0]] > xn[ind[0] + 1])
+ if xor(first_is_down_crossing, kind in ('dw', 'tw')):
+ ind = ind[1::]
+
+ n_c = ind.size # number of level v crossings
+ # make sure the number of troughs and crests are according to the
+ # wavedef, i.e., make sure length(ind) is odd if dw or uw
+ # and even if tw or cw
+ is_odd = mod(n_c, 2)
+ if xor(is_odd, kind in ('dw', 'uw')):
+ ind = ind[:-1]
+ else:
+ raise ValueError('Unknown wave/crossing definition!'
+ ' ({})'.format(kind))
+ return ind
+
+def findtc(x_in, v=None, kind=None):
+ """
+ Return indices to troughs and crests of data.
+
+ Parameters
+ ----------
+ x : vector
+ surface elevation.
+ v : real scalar
+ reference level (default v = mean of x).
+
+ kind : string
+ defines the type of wave. Possible options are
+ 'dw', 'uw', 'tw', 'cw' or None.
+ If None indices to all troughs and crests will be returned,
+ otherwise only the paired ones will be returned
+ according to the wavedefinition.
+
+ Returns
+ --------
+ tc_ind : vector of ints
+ indices to the trough and crest turningpoints of sequence x.
+ v_ind : vector of ints
+ indices to the level v crossings of the original
+ sequence x. (d,u)
+
+ Example:
+ --------
+ >>> import matplotlib.pyplot as plt
+ >>> import wafo.misc as wm
+ >>> t = np.linspace(0,30,500).reshape((-1,1))
+ >>> x = np.hstack((t, np.cos(t)))
+ >>> x1 = x[0:200,:]
+ >>> itc, iv = wm.findtc(x1[:,1],0,'dw')
+ >>> tc = x1[itc,:]
+ >>> np.allclose(itc, [ 52, 105])
+ True
+ >>> itc, iv = wm.findtc(x1[:,1],0,'uw')
+ >>> np.allclose(itc, [ 105, 157])
+ True
+
+ a = plt.plot(x1[:,0],x1[:,1],tc[:,0],tc[:,1],'ro')
+ plt.close('all')
+
+ See also
+ --------
+ findtp
+ findcross,
+ wavedef
+ """
+
+ x = atleast_1d(x_in)
+ if v is None:
+ v = x.mean()
+
+ v_ind = findcross(x, v, kind)
+ n_c = v_ind.size
+ if n_c <= 2:
+ warnings.warn('There are no waves!')
+ return zeros(0, dtype=np.int), zeros(0, dtype=np.int)
+
+ # determine the number of trough2crest (or crest2trough) cycles
+ is_even = mod(n_c + 1, 2)
+ n_tc = int((n_c - 1 - is_even) / 2)
+
+ # allocate variables before the loop increases the speed
+ ind = zeros(n_c - 1, dtype=np.int)
+
+ first_is_down_crossing = (x[v_ind[0]] > x[v_ind[0] + 1])
+ if first_is_down_crossing:
+ f1, f2 = np.argmin, np.argmax
+ else:
+ f1, f2 = np.argmax, np.argmin
+
+ for i in range(n_tc):
+ # trough or crest
+ j = 2 * i
+ ind[j] = f1(x[v_ind[j] + 1:v_ind[j + 1] + 1])
+ # crest or trough
+ ind[j + 1] = f2(x[v_ind[j + 1] + 1:v_ind[j + 2] + 1])
+
+ if (2 * n_tc + 1 < n_c) and (kind in (None, 'tw', 'cw')):
+ # trough or crest
+ ind[n_c - 2] = f1(x[v_ind[n_c - 2] + 1:v_ind[n_c - 1] + 1])
+
+ return v_ind[:n_c - 1] + ind + 1, v_ind
+
+def trough_crest(time, eta, v=None, wavetype=None):
+ """
+ Return trough and crest turning points
+
+ Parameters
+ -----------
+ v : scalar
+ reference level (default v = mean of x).
+
+ wavetype : string
+ defines the type of wave. Possible options are
+ 'dw', 'uw', 'tw', 'cw' or None.
+ If None indices to all troughs and crests will be returned,
+ otherwise only the paired ones will be returned
+ according to the wavedefinition.
+
+ Returns
+ --------
+ tc : TurningPoints object
+ with trough and crest turningpoints
+ """
+ ind = findtc(eta, v, wavetype)[0]
+ eta_tc = eta[ind]
+
+ try:
+ t = time
+ except Exception:
+ t = ind
+ t_tc = t[ind]
+
+ return t_tc, eta_tc
+
+def wave_parameters(time, eta):
+
+ """
+ Returns several wave parameters from data.
+ Parameters
+ ----------
+ rate : scalar integer
+ interpolation rate. Interpolates with spline if greater than one.
+ Returns
+ -------
+ parameters : dict
+ wave parameters such as
+ Ac, At : Crest and trough amplitude, respectively
+ Tc, Tt : crest time and trough time
+ Hu, Hd : zero-up- and down-crossing wave height, respectively.
+ Tu, Td : zero-up- and down-crossing wave period, respectively.
+ """
+ tc_ind, z_ind = findtc(eta, v=None, kind='tw')
+
+ tc_a = eta[tc_ind]
+ tc_t = time[tc_ind]
+
+ Ac = tc_a[1::2]
+ At = -tc_a[0::2]
+
+ Tc = tc_t[1::2]
+ Tt = tc_t[0::2]
+
+ Hu = Ac + At[1:]
+ Hd = Ac + At[:-1]
+
+ tz = time[z_ind]
+ tu = tz[1::2]
+# td = tz[0::2]
+ Tu = np.diff(tu)
+ T_crest = Tc - tu[:-1]
+ return Ac, At, Tc, Tt, Hu, Hd, Tu, tu, T_crest
+
\ No newline at end of file
--
GitLab