%求解Lame波對(duì)稱部分的基本模態(tài)So
format long
syms c;   %定義速度c為符號(hào)函數(shù)
xx=[];err=[];      %xx存貯不同KH值的波速，err存貯誤差
eps=1e-12;         %精度
u=0.3; %泊松比u，越大對(duì)應(yīng)速度越小
E=2.06e11; %彈性模量E,越大對(duì)應(yīng)速度越大
p=7850;%密度越大對(duì)應(yīng)速度越小
cd=((1-u)*E/(1+u)/(1-2*u)/p)^0.5;%cd為縱波波速，cs為橫波波速
cs=(E/2/(1+u)/p)^0.5; 
for kh=0:0.05:8       %k為波數(shù)，h為總的板厚
  if kh<7.65          %利用lamea0作圖法先估計(jì)出c值范圍
     y1=4*tanh(kh/2*(1-(c/cd)^2)^(1/2))*(1-(c/cd)^2)^(1/2)*((c/cs)^2-1)^(1/2);
     y2=(2-(c/cs)^2)^2*tan(kh/2*((c/cs)^2-1)^(1/2));
     if kh<0.05
        c=[5365;5375];
     elseif kh<2.55  c=[cs;5400];
     elseif kh<3  c=[cs;xx(fix(20*kh))];
     else   c=[cs;xx(fix(20*kh))];    
     end
     for i=2:2
         a=c(i-1);b=c(i);
         c=(a+b)/2;
        while (subs(y1-y2)~= 0) & ((b-a)>eps)
              c=(a+b)/2;
              fa=4*tanh(kh/2*(1-(a/cd)^2)^(1/2))*(1-(a/cd)^2)^(1/2)*((a/cs)^2-1)^(1/2)-(2-(a/cs)^2)^2*tan(kh/2*((a/cs)^2-1)^(1/2));
              fb=4*tanh(kh/2*(1-(b/cd)^2)^(1/2))*(1-(b/cd)^2)^(1/2)*((b/cs)^2-1)^(1/2)-(2-(b/cs)^2)^2*tan(kh/2*((b/cs)^2-1)^(1/2));
              fy=4*tanh(kh/2*(1-(c/cd)^2)^(1/2))*(1-(c/cd)^2)^(1/2)*((c/cs)^2-1)^(1/2)-(2-(c/cs)^2)^2*tan(kh/2*((c/cs)^2-1)^(1/2));
            if fa*fy < 0
               b=c;
            else
               a=c;
            end
        end
        xx=[xx;c];
        m=subs(y1-y2);
        err=[err;m];
     end
  else
     y1=4*tanh(kh/2*(1-(c/cd)^2)^(1/2))*(1-(c/cd)^2)^(1/2)*(1-(c/cs)^2)^(1/2);
     y2=(2-(c/cs)^2)^2*tanh(kh/2*(1-(c/cs)^2)^(1/2));
     c=[0;cs];
      for i=2:2     %利用二分法求根
          a=c(i-1);b=c(i);
          c=(a+b)/2;
          while (subs(y1-y2)~= 0) & ((b-a)>eps)
                c=(a+b)/2;
                fa=4*tanh(kh/2*(1-(a/cd)^2)^(1/2))*(1-(a/cd)^2)^(1/2)*(1-(a/cs)^2)^(1/2)-(2-(a/cs)^2)^2*tanh(kh/2*(1-(a/cs)^2)^(1/2));
                fb=4*tanh(kh/2*(1-(b/cd)^2)^(1/2))*(1-(b/cd)^2)^(1/2)*(1-(b/cs)^2)^(1/2)-(2-(b/cs)^2)^2*tanh(kh/2*(1-(b/cs)^2)^(1/2));
                fy=4*tanh(kh/2*(1-(c/cd)^2)^(1/2))*(1-(c/cd)^2)^(1/2)*(1-(c/cs)^2)^(1/2)-(2-(c/cs)^2)^2*tanh(kh/2*(1-(c/cs)^2)^(1/2));
              if fa*fy < 0
                 b=c;
              else
                 a=c;
              end
          end
          xx=[xx;c];
          m=subs(y1-y2);
          err=[err;m];
      end
  end
end
cp=xx;        %計(jì)算So群速度
kh=[0:0.05:8]; 
for i=2:1:161 
    if cp(i-1) ~=cp(i)
    cg(i)=cp(i)+kh(i)*(cp(i)-cp(i-1))/(kh(i)-kh(i-1)); %波的相速度cp與群速度cg的轉(zhuǎn)換公式
    else cg(i)=cg(i-1)-50;
    end
end
cg(1)=0;
%作圖
fh=[];
for kh=0:0.05:8
    f=kh*cp(fix(20*kh)+1)/2/pi; %將kd與cp、cg的關(guān)系轉(zhuǎn)換為fh與cp、cg的關(guān)系
    fh=[fh;f];
end    
figure(1)     %f與cp的圖
plot(fh/2,cp);%板厚為2mm
xlabel('頻率f (KHz)');ylabel('相速度Cp (m/s)');
grid on;
hold on;
figure(2)     %f與cg的圖
plot(fh/2,cg);%板厚為2mm
xlabel('頻率f (KHz)');ylabel('群速度Cg (m/s)');
grid on;
hold on;