Black + improvements on pa reflex/puv
This commit is contained in:
parent
2f16aef816
commit
58ed9a4716
GPG key ID:
9833D3C5A25BD227
2 changed files with 313 additions and 236 deletions
|
|
@ -1,4 +1,4 @@
|
|||
import math
|
||||
import math
|
||||
import os
|
||||
import shutil
|
||||
import subprocess
|
||||
|
|
@ -16,83 +16,99 @@ from scipy.fftpack import fft, ifft
|
|||
from scipy.signal import detrend
|
||||
from numpy import angle
|
||||
|
||||
# fonction moyenne glissante
|
||||
def lissage(Lx,Ly,p):
|
||||
'''Fonction qui débruite une courbe par une moyenne glissante
|
||||
sur 2P+1 points'''
|
||||
Lxout=[]
|
||||
Lyout=[]
|
||||
Lxout = Lx[p: -p]
|
||||
for index in range(p, len(Ly)-p):
|
||||
average = np.mean(Ly[index-p : index+p+1])
|
||||
Lyout.append(average)
|
||||
return Lxout,Lyout
|
||||
g = 9.81
|
||||
|
||||
# h = profondeur locale
|
||||
#fs = fréquence echantillonage
|
||||
# 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)
|
||||
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
|
||||
# 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
|
||||
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 = []
|
||||
ii = 1
|
||||
while ii<iicutoff_xb :
|
||||
k_xb[ii] = newtonpplus(xf[ii],h,0);
|
||||
a1=Ampu[ii];
|
||||
a3=Ampp[ii];
|
||||
phi1=ThetaU[ii];
|
||||
phi3=ThetaP[ii];
|
||||
omega[ii]=2*pi*xf[ii];
|
||||
cc=omega[ii]*cosh(xf[ii]*(zmesU+h))/sinh(xf[ii]*h);
|
||||
ccp=omega[ii]*cosh(xf[ii]*(zmesP+h))/sinh(xf[ii]*h);
|
||||
cs=omega[ii]*sinh(xf[ii]*(zmesU+h))/sinh(xf[ii]*h);
|
||||
#Procedure de calcul des ai et ar sans courant
|
||||
ccc[ii]=cc;
|
||||
ccs[ii]=cs;
|
||||
cccp[ii]=ccp;
|
||||
#calcul des amplitudes des ondes incidentes a partir de uhoriz et p
|
||||
aincident13[ii]=0.5*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]*cos(phi3-phi1)/(cc*ccp*omega[ii]))
|
||||
areflechi13[ii]=0.5*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]*cos(phi3-phi1)/(cc*ccp*omega[ii]))
|
||||
cv=g*xf[ii]/(omega[ii]*ccp)
|
||||
cp=1/cc
|
||||
aprog[ii]= 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[ii]=0
|
||||
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[ii]=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'
|
||||
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
|
||||
|
|
@ -100,19 +116,29 @@ def PUV_dam(u,p,h,fs,zmesP,zmesU):
|
|||
# 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) :
|
||||
def newtonpplus(f, h, u):
|
||||
# calcul de k:
|
||||
g = 9.81
|
||||
kh = 0.5
|
||||
x = 0.001
|
||||
u=-u
|
||||
while (abs((kh - x)/x) > 0.00000001) :
|
||||
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
|
||||
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
|
||||
|
||||
# fin du calcul de k
|
||||
|
|
|
|||
Reference in a new issue