clear all
%close all
clc
tic
wnx=600*2*pi;
kx=10e6;
zetax=0.04;
wny=600*2*pi;
ky=10e6;
zetay=0.04;
Kt=2.2e9;
Kr=0.5e9;
Kn=0.227;
Nt=4;
phis=0*pi/180;
phie=180*pi/180;
a_xx=0.5*((cos(2*phie)-2*Kn*phie+Kn*sin(2*phie))-(cos(2*phis)-2*Kn*phis+Kn*sin(2*phis)));
a_xy=0.5*((-sin(2*phie)-2*phie+Kn*cos(2*phie))-(-sin(2*phis)-2*phis+Kn*cos(2*phis)));
a_yx=0.5*((-sin(2*phie)+2*phie+Kn*cos(2*phie))-(-sin(2*phis)+2*phis+Kn*cos(2*phis)));
a_yy=0.5*((-cos(2*phie)-2*Kn*phie-Kn*sin(2*phie))-(-cos(2*phis)-2*Kn*phis-Kn*sin(2*phis)));
wnmax = max([wnx wny]);
w=(0:1:2*wnmax/2/pi)*2*pi;
G11=(wnx^2/kx)./(wnx^2 - w.^2 + i*2*zetax*wnx.*w);
G22=(wny^2/ky)./(wny^2 - w.^2 + i*2*zetay*wny.*w);
for cnt=1:length(w)
FRF_or = [a_xx*G11(cnt) a_xy*G22(cnt);a_yx*G11(cnt) a_yy.*G22(cnt)];
E = eig(FRF_or);
temp=E(1);
lambda1(cnt)=temp;
temp=E(2);
lambda2(cnt)=temp;
if (cnt>1)
dot_prod1=real(lambda2(cnt)).*real(lambda2(cnt-1))+ imag(lambda2(cnt)).*imag(lambda2(cnt-1));
dot_prod2=real(lambda2(cnt)).*real(lambda1(cnt-1)) + imag(lambda2(cnt)).*imag(lambda1(cnt-1));
if ( dot_prod2 > dot_prod1)
temp = lambda2(cnt);
lambda2(cnt) =lambda1(cnt);
lambda1(cnt) = temp;
end
end
end
lambda1 = lambda1';
bliml = (2*pi/Nt/Kt)./((real(lambda1)).^2 + (imag(lambda1)).^2).*(real(lambda1).*(1+(imag(lambda1)./real(lambda1)).^2));
[index1]=find(bliml >= 0);
bliml = bliml(index1);
bliml = bliml*1e6;
w1 = w(index1);
psi1 = atan2(imag(lambda1),real(lambda1));
psi1 = psi1(index1);
epsilonl = pi - 2*psi1;
lobes = 8;
w=w1;
w=w';
Omega=zeros(lobes,length(w));
for cnt=1:lobes
Omega(cnt,:)=(60/Nt)*w./((cnt-1)*2*pi + epsilonl);
end
x = cell(lobes,1);
y = cell(lobes,1);
z = cell(lobes,1);
CZ = max(Omega);
for cnt = 1:lobes
% x{cnt}=ceil(Omega(cnt,1)):1:floor(Omega(cnt,length(Omega(cnt,:)));
x{cnt}=1:10:CZ;
x{cnt}=x{cnt}';
y{cnt}=interp1 (Omega(cnt,:),bliml,x{cnt,:});
z{cnt}=interp1 (Omega(cnt,:),w,x{cnt,:});
end
A=[];
for cnt=1:lobes
A=cat(1,A,x{cnt});
end
B=unique(A);
C=zeros(length(B),3);
C(:,1)=B;
C(:,2)=1/eps;
C(:,3)=1/eps;
for cnt = lobes:-1:1
for i=1:length(x{cnt})
Nt = find(C(:,1) == x{cnt}(i,1));
if y{cnt}(i,1)<C(Nt,2)
C(Nt,2)=y{cnt}(i,1);
C(Nt,3)=z{cnt}(i,1);
end
end
end
plot(C(:,1),C(:,2),'b-')
hold on
axis([1000 6000 0 500])
set(gca,'FontSize',10)
legend('C1(:,2)')
hold on
wnx1=4035*2*pi;
kx1=2.1425e6;
zetax1=0.02;
wny1=4035*2*pi;
ky1=2.1425e6;
zetay1=0.02;
Kt1=1.482e9;
Kr1=0.889e9;
Kn1=1550/2550;
Nt1=4;
phis1=0*pi/180;
phie1=180*pi/180;
a_xx1=0.5*((cos(2*phie1)-2*Kn1*phie1+Kn1*sin(2*phie1))-(cos(2*phis1)-2*Kn1*phis1+Kn1*sin(2*phis1)));
a_xy1=0.5*((-sin(2*phie1)-2*phie1+Kn1*cos(2*phie1))-(-sin(2*phis1)-2*phis1+Kn1*cos(2*phis1)));
a_yx1=0.5*((-sin(2*phie1)+2*phie1+Kn1*cos(2*phie1))-(-sin(2*phis1)+2*phis1+Kn1*cos(2*phis1)));
a_yy1=0.5*((-cos(2*phie1)-2*Kn1*phie1-Kn1*sin(2*phie1))-(-cos(2*phis1)-2*Kn1*phis1-Kn1*sin(2*phis1)));
wnmax1 = max([wnx1 wny1]);
w2=(0:1:2*wnmax1/2/pi)*2*pi;
G111=(wnx1^2/kx1)./(wnx1^2 - w2.^2 + i*2*zetax1*wnx1.*w2);
G221=(wny1^2/ky1)./(wny1^2 - w2.^2 + i*2*zetay1*wny1.*w2);
for t=1:length(w2)
FRF_or1 = [a_xx1*G111(t) a_xy1*G221(t);a_yx1*G111(t) a_yy1.*G221(t)];
Ee = eig(FRF_or1);
temp1=Ee(1);
lambda11(t)=temp1;
temp1=Ee(2);
lambda21(t)=temp1;
if (t>1)
dot_prod1=real(lambda21(t)).*real(lambda21(t-1))+ imag(lambda21(t)).*imag(lambda21(t-1));
dot_prod2=real(lambda21(t)).*real(lambda11(t-1)) + imag(lambda21(t)).*imag(lambda11(t-1));
if ( dot_prod2 > dot_prod1)
temp1 = lambda21(t);
lambda21(t) =lambda11(t);
lambda11(t) = temp1;
end
end
end
lambda11 = lambda11';
bliml1 = (2*pi/Nt1/Kt1)./((real(lambda11)).^2 + (imag(lambda11)).^2).*(real(lambda11).*(1+(imag(lambda11)./real(lambda11)).^2));
[index11]=find(bliml1 >= 0);
bliml1 = bliml1(index11);
bliml1 = bliml1*1e6;
w3 = w2(index11);
psi11 = atan2(imag(lambda11),real(lambda11));
psi11 = psi11(index11);
epsilonl1 = pi - 2*psi11;
lobess = 6;
w2=w3;
w2=w2';
Omega1=zeros(lobess,length(w2));
for t=1:lobess
Omega1(t,:)=(60/Nt1)*w2./((t-1)*2*pi + epsilonl1);
end
x1 = cell(lobess,1);
y1 = cell(lobess,1);
z1 = cell(lobess,1);
CZ1 = max(Omega1);
for t = 1:lobess
% x{cnt}=ceil(Omega(cnt,1)):1:floor(Omega(cnt,length(Omega(cnt,:)));
x1{t}=1:10:CZ1;
x1{t}=x1{t}';
y1{t}=interp1 (Omega1(t,:),bliml1,x1{t,:});
z1{t}=interp1 (Omega1(t,:),w2,x1{t,:});
end
A1=[];
for t=1:lobess
A1=cat(1,A1,x1{t});
end
B1=unique(A1);
C1=zeros(length(B1),3);
C1(:,1)=B1;
C1(:,2)=1/eps;
C1(:,3)=1/eps;
for t = lobess:-1:1
for i=1:length(x1{t})
Nt1 = find(C1(:,1) == x1{t}(i,1));
if y1{t}(i,1)<C1(Nt1,2)
C1(Nt1,2)=y1{t}(i,1);
C1(Nt1,3)=z1{t}(i,1);
end
end
end
plot(C1(:,1),C1(:,2),'g-')
hold on
axis([10000 90000 0 500])
set(gca,'FontSize',10)
hold on
wnx2=400*2*pi;
kx2=2e6;
zetax2=0.04;
wny2=400*2*pi;
ky2=2e6;
zetay2=0.04;
Kt2=1.482e5;
Kr2=1.68e5;
Kn2=1.5e5;
Nt2=4;
phis2=0*pi/180;
phie2=180*pi/180;
a_xx2=0.5*((cos(2*phie2)-2*Kn2*phie2+Kn2*sin(2*phie2))-(cos(2*phis2)-2*Kn2*phis2+Kn2*sin(2*phis2)));
a_xy2=0.5*((-sin(2*phie2)-2*phie2+Kn2*cos(2*phie2))-(-sin(2*phis2)-2*phis2+Kn2*cos(2*phis2)));
a_yx2=0.5*((-sin(2*phie2)+2*phie2+Kn2*cos(2*phie2))-(-sin(2*phis2)+2*phis2+Kn2*cos(2*phis2)));
a_yy2=0.5*((-cos(2*phie2)-2*Kn2*phie2-Kn2*sin(2*phie2))-(-cos(2*phis2)-2*Kn2*phis2-Kn2*sin(2*phis2)));
wnmax2 = max([wnx2 wny2]);
w4=(0:1:2*wnmax2/2/pi)*2*pi;
G112=(wnx2^2/kx2)./(wnx2^2 - w4.^2 + i*2*zetax2*wnx2.*w4);
G222=(wny2^2/ky2)./(wny2^2 - w4.^2 + i*2*zetay2*wny2.*w4);
for cnt2=1:length(w4);
FRF_or = [a_xx2*G112(cnt2) a_xy2*G222(cnt2);a_yx2*G112(cnt2) a_yy2.*G222(cnt2)];
Ee1 = eig(FRF_or);
temp2=Ee1(1);
lambda12(cnt2)=temp2;
temp2=Ee1(2);
lambda22(cnt2)=temp2;
if (cnt2>1)
dot_prod12=real(lambda22(cnt2)).*real(lambda22(cnt2-1))+ imag(lambda22(cnt2)).*imag(lambda22(cnt2-1));
dot_prod22=real(lambda22(cnt2)).*real(lambda12(cnt2-1)) + imag(lambda22(cnt2)).*imag(lambda12(cnt2-1));
if ( dot_prod22 > dot_prod12)
temp2 = lambda22(cnt2);
lambda22(cnt2) =lambda12(cnt2);
lambda12(cnt2) = temp2;
end
end
end
lambda12 = lambda12';
bliml2 = (2*pi/Nt2/Kt2)./((real(lambda12)).^2 + (imag(lambda12)).^2).*(real(lambda12).*(1+(imag(lambda12)./real(lambda12)).^2));
[index12]=find(bliml2 >= 0);
bliml2 = bliml2(index12);
bliml2 = bliml2*1e6;
w11 = w4(index12);
psi12 = atan2(imag(lambda12),real(lambda12));
psi12 = psi12(index12);
epsilonl2 = pi - 2*psi12;
lobes2 = 6;
w4=w11;
w4=w4';
Omega2=zeros(lobes2,length(w4));
for cnt2=1:lobes2
Omega2(cnt2,:)=(60/Nt2)*w4./((cnt2-1)*2*pi + epsilonl2);
end
x2 = cell(lobes2,1);
y2 = cell(lobes2,1);
z2 = cell(lobes2,1);
CZ2 = max(Omega2);
for cnt2 = 1:lobes2
% x{cnt}=ceil(Omega(cnt,1)):1:floor(Omega(cnt,length(Omega(cnt,:)));
x2{cnt2}=1:10:CZ2;
x2{cnt2}=x2{cnt2}';
y2{cnt2}=interp1 (Omega2(cnt2,:),bliml2,x2{cnt2,:});
z2{cnt2}=interp1 (Omega2(cnt2,:),w4,x2{cnt2,:});
end
A2=[];
for cnt2=1:lobes2
A2=cat(1,A2,x2{cnt2});
end
B2=unique(A2);
C2=zeros(length(B2),3);
C2(:,1)=B2;
C2(:,2)=1/eps;
C2(:,3)=1/eps;
for cnt2 = lobes2:-1:1
for i=1:length(x2{cnt2})
Nt2 = find(C2(:,1) == x2{cnt2}(i,1));
if y2{cnt2}(i,1)<C2(Nt2,2)
C2(Nt2,2)=y2{cnt2}(i,1);
C2(Nt2,3)=z2{cnt2}(i,1);
end
end
end
plot(C2(:,1),C2(:,2),'b-')
hold on
axis([1
- 1
- 2
- 3
- 4
- 5
前往页