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