/// (S. Ayrinhac, IMPMC, Sorbonne Université, 2024
// 
// Calcule la vitesse de l'onde de Rayleigh à partir des vitesses longi. et transverses
// Exemple : "vR=Func_Rayleigh_velocity(5960,3400,2);"

   // option 1 : Viktorov equation
   // option 2 : Rayleigh equation
   
function vR=Func_Rayleigh_velocity(vL,vT,option)

     r=vT/vL;
     vinit=vT*(0.718-r^2)/(0.75-r^2); // D. Royer & E. Dieulesaint, Tome 1, p.251
      
   if option==1
     vR=vinit;
   else
     // Attention il faut choisir une bonne valeur de départ 
     // Ici cette valeur de départ est donnée par l'approximation de Viktorov
     [vR,v,info]=fsolve(vinit,list(Rayleigh_equation,vL,vT));
   end

endfunction


function y=Rayleigh_equation(x,vL,vT)
    
    vR2=x.^2;
    vL2=vL.^2;
    vT2=vT.^2;
    y=16.*(1-vR2/vL2).*(1-vR2/vT2)-(2-vR2/vT2).^4; // D. Royer & E. Dieulesaint, Tome 1, p.250
    
endfunction