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