# -*- coding: utf-8 -*-
"""
Created on Tue May 28 14:13:10 2024

@author: M3051168
"""

import numpy as np
import matplotlib.pyplot as plt

from pack_marc.p_fct_etude import *
from pack_marc.p_fourier import *

cmap=plt.get_cmap('winter',11)
L=400

def flutter(h,L=400,n_jet=0):
    return np.sqrt(2*n_jet+1)/(2*(np.pi))*np.sqrt(9.81/(0.001*h))

def airy(h,L=400,n=0):
    """Calcul la fréquence de résonance théorique de la cuve, les données en mm"""
    g=9.806
    d_b=22
    return np.sqrt((2*n+1)*g/(4*np.pi*L*0.001)*np.tanh((2*n+1)*np.pi*h/L))


def amortissement_red(h,L_1,n,L_2=100):
    nu=1.002*0.000001 #viscosité USI
    h=0.001*h
    L_1=0.001*L_1
    L_2=0.001*L_2
    omega=airy(h,L_1,n)
    tau_l_inv=np.sqrt(nu*omega)*(1/L_1+1/L_2) #1/tau avec tau le taux d'amortissement_red
    lambda_s=1/(2*omega*np.sqrt(1+omega**2*1/tau_l_inv**2))
    return lambda_s

lambda_2=0.1
liste_a=np.linspace(0.01,5,5000)
liste_h=np.linspace(30,400,100)


plt.figure(dpi=500) 
##  
n_surf=0
n_jet=0
i=0
seuil=[]
omega_im=[]
                                                  
liste_omega_1=[airy(liste_h[i],400,n_surf) for i in range(len(liste_h))]
liste_omega_2=[flutter(liste_h[i],400,n_jet) for i in range(len(liste_h))]
liste_lambda_1=[amortissement_red(liste_h[i],400,n_surf) for i in range(len(liste_h))]
ind_a=10
for i in range(len(liste_omega_1)):
    lambda_1=liste_lambda_1[i]
    res_1=[]
    res_2=[]
    for a in liste_a:
        mat_4=np.array([[0, 1, 0, 0],
                    [-liste_omega_1[i]**2, -lambda_1, -a, 0],
                    [0, 0, 0, 1],
                    [a, 0, -liste_omega_2[i]**2, -lambda_2]])
        lamb,v=np.linalg.eig(mat_4)
        idx = lamb.argsort()[::-1]   
        lamb = lamb[idx]
        res_1.append(np.real(lamb)[0])
        res_2.append(np.imag(lamb)[0])
    seuil.append(liste_a[int(np.argwhere(np.array(res_1)>0.01)[0][0])])
    omega_im.append(res_2[ind_a])    
    i=i+1

omega_min=liste_omega_1[np.argmin(seuil)]

plt.plot(liste_h,seuil)
##
n_surf=1
n_jet=0
i=0
seuil=[]
omega_im=[]
                                                  
liste_omega_1=[airy(liste_h[i],400,n_surf) for i in range(len(liste_h))]
liste_omega_2=[flutter(liste_h[i],400,n_jet) for i in range(len(liste_h))]
liste_lambda_1=[amortissement_red(liste_h[i],400,n_surf) for i in range(len(liste_h))]
ind_a=10
for i in range(len(liste_omega_1)):
    lambda_1=liste_lambda_1[i]
    res_1=[]
    res_2=[]
    for a in liste_a:
        mat_4=np.array([[0, 1, 0, 0],
                    [-liste_omega_1[i]**2, -lambda_1, -a, 0],
                    [0, 0, 0, 1],
                    [a, 0, -liste_omega_2[i]**2, -lambda_2]])
        lamb,v=np.linalg.eig(mat_4)
        idx = lamb.argsort()[::-1]   
        lamb = lamb[idx]
        res_1.append(np.real(lamb)[0])
        res_2.append(np.imag(lamb)[0])
    seuil.append(liste_a[int(np.argwhere(np.array(res_1)>0.01)[0][0])])
    omega_im.append(res_2[ind_a])    
    i=i+1

omega_min=liste_omega_1[np.argmin(seuil)]

plt.plot(liste_h,seuil)
##
n_surf=0
n_jet=1
i=0
seuil=[]
omega_im=[]
liste_h=np.linspace(100,400,100)

lambda_2=0.1
                                                  
liste_omega_1=[airy(liste_h[i],400,n_surf) for i in range(len(liste_h))]
liste_omega_2=[flutter(liste_h[i],400,n_jet) for i in range(len(liste_h))]
liste_lambda_1=[amortissement_red(liste_h[i],400,n_surf) for i in range(len(liste_h))]
ind_a=10
for i in range(len(liste_omega_1)):
    res_1=[]
    res_2=[]
    lambda_1=liste_lambda_1[i]
    for a in liste_a:
        mat_4=np.array([[0, 1, 0, 0],
                    [-liste_omega_1[i]**2, -lambda_1, -a, 0],
                    [0, 0, 0, 1],
                    [a, 0, -liste_omega_2[i]**2, -lambda_2]])
        lamb,v=np.linalg.eig(mat_4)
        idx = lamb.argsort()[::-1]   
        lamb = lamb[idx]
        res_1.append(np.real(lamb)[0])
        res_2.append(np.imag(lamb)[0])
    seuil.append(liste_a[int(np.argwhere(np.array(res_1)>0.01)[0][0])])
    omega_im.append(res_2[ind_a])    
    i=i+1

omega_min=liste_omega_1[np.argmin(seuil)]

plt.plot(liste_h,seuil)




plt.xlabel(r'$h (mm)$')
plt.ylabel(r'$\alpha_c$')
plt.ylim([0,2.5])
plt.grid()