tm.m

TDFD,双脊金属加载矩形波导

%-------------------TM波的基模与第一高阶模------------------------%
clear
a=10.16;
b=5.588;
c=1;
d=5.08;
h=0.3;
xmax=fix(a/h+1);
ymax=fix(b/h+1);
kx=xmax-2;
ky=ymax-2;
x1min=fix((a/2-d/2)/h);            %金属脊的位置
x1max=kx-x1min+1;
y1min=fix(c/h);
y1max=ky-y1min+1
%--------------------TM波----------------%
kn=kx*ky;
km=zeros(kn);                       %K矩阵
for i=0:ky-1
    for j=1:kx                      % D矩阵
        km(j+i*kx,j+i*kx)=-4;
        if j<kx
            km(j+i*kx,j+i*kx+1)=1;
            km(j+i*kx+1,j+i*kx)=1;
        end
    end
    if i<ky-1                       %D矩阵两旁的I矩阵
        for j=1:kx
            km(j+i*kx,j+(i+1)*kx)=1;
            km(j+(i+1)*kx,j+i*kx)=1;
        end
    end
end
for i=x1min:x1max                   %金属脊处场为零，k矩阵对应行全赋零
    for j=1:y1min
        km(i+(j-1)*kx,:)=0;
        km(i+(ky-j)*kx,:)=0;
    end
end
%---对k矩阵得处理
[fai,v]=eig(km);
[m,n]=find(v<0);                    %找出所有负特征值的坐标
l=length(m);
% fu(l)=0;
for i=1:l
    fu(i)=v(m(i),n(i));             %将所有负特征值放入数组fu中。
end
beta1=abs(fu);
kc1=sqrt(beta1)*1000/h;
[kc1,index1]=sort(kc1);              %将kc按绝对值从小到大排序，
k1=kc1(1);
hz1=fai(:,m(index1(1)));             %取出前两个特征值，及对应得特征向量
k2=kc1(2);
hz2=fai(:,m(index1(2)));
beta=0:0.1:100;                    %最前两个为绝对值最小得特征值。
k11=sqrt(k1^2+beta.^2);
k22=sqrt(k2^2+beta.^2);
plot(k11,beta,'b');
figure(2)
plot(k22,beta,'r--');

fai1=hz1;
fai2=hz2;
m=length(fai1);                  %基模的特征值
if m>(xmax-2)*(ymax-2)
    nx=xmax;
    ny=ymax;
else
    nx=xmax-2;
    ny=ymax-2;
end
fai11(ny,nx)=0;
for i=1:ny
    for j=1:nx
        fai11(i,j)=fai1(j+(i-1)*(nx));
    end
end
figure(3)
subplot(2,2,1)
[px,py] = gradient(fai11);
contour(fai11);
hold on;
quiver(px,py);
subplot(2,2,2)
mesh(fai11);


m=length(fai2);                  %第一高阶模的特征值
if m>(xmax-2)*(ymax-2)
    nx=xmax;
    ny=ymax;
else
    nx=xmax-2;
    ny=ymax-2;
end
fai22(ny,nx)=0;
for i=1:ny
    for j=1:nx
        fai22(i,j)=fai2(j+(i-1)*(nx));
    end
end
hold on;
subplot(2,2,3)
[px,py] = gradient(fai22);
contour(fai22);
hold on;
quiver(px,py);
subplot(2,2,4)
mesh(fai22);