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