%雙脊金屬加載矩形波導(dǎo)，的基模和第一個(gè)高階模的色散曲線計(jì)算。
%截止波數(shù)與步長的關(guān)系
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矩陣包含邊界點(diǎn)
    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               %金屬脊內(nèi)部場為零，k矩陣對(duì)應(yīng)行全賦零
        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                   %金屬脊的下上兩個(gè)邊界得邊界條件
        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
%---對(duì)k矩陣得處理
    [fai,v]=eig(ke);
    [m,n]=find(v<0);                    %找出所有負(fù)特征值的坐標(biāo)
    l=length(m);
    for i=1:l
        fu2(i)=v(m(i),n(i));             %將所有負(fù)特征值放入數(shù)組fu中。
    end
    beta2=abs(fu2);
    kc2=sqrt(beta2)*1000/h;
    [kc2,index2]=sort(kc2);                 %將fu按絕對(duì)值從小到大排序，
    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的關(guān)系')
xlabel('h(mm)')
ylabel('kc')
figure(2)
plot(hh,k3);
hold on
plot(hh,k4);
%%----------------kc的相對(duì)誤差--------------%%
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('相對(duì)誤差');
title('基模和第一個(gè)高階模的截止波數(shù)相對(duì)誤差')