From 8481dd456dea0121d4b5632110092cee47515a6d Mon Sep 17 00:00:00 2001 From: "Edgar P. Burkhart" Date: Thu, 31 Mar 2022 14:40:19 +0200 Subject: [PATCH 1/2] Functions from PAPoncet for reflection calculation (array and PUV) --- swash/processing/pa/PUV.py | 302 ++++++++++++++++++++++++++++++++ swash/processing/pa/reflex3S.py | 220 +++++++++++++++++++++++ 2 files changed, 522 insertions(+) create mode 100644 swash/processing/pa/PUV.py create mode 100644 swash/processing/pa/reflex3S.py diff --git a/swash/processing/pa/PUV.py b/swash/processing/pa/PUV.py new file mode 100644 index 0000000..38afe6b --- /dev/null +++ b/swash/processing/pa/PUV.py @@ -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= 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 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 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 Date: Thu, 31 Mar 2022 15:23:46 +0200 Subject: [PATCH 2/2] Fixed some issues in PA code --- swash/processing/pa/PUV.py | 208 ++------------------------------ swash/processing/pa/reflex3S.py | 2 +- 2 files changed, 13 insertions(+), 197 deletions(-) diff --git a/swash/processing/pa/PUV.py b/swash/processing/pa/PUV.py index 38afe6b..1c38113 100644 --- a/swash/processing/pa/PUV.py +++ b/swash/processing/pa/PUV.py @@ -10,10 +10,11 @@ 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 +from scipy.signal import detrend +from numpy import angle # fonction moyenne glissante def lissage(Lx,Ly,p): @@ -27,143 +28,31 @@ def lissage(Lx,Ly,p): 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= 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) + #u = detrend(u) + #p = detrend(p) ampliseuil = 0.001 N = len(u) - delta = 1.0/T + 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*T), N//2) - Ampu = np.abs(Yu[1:N/2])/(N/2.0) - Ampp = np.abs(Yp[1:N/2])/(N/2.0) + 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)/deltat*nyquist; + deltaf=1/(N/2)/delta*nyquist; #phase - ThetaU=angle(Yu[1:N/2]); - ThetaP=angle(Yp[1:N/2]); + ThetaU=angle(Yu[1:N//2]); + ThetaP=angle(Yp[1:N//2]); #calcul du coefficient de reflexion - nbf=length(xf); + 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 @@ -205,8 +94,6 @@ def PUV_dam(u,p,h,fs,zmesP,zmesU): 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 @@ -229,74 +116,3 @@ def newtonpplus(f,h,u) : # 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