waveguide.m

计算无限大平板波导的截止频率的matlab源代码。具体geometry信息可以查看m源文件内的说明。

% ece597II HW#1
% Problem 2: Waveguide Cutoff Frequency
% by Wei Wang, ECE, UMass Amherst
% Sep 28, 2007

% geometry set up
% unit: mm
length = 22.86;
width = 10.16;
hx = 0.381;             % grid size
hy = 0.254;
nx = length/hx + 1;     % number of node along x-axis
ny = width/hy + 1;      % number of node along y-axis
node = nx*ny-2*(7.62/hx-1)*(2.54/hy);       % number of unknowns
xval = zeros(node,1);
yval = zeros(node,1);
A = sparse(node,node);  % create the FDM matrix A
%********************************************************

% nodes coordinate value initialization
for i=1:(7.62/hx+1)
    for j=1:ny
        xval(j+(i-1)*ny,1) = (i-1)*hx;
        yval(j+(i-1)*ny,1) = (j-1)*hy;
    end
end
xcurr = j+(i-1)*ny;
ycurr = xcurr;

for i=1:(7.62/hx-1)
    for j=1:(5.08/hy+1)
        xval(xcurr+j+(i-1)*(5.08/hy+1),1) = (i-1)*hx + 7.62 + hx;
        yval(xcurr+j+(i-1)*(5.08/hy+1),1) = (j-1)*hy + 2.54;
    end
end
xcurr = xcurr+j+(i-1)*(5.08/hy+1);
ycurr = xcurr;

for i=1:(7.62/hx+1)
    for j=1:ny
        xval(xcurr+j+(i-1)*ny,1) = (i-1)*hx + 15.24;
        yval(xcurr+j+(i-1)*ny,1) = (j-1)*hy;
    end
end
%******************************************************

t0 = clock; 

% using FDM to create matrix A
% for convenience, define some constants
const1 = -2*(hx^2+hy^2);
const2 = hx^2;
const3 = hy^2;

% for boundary: x=0, x=length.
A(1,1)=const1; A(1,1+ny)=2*const3; A(1,2)=2*const2;
for i=2:(ny-1)
    A(i,i)=const1; A(i,i+ny)=2*const3; A(i,i-1)=const2; A(i,i+1)=const2;
end
A(ny,ny)=const1; A(ny,ny*2)=2*const3; A(ny,ny-1)=2*const2;

A(node-ny+1,node-ny+1)=const1; A(node-ny+1,node-2*ny+1)=2*const3; A(node-ny+1,node-ny+2)=2*const2;
for i=(node-ny+2):(node-1)
    A(i,i)=const1; A(i,i-ny)=2*const3; A(i,i-1)=const2; A(i,i+1)=const2;
end
A(node,node)=const1; A(node,node-1)=2*const2; A(node,node-ny)=2*const3;

% for boundary: x=7.62, x=15.24
tmp1 = (length/3/hx)*ny;
tmp2 = 2.54/hy;
tmp3 = (7.62/hx+1)*ny;
A(tmp1+1,tmp1+1)=const1; A(tmp1+1,tmp1+1-ny)=2*const3; A(tmp1+1,tmp1+2)=2*const2;
for i=(tmp1+2):(tmp1+tmp2)
    A(i,i)=const1; A(i,i-ny)=2*const3; A(i,i-1)=const2; A(i,i+1)=const2;
    A((i+3*tmp2),(i+3*tmp2))=const1; A((i+3*tmp2),(i+3*tmp2-ny))=2*const3; A((i+3*tmp2),(i+3*tmp2-1))=const2;
    A((i+3*tmp2),(i+3*tmp2+1))=const2;
end
A(tmp3,tmp3)=const1; A(tmp3,tmp3-1)=2*const2; A(tmp3,tmp3-ny)=2*const3;

tmp4 = node-tmp3+1;
A(tmp4,tmp4)=const1; A(tmp4,tmp4+ny)=2*const3; A(tmp4,tmp4+1)=2*const2;
for i=(tmp4+1):(tmp4+tmp2-1)
    A(i,i)=const1; A(i,i+ny)=2*const3; A(i,i-1)=const2; A(i,i+1)=const2;
    A((i+3*tmp2),(i+3*tmp2))=const1; A((i+3*tmp2),(i+3*tmp2+ny))=2*const3; A((i+3*tmp2),(i+3*tmp2-1))=const2;
    A((i+3*tmp2),(i+3*tmp2+1))=const2;
end
A((tmp4+ny-1),(tmp4+ny-1))=const1; A((tmp4+ny-1),(tmp4+2*ny-1))=2*const3; A((tmp4+ny-1),(tmp4+ny-2))=2*const2;

tmp5 = tmp1+tmp2+1;
for i=tmp5:(tmp5+2*tmp2)
    A(i,i)=const1; A(i,i-1)=const2; A(i,i+1)=const2; A(i,i-ny)=const3; A(i,i+tmp3+1-tmp5)=const3;
end
tmp6 = node-tmp3+tmp2+1;
for i=tmp6:(tmp6+2*tmp2)
    A(i,i)=const1; A(i,i-1)=const2; A(i,i+1)=const2; A(i,i-(3*tmp2+1))=const3; A(i,i+ny)=const3;
end

% for area (7.62+hx)<=x<=(7.62*2-hx), 2.54<=y<=7.62
i=tmp3+1;
A(i,i)=const1; A(i,i+1)=2*const2; A(i,i-(tmp3-tmp1-tmp2))=const3; A(i,i+2*tmp2+1)=const3;
i=tmp3+1+2*tmp2;
A(i,i)=const1; A(i,i-1)=2*const2; A(i,i-(tmp3-tmp1-tmp2))=const3; A(i,i+2*tmp2+1)=const3;
for i=(tmp3+2):(tmp3+2*tmp2)
    A(i,i)=const1; A(i,i-1)=const2; A(i,i+1)=const2; A(i,i-(tmp3-tmp1-tmp2))=const3; A(i,i+2*tmp2+1)=const3;
end
i=node-tmp3-2*tmp2;
A(i,i)=const1; A(i,i+1)=2*const2; A(i,i-(2*tmp2+1))=const3; A(i,i+3*tmp2+1)=const3;
i=node-tmp3;
A(i,i)=const1; A(i,i-1)=2*const2; A(i,i-(2*tmp2+1))=const3; A(i,i+3*tmp2+1)=const3;
for i=(node-tmp3-2*tmp2+1):(node-tmp3-1)
    A(i,i)=const1; A(i,i-1)=const2; A(i,i+1)=const2; A(i,i-(2*tmp2+1))=const3; A(i,i+3*tmp2+1)=const3;
end

for i=(tmp3+2*tmp2+2):(2*tmp2+1):(node-tmp3-4*tmp2-1)
    A(i,i)=const1; A(i,i+1)=2*const2; A(i,i-(2*tmp2+1))=const3; A(i,i+(2*tmp2+1))=const3;
end
for i=(tmp3+2+4*tmp2):(2*tmp2+1):(node-tmp3-2*tmp2-1)
    A(i,i)=const1; A(i,i-1)=2*const2; A(i,i-(2*tmp2+1))=const3; A(i,i+(2*tmp2+1))=const3;
end

i=tmp3+2*tmp2+3;
for ii=1:(length/3/hx-3)
    for jj=1:(2*tmp2-1)
        A(i+(jj-1)+(ii-1)*(2*tmp2+1),i+(jj-1)+(ii-1)*(2*tmp2+1))=const1; 
        A(i+(jj-1)+(ii-1)*(2*tmp2+1),i+(jj-1)+(ii-1)*(2*tmp2+1)-2*tmp2-1)=const3;
        A(i+(jj-1)+(ii-1)*(2*tmp2+1),i+(jj-1)+(ii-1)*(2*tmp2+1)+2*tmp2+1)=const3; 
        A(i+(jj-1)+(ii-1)*(2*tmp2+1),i+(jj-1)+(ii-1)*(2*tmp2+1)-1)=const2;
        A(i+(jj-1)+(ii-1)*(2*tmp2+1),i+(jj-1)+(ii-1)*(2*tmp2+1)+1)=const2;
    end
end

% for area: hx<=x<=7.62-hx, 0<=y<=width
for i=(ny+1):ny:(tmp1+1-ny)
    A(i,i)=const1; A(i,i+1)=2*const2; A(i,i-ny)=const3; A(i,i+ny)=const3;
end
for i=(2*ny):ny:(tmp3-ny)
    A(i,i)=const1; A(i,i-1)=2*const2; A(i,i-ny)=const3; A(i,i+ny)=const3;
end
i=ny+2;
for ii=1:(length/3/hx-1)
    for jj=1:(ny-2)
        A(i+(jj-1)+(ii-1)*ny,i+(jj-1)+(ii-1)*ny)=const1; 
        A(i+(jj-1)+(ii-1)*ny,i+(jj-1)+(ii-1)*ny-ny)=const3;
        A(i+(jj-1)+(ii-1)*ny,i+(jj-1)+(ii-1)*ny+ny)=const3;
        A(i+(jj-1)+(ii-1)*ny,i+(jj-1)+(ii-1)*ny-1)=const2;
        A(i+(jj-1)+(ii-1)*ny,i+(jj-1)+(ii-1)*ny+1)=const2;
    end
end

% for area: 15.24+hx<=x<=length-hx, 0<=y<=width
for i=(node-tmp3+1+ny):ny:(node-2*ny+1)
    A(i,i)=const1; A(i,i+1)=2*const2; A(i,i-ny)=const3; A(i,i+ny)=const3;
end
for i=(node-tmp3+2*ny):ny:(node-ny)
    A(i,i)=const1; A(i,i-1)=2*const2; A(i,i-ny)=const3; A(i,i+ny)=const3;
end
i=node-tmp3+2+ny;
for ii=1:(length/3/hx-1)
    for jj=1:(ny-2)
        A(i+(jj-1)+(ii-1)*ny,i+(jj-1)+(ii-1)*ny)=const1; 
        A(i+(jj-1)+(ii-1)*ny,i+(jj-1)+(ii-1)*ny-ny)=const3;
        A(i+(jj-1)+(ii-1)*ny,i+(jj-1)+(ii-1)*ny+ny)=const3;
        A(i+(jj-1)+(ii-1)*ny,i+(jj-1)+(ii-1)*ny-1)=const2;
        A(i+(jj-1)+(ii-1)*ny,i+(jj-1)+(ii-1)*ny+1)=const2;
    end
end

%******************************************************************

% solve eigenvalue problem.
[v,d]=eigs(A,3,'SM');           % iterative method to solve eigenvalues.
lamda1=d(2,2);                  % get eigenvalues, excluding zero eigenvalue
lamda2=d(3,3);
f1=10e6*sqrt(-lamda1/const2/const3)/(2*pi*sqrt(4*(10e-7)*(10e-9)/36));
f2=10e6*sqrt(-lamda2/const2/const3)/(2*pi*sqrt(4*(10e-7)*(10e-9)/36));
hz1=v(:,2);     % first mode
hz2=v(:,3);     % second mode
hz1_x=zeros(node,1);
hz1_y=zeros(node,1);
hz2_x=zeros(node,1);
hz2_y=zeros(node,1);
%for i=1:ny
%    hx1(i,1)=0;
%end

%*******************************************************************
% calculate hz/delta-x, hz/delta-y
% choose first two modes.
for i=(ny+1):(ny*(length/3/hx+1)-ny)
    hz1_x(i) = (hz1(i+ny)-hz1(i-ny))/(2*hx);
    hz2_x(i) = (hz2(i+ny)-hz2(i-ny))/(2*hx);
end
for i=(tmp1+tmp2+1):(tmp1+tmp2+1+5.08/hy)
    hz1_x(i) = (hz1(i+ny-2.54/hy)-hz1(i-ny))/(2*hx);
    hz2_x(i) = (hz2(i+ny-2.54/hy)-hz2(i-ny))/(2*hx);
end
for i=(tmp3+1):(tmp3+1+5.08/hy)
    hz1_x(i) = (hz1(i+5.08/hy+1)-hz1(i-(ny-2.54/hy)))/(2*hx);
    hz2_x(i) = (hz2(i+5.08/hy+1)-hz2(i-(ny-2.54/hy)))/(2*hx);
end
for i=(tmp3+1+5.08/hy+1):(node-tmp3-2*tmp2-1)
    hz1_x(i) = (hz1(i+5.08/hy+1)-hz1(i-(5.08/hy+1)))/(2*hx);
    hz2_x(i) = (hz2(i+5.08/hy+1)-hz2(i-(5.08/hy+1)))/(2*hx);
end
for i=(node-tmp3-2*tmp2):(node-tmp3)
    hz1_x(i) = (hz1(i+5.08/hy+1+2.54/hy)-hz1(i-(5.08/hy+1)))/(2*hx);
    hz2_x(i) = (hz2(i+5.08/hy+1+2.54/hy)-hz2(i-(5.08/hy+1)))/(2*hx);
end
for i=tmp6:(tmp6+5.08/hy)
    hz1_x(i) = (hz1(i+ny)-hz1(i-(5.08/hy+1+2.54/hy)))/(2*hx);
    hz2_x(i) = (hz2(i+ny)-hz2(i-(5.08/hy+1+2.54/hy)))/(2*hx);
end
for i=(node-tmp3+1+ny):(node-ny)
    hz1_x(i) = (hz1(i+ny)-hz1(i-ny))/(2*hx);
    hz2_x(i) = (hz2(i+ny)-hz2(i-ny))/(2*hx);
end

for i=1:length/3/hx+1
    for j=2:ny-1
        hz1_y(j+ny*(i-1)) = (hz1(j+ny*(i-1)+1)-hz1(j+ny*(i-1)-1))/(2*hy);
        hz2_y(j+ny*(i-1)) = (hz2(j+ny*(i-1)+1)-hz2(j+ny*(i-1)-1))/(2*hy);
    end
end
for i=1:length/3/hx-1
    for j=tmp3+2:tmp3+5.08/hy
        hz1_y(j+(5.08/hy+1)*(i-1)) = (hz1(j+(5.08/hy+1)*(i-1)+1)-hz1(j+(5.08/hy+1)*(i-1)-1))/(2*hy);
        hz2_y(j+(5.08/hy+1)*(i-1)) = (hz2(j+(5.08/hy+1)*(i-1)+1)-hz2(j+(5.08/hy+1)*(i-1)-1))/(2*hy);
    end
end
for i=1:length/3/hx+1
    for j=node-tmp3+2:node-tmp3-1+ny
        hz1_y(j+ny*(i-1)) = (hz1(j+ny*(i-1)+1)-hz1(j+ny*(i-1)-1))/(2*hy);
        hz2_y(j+ny*(i-1)) = (hz2(j+ny*(i-1)+1)-hz2(j+ny*(i-1)-1))/(2*hy);
    end
end
hz1_y = -hz1_y; %E_x component
hz2_y = -hz2_y; 

% plot the geometry
xx=[0;0;7.62;7.62;15.24;15.24;22.86;22.86;15.24;15.24;7.62;7.62;0];
yy=[0;10.16;10.16;7.62;7.62;10.16;10.16;0;0;2.54;2.54;0;0];
axis([-1 24 -1 11]);
hold on;
line(xx,yy,'LineWidth',3);

% plot 1st mode
title('1st mode');
quiver(xval,yval,hz1_y,hz1_x);

% plot 2nd mode
%title('2nd mode');
%quiver(xval,yval,hz2_y,hz2_x);

time = etime(clock,t0);

% plot hz distribution using surf()
hz = hz1;       % you can change hz=hz2 for 2nd mode.
nxx = 7.62/hx+1;
nyy = ny;
xx=zeros(nyy,nxx);
yy=zeros(nyy,nxx);
zz=zeros(nyy,nxx);
for j=1:nyy
    for i=1:nxx
        xx(j,i) = (i-1)*hx;
        yy(j,i) = (j-1)*hy;
    end
end
for i=1:nxx
    for j=1:nyy
        zz(j,i) = hz(j+(i-1)*nyy);
    end
end
%surf(xx,yy,zz);
hold on;
nxx = 7.62/hx+1;
nyy = 5.08/hy+1;
xx=zeros(nyy,nxx);
yy=zeros(nyy,nxx);
zz=zeros(nyy,nxx);
for j=1:nyy
    for i=1:nxx
        xx(j,i) = 7.62+(i-1)*hx;
        yy(j,i) = 2.54+(j-1)*hy;
    end
end
for i=2:nxx-1
    for j=1:nyy
        zz(j,i) = hz(tmp3+j+(i-2)*nyy);
    end
end
for j=1:nyy
    zz(j,1) = hz(tmp1+tmp2+j);
    zz(j,nxx) = hz(tmp6+j-1);
end
%surf(xx,yy,zz);
hold on;
nxx = 7.62/hx+1;
nyy = ny;
xx=zeros(nyy,nxx);
yy=zeros(nyy,nxx);
zz=zeros(nyy,nxx);
for j=1:nyy
    for i=1:nxx
        xx(j,i) = 15.24+(i-1)*hx;
        yy(j,i) = (j-1)*hy;
    end
end
for i=1:nxx
    for j=1:nyy
        zz(j,i) = hz(node-tmp3+j+(i-1)*nyy);
    end
end
%surf(xx,yy,zz);