163 lines
4.9 KiB
Python
163 lines
4.9 KiB
Python
import math
|
|
import os
|
|
import shutil
|
|
import subprocess
|
|
import numpy as np
|
|
import csv
|
|
import matplotlib.pyplot as plt
|
|
import sys
|
|
from scipy import signal
|
|
import re
|
|
from scipy.interpolate import griddata
|
|
import scipy.io as sio
|
|
from mpl_toolkits import mplot3d
|
|
from scipy import signal
|
|
from scipy.fftpack import fft, ifft
|
|
from scipy.signal import detrend
|
|
from numpy import angle
|
|
|
|
g = 9.81
|
|
|
|
# h = profondeur locale
|
|
# fs = fréquence echantillonage
|
|
# zmesP = profondeur de mesure de la pression
|
|
# zmesU = profondeur de mesure de la vitesse
|
|
def PUV_dam(u, p, h, fs, zmesP, zmesU):
|
|
# u = detrend(u)
|
|
# p = detrend(p)
|
|
ampliseuil = 0.001
|
|
N = len(u)
|
|
delta = fs
|
|
# transformée de fourrier
|
|
Yu = fft(u)
|
|
Yp = fft(p)
|
|
nyquist = 1 / 2
|
|
# Fu_xb = ((1:N/2)-1)/(N/2)/deltat*nyquist
|
|
xf = np.linspace(0.0, 1.0 / (2.0 / fs), N // 2)
|
|
Ampu = np.abs(Yu[1 : N // 2]) / (N / 2.0)
|
|
Ampp = np.abs(Yp[1 : N // 2]) / (N / 2.0)
|
|
# pas de frequence deltaf
|
|
deltaf = 1 / (N / 2) / delta * nyquist
|
|
# phase
|
|
ThetaU = angle(Yu[1 : N // 2])
|
|
ThetaP = angle(Yp[1 : N // 2])
|
|
# calcul du coefficient de reflexion
|
|
nbf = len(xf)
|
|
if max(p) > 0:
|
|
# si pas de données de pression, jsute Sheremet
|
|
# attention : on commence par le deuxieme point, le premier correspondant a frequence nulle
|
|
iicutoff_xb = max(min(np.where(xf > 0.5))) # length(Fu_xb)) ?
|
|
# calcul de l'amplitude en deca de la frequence de coupure
|
|
k_xb = []
|
|
omega = []
|
|
ccc, ccs, cccp = [], [], []
|
|
aincident13, areflechi13 = [], []
|
|
aprog = []
|
|
r13 = []
|
|
ii = 0
|
|
while ii < iicutoff_xb:
|
|
k_xb.append(newtonpplus(xf[ii], h, 0))
|
|
a1 = Ampu[ii]
|
|
a3 = Ampp[ii]
|
|
phi1 = ThetaU[ii]
|
|
phi3 = ThetaP[ii]
|
|
omega.append(2 * np.pi * xf[ii])
|
|
cc = omega[ii] * np.cosh(xf[ii] * (zmesU + h)) / np.sinh(xf[ii] * h)
|
|
ccp = omega[ii] * np.cosh(xf[ii] * (zmesP + h)) / np.sinh(xf[ii] * h)
|
|
cs = omega[ii] * np.sinh(xf[ii] * (zmesU + h)) / np.sinh(xf[ii] * h)
|
|
# Procedure de calcul des ai et ar sans courant
|
|
ccc.append(cc)
|
|
ccs.append(cs)
|
|
cccp.append(ccp)
|
|
# calcul des amplitudes des ondes incidentes a partir de uhoriz et p
|
|
aincident13.append(
|
|
0.5
|
|
* np.sqrt(
|
|
a1 * a1 / (cc * cc)
|
|
+ a3
|
|
* a3
|
|
* g
|
|
* g
|
|
* xf[ii]
|
|
* xf[ii]
|
|
/ (omega[ii] * omega[ii] * ccp * ccp)
|
|
+ 2
|
|
* a1
|
|
* a3
|
|
* g
|
|
* k_xb[ii]
|
|
* np.cos(phi3 - phi1)
|
|
/ (cc * ccp * omega[ii])
|
|
)
|
|
)
|
|
areflechi13.append(
|
|
0.5
|
|
* np.sqrt(
|
|
a1 * a1 / (cc * cc)
|
|
+ a3
|
|
* a3
|
|
* g
|
|
* g
|
|
* k_xb[ii]
|
|
* k_xb[ii]
|
|
/ (omega[ii] * omega[ii] * ccp * ccp)
|
|
- 2
|
|
* a1
|
|
* a3
|
|
* g
|
|
* xf[ii]
|
|
* np.cos(phi3 - phi1)
|
|
/ (cc * ccp * omega[ii])
|
|
)
|
|
)
|
|
cv = g * xf[ii] / (omega[ii] * ccp)
|
|
cp = 1 / cc
|
|
aprog.append(a3 / (g * xf[ii] / (omega[ii] * ccp)))
|
|
# hypothese progressive Drevard et al |u|/Cv
|
|
# aprog(ii)= a1/cp;
|
|
# |p|/Cp
|
|
if aincident13[ii] < 0.18 * ampliseuil:
|
|
r13.append(0)
|
|
else:
|
|
r13.append(areflechi13[ii] / aincident13[ii])
|
|
# calcul des energies incidente et reflechie sans ponderation avec les voisins
|
|
Eprog = 0.5 * aprog**2 / deltaf
|
|
Eixb = 0.5 * aincident13**2 / deltaf
|
|
Erxb = 0.5 * areflechi13**2 / deltaf
|
|
F = Fu_xb[1:iicutoff_xb]
|
|
return Eixb, Erxb, Eprog, F
|
|
"test"
|
|
|
|
|
|
# Calcul du vecteur d'onde en prÈsence d'un courant
|
|
# constant u de sens opposÈ ‡ celui de propagation de l'onde
|
|
# a partir de la frÈquence f, la profondeur d'eau h et
|
|
# la vitesse u du courant
|
|
# kh : vecteur d'onde * profondeur d'eau
|
|
def newtonpplus(f, h, u):
|
|
# calcul de k:
|
|
g = 9.81
|
|
kh = 0.5
|
|
x = 0.001
|
|
u = -u
|
|
i = 0
|
|
while abs((kh - x) / x) > 0.00001:
|
|
i += 1
|
|
if i > 10**5:
|
|
sys.exit(1)
|
|
x = kh
|
|
fx = x * math.tanh(x) - (h / g) * (2 * math.pi * f - (u / h) * x) * (
|
|
2 * math.pi * f - (u / h) * x
|
|
)
|
|
fprimx = (
|
|
math.tanh(x)
|
|
+ x * (1 - (math.tanh(x)) ** 2)
|
|
+ (2 * u / g) * (2 * math.pi * f - (u / h) * x)
|
|
)
|
|
kh = x - (fx / fprimx)
|
|
k = kh / h
|
|
return k
|
|
|
|
|
|
# fin du calcul de k
|