internship/swash/processing/pa/reflex3S.py

import csv
import matplotlib.pyplot as plt
import os
import math
import sys
import numpy as np
from matplotlib import animation
import cmath
from scipy.fft import fft

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