internship/swash/processing/pa/PUV.py

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
from mpl_toolkits import mplot3d
from scipy import signal
from scipy.fftpack import fft, ifft
from scipy.signal import detrend
from numpy import angle

g = 9.81

# 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 = 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 / 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 = []
        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.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
# 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
    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
    return k


# fin du calcul de k