clc

clear

close all

n1=0;

syms c;   %定义速度c为符号函数

err=[];      %xx，err存贮误差

eps=1e-12;         %精度

ct=3.130;cl=6.350;%铝板的横波速度3130m/s，纵波速度6350m/s   这边的速度最好是给出的个位数，如果太大，完蛋！那tanh会趋近于0，所以图形就不对

d=1;

m=0.5:0.1:3.1;

A=[];

cp1 = [];

n1 = 0;

k = 1;

yy = [];

for j=1:length(m)

xx=[];     %%存贮当前cp下频率f的解

cp=m(j);

f1=0.1:0.01:20;

f=0.1;

y1=((cp/ct)^2-2)^2*tanh((d*pi*f/cp)*sqrt(1-(cp/cl)^2));

y2=4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*f/cp)*sqrt(1-(cp/ct)^2));

%plot(f,y1,f,y2)

b1=subs(y1-y2);

y11=[y1];

y22=[y2];

for i=2:length(f1)

f=f1(i);

y1=((cp/ct)^2-2)^2*tanh((d*pi*f/cp)*sqrt(1-(cp/cl)^2));

y2=4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*f/cp)*sqrt(1-(cp/ct)^2));

y11=[y11;y1];

y22=[y22;y2];

b2=subs(y1-y2);

if  b2==0

xx=[xx;f];

else if b1*b2<0     %在b1*b2<0是可以改进，取(f1+f2)/2

%  xx=[xx;(f1(i)+f1(i-1))/2];

a = f1(i-1);b = f1(i);

while ((b-a)>eps)

c = (b+a)/2;

fa = ((cp/ct)^2-2)^2*tanh((d*pi*a/cp)*sqrt(1-(cp/cl)^2))-4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*a/cp)*sqrt(1-(cp/ct)^2));

fb = ((cp/ct)^2-2)^2*tanh((d*pi*b/cp)*sqrt(1-(cp/cl)^2))-4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*b/cp)*sqrt(1-(cp/ct)^2));

fc = ((cp/ct)^2-2)^2*tanh((d*pi*c/cp)*sqrt(1-(cp/cl)^2))-4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*c/cp)*sqrt(1-(cp/ct)^2));

if fa*fc < 0

b=c;

else

a=c;

end

end

xx = [xx;c];

end

end

b1=b2;

end

n=length(xx);

%if n>n1

%   n1=n;

%end

%yy=cp*ones(n,1);

%plot(xx,yy,'*b')  %%%最后结合的时候只需要这一个图形即可

%hold on

%end

if n>0

A(k,1:length(xx))=xx;    %%%A表示群速度的频率矩阵，其中的每一行都代表某一群速度下的各个频率的排列

yy =[yy;cp];

k = k + 1 ;

n1=n;

if n>n1

n1 = n;

end

end

end

for i=1:n1

plot(A(:,i),yy)

hold on

end

%%%%%%%%%%%第二段

m1 = 0;

m2 = 0;

n1=0;                     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m=3.2:0.1:6.3;

cp1 = [];

A = [];

l =1;

yy = [];

for k=1:length(m)

cp=m(k);

xx=[];  %xx1=[];

%  m2=0;

ff=[0];

%(d*pi*f/cp)*sqrt((cp/ct)^2-1)=(2*k+1)*pi/2   %%%其中y相当于频率f，x相当于2k+1中的k

%将上式换成y=……x的关系f=(2*k+1)*((cp/ct)^2-1)^-0.5*cp/(2*d);

b=0:1:20;

for i=1:length(b)

k1=b(i);

fy=(2*k1+1)*((cp/ct)^2-1)^-0.5*cp/(2*d);

if fy<=20

ff=[ff;fy];          %%%断点

else

i=length(b);

end

end

%%%那几个断点f1是不能取的，取得话，会出现tan无穷的情况

for i=1:length(ff)

if i~=length(ff)

f1=ff(i)+0.01:0.01:ff(i+1)-0.01;

f=f1(1);

y1=((cp/ct)^2-2)^2*tanh((d*pi*f/cp)*sqrt(1-(cp/cl)^2));

y2=-4*sqrt((cp/ct)^2-1)*sqrt((1-cp/cl)^2)*tan((d*pi*f/cp)*sqrt((cp/ct)^2-1));

b1=subs(y1-y2);                                 %%%b1=0,没有结果的原因之一！！！

if b1 == 0

xx=[xx;f];

end

for j=2:length(f1)

f=f1(j);

y1=((cp/ct)^2-2)^2*tanh((d*pi*f/cp)*sqrt(1-(cp/cl)^2));

y2=-4*sqrt((cp/ct)^2-1)*sqrt((1-cp/cl)^2)*tan((d*pi*f/cp)*sqrt((cp/ct)^2-1));

b2=subs(y1-y2);

if  b2==0

xx=[xx;f];

else if b1*b2<0     %在b1*b2<0是可以改进，取(f1+f2)/2

% xx=[xx;(f1(j)+f1(j-1))/2];

a = f1(j-1);b = f1(j);

while ((b-a)>eps)

c = (b+a)/2;

fa = ((cp/ct)^2-2)^2*tanh((d*pi*a/cp)*sqrt(1-(cp/cl)^2))-4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*a/cp)*sqrt(1-(cp/ct)^2));

fb = ((cp/ct)^2-2)^2*tanh((d*pi*b/cp)*sqrt(1-(cp/cl)^2))-4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*b/cp)*sqrt(1-(cp/ct)^2));

fc = ((cp/ct)^2-2)^2*tanh((d*pi*c/cp)*sqrt(1-(cp/cl)^2))-4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*c/cp)*sqrt(1-(cp/ct)^2));

if fa*fc < 0

b=c;

else

a=c;

end

end

xx = [xx;c];

end

end

b1=b2;

end

end

end

n=length(xx);

%   if n>n1

%      n1=n;

%   end

%   yy=cp*ones(n,1);

% plot(xx,yy,'*b')  %%%最后结合的时候只需要这一个图形即可

%  hold on

%end

if n>0

A(l,1:length(xx))=xx;    %%%A表示群速度的频率矩阵，其中的每一行都代表某一群速度下的各个频率的排列

yy =[yy;cp];

l = l + 1 ;

if n>n1

n1 = n;

end

end

end

for i=1:n1             %%%%%%for 循环中的初始值，步长，以及终止值是在第一次循环的时候确定的，不会在运行的过程在改变，成为变量一样的存在,所以这里用while代替

yy1 = yy;

xx1 = A(:,i);

n2 = size(xx1);

j=1;

while (j <= n2)     %%%是为了去掉列向量中 的0元素

if xx1(j) == 0

yy1(j) = [];

xx1(j) = [];

n2 = size(xx1);

j=0;

end

j=j+1;

end

plot(xx1,yy1)

hold on

end

%%%%下面的应该是不变的，只要给出A矩阵就可以了

%%%%%%%%%%%%%%%%%%%%第三段

A=[];

ff=[0];

m=6.4:0.1:10;

n1=0;

cp1 = [];

l = 1;

yy = [];

for k=1:length(m)

cp=m(k);

xx=[];

b=0:1:20;

for i = 1:length(b)

k1 = b(i);

fx=((2*k1+1)*cp*((cp/cl)^2-1)^-0.5)/(d*2);

fy=((2*k1+1)*cp*((cp/ct)^2-1)^-0.5)/(d*2);

if fx<=14 && fy<=14

ff=[ff;fx];

ff=[ff;fy];

else if fx<=14 && fy>14

ff=[ff;fx];

else if fx>14 && fy<=14

ff=[ff;fy];

else if fx>14 && fy >14

i=length(b);

end

end

end

end

end         %%%求出各个tan为(2k+1)pi的点

ff=sort(ff);

for i=1:length(ff)-1

if ff(i+1)-ff(i)>=0.1       %%%%为了防止在上面f1出现空区间；本身f1(i)与f(i+1)相差很小，导致出现空区间

f1=ff(i)+0.01:0.01:ff(i+1)-0.01;

f=f1(1);

y1=((cp/ct)^2-2)^2*tan((d*pi*f/cp)*sqrt((cp/cl)^2-1));

y2=-4*sqrt((cp/ct)^2-1)*sqrt((cp/cl)^2-1)*tan((d*pi*f/cp)*sqrt((cp/ct)^2-1));

b1=subs(y1-y2);                                 %%%b1=0,没有结果的原因之一！！！

if b1 == 0

xx=[xx;b1];

end

for j=2:length(f1)

f=f1(j);

y1=((cp/ct)^2-2)^2*tan((d*pi*f/cp)*sqrt((cp/cl)^2-1));

y2=-4*sqrt((cp/ct)^2-1)*sqrt((cp/cl)^2-1)*tan((d*pi*f/cp)*sqrt((cp/ct)^2-1));

b2=subs(y1-y2);

if  b2==0

xx=[xx;f];

else if b1*b2<0     %在b1*b2<0是可以改进，取(f1+f2)/2

% xx=[xx;(f1(j)+f1(j-1))/2];

a = f1(j-1);b = f1(j);

while ((b-a)>eps)

c = (b+a)/2;

fa = ((cp/ct)^2-2)^2*tanh((d*pi*a/cp)*sqrt(1-(cp/cl)^2))-4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*a/cp)*sqrt(1-(cp/ct)^2));

fb = ((cp/ct)^2-2)^2*tanh((d*pi*b/cp)*sqrt(1-(cp/cl)^2))-4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*b/cp)*sqrt(1-(cp/ct)^2));

fc = ((cp/ct)^2-2)^2*tanh((d*pi*c/cp)*sqrt(1-(cp/cl)^2))-4*sqrt(1-(cp/ct)^2)*sqrt(1-(cp/cl)^2)*tanh((d*pi*c/cp)*sqrt(1-(cp/ct)^2));

if fa*fc < 0

b=c;

else

a=c;

end

end

xx = [xx;c];

end

end

b1=b2;

end

% plot(f,y11,f,y22)

%hold on

end

end

n=length(xx);

%   if n>n1

%       n1=n;

%  end

%   yy=cp*ones(n,1);

%   plot(xx,yy,'*b')  %%%最后结合的时候只需要这一个图形即可

%   hold on

%end

if n>0

A(l,1:length(xx))=xx;    %%%A表示群速度的频率矩阵，其中的每一行都代表某一群速度下的各个频率的排列

yy =[yy;cp];

l = l + 1 ;

if n>n1

n1 = n;

end

end

end

%%%%%%%%%%%%对A矩阵进行重新排列，排列成每个模式的点都在每一列上

b=zeros(size(A,1),1);

B = [A b];

for j = 1:n1+1

i1 = [];                 %%%%%

i2 = [];                  %%%% 记录的是本应为0的横坐标

for i = 1:size(B,1)-1

if  B(i,j) ~= 0 && B(i+1,j) ~= 0       %%%其中有00的情况

if  abs(B(i,j)-B(i+1,j))>0.7

for j1 = n1+1:-1:j+1

B(i+1,j1) = B(i+1,j1-1);

end

B(i+1,j) = B(i,j);

i2 = [i2;i+1];

end

else if B(i,j) ~=0 && B(i+1,j) == 0

B(i+1,j) = B(i,j);

i1 = [i1;i+1];

end

end

end

for n = 1:length(i1)

B(i1(n),j) = 0;

end

for n = 1:length(i2)

B(i2(n),j) = 0;

end

end

for i=1:n1             %%%%%%for 循环中的初始值，步长，以及终止值是在第一次循环的时候确定的，不会在运行的过程在改变，成为变量一样的存在,所以这里用while代替

yy1 = yy;

xx1 = B(:,i);

n2 = size(xx1,1);

j=1;

while (j <= n2)     %%%是为了去掉列向量中 的0元素

if xx1(j) == 0

yy1(j) = [];

xx1(j) = [];

n2 = size(xx1,1);

j=0;

end

j=j+1;

end

plot(xx1,yy1)

hold on

end