%notched rectangular waveguide cutoff frequency eigenvalue solver
%
%
%               <---an--->
%               |__       |        ______
%       *********  A      *********     A
%       *       *  |      *       *     |
%       *       *  | bn   *       *     |
%       *       *  V      *       *     |
%       *       ***********       *     b
%       *                         *     |
%       *                         *     |
%       *                         *     |
%       ***************************_____V
%       |                         |
%       <----------- a ----------->
%
%   Comment: In order for the program to work, 'bn' must be greater than
%   'an'
%               bn > an
%            
%   if bn < an, the functions force certain indices to zero and negative
%   numbers. These types of numbers are invalid for accessing adresses
%   within a matrix, and an error message informs the user

%enter dimensions (only built to enter even numbered values)
a=input('enter width of waveguide ');    %width
b=input('enter height of waveguide ');    %height

an=input('enter width of notch ');
while an >= a
    an = input('error! enter notch width smaller than waveguide width ');
end

bn=input('enter height of notch ');
while bn >= b
    bn = input('error! enter notch height smaller than waveguide height ');
end

%enter mesh settings
h=input('enter mesh size ');

%Break Notched wave guide into three areas
%
%   ******      ******
%   *  A *      *  B *
%   *    ********    *
%   *       C        *
%   ******************
%
% In the equal arm fd method, acessing the upper and lower arms will change
% based on differnt jmax due to the fact that the nodes only include the
% area within the wave guide. So the two areas will be:
%   Section I       Section II
%   ***********     **********
%   *  A * B  *     *   C    *
%   ***********     **********

% upper jmax
umax = (a-an)/h+2;

% base jmax
bmax = a/h+1;

%misc node values
imax = b/h+1;           %wg nodal height
inmax = bn/h;           %notch nodal height
ibmax = (b-bn)/h+1;     %lower nodal height
qmax = (a-an)/2/h+1;    %upper ridge width nodal index;
rmax = an/h-1;          %notch nodal width
smax = bn/h;            %notch nodal height

% total number of nodes
Kr = imax*bmax; %notchless waveguide
Kn = rmax*smax; %phantom nodes in notch
Ks = qmax*2*inmax; %Upper nodes
Kb = ((b-bn)/h+1)*bmax; %Lower nodes
K = Kr - Kn;

%Build A Matrix
A = eye(K)*4;
for k = 1:Ks
    if k-1 > 0
        A(k,k-1) = 1;
    end
    if k+1 < Ks
        A(k,k+1) =1;
    end
    if k-umax > 0
        A(k,k-umax) = 1;
    end
    if k+umax < Ks
        A(k,k+umax) = 1;
    end
end

%Enforce Boundary Conditions

%%%%%%%%%%%%%%%%%%%%%%% Section I %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%top
for n = 1:umax
    A(n,n+umax) = 2;
end

%left
for n = 1:inmax
    q = n*umax-(umax-1);
    if q-1 > 0
        A(q,q-1) = 0;
    end
    A(q,q+1) = 2; %doesn't need a conditional statement
end

%right
for n = 1:inmax
    q = n*umax;
    A(q,q+1) = 0;
    A(q,q-1) = 2;
end

%left side notch
for n = 1:inmax
    q = n*umax-qmax-rmax+1;
    A(q,q-1) = 2;
    A(q,q+1) = 0;
end

%right side notch
for n = 1:inmax
    q = n*umax-qmax+1;
    if q-1 > 0
        A(q,q-1) = 0;
    end
    A(q,q+1) = 2;
end

%bottom of notch
bl = umax*inmax+qmax;
for n = 1:rmax
    p = bl+n;
    A(p,p-umax) = 0;
    A(p,p-bmax) = 0;
    A(p,p+bmax) = 2;
end

%%%%%%%%%%%%%%%%%%%%%%%% Section II %%%%%%%%%%%%%%%%%%%%%%%%%%
%bottom
for n = 1:Kb
    k = Ks+n;
    if k-1 > 0
        A(k,k-1) = 1;
    end
    if k+1 < K
        A(k,k+1) = 1;
    end
    if k-bmax > 0
        A(k,k-bmax) = 1;
    end
    if k+bmax < K
        A(k,k+bmax) = 1;
    end
end

%lower left
for n = 1:ibmax
    q = Ks+n*bmax-(bmax-1);
    if q-1 > 0
        A(q,q-1) = 0;
    end
    if q+1 < K
        A(q,q+1) = 2;
    end
end

%lower right
for n = 1:ibmax
    q = Ks + n*bmax;
    if q+1 < K
        A(q,q+1) = 0;
    end
    if q-1 > 0
        A(q,q-1) = 2;
    end
end
%bottom
for n = 1:bmax
    p = K-bmax+n;
    A(p,p-bmax) = 2;
end

%Eigen Values

yc = eig(A)