// Simon Ayrinhac, IMPMC, Sorbonne Université , (2024)
// 
// Vitesse de l'onde de Rayleigh ( de nombre d'onde k)dans un échantillon stratifié (d'épaisseur h)
// 
// D'après l'article de :
//     G.W. Farnell and E.L. Adler, "Elastic wave propagation in thin layers", in physical acoustics. NY: Academic, New York, 1972, Eq. (18), p. 48
// 
//  IMPORTANT, à noter : il y a une coquille dans la 5e ligne et 4e colonne de la matrice, '2b^(5)e^(ib^(6)kh)' doit être remplacé par '2b^(6)e^(ib^(6)kh)'
//
// (Exemple en fin de fichier)


function vR=Func_Farnell_velocity(vinit,kh,vL,vT,rho,vlL,vlT,rhol)
   
   // --- Paramètres d'entrée de la fonction ---------------------------------------
   //   vinit : IMPORTANT, vitesse initiale pour la fonction 'fsolve'
   //   kh : produit nombre d'onde par épaisseur de la couche
   //   vL, vT : vit. longi. et trans. du substrat
   //   vlL, vlT : vit. longi. et trans. de la couche mince ('l' pour 'layer' en anglais)
   //   rho, rhol : masses volumiques du substrat et de la couche ('l' pour 'layer' en anglais)
   
   
   // r : rapport des constantes élastiques transverses (article de Farnell p.48, en-dessous de l'Eq.(18))
     C44=rho*vT^2;
     Cl44=rhol*vlT^2;
     r=C44/Cl44;
     
   // attention il faut choisir une bonne valeur de départ
   [vR0,v,info]=fsolve(vinit,list(Farnell_equation,kh,r,vL,vT,vlL,vlT));
   
   res=Farnell_equation(vR0,kh,r,vL,vT,vlL,vlT);
   if res~=0
       vR=vR0;
   else
       vR=0;
   end
    
endfunction


function y=Farnell_equation(x,kh,r,vL,vT,vlL,vlT)
    
bc=-%i*sqrt(1-(x/vT)^2); // substrat
bd=-%i*sqrt(1-(x/vL)^2);
bc2=bc^2;
bd2=bd^2;

b5=%i*sqrt(1-(x/vlT)^2); // couche ('layer')
b6=%i*sqrt(1-(x/vlL)^2);
b52=b5^2;
b62=b6^2;


// Farnell (1972), Eq. (18), p. 48
M=[b5 -1 -b5 -1 -bc 1 ; 
-1 -b6 -1 b6 1 bd ; 
(1-b52) 2*b6 (1-b52) -2*b6 -r*(1-bc2) -2*r*bd ; 
2*b5 -(1-b52) -2*b5 -(1-b52) -2*r*bc r*(1-bc2) ; 
(1-b52)*exp(%i*b5*kh) 2*b6*exp(%i*b6*kh) (1-b52)*exp(-%i*b5*kh) -2*b6*exp(-%i*b6*kh) 0 0 ; 
2*b5*exp(%i*b5*kh) -(1-b52)*exp(%i*b6*kh) -2*b5*exp(-%i*b5*kh) -(1-b52)*exp(-%i*b6*kh) 0 0 ];

[e,m]=det(M);
y0=m*10^e;
y=abs(y0);
    
endfunction



// Exemple ----------------------------------------------------------
// substrat : verre sodocalcique, couche mince : aluminium (A1050)
//q=1;
//for kh0=0:0.1:5
//  kh(q)=kh0;
//  vR(q)=Func_Farnell_velocity(3000,kh0,5840,3460,2490,6470,3150,2700);
//  q=q+1
//end
//  figure(1)
//    plot(kh,vR,'or')