From 8481dd456dea0121d4b5632110092cee47515a6d Mon Sep 17 00:00:00 2001
From: "Edgar P. Burkhart" <git@edgarpierre.fr>
Date: Thu, 31 Mar 2022 14:40:19 +0200
Subject: [PATCH] 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<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
diff --git a/swash/processing/pa/reflex3S.py b/swash/processing/pa/reflex3S.py
new file mode 100644
index 0000000..9417d29
--- /dev/null
+++ b/swash/processing/pa/reflex3S.py
@@ -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