// S. Ayrinhac (04/2025) - IMPMC
// Modes propres de vibrations d'une lame 
//  - matériau homogène et isotrope
//  - les deux extrémités libres

function [fnb,fni,fnL,fntor]=Func_vibs_lames(dim,mat)
    
    // INPUT -----------------------------
      // dim : vecteur dimensions [longueur largeur épaisseur]
      // Ex. : dim=[300e-3 40e-3 10e-3]; // barre alu IMPMC
      // mat : vecteur élasticité du matériau [masse volumique vit_longi vi_transverse]
      // Ex. mat=[2700 6421.93 3127.5]; // vitesses alu from ultrasons
    // OUTPUT -----------------------------
      // fnb : modes de 'bending'
      // fni : modes de 'in-plane bending'
      // fnL : modes de compression
      // fntor : modes de torsion
    
    affichage=0;
    
    

    
    
L=dim(1); // longueur
l=dim(2); // largeur
h=dim(3); // hauteur
  K=h/sqrt(12);  // radius of gyration
  Ki=l/sqrt(12); // radius of gyration for in-plane bending
  // --- https://www.acs.psu.edu/drussell/Demos/Torsional/torsional.html ------------
  AJ=1-(h^4)/(12*l^4);
  Jexact=((l*h^3)/16)*(16/3-3.36*(h/l)*AJ); // torsional rigidity
  Iexact=((l*h)/12)*(l^2+h^2);  // moment of inertia
rho=mat(1); // masse volumique
cL=mat(2);  // vitesse longi
cT=mat(3);  // vitesse transverse

  
// Grandeurs dérivées

cL2=cL^2;
cT2=cT^2;
G=(rho*cT2);
E=(rho*cT2)*(3*cL2-4*cT2)/(cL2-cT2);
nu=(cL2-2*cT2)/(2*cL2-2*cT2);
M=E*(1-nu)/((1+nu)*(1-2*nu)); // P-wave modulus
cE=sqrt(E/rho);
ctor=(2*h/l)*sqrt(G/rho);
ctorexact=sqrt((G*Jexact)/(rho*Iexact)); // voir site de Russell et Garrett book p.246
   

if affichage==1
  disp('cL = '+string(cL)+' m/s');
  disp('cT = '+string(cT)+' m/s');
    
disp('K ='+string(K));
disp('Jexact ='+string(Jexact));
disp('Iexact ='+string(Iexact));
    disp('G ='+string(G/1e9)+' GPa');
    disp('E ='+string(E/1e9)+' GPa');
    disp('c_E ='+string(cE)+' m/s');
    disp('c_torsion ='+string(ctor)+' m/s');
    disp('c_torsion_exact ='+string(ctorexact)+' m/s');
end

/// Modes en flexion = bending ----------------------------------
nb=1:10;
fnb=zeros(nb);  // initialisation
sb=ones(nb).*1;    // type de mode
Ab=ones(nb).*1; // amplitude des pics

for i=1:length(nb)
  
        if nb(i)==1
           fnb(i)=((%pi.*K)./(8.*L.^2)).*sqrt(E./rho).*(3.01124).^2;  
        end
        if nb(i)==2
           fnb(i)=((%pi.*K)./(8.*L.^2)).*sqrt(E./rho).*(4.99951).^2;  
        end
        if nb(i)==3
           fnb(i)=((%pi.*K)./(8.*L.^2)).*sqrt(E./rho).*(7.00002).^2;  
        end
    if nb(i)>3
       fnb(i)=((%pi.*K)./(8.*L.^2)).*sqrt(E./rho).*(2.*nb(i)+1).^2; 
    end
end

/// Modes in-plane bending ----------------------------------
ni=1:10;
fni=zeros(ni); // initialisation
si=ones(ni).*2;    // type de mode
Ai=ones(ni).*0.3; // amplitude des pics

for i=1:length(ni)
        if ni(i)==1
           fni(i)=((%pi.*Ki)./(8.*L.^2)).*sqrt(E./rho).*(3.01124).^2;  
        end
        if ni(i)==2
           fni(i)=((%pi.*Ki)./(8.*L.^2)).*sqrt(E./rho).*(4.99951).^2;  
        end
        if ni(i)==3
           fni(i)=((%pi.*Ki)./(8.*L.^2)).*sqrt(E./rho).*(7.00002).^2;  
        end
    if ni(i)>3
       fni(i)=((%pi.*Ki)./(8.*L.^2)).*sqrt(E./rho).*(2.*ni(i)+1).^2; 
    end
end

// Modes de compression suivant la longueur
nL=1:10;
sL=ones(nL).*3;    // type de mode
AL=ones(nL).*0.02; // amplitude des pics
fnL=(nL.*cE)./(2.*L);

// Modes de torsion -------------------------------
ntor=1:10;
stor=ones(ntor).*4;   // type de mode
Ator=ones(ntor).*0.2; // amplitude des pics
fntor=(ntor.*ctorexact)./(2.*L);

//fnb=fnb';
//fni=fni';
//fnL=fnL';
//fntor=fntor';

// f=[fnb fni fnL fntor];
// n=[nb ni nL ntor];
// s=[sb si sL stor];
// A=[Ab Ai AL Ator];


endfunction



//[fnb,fni,fnL,fntor]=Func_vibs_lames([300e-3 40e-3 10e-3],[2700 6421.93 3127.5]);
//disp([fnb fni fnL fntor])
//
//abort