Merge branch 'reflection' into swash
This commit is contained in:
commit
6ae1ea9899
GPG key ID:
9833D3C5A25BD227
2 changed files with 338 additions and 0 deletions
118
swash/processing/pa/PUV.py
Normal file
118
swash/processing/pa/PUV.py
Normal file
|
|
@ -0,0 +1,118 @@
|
||||||
|
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
|
||||||
|
|
||||||
|
# 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
|
||||||
|
|
||||||
|
# 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 = []
|
||||||
|
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
|
||||||
|
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'
|
||||||
|
|
||||||
|
|
||||||
|
# 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
|
||||||
|
while (abs((kh - x)/x) > 0.00000001) :
|
||||||
|
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
|
||||||
|
|
||||||
220
swash/processing/pa/reflex3S.py
Normal file
220
swash/processing/pa/reflex3S.py
Normal file
|
|
@ -0,0 +1,220 @@
|
||||||
|
import csv
|
||||||
|
import matplotlib.pyplot as plt
|
||||||
|
import os
|
||||||
|
import math
|
||||||
|
import sys
|
||||||
|
import numpy as np
|
||||||
|
from matplotlib import animation
|
||||||
|
import cmath
|
||||||
|
from scipy.fft import fft
|
||||||
|
|
||||||
|
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) :
|
||||||
|
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
|
||||||
|
|
||||||
|
def newtonpmoins(f,h,u0) :
|
||||||
|
# calcul de k:
|
||||||
|
g = 9.81;
|
||||||
|
kh = 0.5;
|
||||||
|
x = 0.01;
|
||||||
|
x = 6*h;
|
||||||
|
|
||||||
|
while (np.abs((kh - x)/x) > 0.00000001):
|
||||||
|
x = kh;
|
||||||
|
fx = x*math.tanh(x) - (h/g)*(2*math.pi*f-(u0/h)*x)*(2*math.pi*f-(u0/h)*x);
|
||||||
|
fprimx = math.tanh(x) + x*(1- (math.tanh(x))**2)+(2*u0/g)*(2*math.pi*f-(u0/h)*x);
|
||||||
|
kh = x - (fx/fprimx);
|
||||||
|
k = kh/h
|
||||||
|
return k
|
||||||
|
|
||||||
|
#Calcul du vecteur d'onde a partir de la frÈquence
|
||||||
|
#kh : vecteur d'onde * profondeur d'eau
|
||||||
|
def newtonpropa(hi,f):
|
||||||
|
# calcul de k:
|
||||||
|
g=9.81;
|
||||||
|
si = (2*math.pi*f)**2/g;
|
||||||
|
kh = 0.5;
|
||||||
|
x = 0.001;
|
||||||
|
while (np.abs((kh - x)/x) > 0.00000001) :
|
||||||
|
x = kh;
|
||||||
|
fx = x*math.tanh(x) - si*hi;
|
||||||
|
fprimx = math.tanh(x) + x*(1- (math.tanh(x))**2);
|
||||||
|
kh = x - (fx/fprimx);
|
||||||
|
kpropa = kh/hi;
|
||||||
|
return kpropa
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
def reflex3S(x1,x2,x3,xs1,xs2,xs3,h,mean_freq,fmin,fmax) :
|
||||||
|
# Analyse avec transformee de fourier d un signal en sinus
|
||||||
|
# calcul du coefficient de reflexion en presence d un courant u
|
||||||
|
# tinit : temps initial
|
||||||
|
# tfinal : temps final
|
||||||
|
# deltat : pas de temps
|
||||||
|
# fech= 1/deltat : frequence d echantillonnage
|
||||||
|
# T : periode de la houle
|
||||||
|
# nbrepoints : nombre de points
|
||||||
|
# fmin : frequence minimale de recherche des maxima du signal de fourier
|
||||||
|
# fmax : frequence maximale de recherche des maxima du signal de fourier
|
||||||
|
# ampliseuil : amplitude minimale des maxima recherches
|
||||||
|
# valeur du courant (u > 0 correspond a un courant dans le sens
|
||||||
|
# des x croissants)
|
||||||
|
#u=0;
|
||||||
|
|
||||||
|
# h profondeur d'eau
|
||||||
|
|
||||||
|
#hold on;
|
||||||
|
#fech=16;
|
||||||
|
#fmin=0.1;fmax=4;ampliseuil=0.005;
|
||||||
|
#nbrepoints=fech*T*1000;
|
||||||
|
#deltat=1./fech;
|
||||||
|
#tinit=0;tfinal=tinit+deltat*(nbrepoints-1);
|
||||||
|
|
||||||
|
#aitheo=1;artheo=0.4;
|
||||||
|
#h=3;
|
||||||
|
#T=1.33;
|
||||||
|
#ftheo=1/T;
|
||||||
|
#fech=16;
|
||||||
|
#fmin=0.1;fmax=4;ampliseuil=0.005;
|
||||||
|
#nbrepoints=fech*T*1000;
|
||||||
|
#deltat=1./fech;
|
||||||
|
#tinit=0;tfinal=tinit+deltat*(nbrepoints-1);
|
||||||
|
#t = [tinit:deltat:tfinal];
|
||||||
|
#ktheo = newtonpropa(ftheo,h);
|
||||||
|
|
||||||
|
#Positions respectives sondes amonts entr'elles et sondes aval
|
||||||
|
# entr'elles
|
||||||
|
|
||||||
|
# xs1=0;xs2=0.80;xs3=1.30;
|
||||||
|
|
||||||
|
#ENTREES DONNEES DES 3 SONDES AMONT et des 2 SONDES AVAL
|
||||||
|
ampliseuil=0.005;
|
||||||
|
#'check'
|
||||||
|
#pause
|
||||||
|
#PAS DE TEMPS
|
||||||
|
deltat1=1/mean_freq;
|
||||||
|
deltat2=1/mean_freq;
|
||||||
|
deltat3=1/mean_freq;
|
||||||
|
#transformees de Fourier
|
||||||
|
Y1 = fft(x1,len(x1));
|
||||||
|
N1 = len(Y1);
|
||||||
|
Y2 = fft(x2,len(x2));
|
||||||
|
N2 = len(Y2);
|
||||||
|
Y3 = fft(x3,len(x3));
|
||||||
|
N3 = len(Y3);
|
||||||
|
#amplitudes normalisees, soit coef de fourier
|
||||||
|
amplitude1=np.abs(Y1[1:N1//2])/(N1//2);
|
||||||
|
nyquist = 1/2;
|
||||||
|
freq1 = (np.arange(1, (N1//2)+1, 1)-1)/(N1//2)/deltat1*nyquist;
|
||||||
|
amplitude2=np.abs(Y2[1:N2//2])/(N2//2);
|
||||||
|
nyquist = 1/2;
|
||||||
|
freq2 = (np.arange(1, (N2//2)+1, 1)-1)/(N2//2)/deltat2*nyquist;
|
||||||
|
amplitude3=np.abs(Y3[1:N3//2])/(N3//2);
|
||||||
|
nyquist = 1/2;
|
||||||
|
freq3 = (np.arange(1, (N3//2)+1, 1)-1)/(N3//2)/deltat3*nyquist;
|
||||||
|
#recherche de la phase
|
||||||
|
theta1=np.angle(Y1[1:N1//2]);
|
||||||
|
theta2=np.angle(Y2[1:N2//2]);
|
||||||
|
theta3=np.angle(Y3[1:N3//2]);
|
||||||
|
#pas de frequence deltaf
|
||||||
|
deltaf=1/(N1//2)/deltat1*nyquist;
|
||||||
|
nbrefreq=len(freq1);
|
||||||
|
#Caracteristiques fondamentaux,sondes canaux 1 et 3
|
||||||
|
#distances entre les sondes
|
||||||
|
x12=xs2-xs1;
|
||||||
|
x13=xs3-xs1;
|
||||||
|
x23=xs3-xs2;
|
||||||
|
#Debut calcul des coefficients de reflexion
|
||||||
|
indmin=np.min(np.where(freq1>0.02));
|
||||||
|
indfmin=np.min(np.where(freq1>fmin));
|
||||||
|
indfmax=np.max(np.where(freq1<fmax));
|
||||||
|
|
||||||
|
|
||||||
|
T = []
|
||||||
|
fre = []
|
||||||
|
aincident12 = []
|
||||||
|
aincident23 = []
|
||||||
|
aincident13 = []
|
||||||
|
areflechi12 = []
|
||||||
|
areflechi23 = []
|
||||||
|
areflechi13 = []
|
||||||
|
r13 = []
|
||||||
|
ai = []
|
||||||
|
ar = []
|
||||||
|
Eincident123 = []
|
||||||
|
Ereflechi123 = []
|
||||||
|
count = 0
|
||||||
|
for jj in np.arange(indfmin, indfmax, 1):
|
||||||
|
f=freq1[jj]
|
||||||
|
#periode
|
||||||
|
T.append(1/f)
|
||||||
|
fre.append(f)
|
||||||
|
#calcul des vecteurs d'onde
|
||||||
|
kplus = newtonpplus(f,h,0)
|
||||||
|
kmoins = newtonpmoins(f,h,0)
|
||||||
|
k = newtonpropa(h,f)
|
||||||
|
deltaku=k-(kmoins+kplus)/2
|
||||||
|
#amplitude des signaux pour la frequence f:
|
||||||
|
a1=amplitude1[jj]
|
||||||
|
a2=amplitude2[jj]
|
||||||
|
a3=amplitude3[jj]
|
||||||
|
#dephasages entre les signaux experimentaux des 3 sondes amont
|
||||||
|
phi1=theta1[jj]
|
||||||
|
phi2=theta2[jj]
|
||||||
|
phi3=theta3[jj]
|
||||||
|
phi12=phi2-phi1
|
||||||
|
phi13=phi3-phi1
|
||||||
|
phi23=phi3-phi2
|
||||||
|
#evolution theorique entre les sondes de la phase pour une onde progressive
|
||||||
|
delta12p= -kplus*x12
|
||||||
|
delta13p= -kplus*x13
|
||||||
|
delta23p= -kplus*x23
|
||||||
|
delta12m= -kmoins*x12
|
||||||
|
delta13m= -kmoins*x13
|
||||||
|
delta23m= -kmoins*x23
|
||||||
|
#calcul du coefficient de reflexion a partir des sondes 1 et 2
|
||||||
|
aincident12.append(math.sqrt(a1*a1+a2*a2-2*a1*a2*math.cos(phi12+delta12p))/(2*np.abs(math.sin((delta12p+delta12m)/2))))
|
||||||
|
areflechi12.append(math.sqrt(a1*a1+a2*a2-2*a1*a2*math.cos(phi12-delta12m))/(2*np.abs(math.sin((delta12p+delta12m)/2))))
|
||||||
|
#r12(jj)=areflechi12(jj)/aincident12(jj);
|
||||||
|
#calcul du coefficient de reflexion a partir des sondes 2 et 3
|
||||||
|
aincident23.append(math.sqrt(a2*a2+a3*a3-2*a2*a3*math.cos(phi23+delta23p))/(2*np.abs(math.sin((delta23p+delta23m)/2))))
|
||||||
|
areflechi23.append(math.sqrt(a2*a2+a3*a3-2*a2*a3*math.cos(phi23-delta23m))/(2*np.abs(math.sin((delta23p+delta23m)/2))))
|
||||||
|
#r23(jj)=areflechi23(jj)/aincident23(jj);
|
||||||
|
#calcul du coefficient de reflexion a partir des sondes 1 et 3
|
||||||
|
aincident13.append(math.sqrt(a1*a1+a3*a3-2*a1*a3*math.cos(phi13+delta13p))/(2*np.abs(math.sin((delta13p+delta13m)/2))))
|
||||||
|
areflechi13.append(math.sqrt(a1*a1+a3*a3-2*a1*a3*math.cos(phi13-delta13m))/(2*np.abs(math.sin((delta13p+delta13m)/2))))
|
||||||
|
#r13.append(areflechi13[jj]/aincident13[jj])
|
||||||
|
#calcul du coefficient de reflexion par methode des 3 sondesavec moindres carres
|
||||||
|
delta1m=0
|
||||||
|
delta2m=delta12m
|
||||||
|
delta3m=delta13m
|
||||||
|
delta1p=0
|
||||||
|
delta2p=delta12p
|
||||||
|
delta3p=delta13p
|
||||||
|
s1=cmath.exp(-1j*2*delta1m)+cmath.exp(-1j*2*delta2m)+cmath.exp(-1j*2*delta3m)
|
||||||
|
s2=cmath.exp(+1j*2*delta1p)+cmath.exp(+1j*2*delta2p)+cmath.exp(+1j*2*delta3p)
|
||||||
|
s12=cmath.exp(1j*(delta1p-delta1m))+cmath.exp(1j*(delta2p-delta2m))+cmath.exp(1j*(delta3p-delta3m))
|
||||||
|
s3=a1*cmath.exp(-1j*(phi1+delta1m))+a2*cmath.exp(-1j*(phi2+delta2m))+a3*cmath.exp(-1j*(phi3+delta3m))
|
||||||
|
s4=a1*cmath.exp(-1j*(phi1-delta1p))+a2*cmath.exp(-1j*(phi2-delta2p))+a3*cmath.exp(-1j*(phi3-delta3p))
|
||||||
|
s5=s1*s2-s12*s12
|
||||||
|
ai.append(abs((s2*s3-s12*s4)/s5))
|
||||||
|
ar.append(abs((s1*s4-s12*s3)/s5))
|
||||||
|
#refl[jj]=ar[jj]/ai[jj];
|
||||||
|
#Calcul de l'energie, on divise par le pas de frequence deltaf
|
||||||
|
#calcul de l'energie incidente sans ponderation avec les voisins
|
||||||
|
Eincident123.append(0.5*ai[count]*ai[count]/deltaf)
|
||||||
|
Ereflechi123.append(0.5*ar[count]*ar[count]/deltaf)
|
||||||
|
count+=1
|
||||||
|
|
||||||
|
|
||||||
|
return ai,ar,Eincident123, Ereflechi123, indfmin, indfmax, fre
|
||||||
Loading…
Reference in a new issue