fdtd.m

FDTD（时域有限差分法）写的一个双脊金属加载矩形波导

%双脊金属加载矩形波导，的基模和第一个高阶模的色散曲线计算。
% function fdtd()
clear
a=10.16;                    %波导的长度
b=5.588;                    %波导宽
c=1;                        %金属脊高
d=5.08;                     %金属脊位置
h=0.2;                     %网格步长
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)+1;
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);
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)));

%-----------------------------------TE波--------------------------------%
kx=xmax;                            %K矩阵包含边界点
ky=ymax;
x1min=fix((a/2-d/2)/h)+1;            %金属脊的位置
x1max=xmax-x1min+1;
y1min=fix(c/h)+1;
y1max=ymax-y1min+1;
kn=xmax*ymax;
ke=zeros(kn);                       %K矩阵
for i=0:ky-1
    for j=1:kx                      %D矩阵
        ke(j+i*kx,j+i*kx)=-4;
        if j<kx
            ke(j+i*kx,j+i*kx+1)=1;
            ke(j+i*kx+1,j+i*kx)=1;
        end
    end
    if i<ky-1                       %D矩阵两旁的I矩阵
        for j=1:kx
            ke(j+i*kx,j+(i+1)*kx)=1;
            ke(j+(i+1)*kx,j+i*kx)=1;
        end
    end
end
for i=1:kx                          %传输TE波的边界条件
    ke(i,i+kx)=2;                   %下边界
    ke((ky-1)*kx+i,(ky-2)*kx+i)=2;  %上边界
end
for j=1:ky
    ke((j-1)*kx+1,(j-1)*kx+2)=2;    %左边界
    ke(j*kx,j*kx-1)=2;              %右边界
end

for i=x1min+1:x1max-1               %金属脊内部场为零，k矩阵对应行全赋零
    for j=1:y1min-1
        ke(i+(j-1)*kx,:)=0;
    end
    for j=y1max+1:ky
        ke(i+(j-1)*kx,:)=0;
    end
end
for i=x1min+1:x1max-1                   %金属脊的下上两个边界得边界条件
    ke(kx*(y1min-1)+i,kx*(y1min-2)+i)=0;
    ke((y1min-1)*kx+i,y1min*kx+i)=2;
    ke(kx*(y1max-1)+i,kx*y1max+i)=0;
    ke((y1max-1)*kx+i,(y1max-2)*kx+i)=2;
end
for i=1:y1min-1                      %下金属脊的左右两边界
    ke((i-1)*kx+x1min,(i-1)*kx+x1min+1)=0;  
    ke((i-1)*kx+x1min,(i-1)*kx+x1min-1)=2;
    ke((i-1)*kx+x1max,(i-1)*kx+x1max-1)=0;  
    ke((i-1)*kx+x1max,(i-1)*kx+x1max+1)=2;
end
for i=y1max+1:ky                      %上金属脊的左右两边界
    ke((i-1)*kx+x1min,(i-1)*kx+x1min+1)=0;  
    ke((i-1)*kx+x1min,(i-1)*kx+x1min-1)=2;
    ke((i-1)*kx+x1max,(i-1)*kx+x1max-1)=0;  
    ke((i-1)*kx+x1max,(i-1)*kx+x1max+1)=2;
end
%---对k矩阵得处理
[fai,v]=eig(ke);
[m,n]=find(v<0);                    %找出所有负特征值的坐标
l=length(m);
for i=1:l
    fu2(i)=v(m(i),n(i));             %将所有负特征值放入数组fu中。
end
beta2=abs(fu2);
kc2=sqrt(beta2)*1000/h;
[kc2,index2]=sort(kc2);                 %将fu按绝对值从小到大排序，
if kc2(1)>1e-2
    k3=kc2(1);
    ez1=fai(:,m(index2(1)));                %最小TE模的特征向量
    k4=kc2(2);
    ez2=fai(:,m(index2(2)));                %次小TE模的特征向量
else
    k3=kc2(2);
    ez1=fai(:,m(index2(2)));                %最小TE模的特征向量
    k4=kc2(3);
    ez2=fai(:,m(index2(3)));                %次小TE模的特征向量
end

kc=sort([k1 k2 k3 k4]);
%%-------------------画图-----%
beta=0:0.1:100;                    %最前两个为绝对值最小得特征值。
k11=sqrt(kc(1)^2+beta.^2);
k22=sqrt(kc(2)^2+beta.^2);
plot(k11,beta,'b');
figure(2)
plot(k22,beta,'r--');
hold on
ylabel('beta');
xlabel('k')
if k1>kc(2)                          %分辨基模与第一高阶模的类型
    kc11=k3;
    fai1=ez1;
    kc22=k4;
    fai2=ez2;
    text(0.2,0.5,'蓝色-基模,为TE模')
    text(0.5,0.2,'红色-第一高次模，为TE模')
else
    if k1>kc(1)
        kc11=k3;
        fai1=ez1;
        kc22=k1;
        fai2=hz1;
        text(0.2,0.5,'蓝色-基模,为TE模')
        text(0.5,0.2,'红色-第一高次模，为TM模') 
    else
        if k2>k3
            kc11=k1;
            fai1=hz1;
            kc22=k3;
            fai2=ez1;
            text(0.2,0.5,'蓝色-基模,为TM模')
            text(0.5,0.2,'红色-第一高次模，为TE模')
        else
            kc11=k1;
            fai1=hz1
            kc22=k2;
            fai2=hz2;
            text(0.2,0.5,'蓝色-基模,为TM模')
            text(0.5,0.2,'红色-第一高次模，为TM模') 
        end
    end
end
%-----------------------画出横界面模式场分布------------------%

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);
xlabel('基模的横截面场分布（此处为磁场）')
% 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
subplot(2,2,2)
[px,py] = gradient(fai22);
contour(fai22);
hold on;
quiver(px,py);
xlabel('第一高阶模的横截面场分布（此处为磁场）')
% mesh(fai22);

%%-------------------------TM模的两个最小模式场分布----------%%

fai3=hz1;
fai4=hz2;
m=length(fai3);                  %TM模的基模的特征值
if m>(xmax-2)*(ymax-2)
    nx=xmax;
    ny=ymax;
else
    nx=xmax-2;
    ny=ymax-2;
end
fai33(ny,nx)=0;
for i=1:ny
    for j=1:nx
        fai33(i,j)=fai3(j+(i-1)*(nx));
    end
end
subplot(2,2,3)
[px,py] = gradient(fai33);
contour(fai33);
hold on;
quiver(px,py);
xlabel('TM波基模的横截面场分布（此处为电场）')
% mesh(fai33);
m=length(fai4);                  %TM模的第一高阶模的特征值
if m>(xmax-2)*(ymax-2)
    nx=xmax;
    ny=ymax;
else
    nx=xmax-2;
    ny=ymax-2;
end
fai44(ny,nx)=0;
for i=1:ny
    for j=1:nx
        fai44(i,j)=fai4(j+(i-1)*(nx));
    end
end
hold on;
subplot(2,2,4)
[px,py] = gradient(fai44);
contour(fai44);
hold on;
quiver(px,py);
xlabel('TM波第一高阶模的横截面场分布（此处为电场）')
% mesh(fai44);