Matlab实现双脊波导基模色散特性曲线计算

%双脊金属加载矩形波导，的基模和第一个高阶模的色散曲线计算%
clear 
close all
clc
Mu=1/(36*pi*1e9);  %真空中介电常数
Epsilon=4*pi*1e-7;  %真空中磁导率
GFc=zeros(7,1);
a=10e-3;                    %波导的长度
b=5e-3;                    %波导宽
C=2e-3;                        %金属脊高
d=1.5e-3;                     %金属脊位置
h=0.3e-3                      %网格步长 
D=3e-3; 
c=3e8;
xmax=fix(a/h+1);
ymax=fix(b/h+1);
kx=xmax-2;
ky=ymax-2;
x1min=fix((a/2-d-D/2)/h+1);  %第一个金属脊的位置
x1max=x1min+fix(d/h);
y=fix(C/h+1);
x2min=fix(a/2/h+D/2/h+1);  %第二个金属脊的位置
x2max=x2min+fix(d/h);
%-----------------------------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-1:x1max-1                   %金属脊处场为零，k矩阵对应行全赋零
    for j=1:y
        km(i+(j-1)*kx,:)=0;
    end
end
for i=x2min-1:x2max-1                   %金属脊处场为零，k矩阵对应行全赋零
    for j=1:y
        km(i+(j-1)*kx,:)=0;
    end
end
%---对k矩阵得处理
[fai,v]=eig(km);
[m,n]=find(v<0);                    %找出所有负特征值的坐标
l=length(m);
for i=1:l
    fu(i)=v(m(i),n(i));             %将所有负特征值放入数组fu中。
end
beta1=abs(fu);
kc1=sqrt(beta1)/h;