%function y=Lamb()
%BONNE VERSION

% This function compute the dispersion curve omega(k) for S0 and A0 Lamb
% modes using a mean square fit to minimize the equation dispersion given
% in Royer page ?
% X is wavevector x h / pi
% Y is omega x h / pi
% Need to give pin=Y0 for initial fit at X=2 (where S0 and A0 are almost equal
% to the Rayleigh dispersion)
% Than, X decreases to zero, using pout as the new pin at each step.
% (this is because it is hard to find pin for small values of X)
% (To calculate higher modes, must know cutoff frequency and go increasing
% with bounds pinmin & pinmax (fminsearchbnd) to avoid jumps from on mode
% to the other).
% Patrick July 2007
%figure
h=.352/2; % half-thickness in mm
%%% 100
Vt=4630;  % m/s
% Vt=5840;
Vl=9134;  % m/s
a=5.2e-2;
b=3.6e-2;
% autre : APL 88,251907(2006)
% Vt=5360;
% Vl=8950;
% %%% 110
% Vt=4673;  % m/s
% Vl=9134;  % m/s
% %%% 111
% Vt=5085;  % m/s
% Vl=9360;  % m/s
% ajustement
%Vt=4000;  % m/s
% Vl=6360; %Alu
% Vt=3140;
X=0:0.001:0.3;

%alpha=0;    % 0 si mode symétrique / pi/2 si mode antisymétrique
Z=[];

options = optimset('Display','notify','MaxFunEvals',100000,'TolFun',1e-8,'MaxIter',20000,'TolX',1e-12*10000);
% Display : Level of display.
    % 'off' displays no output;
    % 'iter' displays output at each iteration;
    % 'final' displays just the final output;
    % 'notify' (default) displays output only if the function does not converge.
    % FunValCheck: Check whether objective function values are valid.
    % 'on' displays a warning when the objective function returns a value that is complex, Inf or NaN.
    % 'off' (the default) displays no warning.
% MaxFunEvals: Maximum number of function evaluations allowed
% MaxIter: Maximum number of iterations allowed
% OutputFcn: Specify a user-defined function that the optimization function calls at each iteration.
% TolFun: Termination tolerance on the function value
% TolX: Termination tolerance on x

%figure

%% Modes anti-symétriques
alpha=pi/2;
pin=1200;          % initial value, to be changed for m>0
Y=[];
for x=fliplr(X)
    [pout,fval,exitflag]=    fminsearch(@(y) dispersion(x,y),pin,options);
    pin=pout;
    Y=[Y pout];
end
Z=[Z;fliplr(Y)];

figure,
%  plot(pi*X/h,Z*0.001/(2*h))
%k=sqrt((n*pi/a)^2+(m*pi/b)^2);
plot(pi*X/h*1000,1e6*Z*0.001/(2*h),'r')
ylabel('f(Hz)')

%plot(pi*X/h*1000,1e6*2*pi*Z*0.001/(2*h))

%xlim([0,max(X)])
% xlim([0 1200])
%ylim([0 2.5])
title(['dispersion: variation de v en fonction de k epaisseur:h=',num2str(2*h),'mm'])
 xlabel('k(m^{-1})')
%xlabel('\lambda(mm)')
ylabel('\omega(rad/s)')
%ylim([0 1.6])
%xlim([0 5])
%legend('S0','A0','Location','NorthWest');
y=Z;    % output from function

hold on

% Vp=2*Vt*(1-(Vt^2/Vl^2));
% ohmega=1.15*(Vp./sqrt(3)).*(pi*X/h*1000).^2*h*1e-3;
%
%
% plot(pi*X/h*1000,ohmega/(2*pi))
%
% h = legend('methode numerique','equation theorique',2);
% set(h,'Interpreter','none')
% ylim([0 1e5])

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     function sse=dispersion(x,y)
%         p2=y^2/Vl^2-x^2;     % (ph/pi)^2
%         q2=y^2/Vt^2-x^2;     % (qh/pi)^2
%         sse = abs(y^4/Vt^4 - (4*x^2*q2*(1-sqrt(p2)/sqrt(q2)*tan(pi*sqrt(p2)+alpha)/tan(pi*sqrt(q2)+alpha))))^2;
%        sse = abs(y^4/Vt^4 - (4*x^2*(y^2/Vt^2-x^2)*(1-sqrt((y^2/Vl^2-x^2)/(y^2/Vt^2-x^2))*tan(pi*sqrt((y^2/Vl^2-x^2))+alpha)/tan(pi*sqrt((y^2/Vt^2-x^2))+alpha))))^2;
%        sse = abs(y^4/Vt^4 - (4*x^2*(y^2/Vl^2-x^2)*(1-sqrt((y^2/Vt^2-x^2)/(y^2/Vl^2-x^2))*tan(pi*sqrt((y^2/Vt^2-x^2))+alpha)/tan(pi*sqrt((y^2/Vl^2-x^2))+alpha))))^2;
% Check in Royer the correct expression
%     end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end