kc_h.m

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

%双脊金属加载矩形波导，的基模和第一个高阶模的色散曲线计算。
%截止波数与步长的关系
clear
a=10.16;
b=5.588;
c=1;
d=5.08;
for p=1:11
    h=0.2+(p-1)*0.02;
    xmax=fix(a/h+1);
    ymax=fix(b/h+1);
    x1min=fix((a/2-d/2)/h);            %金属脊的位置
    x1max=xmax-x1min-1;
    y1min=fix(c/h)+1;
    y1max=ymax-y1min+1;
%-----------------------------------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;
%         ke(i+(ky-j)*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(p)=kc2(1);
        ez1=fai(:,m(index2(1)));                %最小TE模的特征向量
        k4(p)=kc2(2);
        ez2=fai(:,m(index2(2)));                %次小TE模的特征向量
    else
        k3(p)=kc2(2);
        ez1=fai(:,m(index2(2)));                %最小TE模的特征向量
        k4(p)=kc2(3);
        ez2=fai(:,m(index2(3)));                %次小TE模的特征向量
    end
end
hh=0.2:0.02:0.4;
figure(1)
subplot(1,2,1)
plot(hh,k3);
hold on
subplot(1,2,2)
plot(hh,k4);
title('kc与h的关系')
xlabel('h(mm)')
ylabel('kc')
figure(2)
plot(hh,k3);
hold on
plot(hh,k4);
%%----------------kc的相对误差--------------%%
kc_te10=pi*1000/(a+2*c);
kc_te20=2*kc_te10;
k3=(k3-kc_te10)/kc_te10;
k4=(k4-kc_te20)/kc_te20;
figure(3)
plot(hh,k3,'g',hh,k4,'r--')
xlabel('h');
ylabel('相对误差');
title('基模和第一个高阶模的截止波数相对误差')