Fixed some issues in PA code

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<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)
+    #u = detrend(u)
-    p = detrend(p)
+    #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)
+    xf = np.linspace(0.0, 1.0/(2.0/fs), N//2)
-    Ampu = np.abs(Yu[1:N/2])/(N/2.0)
+    Ampu = np.abs(Yu[1:N//2])/(N/2.0)
-    Ampp = np.abs(Yp[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]);
+    ThetaU=angle(Yu[1:N//2]);
-    ThetaP=angle(Yp[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<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

swash/processing/pa/reflex3S.py

@@ -5,8 +5,8 @@ import math
 import sys
 import numpy as np
 from matplotlib import animation
-import imageio
 import cmath
+from scipy.fft import fft
 
 def newtonpplus(f,h,u) :
     # calcul de k: