% Source code waveguide.m Courtesy: Sharmistha Dey
% this software computes the ratio of the cut-off wave length in a ridge wave guide
% with respect to the cut-off wave length of a rectangular wave guide
% the waveguide dimensions are a1Xb1 with b1=a1/2.
% the ridge dimensions are a2Xb2
% the tolerance in k is specified to be 0.1%.
% the function ridge_wg is used in the software.
clear all;
close all;
clf;
clc;
tic;
a1=0.02; % width of waveguide
b1=a1/2; % height of waveguide
lambda_c = 2*a1; % cut-off wavelength in the waveguide
lambda = [;];
i =1;
for beta=0.8:0.05:0.95 %beta is the ratio b2/b1
j=0;
for alpha=0.1:0.1:0.9 %alpha is the ratio a2/a1
j=j+1;
kc=ridge_wg(a1,b1,alpha,beta);
lambda1(i,j)=(2*pi)/kc;
a2(j)=alpha;
end
i=i+1;
end
lp_l=lambda1/lambda_c;
% plotting normalized cut-off wavenumber
plot(a2,lp_l(1,1:j),'bo-',a2,lp_l(2,1:j),'ro-',a2,lp_l(3,1:j),'mo-',a2,lp_l(4,(1:j)),'ko-');
legend('b_2/b_1=0.80','b_2/b_1=0.85','b_2/b_1=0.90','b_2/b_1=0.95');
title('Plot of \lambda_c^`/\lambda_c with a_2/a_1 and b_2/b_1');
xlabel('a_2/a_1');
ylabel('\lambda_c^`/\lambda_c');
hold off;
toc;
function kk=ridge_wg(a1,b1,alpha,beta)
% this function calculates the cut-off wavenumber k for the dominant mode in a ridge
% waveguide. The inputs to this function are the dimensions of the
% waveguide (a1 and b1), alpha=a2/a1; and beta=b2/b1.
% magnetic wall symmetry at a1/2 is used to reduce the size of the problem.
% square grid is employed.
% conventionally, the top left corner of the waveguide is assigned row
% number one and column number one
%SOR is not used
a2=alpha*a1;
b2=beta*b1;
Nx=200; %number of divisions along the x-axis
Ny=100; %number of divisions along the y-axis
nx=Nx+1; % number of nodes along the x-axis
ny=Ny+1; %number of nodes along the y-axis
nx_c=(nx+1)/2; %defines the symmetry plane
xridge=nx_c-(0.5*alpha*Nx); %defines the starting position of ridge from the first column
xridge=round(xridge);
yridge =ny-(beta*Ny); %defines the starting position of the ridge from the first row
yridge=round(yridge);
k=2*pi/a1; % cut-off wavenumer for TE10 mode in the rectangular waveguide
h=a1/Nx; % size of the square grid
hk2=(h*k)^2;
kc=0; % initialize kc
%symmetry about the mid-plane employed
H=zeros(ny,Nx/2+1); %initialize Hz over the nodes
H0=ones(ny,Nx/2+1); % initialization of Hz0 over the nodes
% =zeros(ny,nx); %initialize Hz over the nodes
% H0=ones(ny,nx); % initialization of Hz0 over the nodes
H(1:ny,1)=1; % boundary condition at the first column
while abs(k-kc)>=1e-4*abs(kc) %specifies tolerance in k
kc=k;
L=4-((h*k)^2);
phi_nu=0;
phi_dw=0;
nu=[];
dw=[];
R = [];
for x=1:5 % Computation of the node potential for the first 5 iterations
H0=H;
% corner nodes
H(1,1)=(2/L)*(H0(1,2)+H0(2,1));
H(ny,1)=(2/L)*(H0(ny-1,1)+H0(ny,2));
H(ny,xridge)=(2/L)*(H0(ny-1,xridge)+H0(ny,xridge-1));
H(yridge,xridge)=(1/L)*(H0(yridge-1,xridge)+H0(yridge+1,xridge)+H0(yridge,xridge-
1)+H0(yridge,xridge+1));
% boundary nodes
for j=2:nx/2
H(1,j)=(1/L)*(H0(1,j-1)+H0(1,j+1)+(2*H0(2,j))); % 1st row
if j H(ny,j)=(1/L)*(H0(ny,j-1)+H0(ny,j+1)+(2*H0(ny-1,j))); % last row
end
end
for i=2:ny-1 % 1st column
H(i,1)=(1/L)*(H0(i-1,1)+H0(i+1,1)+(2*H(i,2)));
end
for i=(yridge+1):(ny-1)
H(i,xridge)=(1/L)*(H0(i-1,xridge)+H0(i+1,xridge)+(2*H0(i,xridge-1)));
end
for j=xridge+1:nx/2
H(yridge,j)=(1/L)*(H0(yridge,j-1)+H0(yridge,j+1)+(2*H0(yridge-1,j)));
end
% other nodes
for i=2:(yridge-1)
for j=2:nx/2
H(i,j)=(1/L)*(H0(i-1,j)+H0(i+1,j)+H0(i,j-1)+H0(i,j+1)); % other columns
end
end
for i=yridge:(ny-1)
for j=2:nx/2
if j H(i,j)=(1/L)*(H0(i-1,j)+H0(i+1,j)+H0(i,j-1)+H0(i,j+1));
end
end
end
end % end of 5 iteration loop
nu(1:ny,101)=0;
% corner nodes
nu(1,1)=(2*H(1,2))+(2*H(2,1))-(4*H(1,1));
nu(ny,1)=(2*H(ny-1,1))+(2*H(ny,2))-(4*H(ny,1));
nu(ny,xridge)=(2*H(ny-1,xridge))+(2*H(ny,xridge-1))-(4*H(ny,xridge));
nu(yridge,xridge)=H(yridge-1,xridge)+H(yridge+1,xridge)+H(yridge,xridge-
1)+H(yridge,xridge+1)-(4*H(yridge,xridge));
% boundary nodes
for j=2:nx/2
if j nu(ny,j)=H(ny,j-1)+H(ny,j+1)+(2*H(ny-1,j))-(4*H(ny,j)); %last row
else
nu(1,j)=H(1,j-1)+H(1,j+1)+(2*H(2,j))-(4*H(1,j)); %first row
end
end
for i=2:ny-1
nu(i,1)=H(i-1,1)+H(i+1,1)+(2*H(i,2))-(4*H(i,1)); %first column
end
for i=yridge+1:ny-1
nu(i,xridge)=H(i-1,xridge)+H(i+1,xridge)+(2*H(i,xridge-1))-(4*H(i,xridge)); %side
wall of the ridge
end
for j=xridge+1:nx/2
nu(yridge,j)=H(yridge,j-1)+H(yridge,j+1)+(2*H(yridge-1,j))-(4*H(yridge,j)); %top wall
of the ridge
end
for i=2:yridge-1
for j=2:nx/2
nu(i,j)=H(i-1,j)+H(i+1,j)+H(i,j-1)+H(i,j+1)-(4*H(i,j)); %nodes between w/g top and
upto top of ridge
end
end
%nodes between the ridge and side wall
for i=yridge:ny-1
for j = 2:xridge-1
nu(i,j) =(H(i-1,j)+H(i+1,j)+H(i,j-1)+H(i,j+1)-4*H(i,j));
end
end
nu=H.*nu; % calculation of numerator
dw=H.*H; %calculation of denominator
phi_nu=sum(sum(nu)); %double summation for the numerator
phi_dw=sum(sum(dw)); %double summation for the denominator
hk2=(-1)*(phi_nu/phi_dw);
k=sqrt((hk2/(h^2)));
end
kk=k;
% Plot H field contours
% figure;
% [C,h] = contour(H);
%set(h,'ShowText','on','TextStep',get(h,'LevelStep')*10)