1
Fork 0

Functions from PAPoncet for reflection calculation (array and PUV)

This commit is contained in:
Edgar P. Burkhart 2022-03-31 14:40:19 +02:00
parent e95d679def
commit 8481dd456d
Signed by: edpibu
GPG key ID: 9833D3C5A25BD227
2 changed files with 522 additions and 0 deletions

302
swash/processing/pa/PUV.py Normal file
View file

@ -0,0 +1,302 @@
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
import imageio
from mpl_toolkits import mplot3d
from scipy import signal
from scipy.fftpack import fft, ifft
# 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
# import orientation data
with open('ADV/ARTHA_aut18.sen','r') as n:
content = n.readlines()
heading = []
pitch = []
roll = []
i = 0
for i in range(len(content)-1):
row = content[i].split()
heading.append(float(row[10]))
pitch.append(float(row[11]))
roll.append(float(row[12]))
i+=1
# import Pressure and velocity measurement
with open('ADV/ARTHA_aut18.dat','r') as n:
content = n.readlines()
E = []
N = []
Up = []
P = []
burstnb = []
position = []
i = 0
for i in range(len(content)-1):
row = content[i].split()
burstnb.append(int(row[0]))
position.append(int(row[1]))
E.append(float(row[2]))
N.append(float(row[3]))
Up.append(float(row[4]))
P.append(float(row[14]))
i+=1
# import seal level pressure and reshape vector to fit the ADV data
with open('press_socoa_0110_2311_18.data','r') as n:
content = n.readlines()[1:]
burstind = []
date = []
slp = []
ind = 1
i = 0
for i in range(len(content)-1):
row = content[i].split(";")
date.append(int(row[1]))
slp.append(float(row[3]))
if date[i] >= 2018101004 and date[i]%2 == 0:
burstind.append(int(ind))
ind+=1
else :
burstind.append("nan")
i+=1
slpind = []
i=0
while i<len(slp):
if burstind[i] != 'nan' :
slpind.append(slp[i])
i+=1
slpind = slpind[:304]
# height from ground to ADV
delta = 0.4
# gravity
g = 9.81
#define time vector
t = np.linspace(0, 1200, 2400)
burstP = np.zeros((304,2400))
burstE = np.zeros((304,2400))
burstN = np.zeros((304,2400))
i = 0
j = 0
n = 1
while i < len(P) :
burstP[n,j] = P[i] - slpind[n]/100*0.4/10.13
burstE[n,j] = E[i]
burstN[n,j] = N[i]
if i >= n*2400 - 1:
n+=1
j=-1
i+=1
j+=1
burstP = np.delete(burstP,0,0)
burstE = np.delete(burstE,0,0)
burstN = np.delete(burstN,0,0)
# number of point in the burst
N = 2400
# time step
T = 0.5
#frequency vector
xf = np.linspace(0.0, 1.0/(2.0*T), N//2)
# hs calculation interval
x1 = np.max(np.where(xf < 0.05))
x2 = np.min(np.where(xf > 0.17))
# u = vitesse hori
# p = pression
# 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 = 1.0/T
#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*T), 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)/deltat*nyquist;
#phase
ThetaU=angle(Yu[1:N/2]);
ThetaP=angle(Yp[1:N/2]);
#calcul du coefficient de reflexion
nbf=length(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'
fs =2
# 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
u = signal.detrend(U[1])
p = signal.detrend(burstP[1])
ampliseuil = 0.001
N = len(u)
deltat = 1.0/fs
#transformée de fourrier
Yu = fft(u)
Yp = fft(p)
nyquist = 1/2.0
# h = profondeur locale
h = np.mean(burstP[1])
#Fu_xb = ((1:N/2)-1)/(N/2)/deltat*nyquist
xf = np.linspace(0.0, 1.0/(2.0*T), 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.0/(N/2.0)/deltat*nyquist
#phase
ThetaU=np.angle(Yu[1:N/2])
ThetaP=np.angle(Yp[1:N/2])
#calcul du coefficient de reflexion
#fs = fréquence echantillonage
fs = 2
# zmesP = profondeur de mesure de la pression
zmesP = h - 0.4
# zmesU = profondeur de mesure de la vitesse
zmesU = zmesP
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 =[]
r13 =[]
aprog = []
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*math.pi*xf[ii])
cc=omega[ii]*math.cosh(xf[ii]*(zmesU+h))/math.sinh(xf[ii]*h);
ccp=omega[ii]*math.cosh(xf[ii]*(zmesP+h))/math.sinh(xf[ii]*h);
cs=omega[ii]*math.sinh(xf[ii]*(zmesU+h))/math.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*math.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]*math.cos(phi3-phi1)/(cc*ccp*omega[ii])))
areflechi13.append(0.5*math.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]*math.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])
ii+=1

View 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 imageio
import cmath
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