// Simulation Acoustique - Modele Limp Frame
// -------------------------------------------------------------------

// Color : red=21, bleu=11, green=13, purple=24, yellow=32
style = 13;    

// Air
rho0 = 1.21;
c0   = 343;
Z0   = rho0 * c0;

// Paramètres "limp frame" 
// RockFeu 85 kg/m3 AFR 40000 best 12000
// Pavacell 40 kg/M3 AFR 10000 best 8000
phi_t = 0.98;    // Porosité
sig_t = 8000 ;   // resistivite au passage de l'air
rho_t = 40 ;     // densité du materiau
e_t = 0.1 ;      // épaisseur

// Couche sd'air
e_air = 0.2 ; // 0.4;       // epaisseur coté mur

// Laine de verre plenum
e_p = 0.8 ; // 0.6      // epaisseur
sig_p = 4000;  // resistivite au passege de l'air

// --- FONCTION MODELE LIMP (Masse + Poreux) ---
function [M,rho_l] = matrice_t(L, sig, rho_mat)
    omega = 2 * %pi * f;

    // 1) Calcul du fluide equivalent (Miki) 
    f_s = f / sig;
    Z_p = Z0 * (1 + 0.070*f_s^-0.632 - %i*0.107*f_s^-0.632);
    k_p = (omega/c0) * (1 + 0.109*f_s^-0.618 - %i*0.160*f_s^-0.618);
    gamma_l = k_p ;
    Z_l     = Z_p ;    
    
    // Test si modele limp
    if (rho_mat <> 0) then

      // Densité effective associée au couple (Z_p, k_p) :
      rho_p = (Z_p * k_p) / omega;
    
      // Prise en compte LIMP
      rho_l = (rho_p * rho_mat - rho0^2) / (rho_p + rho_mat - 2*rho0);

      // Nouvelles proprietes du milieu fusionne 
      ratio   = sqrt(rho_l / rho_p);
      gamma_l = k_p * ratio;
      Z_l     = Z_p * ratio;
    end

    // Matrice de transfert
    M = [cos(gamma_l*L),  %i*Z_l*sin(gamma_l*L);
         (%i/Z_l)*sin(gamma_l*L),  cos(gamma_l*L)];
endfunction


// --- BOUCLE DE CALCUL ---
freqs = logspace(log10(10), log10(1000), 500);
alpha = zeros(1, length(freqs));

for i = 1:length(freqs)
    f = freqs(i);
    T = [1, 0; 0, 1];

    // Couche 1 limp
    [Mat,rho_tab(i)] = matrice_t(e_t, sig_t, rho_t);
    T = T * Mat;
    
    // Couche 2 : Air
    k_air = (2*%pi*f) / c0;
    T = T * [cos(k_air*e_air),  %i*Z0*sin(k_air*e_air);
            (%i/Z0)*sin(k_air*e_air),  cos(k_air*e_air)];

    // Couche 3 : LDV 
    T = T * matrice_t(e_p, sig_p, 0);

    // Impédance de surface (terminaison rigide implicite via T)
    Zs = T(1,1) / T(2,1);

    // Absorption incidence normale
    alpha(i) = 1 - abs((Zs - Z0) / (Zs + Z0))^2;
end

// Calcul fréquende de transition
f_T = sig_t*phi_t/(2*%pi*rho_t) ;
disp(f_T) ;

// Trace coefficient absorption
scf(0);
plot2d("ln", freqs, alpha, style=style, rect=[10,0,1000,1]);
h=gca();
h.grid=[1,1]
h.log_flags="lnn"; // x:log,y=log,z=lin
h.sub_ticks=[0,0];

// X ticks
ticks = h.x_ticks
ticks.labels = ["10";"20";"50";"100";"200";"500";"1000"];
ticks.locations = [10;20;50;100;200;500;1000];
h.x_ticks = ticks;
xtitle("Absorption - LIMP FRAME 40 kg/m3", "Hz", "Alpha");
legends(['AFR10';'AFR8'],[21,13],opt="lr")
//legends(['AFR'],[13],opt="lr")

// Trace densité effecive
scf(1);
subplot(2,1,1) ;
plot2d(freqs, real(rho_tab/rho0), style=3,rect=[0,0,1000,30]);
subplot(2,1,2) ;
plot2d(freqs, imag(rho_tab/rho0), style=3,rect=[0,-30,1000,0]);
