diff --git a/swash/processing/pa/PUV.py b/swash/processing/pa/PUV.py
index 1c38113..858bb0f 100644
--- a/swash/processing/pa/PUV.py
+++ b/swash/processing/pa/PUV.py
@@ -1,4 +1,4 @@
-import math 
+import math
 import os
 import shutil
 import subprocess
@@ -16,83 +16,99 @@ from scipy.fftpack import fft, ifft
 from scipy.signal import detrend
 from numpy import angle
 
-# 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
+g = 9.81
 
 # h = profondeur locale
-#fs = fréquence echantillonage
+# 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)
+def PUV_dam(u, p, h, fs, zmesP, zmesU):
+    # u = detrend(u)
+    # p = detrend(p)
     ampliseuil = 0.001
     N = len(u)
     delta = fs
-    #transformée de fourrier
+    # 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/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)/delta*nyquist;
-    #phase
-    ThetaU=angle(Yu[1:N//2]);
-    ThetaP=angle(Yp[1:N//2]);
-    #calcul du coefficient de reflexion
-    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
+    nyquist = 1 / 2
+    # Fu_xb = ((1:N/2)-1)/(N/2)/deltat*nyquist
+    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) / delta * nyquist
+    # phase
+    ThetaU = angle(Yu[1 : N // 2])
+    ThetaP = angle(Yp[1 : N // 2])
+    # calcul du coefficient de reflexion
+    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 = []
-        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
+        omega = []
+        ccc, ccs, cccp = [], [], []
+        aincident13, areflechi13 = [], []
+        aprog = []
+        r13 = []
+        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 * np.pi * xf[ii])
+            cc = omega[ii] * np.cosh(xf[ii] * (zmesU + h)) / np.sinh(xf[ii] * h)
+            ccp = omega[ii] * np.cosh(xf[ii] * (zmesP + h)) / np.sinh(xf[ii] * h)
+            cs = omega[ii] * np.sinh(xf[ii] * (zmesU + h)) / np.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 * np.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] * np.cos(phi3 - phi1) / (cc * ccp * omega[ii])
+            ))
+            areflechi13.append(0.5 * np.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] * np.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[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'
+                r13.append(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"
 
 
 # Calcul du vecteur d'onde en prÈsence d'un courant
@@ -100,19 +116,29 @@ def PUV_dam(u,p,h,fs,zmesP,zmesU):
 # 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) :
+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) :
+    u = -u
+    i = 0
+    while abs((kh - x) / x) > 0.00001:
+        i += 1
+        if i > 10**5:
+            sys.exit(1)
         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
+        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
 
+# fin du calcul de k
diff --git a/swash/processing/pa/reflex3S.py b/swash/processing/pa/reflex3S.py
index 3a6496c..84b695a 100644
--- a/swash/processing/pa/reflex3S.py
+++ b/swash/processing/pa/reflex3S.py
@@ -8,54 +8,68 @@ from matplotlib import animation
 import cmath
 from scipy.fft import fft
 
-def newtonpplus(f,h,u) :
+
+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) :
+    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
+        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) :
+
+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
+    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 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;
+    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) :
+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
@@ -69,77 +83,76 @@ def reflex3S(x1,x2,x3,xs1,xs2,xs3,h,mean_freq,fmin,fmax) :
     # ampliseuil : amplitude minimale des maxima recherches
     # valeur du courant (u > 0 correspond a un courant dans le sens
     # des x croissants)
-    #u=0;
-    
+    # 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
+
+    # 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;
+
+    # 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));
-    
-    
+    # 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 = []
@@ -155,66 +168,104 @@ def reflex3S(x1,x2,x3,xs1,xs2,xs3,h,mean_freq,fmin,fmax) :
     Ereflechi123 = []
     count = 0
     for jj in np.arange(indfmin, indfmax, 1):
-        f=freq1[jj]
-        #periode
-        T.append(1/f)
+        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
+        # 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
+    return ai, ar, Eincident123, Ereflechi123, indfmin, indfmax, fre
diff --git a/swash/processing/post_reflex.py b/swash/processing/post_reflex.py
new file mode 100644
index 0000000..7833f00
--- /dev/null
+++ b/swash/processing/post_reflex.py
@@ -0,0 +1,83 @@
+import argparse
+import configparser
+import logging
+import pathlib
+
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy.signal as sgl
+
+from .reflex3s import reflex3s
+from .pa.reflex3S import reflex3S
+from .pa.PUV import PUV_dam
+
+parser = argparse.ArgumentParser(description="Post-process swash output")
+parser.add_argument("-v", "--verbose", action="count", default=0)
+parser.add_argument("-c", "--config", default="config.ini")
+args = parser.parse_args()
+
+logging.basicConfig(level=max((10, 20 - 10 * args.verbose)))
+log = logging.getLogger("post")
+
+log.info("Starting post-processing")
+config = configparser.ConfigParser()
+config.read(args.config)
+
+inp = pathlib.Path(config.get("post", "inp"))
+root = pathlib.Path(config.get("swash", "out"))
+
+log.info(f"Reading data from '{inp}'")
+x = np.load(inp.joinpath("xp.npy"))
+t = np.load(inp.joinpath("tsec.npy"))
+
+botl = np.load(inp.joinpath("botl.npy"))
+watl = np.load(inp.joinpath("watl.npy"))
+vel = np.load(inp.joinpath("vel.npy"))
+
+x0 = config.getint("post", "x0")
+arg_x0 = np.abs(x - x0).argmin()
+t0 = config.getfloat("post", "t0")
+arg_t0 = np.abs(t - t0).argmin()
+dt = np.diff(t).mean()
+f = 1 / dt
+
+idx = [arg_x0 - 5, arg_x0, arg_x0 + 7]
+wl = watl[arg_t0:, idx]
+
+# Reflex3S
+res = reflex3S(*wl.T, *x[idx], botl[idx].mean(), f, 0.02, 0.2)
+f_, ai_, ar_ = reflex3s(wl.T, x[idx], botl[idx].mean(), f)
+ai, ar, _, _, _, _, fre = res
+
+# Custom PUV
+eta = watl[t > t0, arg_x0]
+u = vel[t > t0, 0, arg_x0]
+
+phi_eta = sgl.welch(eta, f)
+phi_u = sgl.welch(u, f)
+phi_eta_u = sgl.csd(eta, u, f)
+
+R = np.sqrt(
+    (np.abs(phi_eta[1]) + np.abs(phi_u[1]) - 2 * phi_eta_u[1].real)
+    / (np.abs(phi_eta[1]) + np.abs(phi_u[1]) + 2 * phi_eta_u[1].real)
+)
+
+out = pathlib.Path(config.get("post", "out")).joinpath(f"t{t0}x{x0}")
+log.info(f"Saving plots in '{out}'")
+out.mkdir(parents=True, exist_ok=True)
+
+plt.plot(fre, np.asarray(ar) / np.asarray(ai))
+plt.xlim((0, 0.3))
+plt.ylim((0, 1))
+plt.savefig(out.joinpath("reflex3s_pa.pdf"))
+plt.clf()
+plt.plot(f_, ar_ / ai_)
+plt.plot(phi_eta[0], R)
+plt.xlim((0, 0.3))
+plt.ylim((0, 1))
+plt.savefig(out.joinpath("reflex3s.pdf"))
+
+# pressk = np.load(inp.joinpath("pressk.npy"))
+# press = pressk[:, 1, :]
+# res = PUV_dam(vel[arg_t0:, 0, 500], press[arg_t0:, 500], botl[500], f, botl[500]/2, botl[500]/2)
+# print(res)
diff --git a/swash/processing/reflex3s.py b/swash/processing/reflex3s.py
new file mode 100644
index 0000000..07fce52
--- /dev/null
+++ b/swash/processing/reflex3s.py
@@ -0,0 +1,45 @@
+import numpy as np
+from scipy import fft
+from scipy import signal as sgl
+from scipy import optimize as opti
+
+
+def reflex3s(eta, x, h, f):
+    eta_psd = np.stack([fft.rfft(z) for z in eta])
+    f_psd = fft.rfftfreq(eta.shape[1], 1 / f)
+
+    eta_amp = np.abs(eta_psd) / (f_psd.size - 1)
+    eta_phase = np.angle(eta_psd)
+
+    g = 9.81
+    k = np.asarray(
+        [
+            opti.root_scalar(
+                f=lambda k: k * np.tanh(k) - (2 * np.pi * f) ** 2 / g * h,
+                fprime=lambda k: np.tanh(k) + k * (1 - np.tanh(k) ** 2),
+                x0=0.5,
+            ).root
+            / h
+            for f in f_psd
+        ]
+    )
+
+    s1 = 1 + np.exp(2j * k * (x[1] - x[0])) + np.exp(2j * k * (x[2] - x[0]))
+    s2 = 1 + np.exp(-2j * k * (x[1] - x[0])) + np.exp(-2j * k * (x[2] - x[0]))
+    s12 = 3
+    s3 = (
+        eta_amp[0] * np.exp(-1j * (eta_phase[0]))
+        + eta_amp[1] * np.exp(-1j * (eta_phase[1] - k * (x[1] - x[0])))
+        + eta_amp[2] * np.exp(-1j * (eta_phase[2] - k * (x[2] - x[0])))
+    )
+    s4 = (
+        eta_amp[0] * np.exp(-1j * (eta_phase[0]))
+        + eta_amp[1] * np.exp(-1j * (eta_phase[1] + k * (x[1] - x[0])))
+        + eta_amp[2] * np.exp(-1j * (eta_phase[2] + k * (x[2] - x[0])))
+    )
+    s5 = s1 * s2 - s12**2
+
+    ai = np.abs((s2 * s3 - s12 * s4) / s5)
+    ar = np.abs((s1 * s4 - s12 * s3) / s5)
+
+    return f_psd, ai, ar